text:
function rts_initRTSPropPanel(hRtsMainGUI, hPropFig, hClockLabel, connStatusLabel, vesselGUID, hResUITable1, hResSumLabel1, hResUITable2, hResSumLabel2, hResourceCombo1, hResourceCombo2, hMainEngineUITable, hRcsEngineUITable, resourcePanel1, resourcePanel2)
%rts_initRTSPropPanel Summary of this function goes here
%   Detailed explanation goes here

    hResSumProgBar1 = javacomponent('javax.swing.JProgressBar',get(hResSumLabel1,'Position'), resourcePanel1);
    hResSumProgBar2 = javacomponent('javax.swing.JProgressBar',get(hResSumLabel2,'Position'), resourcePanel2);

    hResSumProgBar1.setMinimum(0);
    hResSumProgBar1.setMaximum(100);
    hResSumProgBar2.setMinimum(0);
    hResSumProgBar2.setMaximum(100);

    delete(hResSumLabel1);
    delete(hResSumLabel2);

    timers = getappdata(hRtsMainGUI,'timers');
    utTcpipClient = instrfind('Name','TM_UT');
    resourcesTcpipClient = instrfind('Name','TM_VESSEL_RESOURCES');
    tcpipClientMainEng = instrfind('Name','TM_VESSEL_MAIN_ENGINES');
    tcpipClientRcsEng = instrfind('Name','TM_VESSEL_RCS_ENGINES');

    rtsOpts = getappdata(hRtsMainGUI,'RTSOptions');
    procExecPeriod = 1/str2double(rtsOpts.processFreq);
    procExecMode = rtsOpts.procExecMode;

    clockUpdtFunc = @(obj, event) rts_mainGUIClockUpdateFunc(obj, event, hClockLabel, utTcpipClient);
   	propClockTimer = timer('Name', 'RTSPropClockUpdate',...
                                    'ExecutionMode', 'fixedRate', ...   % Run timer repeatedly
                                    'Period', 1, ...                % Initial period is 1 sec.
                                    'StartDelay',0.1, ...
                                    'TimerFcn', clockUpdtFunc); % Specify callback
	set(propClockTimer,'UserData',{false}); %disallow user edits

    connUpdtFunc = @(obj, event) rts_pageGetConnStatus(obj, event, vesselGUID, hRtsMainGUI, connStatusLabel);
    connStopFcn = @(obj, event) rts_pageGetConnStatus(obj, event, vesselGUID, hRtsMainGUI, connStatusLabel, false);
   	propConnStatusTimer = timer('Name', 'RTSPropConnStatusUpdate',...
                                    'ExecutionMode', 'fixedRate', ...   % Run timer repeatedly
                                    'Period', 1, ...                % Initial period is 1 sec.
                                    'StartDelay',0.3, ...
                                    'StopFcn', connStopFcn, ...
                                    'TimerFcn', connUpdtFunc); % Specify callback
	set(propConnStatusTimer,'UserData',{false}); %disallow user edits

    resUpdt1Func = @(obj, event) rts_propUpdateFunc(obj, event, hPropFig, resourcesTcpipClient, tcpipClientMainEng, tcpipClientRcsEng, hResUITable1, hResSumProgBar1, hResUITable2, hResSumProgBar2, hResourceCombo1, hResourceCombo2, hMainEngineUITable, hRcsEngineUITable);
   	resUpt1Timer = timer('Name', 'RTSPowerThermResourceUpdate',...
                                    'ExecutionMode', procExecMode, ...   % Run timer repeatedly
                                    'Period', procExecPeriod, ...                % Initial period is 1 sec.
                                    'StartDelay',0.0, ...
                                    'TimerFcn', resUpdt1Func); % Specify callback
    set(resUpt1Timer,'UserData',{true}); %allow user edits

    start(propClockTimer);
    start(propConnStatusTimer);
    start(resUpt1Timer);

    timers(end+1) = propClockTimer;
    timers(end+1) = propConnStatusTimer;
    timers(end+1) = resUpt1Timer;
    setappdata(hRtsMainGUI,'timers', timers);

    timersProp = timer.empty(1,0);
    timersProp(end+1) = propClockTimer;
    timersProp(end+1) = propConnStatusTimer;
    timersProp(end+1) = resUpt1Timer;
    setappdata(hPropFig,'timers', timersProp);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [spec,idctm] = cep2spec(cep, nfreq, type)
% spec = cep2spec(cep, nfreq, type)
%   Reverse the cepstrum to recover a spectrum.
%   i.e. converse of spec2cep
%   nfreq is how many points to reconstruct in spec
% 2005-05-15 dpwe@ee.columbia.edu

if nargin < 2;   nfreq = 21;   end
if nargin < 3;   type = 2;   end   % type of DCT

[ncep,ncol] = size(cep);

% Make the DCT matrix
dctm = zeros(ncep, nfreq);
idctm = zeros(nfreq, ncep);
if type == 2 || type == 3
  % this is the orthogonal one, so inv matrix is same as fwd matrix
  for i = 1:ncep
    dctm(i,:) = cos((i-1)*[1:2:(2*nfreq-1)]/(2*nfreq)*pi) * sqrt(2/nfreq);
  end
  if type == 2
    % make it unitary! (but not for HTK type 3)
    dctm(1,:) = dctm(1,:)/sqrt(2);
  else
    dctm(1,:) = dctm(1,:)/2;
  end
  idctm = dctm';
elseif type == 4 % type 1 with implicit repetition of first, last bins
  % so all we do is reconstruct the middle nfreq rows of an nfreq+2 row idctm
  for i = 1:ncep
    % 2x to compensate for fact that only getting +ve freq half
    idctm(:,i) = 2*cos((i-1)*[1:nfreq]'/(nfreq+1)*pi);
  end
  % fixup 'non-repeated' basis fns
  idctm(:, [1 ncep]) = idctm(:, [1 ncep])/2;
else % dpwe type 1 - idft of cosine terms
  for i = 1:ncep
    % 2x to compensate for fact that only getting +ve freq half
    idctm(:,i) = 2*cos((i-1)*[0:(nfreq-1)]'/(nfreq-1)*pi);
  end
  % fixup 'non-repeated' basis fns
  idctm(:, [1 ncep]) = 0.5* idctm(:, [1 ncep]);
end

spec = exp(idctm*cep);

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
%% Arrayfun Example
clear all; close all; clc;

%% Example
% Run a function using the GPU by asking matlab to evaluate this function
% with data created inside our Gpu memory (gpuarray type data)
D = gpuDevice(1); % because, I have 2 GPUs in my PC.

% Create gpuarrays data,
s1 = gpuArray(rand(400));
s2 = gpuArray(rand(400));
s3 = gpuArray(rand(400));

% Evaluate aGpuFunction with gpuarray arguments,
[o1,o2] = arrayfun(@aGpuFunction,s1,s2,s3);

% We get a gpurray solution, therefore we need to gather it from gpuDevice
d = gather(o2);

% See variables
whos

% Clear GPU Memory
reset(D);

repo: Aero-matlab


text:
function [Psi, algParams, outParams] = initBPHMMSeq( data, model, initParams, algParams, outParams )
% Initialize BP-HMM hidden variables for MCMC sampling
%  via a **sequential** process
% That first fits a crude model to a subset of M sequences alone,
%  and then iteratively adds sequences ONE AT A TIME
%    trying to create unique features when appropriate
% Be warned: there is no guaranteed number of features that exist
%  after this procedure. It depends not only on number of sequences
%  but also on the complexity of those sequences

M = initParams.nSubsetObj;

% ------------------------------- Remember old state to use again afterward
curStream = RandStream.getGlobalStream();
entryState = curStream.State;

% Reset PRNG state to associate with specific *task*
reset( RandStream.getGlobalStream(), outParams.taskID );

% -------------------------------------------- Init BPHMM for first M sequences
if strcmp( class(data), 'SeqData' )
    Dsubset = SeqData();
elseif strcmp( class(data), 'ARSeqData' )
    Dsubset = ARSeqData( data.R );
end
subsetObjIDs = randsample( data.N, min(M,data.N)  );
remObjIDs    = setdiff( 1:data.N, subsetObjIDs );
for ii = subsetObjIDs'
   Dsubset = Dsubset.addSeq( data.seq(ii), data.seqNames{ii} );
end
outP = outParams;
outP.doPrintHeaderInfo = 0;

fprintf( 'Initializing: %d unique on subset of %d sequences\n', initParams.F.nUniquePerObj, M );

% Initialize model on first M sequences
%  using a few unique featuress per sequence
subPsi = initBPHMMFresh( Dsubset, model, initParams, algParams, outP );

if ~isempty( algParams.doAnneal )
   subPsi.invTemp = 0; % forget reversibility... give me a good init!
end

% Run through a few iters of full sampler on this small dataset
for warmiter = 1:5
    subPsi = sampleSharedFeats( subPsi, Dsubset );
    subPsi = sampleStateSeq( subPsi, Dsubset );

    subPsi = sampleSplitMerge_SeqAlloc( subPsi, Dsubset, algParams );

    subPsi.ThetaM = subPsi.ThetaM.sampleAllTheta( Dsubset, subPsi.stateSeq );
    subPsi.TransM = subPsi.TransM.sampleAllEta( subPsi.F, subPsi.stateSeq );

    % Try to delete extra features if possible
    algPdeath = algParams;
    algPdeath.Debug.MoveType = 0; %force death move!
    algPdeath.doAvoidCache   = 1;
    if isfield( subPsi, 'cache' )
        subPsi = rmfield( subPsi, 'cache' );
    end
    subPsi = sampleUniqueFeats( subPsi, Dsubset, algPdeath, 0 );
    subPsi = sampleStateSeq( subPsi, Dsubset );

    checkPsiFeatureConsistency( subPsi );

end

fprintf( 'Sampling features for remaining sequences\n' );
tic;
for rr = 1:length( remObjIDs )
    Dsubset = Dsubset.addSeq( data.seq( remObjIDs(rr) ), data.seqNames{remObjIDs(rr)} );
    curObjID = Dsubset.N;
    % Assume new object possesses all features
    %   and has prior mean on transition params
    subPsi.F(curObjID,:) = 1;
    subPsi.TransM = subPsi.TransM.getAllEta_PriorMean( 1:curObjID, subPsi.F, []);

    %if isfield( subPsi, 'cache' )
    %       subPsi = rmfield( subPsi, 'cache' );
    %end
    %subPsi = sampleStateSeq( subPsi, Dsubset, curObjID );
    %checkPsiFeatureConsistency( subPsi );

    % Perform MH updates to turn shared features on and off
    [subPsi] = sampleSharedFeats( subPsi, Dsubset, curObjID );

    % Attempt to do birth move using *entire sequence* data
    algPbirth = algParams;
    algPbirth.Debug.MoveType = 1;
    algPbirth.Debug.wstart   = 1;
    algPbirth.Debug.wend     = Dsubset.Ts( curObjID );
    [subPsi,~, ~] = sampleUniqueFeats( subPsi, Dsubset, algPbirth, 0, curObjID );
    subPsi.TransM.K = size( subPsi.F, 2 );
    subPsi.TransM.N = size( subPsi.F, 1 );

    % Every so often, resample ALL shared feat assignments + emit params
    if rr < 3 || mod(rr,5)==0 || rr == length(remObjIDs)
        if isfield( subPsi, 'cache' )
            subPsi = rmfield( subPsi, 'cache' );
        end
        subPsi = sampleSharedFeats( subPsi, Dsubset );
        % Try to do feature death moves
        subPsi = sampleUniqueFeats( subPsi, Dsubset, algPdeath, 0 );

        if isfield( subPsi, 'cache' )
            subPsi = rmfield( subPsi, 'cache' );
        end
        subPsi = sampleStateSeq( subPsi, Dsubset );
        checkPsiFeatureConsistency( subPsi );

        subPsi.ThetaM = subPsi.ThetaM.sampleAllTheta( Dsubset, subPsi.stateSeq );
        subPsi.TransM = subPsi.TransM.sampleAllEta( subPsi.F, subPsi.stateSeq );

        checkPsiFeatureConsistency( subPsi );
    end

    if mod(rr,20)==0
        fprintf( '    seq %4d/%4d  after %4.0f sec | nF= %d\n', rr, length(remObjIDs), toc, size( subPsi.F,2) );
    end
end

% Finally,
%   shuffle so that the data sequence is in the right (original) order
if isfield( subPsi, 'cache' )
    subPsi = rmfield( subPsi, 'cache' );
end
Psi = subPsi;
reorderIDs = [subsetObjIDs(:)' remObjIDs(:)'];
Psi.F = subPsi.F( reorderIDs, : );
Psi.stateSeq = subPsi.stateSeq( reorderIDs );

Psi.TransM = Psi.TransM.sampleAllEta( Psi.F, Psi.stateSeq );
Psi.ThetaM = Psi.ThetaM.sampleAllTheta( data, Psi.stateSeq );

checkPsiFeatureConsistency( Psi );
Psi = reallocateFeatIDs( Psi );
checkPsiFeatureConsistency( Psi );

% ---------------------------------------------------------  Reset stream
curStream = RandStream.getDefaultStream();
curStream.State = entryState;

author: michaelchughes

repo: NPBayesHMM


text:
function [p,x_min,upper] = initializesolution(p,x_min,upper)

if p.options.warmstart
    x = p.x0;
    z = evaluate_nonlinear(p,x);
    residual = constraint_residuals(p,z);
    relaxed_feasible = all(residual(1:p.K.f)>=-p.options.bmibnb.eqtol) & all(residual(1+p.K.f:end)>=-p.options.bmibnb.pdtol);
    if relaxed_feasible
        % upper_ = p.f+p.c'*z+z'*p.Q*z;
        upper_ = computecost(p.f,p.c,p.Q,z,p);
        if upper_ <= upper
            upper = upper_;
            x_min = z;
        end
    end
else
    x0 = p.x0;
    p.x0 = zeros(length(p.c),1);
    % Avoid silly warnings
    if ~isempty(p.evalMap)
        for i = 1:length(p.evalMap)
            if (isequal(p.evalMap{i}.fcn,'log') | isequal(p.evalMap{i}.fcn,'log2') | isequal(p.evalMap{i}.fcn,'log10'))
                p.x0(p.evalMap{i}.variableIndex) = (p.lb(p.evalMap{i}.variableIndex) +  p.ub(p.evalMap{i}.variableIndex))/2;
            end
        end
    end
    p = correctEXPConeClosureInitial(p);
    x = p.x0;
    z = evaluate_nonlinear(p,x);
    z = propagateAuxilliary(p,z);

    residual = constraint_residuals(p,z);
    relaxed_feasible = all(residual(1:p.K.f)>=-p.options.bmibnb.eqtol) & all(residual(1+p.K.f:end)>=-p.options.bmibnb.pdtol);
    if relaxed_feasible
        for i = 1:length(p.evalMap)
            if ~isempty(p.evalMap{i}.properties.forbidden)
                if z(p.evalMap{i}.variableIndex) > p.evalMap{i}.properties.forbidden(1) && z(p.evalMap{i}.variableIndex) < p.evalMap{i}.properties.forbidden(2)
                    relaxed_feasible = 0;
                end
            end
        end
    end
    if relaxed_feasible
        infs = isinf(z);
        if ~any(infs)
            % upper_ = p.f+p.c'*z+z'*p.Q*z;
            upper_ = computecost(p.f,p.c,p.Q,z,p);
            if upper_ <= upper
                upper = upper_;
                x_min = z;
                x0 = x_min;
            end
        else
            % Allow inf solutions if variables aren't used in objective
            if all(p.c(infs)==0) & nnz(p.Q(infs,:))==0
                ztemp = z;ztemp(infs)=0;
                %upper_ = p.f+p.c'*ztemp+ztemp'*p.Q*ztemp;
                upper_ = computecost(p.f,p.c,p.Q,ztemp,p);
                if upper_ <= upper
                    upper = upper_;
                    x_min = z;
                    x0 = x_min;
                end
            end
        end
    end
    p.x0 = (p.lb + p.ub)/2;
    p.x0(isinf(p.ub))=p.lb(isinf(p.ub));
    p.x0(isinf(p.lb))=p.ub(isinf(p.lb));
    p.x0(isinf(p.lb) & isinf(p.ub))=0;
    if ~isempty(p.integer_variables)
        p.x0(p.integer_variables) = round(p.x0(p.integer_variables));
    end
    if ~isempty(p.binary_variables)
        p.x0(p.binary_variables) = round(p.x0(p.binary_variables));
    end
    p = correctEXPConeClosureInitial(p);
    x = p.x0;
    x(isinf(x))=eps;
    x(isnan(x))=eps;
    z = evaluate_nonlinear(p,x);
    z = propagateAuxilliary(p,z);
    residual = constraint_residuals(p,z);
    relaxed_feasible = all(residual(1:p.K.f)>=-p.options.bmibnb.eqtol) & all(residual(1+p.K.f:end)>=p.options.bmibnb.pdtol);
    if relaxed_feasible
        tempcost = computecost(p.f,p.c,p.Q,z,p);
        if tempcost < upper
            upper = tempcost;
            x_min = z;
            x0 = x_min;
        end
    end
    p.x0 = x0;
end

function z = propagateAuxilliary(p,z)

try
    % New feature. If we introduce new variables xx = f(x) to be used
    % in a nonlinear operator, we can derive its value when x is chosen

    if ~isempty(p.aux_variables)
        if p.K.f > 1
            A = p.F_struc(1:p.K.f,2:end);
            b = p.F_struc(1:p.K.f,1);
            for i = 1:length(p.aux_variables)
                j = find(A(:,p.aux_variables(i)));
                if length(j)==1
                    if A(j,p.aux_variables(i))==1
                        z(p.aux_variables(i)) = -b(j)-A(j,:)*z;
                    end
                end
            end
            z = evaluate_nonlinear(p,z);
        end
    end
catch
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function mpc = t_case9_opfv2
%T_CASE9_OPFV2   Power flow data for 9 bus, 3 generator case, with OPF data.
%   Please see CASEFORMAT for details on the case file format.

%   MATPOWER

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [
	1	3	0	0	0	0	1	1	0	345	1	1.1	0.9;
	2	2	0	0	0	0	1	1	0	345	1	1.1	0.9;
	30	2	0	0	0	0	1	1	0	345	1	1.1	0.9;
	4	1	0	0	0	0	1	1	0	345	1	1.1	0.9;
	5	1	90	30	0	0	1	1	0	345	1	1.1	0.9;
	6	1	0	0	0	0	1	1	0	345	1	1.1	0.9;
	7	1	100	35	0	0	1	1	0	345	1	1.1	0.9;
	8	1	0	0	0	0	1	1	0	345	1	1.1	0.9;
	9	1	125	50	0	0	1	1	0	345	1	1.1	0.9;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	0	0	300	-300	1	100	1	250	90	0	0	0	0	0	0	0	0	0	0	0;
	30	85	0	300	-300	1	100	1	270	10	0	200	-30	30	-15	15	0	0	0	0	0;
	2	163	0	300	-300	1	100	1	300	10	0	200	-20	20	-10	10	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [
	1	4	0	0.0576	0	0	250	250	0	0	1	-360	2.48;
	4	5	0.017	0.092	0.158	0	250	250	0	0	1	-360	360;
	5	6	0.039	0.17	0.358	150	150	150	0	0	1	-360	360;
	30	6	0	0.0586	0	0	300	300	0	0	1	-360	360;
	6	7	0.0119	0.1008	0.209	40	150	150	0	0	1	-360	360;
	7	8	0.0085	0.072	0.149	250	250	250	0	0	1	-360	360;
	8	2	0	0.0625	0	250	250	250	0	0	1	-360	360;
	8	9	0.032	0.161	0.306	250	250	250	0	0	1	-360	360;
	9	4	0.01	0.085	0.176	250	250	250	0	0	1	-2	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	1	0	0	4	0	0	100	2500	200	5500	250	7250;
	1	0	0	3	0	0	200	3000	300	5000	0	0;
	2	0	0	2	24.035	-403.5	0	0	0	0	0	0;
];

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function dmapErrors = wrapperEvalSIRFS(classes)
    globals
    dmapErrors = cell(length(classes),1);

    map = load(fullfile(cachedir,'test_sets.mat'));
    map = map.map; map = map(~cellfun(@isempty,map));
    map = cellfun(@(x)(strcat(x,'.mat')),map,'UniformOutput',false);
    map = map(:);

    for ii=1:length(classes)
        dmapDir = fullfile(cachedir,classes{ii},'sirfsDir');
        gtDmapDir = fullfile(cachedir,classes{ii},'gtDepthMap');
        evalFileName = fullfile(cachedir,classes{ii},'evalStructSIRFS.mat');
        disp(classes{ii});
        if(exist(evalFileName,'file'))
            tmp = load(evalFileName);
            idx = ismember(tmp.fnames(:),map);
            dmapErrors{ii} = tmp.dmapErrors(idx);
            continue;
        end
        dmapErrors{ii} = evalSIRFS(dmapDir,gtDmapDir,evalFileName);
        disp(size(dmapErrors{ii},1));
    end
    fprintf('Class\t\tDepth Correlation (Mean)\n');
    for ii=1:length(classes)
        fprintf('%s\t\t%1.3f\n',classes{ii},mean(dmapErrors{ii}(~isnan(dmapErrors{ii}))));
    end

    fprintf('Class\t\tDepth Correlation (Median)\n');
    for ii=1:length(classes)
        fprintf('%s\t\t%1.3f\n',classes{ii},median(dmapErrors{ii}(~isnan(dmapErrors{ii})))*100);
    end

end

function dmapErrors = evalSIRFS(dmapDir,gtDmapDir,evalFileName)
    globals;
    fnames = getFileNamesFromDirectory(dmapDir,'types',{'.mat'});
    disp('Depth map evaluation');
    dmapErrors = zeros(length(fnames),1);
    p =  TimedProgressBar( length(fnames), round(length(fnames)/4), ...
    'Depth Map Errors: Remaining ', ', Completed ', 'Total Time: ' );

    parfor i=1:length(fnames)
        GtFile = fullfile(gtDmapDir,fnames{i});
        dmapFile = fullfile(dmapDir,fnames{i});
        try
            gtdmap = load(GtFile,'dmap');gtdmap = gtdmap.dmap;
        catch
            dmapErrors(i) = nan;
            continue;
        end
        dmap = load(dmapFile);
        mask = dmap.state.mask; dmap = dmap.state.dmap;
        %dmapErrors(i) = dmapMetric(dmap,gtdmap,mask);
        dmapErrors(i) = dmapMetricZMAE(dmap,gtdmap,mask);
        p.progress;
    end
    p.stop;
    save(evalFileName,'dmapErrors','fnames');
end

repo: CategoryShapes


text:
function uDot=uDotEllipsoid(u,x,a,f)
%%UDOTELLIPSOID Return the derivative of the basis vectors for the local
%               coordinate system of a moving observer when using an
%               ellipsoidal approximation of the curved Earth. Such an
%               evolving coordinate system has been termed "wander
%               coordinate" or "naturally evolving coordinates". Targets
%               moving locally in non-maneuvering, level flight with
%               a coordinate system evolving as per this coordinate system
%               will follow geodesic curves, not rhumb lines.
%
%INPUTS: u The 3X3 orthonormal basis vectors for the local coordinate
%          system. u(:,3) is the local vertical.
%        x An nX1 vector whose first three elements are Cartesian position
%          in the global ECEF coordinate system and whose next three
%          elements are velocity in the local (flat-Earth) coordinate
%          system. Other elements of x do not matter.
%        a The semi-major axis of the reference ellipsoid. If this argument
%          is omitted, the value in Constants.WGS84SemiMajorAxis is used.
%        f The flattening factor of the reference ellipsoid. If this
%          argument is omitted, the value in Constants.WGS84Flattening is
%          used.
%
%OUTPUTS: uDot The derivative of the basis vector u with respect to time.
%
%The solution is detailed in [1], and has been simplified for an
%ellipsoidal Earth so that numerical differentiation is not necessary.
%
%REFERENCES:
%[1] D. F. Crouse, "Simulating aerial targets in 3D accounting for the
%    Earth's curvature," Journal of Advances in Information Fusion, vol.
%    10, no. 1, Jun. 2015.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<4||isempty(f))
        f=Constants.WGS84Flattening;
    end

    if(nargin<3||isempty(a))
        a=Constants.WGS84SemiMajorAxis;
    end

    %The first numerical eccentricity of the ellipsoid.
    e=sqrt(2*f-f^2);

    %The current position.
    r=x(1:3,1);
    %The local velocity.
    rDotLocal=x(4:6,1);

    %Get the velocity in global coordinates.
    rDot=getGlobalVectors(rDotLocal,u);

    %Convert to ellipsoidal coordinates.
    rEllipse=Cart2Ellipse(r,[],a,f);

    %The ellipsoidal latitude
    phi=rEllipse(1);

    %The height above the reference ellipsoid.
    h=rEllipse(3);

    %Obtain the local East-North-Up axes.
    ENUAxes=getENUAxes(rEllipse,false,a,f);
    EAxis=ENUAxes(:,1);
    NAxis=ENUAxes(:,2);

    %The local North and East components of the velocity vector.
    vE=dot(EAxis,rDot(1:3,1));
    vN=dot(NAxis,rDot(1:3,1));

    %The normal radius of curvature.
    Rp=a/sqrt(1-e^2*sin(phi)^2);

    %The radius of curvature of the meridian line.
    Rm=Rp*(1-e^2)/(1-e^2*sin(phi)^2);

    %The rotation rate North.
    omegaN=vE/(Rp+h);

    %The rotation rate East.
    omegaE=-vN/(Rm+h);

    %The 3D rotation vector.
    Omega=omegaN*NAxis+omegaE*EAxis;

    %Compute the derivatives of the coordinate system basis vectors with
    %respect to time.
    uDot(:,1)=cross(Omega,u(:,1));
    uDot(:,2)=cross(Omega,u(:,2));
    uDot(:,3)=cross(Omega,u(:,3));
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:
%% Lucas-Kanade Tracker
%
% Lucas-Kanade sparse optical flow demo. Uses |cv.goodFeaturesToTrack| for
% track initialization and back-tracking for match verification between
% frames.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/lk_track.py>
%

%% Options
track_len = 10;       % max number of locations of point to remember
detect_interval = 5;  % detect new corners every X iterations

% params for corner detection and flow computation
shi_params = {'MaxCorners',500, 'QualityLevel',0.3, 'MinDistance',7, ...
    'BlockSize',7};
lk_params = {'WinSize',[15 15], 'MaxLevel',2, ...
    'Criteria',struct('type','Count+EPS', 'maxCount',10, 'epsilon',0.03)};

%% Video
% Prepare video source
if true
    vid = 0;
elseif true
    vid = fullfile(mexopencv.root(), 'test', '768x576.avi');
elseif mexopencv.require('vision')
    vid = fullfile(toolboxdir('vision'), 'visiondata', 'visiontraffic.avi');
end
cap = cv.VideoCapture(vid);
assert(cap.isOpened(), 'Failed to initialize capturing');

%% First frame
% Grab first frame
frame = cap.read();
assert(~isempty(frame), 'Failed to read frame');
gray0 = cv.cvtColor(frame, 'RGB2GRAY');

%% Initialize
% stores history of locations for a set of points (cell array of Nx2 matrices)
tracks = {};
clr = [0 255 0];

%%
% Plot
hImg = imshow(frame);
title('Lucas-Kanade tracker')

%% Main loop
counter = 0;  % iterations counter
while ishghandle(hImg)
    % Grab next frame
    frame = cap.read();
    if isempty(frame), break; end
    gray1 = cv.cvtColor(frame, 'RGB2GRAY');

    if ~isempty(tracks)
        % track last position of points, in forward and backward direction
        p0 = cellfun(@(tr) tr(end,:), tracks, 'UniformOutput',false);
        p1 = cv.calcOpticalFlowPyrLK(gray0, gray1, p0, lk_params{:});
        p0r = cv.calcOpticalFlowPyrLK(gray1, gray0, p1, lk_params{:});
        % keep only good matches
        good = cellfun(@(a,b) max(abs(a - b)), p0, p0r) < 1;
        tracks = tracks(good);
        p1 = p1(good);
        if any(good)
            % append new locations to existing tracked points
            tracks = cellfun(@(tr,p) [tr; p], tracks, p1, 'UniformOutput',false);
            % keep only the last 10 locations in each track (fixed size queue)
            idx = cellfun(@(tr) size(tr,1), tracks) > track_len;
            tracks(idx) = cellfun(@(tr) tr(2:end,:), tracks(idx), 'UniformOutput',false);
            % draw latest points and their tracks (comet-like plot)
            frame = cv.circle(frame, p1, 2, 'Thickness','Filled', 'Color',clr);
            frame = cv.polylines(frame, tracks, 'Closed',false, 'Color',clr);
        end
    end

    % display number of tracked points
    frame = cv.putText(frame, ...
        sprintf('track count: %d', numel(tracks)), [20 20], ...
        'FontScale',0.5, 'Color',[255 255 0], 'LineType','AA');

    if rem(counter, detect_interval) == 0 || isempty(tracks)
        % region of interest mask
        if ~isempty(tracks)
            % try to find new points by masking-out last track positions
            mask = 255 * ones(size(gray1), 'uint8');
            p = cellfun(@(tr) tr(end,:), tracks, 'UniformOutput',false);
            mask = cv.circle(mask, p, 5, 'Thickness','Filled', 'Color',0);
            roi_params = {'Mask',mask};
        else
            roi_params = {};
        end
        % detect corners
        p = cv.goodFeaturesToTrack(gray1, roi_params{:}, shi_params{:});
        if ~isempty(p)
            % append new set of points (p is a cell array of 1x2 vectors)
            tracks = [tracks, p];
        end
    end

    % Display result
    set(hImg, 'CData',frame);
    drawnow;

    % Next iteration
    counter = counter + 1;
    gray0 = gray1;
end
cap.release();

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function c = triu(a,k)
%TRIU         Implements  triu(a,k)  for intervals
%
%   c = triu(a,k)
%
% functionality as Matlab function triu for matrices
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump  improved speed
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if nargin==1
    k = 0;
  end

  c = a;
  if a.complex
    c.mid = triu(a.mid,k);
    c.rad = triu(a.rad,k);
  else
    c.inf = triu(a.inf,k);
    c.sup = triu(a.sup,k);
  end

repo: OpenMAS


text:
% BilinearClPooling is the dagnn wapper of vl_nnbilinearclpool which
% computes outer product of outputs of two layers and pool the features
% across all locations

% Copyright (C) 2015 Tsung-Yu Lin, Aruni RoyChowdhury, Subhransu Maji.
% All rights reserved.
%
% This file is part of the BCNN and is made available under
% the terms of the BSD license (see the COPYING file).

classdef BilinearClPooling < dagnn.Filter
  properties
    method = 'sum'
    normalizeGradients = false;
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnbilinearclpool(inputs{1}, inputs{2});
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      [derInputs{1}, derInputs{2}] = vl_nnbilinearclpool(inputs{1}, inputs{2}, derOutputs{1});
      if obj.normalizeGradients
        gradNorm = sum(abs(derInputs{1}(:))) + 1e-8;
        derInputs{1} = derInputs{1}/gradNorm;

        gradNorm = sum(abs(derInputs{2}(:))) + 1e-8;
        derInputs{2} = derInputs{2}/gradNorm;
      end
      derParams = {} ;
    end

    function rfs = getReceptiveFields(obj)
      rfs(1,1).size = [NaN NaN] ;
      rfs(1,1).stride = [NaN NaN] ;
      rfs(1,1).offset = [NaN NaN] ;
      rfs(2,1) = rfs(1,1) ;
    end

    function obj = BilinearClPooling(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: DBCNN


text:
function [net, info] = cnn_mnist(varargin)
%CNN_MNIST  Demonstrates MatConvNet on MNIST

run(fullfile(fileparts(mfilename('fullpath')),...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

opts.batchNormalization = false ;
opts.networkType = 'simplenn' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.networkType ;
if opts.batchNormalization, sfx = [sfx '-bnorm'] ; end
opts.expDir = fullfile(vl_rootnn, 'data', ['mnist-baseline-' sfx]) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.dataDir = fullfile(vl_rootnn, 'data', 'mnist') ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.train = struct() ;
opts = vl_argparse(opts, varargin) ;
if ~isfield(opts.train, 'gpus'), opts.train.gpus = [2]; end;

% --------------------------------------------------------------------
%                                                         Prepare data
% --------------------------------------------------------------------

net = cnn_mnist_init('batchNormalization', opts.batchNormalization, ...
                     'networkType', opts.networkType) ;

if exist(opts.imdbPath, 'file')
  imdb = load(opts.imdbPath) ;
else
  imdb = getMnistImdb(opts) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

net.meta.classes.name = arrayfun(@(x)sprintf('%d',x),1:10,'UniformOutput',false) ;

% --------------------------------------------------------------------
%                                                                Train
% --------------------------------------------------------------------

switch opts.networkType
  case 'simplenn', trainfn = @cnn_train ;
  case 'dagnn', trainfn = @cnn_train_dag ;
end

[net, info] = cnn_train_binary(net, imdb, getBatch(opts), ...
  'expDir', opts.expDir, ...
  net.meta.trainOpts, ...
  opts.train, ...
  'val', find(imdb.images.set == 3)) ;

% --------------------------------------------------------------------
function fn = getBatch(opts)
% --------------------------------------------------------------------
switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(x,y) ;
  case 'dagnn'
    bopts = struct('numGpus', numel(opts.train.gpus)) ;
    fn = @(x,y) getDagNNBatch(bopts,x,y) ;
end

% --------------------------------------------------------------------
function [images, labels] = getSimpleNNBatch(imdb, batch)
% --------------------------------------------------------------------
images = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;

% --------------------------------------------------------------------
function inputs = getDagNNBatch(opts, imdb, batch)
% --------------------------------------------------------------------
images = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;
if opts.numGpus > 0
  images = gpuArray(images) ;
end
inputs = {'input', images, 'label', labels} ;

% --------------------------------------------------------------------
function imdb = getMnistImdb(opts)
% --------------------------------------------------------------------
% Preapre the imdb structure, returns image data with mean image subtracted
files = {'train-images-idx3-ubyte', ...
         'train-labels-idx1-ubyte', ...
         't10k-images-idx3-ubyte', ...
         't10k-labels-idx1-ubyte'} ;

if ~exist(opts.dataDir, 'dir')
  mkdir(opts.dataDir) ;
end

for i=1:4
  if ~exist(fullfile(opts.dataDir, files{i}), 'file')
    url = sprintf('http://yann.lecun.com/exdb/mnist/%s.gz',files{i}) ;
    fprintf('downloading %s\n', url) ;
    gunzip(url, opts.dataDir) ;
  end
end

f=fopen(fullfile(opts.dataDir, 'train-images-idx3-ubyte'),'r') ;
x1=fread(f,inf,'uint8');
fclose(f) ;
x1=permute(reshape(x1(17:end),28,28,60e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 't10k-images-idx3-ubyte'),'r') ;
x2=fread(f,inf,'uint8');
fclose(f) ;
x2=permute(reshape(x2(17:end),28,28,10e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 'train-labels-idx1-ubyte'),'r') ;
y1=fread(f,inf,'uint8');
fclose(f) ;
y1=double(y1(9:end)')+1 ;

f=fopen(fullfile(opts.dataDir, 't10k-labels-idx1-ubyte'),'r') ;
y2=fread(f,inf,'uint8');
fclose(f) ;
y2=double(y2(9:end)')+1 ;

set = [ones(1,numel(y1)) 3*ones(1,numel(y2))];
data = single(reshape(cat(3, x1, x2),28,28,1,[]));
dataMean = mean(data(:,:,:,set == 1), 4);
data = bsxfun(@minus, data, dataMean) ;

imdb.images.data = data ;
imdb.images.data_mean = dataMean;
imdb.images.labels = cat(2, y1, y2) ;
imdb.images.set = set ;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = arrayfun(@(x)sprintf('%d',x),0:9,'uniformoutput',false) ;

author: jiangqy


text:
function [z,dz,ymu,ys,fmu,fs,fpi] = acqNegGEI(xi,target,gpstruct,optimState,grad)
%ACQNEGEI Acquisition function for (negative) generalized expected improvement (unsupported).

if nargin < 5 || isempty(grad); grad = 0; end

n = size(xi,1);

if grad == 1 && n > 1
    error('acqNegGEI:gradient', ...
        'Gradient of acquisition function is provided only at one test point XI (row vector).');
end

if grad
    [ymu,ys2,fmu,fs2,hypw,dymu,dys2,dfmu,dfs2] = gppred(xi,gpstruct,'central');
else
    [ymu,ys2,fmu,fs2,hypw] = gppred(xi,gpstruct);
end
fs = sqrt(fs2);
ys = sqrt(ys2);

% Probability of improvement
gammaz = real((target - fmu)./fs);
fpi = 0.5*erfc(-gammaz/sqrt(2));

% Squared expected improvement
z1 = -(fs.*(gammaz.*fpi + exp(-0.5*(gammaz.^2))/sqrt(2*pi))).^2;

% Expected squared improvement
z2 = -fs.^2.*((gammaz.^2+1).*fpi + gammaz.*exp(-0.5*(gammaz.^2))/sqrt(2*pi));

% gamma(1,1,:) = linspace(0,1,11);
gamma(1,1,:) = 0.5;
z = bsxfun(@times,gamma,z1) + bsxfun(@times,1-gamma,z2);

z = sum(bsxfun(@times,hypw(~isnan(hypw)),z(~isnan(hypw),:,:)),1);

[zi,idx] = min(z,[],3);
[~,idx2] = min(zi,[],2);
idx = idx(idx2);
z = z(:,:,idx);

%try
%    z = sum(bsxfun(@times,hypw(~isnan(hypw)),z(~isnan(hypw),:)),1);
%catch
%    z = Inf(1,n);
%    dz = NaN(n,size(xi,2));
%    return;
%end

if grad
    % Gradient of probability of improvement
    dfs = 0.5*dfs2./fs;
    dgammaz = -(dfmu.*fs + (target - fmu).*dfs)./fs2;
    dfpi = -0.5*dgammaz/sqrt(2)*(-2*exp(-gammaz.^2/2)/sqrt(pi));

    % Gradient of expected improvement
    %dz = -(dfs.*gammaz.*fpi + dgammaz.*fs.*fpi + dfpi.*gammaz.*fs) ...
    %    + (fs.*gammaz.*dgammaz - dfs).*exp(-0.5*gammaz.^2)/sqrt(2*pi);
    %dz = sum(bsxfun(@times,hypw,dz(~isnan(hypw),:)),1);
else
    dz = NaN(n,size(xi,2));     % Gradient not estimated
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [frames, descrs] = vl_phow(im, varargin)
% VL_PHOW  Extract PHOW features
%   [FRAMES, DESCRS] = VL_PHOW(IM) extracts PHOW features from the
%   image IM. This function is a wrapper around VL_DSIFT() and
%   VL_IMSMOOTH().
%
%   The PHOW descriptors where introduced in [1]. By default,
%   VL_PHOW() computes the gray-scale variant of the descriptor.  The
%   COLOR option can be used to compute the color variant instead.
%
%   Verbose:: [false]
%     Set to true to turn on verbose output.
%
%   Sizes:: [[4 6 8 10]]
%     Scales at which the dense SIFT features are extracted. Each
%     value is used as bin size for the VL_DSIFT() function.
%
%   Fast:: [true]
%     Set to false to turn off the fast SIFT features computation by
%     VL_DSIFT().
%
%   Step:: [2]
%     Step (in pixels) of the grid at which the dense SIFT features
%     are extracted.
%
%   Color:: [GRAY]
%     Choose between GRAY (PHOW-gray), RGB, HSV, and OPPONENT
%     (PHOW-color).
%
%   ContrastThreshold:: [0.005]
%     Contrast threshold below which SIFT features are mapped to
%     zero. The input image is scaled to have intensity range in [0,1]
%     (rather than [0,255]) and this value is compared to the
%     descriptor norm as returned by VL_DSIFT().
%
%   WindowSize:: [1.5]
%     Size of the Gaussian window in units of spatial bins.
%
%   Magnif:: [6]
%     The image is smoothed by a Gaussian kernel of standard deviation
%     SIZE / MAGNIF.
%
%   FloatDescriptors:: [false]
%     If set to TRUE, the descriptors are returned in floating point
%     format.
%
%   See also: VL_DSIFT(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% -------------------------------------------------------------------
%                                                 Parse the arguments
% -------------------------------------------------------------------

  opts.verbose = false ;
  opts.fast = true ;
  opts.sizes = [4 6 8 10] ;
  opts.step = 2 ;
  opts.color = 'gray' ;
  opts.floatdescriptors = false ;
  opts.magnif = 6 ;
  opts.windowsize = 1.5 ;
  opts.contrastthreshold = 0.005 ;
  opts = vl_argparse(opts,varargin) ;

  dsiftOpts = {'norm', 'windowsize', opts.windowsize} ;
  if opts.verbose, dsiftOpts{end+1} = 'verbose' ; end
  if opts.fast, dsiftOpts{end+1} = 'fast' ; end
  if opts.floatdescriptors, dsiftOpts{end+1} = 'floatdescriptors' ; end
  dsiftOpts(end+(1:2)) = {'step', opts.step} ;

% -------------------------------------------------------------------
%                                                Extract the features
% -------------------------------------------------------------------

  % standarize the image
  imageSize = [size(im,2) ; size(im,1)] ;
  if strcmp(lower(opts.color), 'gray')
    numChannels = 1 ;
    if size(im,3) > 1, im = rgb2gray(im) ; end
  else
    numChannels = 3 ;
    if size(im,3) == 1, im = cat(3, im, im, im) ; end
    switch lower(opts.color)
      case 'rgb'
      case 'opponent'
        % Note that the mean differs from the standard definition of opponent
        % space and is the regular intesity (for compatibility with
        % the contrast thresholding).
        %
        % Note also that the mean is added pack to the other two
        % components with a small multipliers for monochromatic
        % regions.
        mu = 0.3*im(:,:,1) + 0.59*im(:,:,2) + 0.11*im(:,:,3) ;
        alpha = 0.01 ;
        im = cat(3, mu, ...
                 (im(:,:,1) - im(:,:,2))/sqrt(2) + alpha*mu, ...
                 (im(:,:,1) + im(:,:,2) - 2*im(:,:,3))/sqrt(6) + alpha*mu) ;
      case 'hsv'
        im = rgb2hsv(im) ;
      otherwise
        opts.color = 'hsv' ;
        warning('Color space not recongized, defaulting to HSV color space.') ;
    end
  end

  if opts.verbose
    fprintf('%s: color space: %s\n', mfilename, opts.color) ;
    fprintf('%s: image size: %d x %d\n', mfilename, imageSize(1), imageSize(2)) ;
    fprintf('%s: sizes: [%s]\n', mfilename, sprintf(' %d', opts.sizes)) ;
  end

  for si = 1:length(opts.sizes)

    % Recall from VL_DSIFT() that the first descriptor for scale SIZE has
    % center located at XC = XMIN + 3/2 SIZE (the Y coordinate is
    % similar). It is convenient to align the descriptors at different
    % scales so that they have the same geometric centers. For the
    % maximum size we pick XMIN = 1 and we get centers starting from
    % XC = 1 + 3/2 MAX(OPTS.SIZES). For any other scale we pick XMIN so
    % that XMIN + 3/2 SIZE = 1 + 3/2 MAX(OPTS.SIZES).
    %
    % In pracrice, the offset must be integer ('bounds'), so the
    % alignment works properly only if all OPTS.SZES are even or odd.

    off = floor(1 + 3/2 * (max(opts.sizes) - opts.sizes(si))) ;

    % smooth the image to the appropriate scale based on the size
    % of the SIFT bins
    sigma = opts.sizes(si) / opts.magnif ;
    ims = vl_imsmooth(im, sigma) ;

    % extract dense SIFT features from all channels
    for k = 1:numChannels
      [f{k}, d{k}] = vl_dsift(...
        ims(:,:,k), ...
        dsiftOpts{:},  ...
        'size', opts.sizes(si), ...
        'bounds', [off off +inf +inf]) ;
    end

    % remove low contrast descriptors
    % note that for color descriptors the V component is
    % thresholded
    switch opts.color
      case {'gray', 'opponent'}
        contrast = f{1}(3,:) ;
      case 'rgb'
        contrast = mean([f{1}(3,:) ; f{2}(3,:) ; f{3}(3,:)],1) ;
      otherwise % hsv
        contrast = f{3}(3,:) ;
    end
    for k = 1:numChannels
      d{k}(:, contrast < opts.contrastthreshold) = 0 ;
    end

    % save only x,y, and the scale
    frames{si} = [f{1}(1:3, :) ; opts.sizes(si) * ones(1,size(f{1},2))] ;
    descrs{si} = cat(1, d{:}) ;
  end
  descrs = cell2mat(descrs) ;
  frames = cell2mat(frames) ;
end

author: Cloud-CV

repo: object-proposals


text:
function cost = costfunction(mpcModel, u, varargin)
%fst.runningcosts, fst.terminalcosts, fst, fst.horizon, fst.xmeasure, u, fst.price, fst.battery
%UNTITLED3 Summary of this function goes here
%	SOME WORD

    cost = 0;
    % x = zeros(fst.horizon+1, length(fst.xmeasure));
    x = computeOpenloopSolution(mpcModel , u );

    for k=1:mpcModel.horizon
        cost = cost+mpcModel.runningcosts(k, x(k,:), u(:,k), mpcModel); %, varargin
    end
        cost = cost+mpcModel.terminalcosts( mpcModel.horizon+1, x(mpcModel.horizon+1,:), mpcModel );
end

author: juchengquan


text:
function [loser] = wholoses(plr,box)
    loser=0; % "0" means no loser yet
    % the player having no piece on table will lose
    for x=1:8
        for y=1:8
            if plr==1 && box(x,y)==10   %king1 is present KING=10
                loser=0;
                return;
            elseif plr==2 && box(x,y)==-10 %king2 is present king=-10
                loser=0;
                return;
            else
                loser=plr;
            end
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function header = egi_read_header(fid)
% EGI_READ_HEADER: Read header from an open EGI .raw file.
%
% USAGE:  header = egi_read_header(fid);

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2010

% ===== READ RAW HEADER =====
header.versionNumber          = fread(fid, 1, 'integer*4');
header.recordingTime.Year     = fread(fid, 1, 'integer*2');
header.recordingTime.Month    = fread(fid, 1, 'integer*2');
header.recordingTime.Day      = fread(fid, 1, 'integer*2');
header.recordingTime.Hour     = fread(fid, 1, 'integer*2');
header.recordingTime.Minute   = fread(fid, 1, 'integer*2');
header.recordingTime.Second   = fread(fid, 1, 'integer*2');
header.recordingTime.Millisec = fread(fid, 1, 'integer*4');
header.samplingRate           = fread(fid, 1, 'integer*2');
header.numChans               = fread(fid, 1, 'integer*2');
header.boardGain              = fread(fid, 1, 'integer*2');
header.numConvBits            = fread(fid, 1, 'integer*2');
header.ampRange               = fread(fid, 1, 'integer*2');

% ===== DATA FORMAT =====
switch header.versionNumber
    case {2,3}
        header.byteformat = 'integer*2';
        header.bytesize   = 2;
    case {4,5}
        header.byteformat = 'real*4';
        header.bytesize   = 4;
    case {6,7}
        header.byteformat = 'real*8';
        header.bytesize   = 8;
    otherwise
        fclose(fid);
        error('Brainstorm:InvalidRawFile', 'Error (Version ID): %d', header.versionNumber);
end
isReadCell = ismember(header.versionNumber, [3 5 7]);
if isReadCell
    fclose(fid);
    error('Support only for continuous EGI RAW files.');
end

% % ===== READ CELL NAMES =====
% if isReadCell
%     header.numCell = fread(fid, 1, 'int16');
%     % Read cells info
%     header.cellCodes = cell(1,header.numCell);
%     for i = 1 : header.numCell
%         codeSize = fread(fid, 1, 'int8');
%         header.cellCodes{i} = fread(fid, [1, codeSize], '*char');
%     end
%     metaData.nSeg = fread(fid, 1, 'int16');
% else
%     header.numCell = 0;
%     header.nSeg = 0;
%     header.cellCodes = {};
% end

% ===== READ EVENTS =====
header.numSamples = fread(fid, 1, 'integer*4');
header.numEvents  = fread(fid, 1, 'integer*2');
% File contains event info.
if (header.numEvents ~= 0)
    header.eventCodes = cell(header.numEvents, 1);
    for i = 1:header.numEvents
        header.eventCodes{i} = fread(fid, [1 4], '*char');
    end
else  % File does not contain event info.
    % warning('Brainstorm:IncompleteRawFile', 'RAW File does not contain event information');
    header.evtCode = [];
end

% Save beginning of the data block
header.datapos = ftell(fid);

author: brainstorm-tools


text:
function [U,Sigma,V,numiter,out]  = SVT(n,Omega,b,tau,delta,maxiter,tol,EPS)
% [U,Sigma,V,numiter,output]  = SVT(n,Omega,b,tau,delta,maxiter,tol,EPS)
%
% Finds the minimum of   tau ||X||_* + .5 || X ||_F^2
%
% subject to P_Omega(X) = P_Omega(M)
%
% using linear Bregman iterations
%
% Usage:  [U,S,V,numiter]  = SVT(n,Omega,b,delta,maxiter,tol)
%
% Inputs:
%
%   n - size of the matrix X assumed n(1) by n(2). If n is a single integer, it
% is understood that n(1) = n(2).
%
%   Omega - set of observed entries.  Should be linearly indexed.
%
%   b - data vector of the form M(Omega)
%
%   tau - parameter defining the objective functional
%
%   delta - step size.  Choose delta less than 2 to be safe but
%       conservative; choose delta closer to n(1)*n(2)/length(Omega)
%       to be riskier (i.e. algorithm may diverge)
%
%   maxiter - maximum number of iterations
%
%   tol - stopping criteria (default: 1e-4)
%
%   EPS - noise constraint.  This relaxes the constraints, so that they
%       are now of the form | X(i,j) - M(i,j) | <= EPS,
%       for all indices (i,j) in omega.  Default: 0
%
% Outputs: matrix X stored in SVD format X = U*diag(S)*V'
%
%   U - n1xr left singular vectors
%
%   S - rx1 singular values
%
%   V - n2xr right singular vectors
%
%   numiter - number of iterations to achieve convergence
%
%   output - a structure with data from each iteration.  Includes:
%       output.nuclearNorm  - nuclear norm of current iterate
%       output.rank         - rank of current iterate
%       output.time         - time taken for one iteraration
%       output.residual     - the relative residual, norm(x-b)/norm(b)
% Description:
% Reference:
%
%    Cai, Candes and Shen
%    A singular value thresholding algorithm for matrix completion
%    Submitted for publication, October 2008.
%
%    See also more general code as part of the TFOCS package,
%    available at tfocs.stanford.edu as of November 2010.
%
% Written by: Emmanuel Candes
% Email: emmanuel@acm.caltech.edu
% Created: October 2008
% Efficient mex-file and PROPACK version: Stephen Becker, Nov 2008
% Modified: Stephen Becker, March 2009
% Modified: Stephen Becker, May 2009  works with complex numbers
% Modified: Farshad Harirchi and Stephen Becker, April 2011, fixing a bug.

global VERBOSE
if isempty(VERBOSE)
    % -- feel free to change these 'verbosity' parameters
    % VERBOSE = false;
    VERBOSE = 1;    % a little bit of output
    % VERBOSE = 2;    % even more output
end

time1 = cputime;
if nargin < 8 || isempty(EPS)
    EPS = false;
end
if nargin < 7 || isempty(tol)
    tol = 1e-4;
end
if nargin < 6 || isempty(maxiter)
    maxiter = 500;
end

if length(n) == 1,
    n1 = n(1); n2 = n1;
elseif length(n) == 2,
    n1 = n(1); n2 = n(2);
end
if n1*n2 < 100*100, SMALLSCALE = true; else SMALLSCALE = false; end

m = length(Omega); [temp,indx] = sort(Omega);
% simpler: sort b also
incre = 5;
normb = norm(b);

[i, j] = ind2sub([n1,n2], Omega);
USE_SLOW_UPDATE     = false;
if EPS
    % with inequality constraints, should take delta = delta/sqrt(2) at
    % least, or delta = delta/2
    delta = delta/sqrt(2);
%     delta = delta/2; tau = 2*tau;
%     y1 = max(b-EPS,0); y2 = max(-b-EPS,0); % doesn't work well
    y1 = max(b,0); y2 = max(-b,0);
    Y = sparse(i,j,y1-y2,n1,n2,m);
    normProjM = normest(Y,1e-2);
    k0 = ceil(tau/(delta*normProjM));
    y1 = k0*delta*y1;
    y2 = k0*delta*y2;
    try
        updateSparse(Y,y1-y2,indx);
    catch
        l = lasterror;
        if strcmpi( l.identifier, 'MATLAB:UndefinedFunction')
            % mex file not installed, so do this instead:
            [indx_i,indx_j,s] = find(Y);
            Y = updateSparse_slow(Y,y1-y2,indx,indx_i,indx_j);
            USE_SLOW_UPDATE     = true;
        else
            % some other error (unexpected)
            rethrow(lasterror)
        end
    end
else
    Y = sparse(i,j,b,n1,n2,m);
    normProjM = normest(Y,1e-2);
    k0 = ceil(tau/(delta*normProjM));
    normb = norm(b);
    y = k0*delta*b; % kicking by k0 steps
    try
        updateSparse(Y,y,indx);
    catch
        l = lasterror;
        if strcmpi( l.identifier, 'MATLAB:UndefinedFunction')
            % mex file not installed, so do this instead:
            [indx_i,indx_j,s] = find(Y);
            Y = updateSparse_slow(Y,y,indx,indx_i,indx_j);
            USE_SLOW_UPDATE     = true;
        else
            % some other error (unexpected)
            rethrow(lasterror)
        end
    end

end
r = 0;

out.residual = zeros(maxiter,1);
out.rank= zeros(maxiter,1);
out.time = zeros(maxiter,1);
out.nuclearNorm = zeros(maxiter,1);

% What the best way to multiply a sparse matrix?
[forwardType, transposeType] = findBestMultiply(Y,.2);

if VERBOSE==1, fprintf('\nIteration:   '); end
for k = 1:maxiter,
    if VERBOSE==1, fprintf('\b\b\b\b%4d',k);  end
    s = r + 1;

    rInc = 4;  % make this larger for more accuracy
    %if tol < 1e-4  && relRes < 1e-1
    %rInc = rInc + max( round(log10( abs(1e-1/relRes) )), 5 );
    %end
    s = min( [r + rInc, n1, n2] );

    if SMALLSCALE
        [U,Sigma,V] = svd(full(Y),'econ');
    else
        % Make routines for multiplying by a sparse matrix
        Yt = Y';
        switch forwardType
            case 1, Yforward = @(x) Y*x;
            case 2, Yforward = @(x) Yt'*x;
            case 3, Yforward = @(x) smvp(Y,x);
        end
        switch transposeType
            case 1, Ytranspose = @(x) Yt*x;
            case 2, Ytranspose = @(x) Y'*x;
            case 3, Ytranspose = @(x) smvp(Yt,x);
        end
        OK = 0;
        while ~OK
            opts = [];
            if ~isreal(b), opts.eta = 1e-16; end
            [U,Sigma,V] = lansvd(Yforward,Ytranspose,n1,n2,s,'L',opts);
            %[U,Sigma,V] = lansvd(Y,s,'L');
            OK = (Sigma(s,s) <= tau) || ( s == min(n1,n2) );
            s = min(s + incre, min(n1,n2));
        end
    end

    sigma = diag(Sigma); r = sum(sigma > tau);
    U = U(:,1:r); V = V(:,1:r); sigma = sigma(1:r) - tau; Sigma = diag(sigma);

    x = XonOmega(U*diag(sigma),V,Omega);
    eTime = cputime - time1;
    if VERBOSE == 2
        fprintf('iteration %4d, rank is %2d, rel. residual is %.1e\n',k,r,norm(x-b)/normb);
    end
    relRes = norm(x-b)/normb;
    out.residual(k) = relRes;
    out.time(k) = eTime;
    out.rank(k) = r;
    out.nuclearNorm(k) = sum(sigma);

    time1 = cputime;

    if (relRes < tol)
        break
    end
    if EPS && norm(x-b,'inf') < 2*EPS
        break
    end
    if (norm(x-b)/normb > 1e5)
        disp('Divergence!');
        break
    end

    if EPS
        y1 = max( y1 + delta*( -(x-b) - EPS), 0 );
        y2 = max( y2 + delta*(  (x-b) - EPS), 0 );
        if USE_SLOW_UPDATE
            % mex file not installed, so do this instead:
            Y = updateSparse_slow(Y,y1-y2,indx,indx_i,indx_j);
        else
            updateSparse(Y,y1-y2,indx);
        end
    else
        y = y + delta*(b-x);
        if USE_SLOW_UPDATE
            % mex file not installed, so do this instead:
            Y = updateSparse_slow(Y,y,indx,indx_i,indx_j);
        else
            updateSparse(Y,y,indx);
        end
    end
end

if VERBOSE==1, fprintf('\n'); end
numiter = k;
out.residual = out.residual(1:k,:);
out.time = out.time(1:k,:);
out.rank= out.rank(1:k,:);
out.nuclearNorm= out.nuclearNorm(1:k,:);

author: andrewssobral

repo: lrslibrary


text:
function MedianCutAux(xMin, xMax, yMin, yMax, iter)
%
%
%        MedianCutAux(xMin, xMax, yMin, yMax, iter)
%
%
%     Copyright (C) 2011-14  Francesco Banterle
%
%     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/>.
%

global L;
global imgWork;
global lights;

lx = xMax - xMin;
ly = yMax - yMin;

if((lx > 2) && (ly > 2) && (iter > 0))
    tot = sum(sum(L(yMin:yMax, xMin:xMax)));
    pivot = -1;

    if(lx > ly)
        %cut on the X-axis
        for i=xMin:(xMax - 1)
            c = sum(sum(L(yMin:yMax, xMin:i)));
            if(c >= (tot - c))
                pivot = i;
                break;
            end
        end

        if(pivot == -1)
            pivot = xMax-1;
        end

        MedianCutAux(xMin,    pivot, yMin, yMax, iter-1);
        MedianCutAux(pivot+1, xMax,  yMin, yMax, iter-1);
    else
        %cut on the Y-axis
        for i=yMin:(yMax - 1)
            c = sum(sum(L(yMin:i, xMin:xMax)));
            if(c >= (tot - c))
                pivot = i;
                break;
            end
        end

        if(pivot == -1)
            pivot = yMax-1;
        end

        MedianCutAux(xMin, xMax, yMin,    pivot, iter-1);
        MedianCutAux(xMin, xMax, pivot+1, yMax,  iter-1);
    end
else
    %Generation of the light source
    lights = [lights, CreateLight(xMin, xMax, yMin, yMax, L, imgWork)];
end

end

author: banterle


text:
 function zwhite
%function zwhite
%	in black background mode, the z axis does not show up properly
%	and this fixes it for current plot
set(gca, 'zcolor', [1 1 1])

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [n, nm, nl, ts, names, m] = nex_marker(filename, varname)
% nex_marker(filename, varname): Read a marker variable from a .nex file
%
% [n, nm, nl, ts, names, m] = nex_marker(filename, varname)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   varname - variable name
%
%           continuous (a/d) data come in fragments. Each fragment has a timestamp
%           and a number of a/d data points. The timestamp corresponds to
%           the time of recording of the first a/d value in this fragment.
%           All the data values stored in the vector d.
% OUTPUT:
%   n - number of markers
%   nm - number of fields in each marker
%   nl - number of characters in each marker field
%   ts - array of marker timestamps (in seconds)
%   names - names of marker fields ([nm 64] character array)
%   m - character array of marker values [n nl nm]

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

n = 0;
nm = 0;
nl = 0;
ts = 0;
m = 0;
names = 0;

if(nargin ~= 2)
   disp('2 input arguments are required')
   return
end

if(ischar(filename) == 0)
   disp('input arguments should be character arrays')
   return
end

if(ischar(varname) == 0)
   disp('input arguments should be character arrays')
   return
end

if(length(filename) == 0)
   [fname, pathname] = uigetfile('*.nex', 'Select a Nex file');
    filename = strcat(pathname, fname);
end

fid = fopen(filename, 'r', 'ieee-le');
if(fid == -1)
    disp('cannot open file');
   return
end

disp(strcat('file = ', filename));
magic = fread(fid, 1, 'int32');
version = fread(fid, 1, 'int32');
comment = fread(fid, 256, 'char');
freq = fread(fid, 1, 'double');
tbeg = fread(fid, 1, 'int32');
tend = fread(fid, 1, 'int32');
nvar = fread(fid, 1, 'int32');
fseek(fid, 260, 'cof');
name = zeros(1, 64);
found = 0;
for i=1:nvar
    type = fread(fid, 1, 'int32');
    var_version = fread(fid, 1, 'int32');
    name = fread(fid, [1 64], 'char');
    offset = fread(fid, 1, 'int32');
    n = fread(fid, 1, 'int32');
    dummy = fread(fid, 32, 'char');
    adfreq = fread(fid, 1, 'double');
    adtomv = fread(fid, 1, 'double');
    npw = fread(fid, 1, 'int32');
    nm = fread(fid, 1, 'int32');
    nl = fread(fid, 1, 'int32');
    dummy = fread(fid, 68, 'char');
    name = char(name);
    name = deblank(name);
    k = strcmp(name, deblank(varname));
    if(k == 1)
        if type ~= 6
            disp(sprintf('%s is not a marker variable', deblank(varname)));
            return;
        end
        found = 1;
        fseek(fid, offset, 'bof');
        ts = fread(fid, [1 n], 'int32');
        names = zeros(1,64);
        m = zeros(n, nl, nm);
        for j=1:nm
            names(j, :) = fread(fid, [1 64], 'char');
            for p = 1:n
                m(p, :, j) = fread(fid, [1 nl], 'char');
            end
        end
        break
    end
end

fclose(fid);

if found == 0
    disp('did not find variable in the file');
else
    names = char(names);
    m = char(m);
    ts = ts/freq;
    disp(strcat('number of markers = ', num2str(n)));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [mu sigma2] = estimateGaussian(X)
%ESTIMATEGAUSSIAN This function estimates the parameters of a
%Gaussian distribution using the data in X
%   [mu sigma2] = estimateGaussian(X),
%   The input X is the dataset with each n-dimensional data point in one row
%   The output is an n-dimensional vector mu, the mean of the data set
%   and the variances sigma^2, an n x 1 vector
%

% Useful variables
[m, n] = size(X);

% You should return these values correctly
mu = zeros(n, 1);
sigma2 = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the mean of the data and the variances
%               In particular, mu(i) should contain the mean of
%               the data for the i-th feature and sigma2(i)
%               should contain variance of the i-th feature.
%

mu = mean(X);
sigma2 = var(X,opt=1);

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

end

author: vugsus

repo: coursera-machine-learning


text:
function [inflMap, colXCoord, rowYCoord, colDividerXCoord, rowDividerYCoord, rowLeafPositions,MLCopeningSize] = getLSInfluenceMapFactorNoBar(LS,leak)
%"getLSInfluenceMap"
%   Gets an image of the influence generated by the beam described in LS.
%   Use getDICOMLeafPositions to generate LS.
%
%JRA&KZ 02/8/05
%
%Usage:
%   function inflMap = getLSInfluenceMap(LS);
%
% 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/>.

%Maximum precision of leaf position, in mm. Varian End and Side Accuracy 1.0 mm at
%isocenter. End and Side Repeatability 0.5 mm
precision = .5;

%Get x max, min and round to precision value.

%If X jaws dosn't exist in DICOM
if ~isfield(LS,'xLimits')
    xMax = ceil(max(vertcat(LS.xLeafPositions{:}),[],1) / precision) * precision;
    xMin = floor(min(vertcat(LS.xLeafPositions{:}),[],1) / precision) * precision;
    LS.xLimits{1}(1) = xMin;
    LS.xLimits{1}(2) = xMax;
end

xMax = ceil(max(vertcat(LS.xLimits{:}),[],1) / precision) * precision;
xMin = floor(min(vertcat(LS.xLimits{:}),[],1) / precision) * precision;
fieldSize.x = max(xMax) - min(xMin);
fieldLim.x  = [max(xMax) min(xMin)];

yMax = ceil(max(vertcat(LS.yLimits{:}),[],1) / precision) * precision;
yMin = floor(min(vertcat(LS.yLimits{:}),[],1) / precision) * precision;
fieldSize.y = max(yMax) - min(yMin);
fieldLim.y  = [max(yMax) min(yMin)];

yRes = precision;
nyElements = ceil(fieldSize.y/yRes);
xRes = precision;
nxElements = ceil(fieldSize.x/xRes);

inflMap=zeros(nyElements, nxElements);
colDividerXCoord = linspace(fieldLim.x(2), fieldLim.x(1), nxElements+1);
rowDividerYCoord = linspace(fieldLim.y(2), fieldLim.y(1), nyElements+1);

if isfield(LS, 'yLeafPositions')
    rowLeafPositions = round(interp1(rowDividerYCoord, 1:nyElements+1, LS.yLeafPositions,'linear', 'extrap'));
    rowLeafPositions = clip(rowLeafPositions, 1, nyElements+1, 'limits');
    leafBoundariesToKeep = [diff(rowLeafPositions)>0;true];
    rowLeafPositions = rowLeafPositions(leafBoundariesToKeep);
    leavesToKeep = leafBoundariesToKeep(1:end-1);
else
    LS.xLeafPositions{1} = [xMin xMax];
    LS.meterSetWeight = {1};
    rowLeafPositions = [1 nyElements+1];
    leavesToKeep = 1;
end

if length(LS.meterSetWeight) == 1
    doses = LS.meterSetWeight{:};
else
    doses = [0 diff([LS.meterSetWeight{:}])];
end

% backupMap = inflMap;
%h = waitbar(0,['Generating Fluence Map From MLC Positions For Beam ',num2str(beamIndex)],'Name','Please wait...');

%HCF - head scatter radiation parametrs for varian 2100CD Zhu, MedPhys,
%Vol. 31, No 9 , Sept 2004
if LS.energy == 6 && strcmpi(LS.manufacturer(1),'V')
    A1 = 0.0013;
    A2 = 0.078;
    K = 1.5;
    LAMDA = 7.69;
elseif LS.energy == 18 && strcmpi(LS.manufacturer(1),'V')
    A1 = 0.0013;
    A2 = 0.082;
    K = 1.5;
    LAMDA = 8.16;
elseif LS.energy == 6 && strcmpi(LS.manufacturer(1),'E')
    A1 = 0.0005;
    A2 = 0.072;
    K = 1.2;
    LAMDA = 9.85;
elseif LS.energy == 18 && strcmpi(LS.manufacturer(1),'E')
    A1 = 0.0004;
    A2 = 0.088;
    K = 1.2;
    LAMDA = 8.57;
elseif LS.energy == 6 && strcmpi(LS.manufacturer(1),'S')
    A1 = 0.0004;
    A2 = 0.094;
    K = 1.5;
    LAMDA = 9.8;
elseif LS.energy == 18 && strcmpi(LS.manufacturer(1),'S')
    A1 = 0.0004;
    A2 = 0.099;
    K = 1.5;
    LAMDA = 9.52;
else
    errordlg('This MLC Manufacturer or Energy is not supported or Missed Info in DICOM');
    return
end

for i=1:length(LS.xLeafPositions)
    %    inflMap = backupMap;
    nLeaves = length(LS.xLeafPositions{i})/2;

    if length(LS.xLimits) > 1
        jpL = LS.xLimits{i}(1);
        jpR = LS.xLimits{i}(2);
    else
        jpL = LS.xLimits{1}(1);
        jpR = LS.xLimits{1}(2);
    end

    lpL = LS.xLeafPositions{i}(1:nLeaves);
    lpR = LS.xLeafPositions{i}(nLeaves+1:end);
    lpLK = lpL(leavesToKeep);
    lpRK = lpR(leavesToKeep);

    MLCopeningSize(:,i) = lpRK - lpLK;

    lpLCols = interp1(colDividerXCoord, 1:nxElements+1, lpLK, 'linear', 'extrap');
    lpRCols = interp1(colDividerXCoord, 1:nxElements+1, lpRK, 'linear', 'extrap');

    %Column divider positions of jaws.
    jpLCol = interp1(colDividerXCoord, 1:nxElements+1, jpL, 'linear', 'extrap');
    jpRCol = interp1(colDividerXCoord, 1:nxElements+1, jpR, 'linear', 'extrap');

    jpLCol = round(jpLCol);
    jpRCol = round(jpRCol);

    lpLCols = clip(lpLCols, jpLCol, jpRCol, 'limits');
    lpRCols = clip(lpRCols, jpLCol, jpRCol, 'limits');

    lpLCols = round(lpLCols);
    lpRCols = round(lpRCols);

   for j=1:length(lpLCols)
        %HCF from output ratio for MLC fields Zhu, MedPhys
        F_X = abs(lpLCols(j) - (lpRCols(j)-1))*precision/10;
        F_Y = abs(rowLeafPositions(j+1) - rowLeafPositions(j))*precision/10;
        F = (1+K)*F_X * F_Y/(K*F_X + F_Y);
        HCF = (1+A1*F)*(1+A2*(erf(F/LAMDA))^2)/((1+A1*10)*(1+A2*(erf(10/LAMDA))^2));
        inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpLCols(j):lpRCols(j)-1) = inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpLCols(j):lpRCols(j)-1) + HCF*doses(i);
        inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, jpLCol:lpLCols(j)-1) = inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, jpLCol:lpLCols(j)-1) + leak*doses(i);
        inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpRCols(j):jpRCol-1) = inflMap(rowLeafPositions(j):rowLeafPositions(j+1)-1, lpRCols(j):jpRCol-1) + leak*doses(i);
    end

    %waitbar(i/length(LS.xLeafPositions));
    %    frame = inflMap;
    %    imagesc(inflMap);
    %    mi(:,:,i) = inflMap;
    %    inflMap(inflMap == 0) = 1;
    %    inflMap(inflMap ~= 0) = 2;
    %    colormap([0 0 0; 1 1 1]);
    %    %mi(:,:,i) = inflMap;
    %    %mi(i) = im2frame(inflMap, [0 0 0; 1 1 1]);
    %    drawnow;
    %    pause(.006);
end
%close(h);
colXCoord = colDividerXCoord(1:end-1) + precision/2;
rowYCoord = rowDividerYCoord(1:end-1) + precision/2;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function movie = avi (a)
% AVI Avi Reading Class
% AVI() creates an avi class with no set values
% AVI(avi class name) copies the contents of the passed avi class
% AVI(filename) opens the file with the passed filename

if nargin == 0
   movie.fid 				= 0;
   movie.filename 		= '';
   movie.lines 			= 0;
   movie.columns 			= 0;
   movie.bytes 			= 0;
   movie.color_depth		= 0;
   movie.num_frames 		= 0;
   movie.frame_length	= 0;
   movie.rate				= 0;
   movie.movi_loc			= 0;
   movie.num_input		= 5000;
   movie.idx1				= 0;
   movie.index				= 0;
   movie						= class (movie, 'avi');
elseif isa (a, 'avi')
   movie = a;
elseif isa (a, 'char')
   movie = avi;
   movie = open (movie, a);
else
   error ('#AVI Error: input arguments unacceptable');
end


text:
function [U, L, mu, mse] = pca(X, m)
% Principal component analysis
% Input:
%   X: d x n data matrix
%   m: target dimension
% Output:
%   U: d x m Projection matrix
%   L: m x 1 Eigen values
%   mu: d x 1 mean
%   mse: mean square error
% Written by Mo Chen (sth4nth@gmail.com).
n = size(X,2);
mu = mean(X,2);
Xo = bsxfun(@minus,X,mu);
S = Xo*Xo'/n;                   % 12.3
[U,L] = eig(S);                         % 12.5
[L,idx] = sort(diag(L),'descend');
mse = sum(L)-sum(L(1:m));
U = U(:,idx(1:m));
L = L(1:m);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
classdef C10MOP2 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'nb101';
            config.args.objs = 'err&params&flops';
            config.args.normalized_objectives = true;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [X, spectrum] = slcmds(D, d, w, ty)
%SLMDS Performs Classical Multidimensional scaling
%
% $ Syntax $
%   - X = slcmds(D, d)
%   - X = slcmds(D, d, w)
%   - X = slcmds(D, d, w, 'sqr')
%   - [X, spectrum] = slcmds(...)
%
% $ Arguments $
%   - D:        The pairwise distance matrix (n x n)
%   - d:        The dimension of the embedding space
%   - w:        The weights of samples (1 x n or [])
%   - X:        The embedded samples (d x n)
%
% $ Description $
%   - X = slcmds(D, d) performs classic multidimensional scaling to
%     pursue an embedding space of d-dimension and the vector
%     representation in that space of the objects, such that the
%     distances are optimally preserved.
%
%   - X = slcmds(D, d, w) If w is not empty, it performs classic
%     multidimensional scaling on weighted samples.
%
%   - X = slcmds(D, d, w, 'sqr') indicates that D contains the square
%     of distances.
%
%   - [X, spectrum] = slcmds(...) additionally outputs the spectrum of
%     the embedded space
%
% $ History $
%   - Created by Dahua Lin, on Sep 8th, 2006
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('slcmds', 2);
end

if ndims(D) ~= 2 || size(D, 1) ~= size(D, 2)
    error('sltoolbox:invalidarg', ...
        'The D should be a square matrix');
end
n = size(D, 1);

if d >= n
    error('sltoolbox:exceedbound', ...
        'The dimension d should be less than the number of samples n');
end

if nargin < 3
    w = [];
else
    if ~isempty(w)
        if ~isequal(size(w), [1, n])
            error('sltoolbox:sizmismatch', ...
                'If w is specified, it should be an 1 x n row vector');
        end
    end
end

if nargin >= 4 && strcmpi(ty, 'sqr')
    is_sqr = true;
else
    is_sqr = false;
end

%% compute

if ~is_sqr
    K = sldists2kernels(D);
else
    K = sldists2kernels(D, 'sqr');
end

[X, spectrum] = slkernelembed(K, d, w);

author: lmthang


text:
function dlap_rhs_print ( n, rhs )

%*****************************************************************************80
%
%% DLAP_RHS_PRINT prints the right hand side vector from a DLAP file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Mark Seager,
%    A SLAP for the Masses,
%    Lawrence Livermore National Laboratory,
%    Technical Report UCRL-100267, December 1988.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real RHS(N), the right hand side vector to be written.
%
   r8vec_print ( n, rhs, '  DLAP right hand side vector:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h=holder(tfr,f,n1,n2,t,pl)
%HOLDER	Estimate the Holder exponent through an affine TFR.
%	H=HOLDER(TFR,F,N1,N2,T) estimates the Holder exponent of a
%	function through an affine time-frequency representation of it,
%	and plots the frequency marginal and the regression line.
%
%	TFR : affine time-frequency representation.
%	F   : frequency values of the spectral analysis.
%	N1  : indice of the  minimum frequency for the linear regression.
%                                         (default : 1).
%	N2  : indice of the  maximum frequency for the linear regression.
%                                         (default : length(F)).
%	T : time vector. If T is omitted, the function returns the
%	    global estimate of the Holder exponent. Otherwise, it
%	    returns the local estimates H(T) at the instants specified
%	    in T.
%	H : output value (if T omitted) or vector (otherwise) containing
%	    the Holder estimate(s).
%
%	Example :
%	 S=altes(128); [TFR,T,F]=tfrscalo(S,1:128,8);
%	 H=holder(TFR,F,1,length(F));

%	P. Goncalves, October 1995
%	Copyright (c) 1995 Rice University
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

clf;
tfr = abs(tfr) ;

if nargin<2,
 error('There must be at least 2 input parameters');
elseif nargin==2,
 n1=1; n2=length(f);
elseif nargin==3,
 n2=length(f);
end
if nargin<=5,
 pl=1;
else
 pl=0;
end

if nargin <= 4
  fmarg = mean(tfr.').' ;
  if pl, plot(log(f),log(fmarg)) , hold on; end;
  p = polyfit(log(f(n1:n2)),log(fmarg(n1:n2)),1) ;
  drt = p(1)*log(f)+p(2) ;
  if pl,
   plot(log(f),drt,'--g') ;
   legend('Frequency marginal','Regression line');
   plot(log(f([n1 n2])),log(fmarg([n1 n2])),'+r') ;
   xticklabels = round(logspace(log10(min(f)),log10(max(f)),4)*1000)./1000 ;
   set(gca,'XTick',log(xticklabels)) ;
   xlabel('frequency (logarithmically spaced)') ;
   hold off ;
  end
  h = (-p(1)-1)/2 ;
else
  if length(t) == 1
    if pl, plot(log(f),log(tfr(:,t))) ; hold on; end
    p = polyfit(log(f(n1:n2)),log(tfr(n1:n2,t)),1) ;
    drt = p(1)*log(f)+p(2) ;
    if pl,
     plot(log(f),drt,'--g') ;
     legend('Frequency marginal','Regression line');
     plot(log(f([n1 n2])),log(tfr([n1 n2],t)),'+r') ;
     xticklabels = round(logspace(log10(min(f)),log10(max(f)),4)*1000)./1000 ;
     set(gca,'XTick',log(xticklabels)) ;
     xlabel('frequency (logarithmically spaced)') ;
     hold off ;
    end
    h = (-p(1)-1)/2 ;
  elseif length(t) > 1
    [yt,xt] = size(t) ; if yt>xt, t = t.' ; end ;
    j = 1 ;
    for k = t,
      p = polyfit(log(f(n1:n2)),log(tfr(n1:n2,k)),1) ;
      h(j) = (-p(1)-1)/2 ; j = j + 1 ;
    end
    if pl,
     plot(t,h) ;grid
     title('Holder estimates at time instants T');
    end
  end,
end

author: HeLiangHIT


text:
% configure parameters for wheel-slam

paras.datapath       = "data\wheel-ins_DR.bin";
paras.gt_path        = "data\ground_truth.bin";
paras.initheading    = 27.58;  %only for comparison with ground truth

paras.sigmaPhi       = (0.05 * pi/180);
paras.sigmaDis_scale = 0.05;

paras.odomdata_scale    = 0.5;%only process half of the data
paras.NPARTICLES        = 100;
paras.NEFFECTIVE        = paras.NPARTICLES * 0.5;
paras.rollSeqdimension  = 100;
paras.gridlen           = 1.5;
paras.conRevisitNumThr  = 5;
paras.corrCoefThr       = 0.4;
paras.corrCoefNumThr    = 2;

repo: Wheel-SLAM


text:
function v = drosInsituPositives(ids, drosinsitu, drosexp),

% DROSINSITUPOSITIVES Filter genes by looking for positive in-situ annotations.
% FORMAT
% DESC Filter a set of genes by retaining ones with positive in-situ
% annotations, keeping the order of the original list.
% ARG ids : Genes to look, expressed as indices in drosexp
% ARG drosinsitu : drosinsitu data set returned by drosLoadData
% ARG drosexp : drosexp data set returned by drosLoadData
% RETURN v : the filtered indices
%
% SEEALSO : drosLoadData
%
% COPYRIGHT : Antti Honkela, 2009

% SHEFFIELDML

isgenes = drosinsitu.genes(find(sum(drosinsitu.data, 2)));
I = drosFindGeneinds(drosexp, isgenes, 0);

[C, IA, IB] = intersect(ids, I);

v = ids(sort(IA));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [x,p] = getLastState(stage,colloc,stageVars)
N = stage.N;
nx = stage.nx;
ni = colloc.num_i;
nu = stage.nu;
np = stage.np;

[X_indizes, ~, ~, P_indizes, ~] = ocl.simultaneous.indizes(N, nx, ni, nu, np);
x = stageVars(X_indizes(:,end));
p = stageVars(P_indizes(:,end));

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function [nnz,ind, val] = pennlp_congrad(i,x)

try
    f = datasaver(5,x,i+1);
    [ind,dummy,val] = find(f(:));
    nnz = length(ind);
catch

end

%[nnz,ind, val] = dg(i,x);

author: yalmip

repo: YALMIP


text:
function value = r8vec_min_pos ( n, a )

%*****************************************************************************80
%
%% R8VEC_MIN_POS returns the minimum positive value of an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, real A(N), the array.
%
%    Output, real VALUE, the smallest positive entry.
%
  value = realmax ( );

  for i = 1 : n
    if ( 0.0 < a(i) )
      value = min ( value, a(i) );
    end
  end

  return
end


text:
function plotExcitationMagnitude(varargin)
% Plots the excitation force magnitude for each hydro structure's bodies in
% the heave, surge and pitch degrees of freedom.
%
% Usage:
% ``plotExcitationMagnitude(hydro, hydro2, hydro3, ...)``
%
% Parameters
% ----------
%     varargin : struct(s)
%         The hydroData structure(s) created by the other BEMIO functions.
%         One or more may be input.
%

if isempty(varargin)
    error(['plotExcitationMagnitude: No arguments passed. Include one or more hydro ' ...
        'structures when calling: plotExcitationMagnitude(hydro1, hydro2, ...)']);
end

B=1;  % Wave heading index
figHandle = figure('Position',[950,500,975,521]);
titleString = ['Normalized Excitation Force Magnitude: ',...
    '$$\bar{X_i}(\omega,\theta) = {\frac{X_i(\omega,\theta)}{{\rho}g}}$$'];
subtitleString = {'Surge','Heave','Pitch'};
xString = {'$$\omega (rad/s)$$','$$\omega (rad/s)$$','$$\omega (rad/s)$$'};
yString = {['$$\bar{X_1}(\omega,\theta$$',' = ',num2str(varargin{1}.theta(B)),'$$^{\circ}$$)'],...
    ['$$\bar{X_3}(\omega,\theta$$',' = ',num2str(varargin{1}.theta(B)),'$$^{\circ}$$)'],...
    ['$$\bar{X_5}(\omega,\theta$$',' = ',num2str(varargin{1}.theta(B)),'$$^{\circ}$$)']};

notes = {''};

numHydro = length(varargin);
for ii = 1:numHydro
    numBod = varargin{ii}.Nb;
    tmp1 = strcat('X',num2str(ii));
    X.(tmp1) = varargin{ii}.w;
    tmp2 = strcat('Y',num2str(ii));
    a = 0;
    for i = 1:numBod
        m = varargin{ii}.dof(i);
        Y.(tmp2)(1,i,:) = squeeze(varargin{ii}.ex_ma(a+1,B,:));
        Y.(tmp2)(2,i,:) = squeeze(varargin{ii}.ex_ma(a+3,B,:));
        Y.(tmp2)(3,i,:) = squeeze(varargin{ii}.ex_ma(a+5,B,:));
        legendStrings{i,ii} = [varargin{ii}.body{i}];
        a = a + m;
    end
end

formatPlot(figHandle,titleString,subtitleString,xString,yString,X,Y,legendStrings,notes);
saveas(figHandle,'Excitation_Magnitude.png');

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function g = linearLogLikeGradients(model)

% LINEARLOGLIKEGRADIENTS Linear model gradients.
% FORMAT
% DESC computes the gradients of the log likelihood of a
% linear model with respect to the parameters.
% ARG model : the model structure for computing the log likelihood.
% RETURN g : the gradients of the model log likelihood.
%
% SEEALSO : modelLogLikeihood, lineargrad
%
% COPYRIGHT : Neil D. Lawrence, 2006

% MLTOOLS

Xo = [model.X ones(size(model.X, 1), 1)];
W = [model.W; model.b];
G = -model.beta*(Xo'*Xo*W - Xo'*model.y);
gW = G(1:end-1, :);
gb = G(end, :);
g = [gW(:)' gb];


text:
function title = p20_title ( option )

%*****************************************************************************80
%
%% P20_TITLE sets the title for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, string TITLE, the title of the problem.
%    TITLE will never be longer than 80 characters.
%
  if ( option == 1 )
    title = 'The Buckling Spring, F(L,Theta,Lambda,Mu).';
  else
    title = '???';
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P20_TITLE - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized option number.\n' );
    error ( 'P20_TITLE - Fatal error!' );
  end

  return
end


text:
% FindInl    find inliers in joint image matrix
%			 by pairwise epipolar geometry
%
% function IdMatIn = findinl(Ws,IdMat,tol)
% Ws ... 3MxN joint image matrix
% IdMat ... MxN ... 0 -> no point detected
%                   1 -> point detected
% tol ... [pixels] tolerance for the epipolar geometry
%         the point are accpted as outliers only if they
%         are closer to the epipolar line than tol

% $Author: svoboda $
% $Revision: 2.0 $
% $Id: findinl.m,v 2.0 2003/06/19 12:07:09 svoboda Exp $
% $State: Exp $

function IdMatIn = findinl(Ws,IdMat,tol)

NoCams = size(IdMat,1);

% fill the array of structures not_used denoted as 0
% allocate the array of structures for used
for i=1:NoCams,
  not_used(i).pts = sum(IdMat(i,:));
  used(i).pts	  = -1;
end

% allocate IdMat for outliers
IdMatIn = zeros(size(IdMat));

while (sum([not_used.pts])>1-NoCams),
  [buff, id.cam_max]  = max([not_used.pts]);
  used	   = add(used, id.cam_max, not_used(id.cam_max).pts);
  not_used = remove(not_used, id.cam_max);
  Mask	   = repmat(IdMat(id.cam_max,:),NoCams,1);
  Corresp  = Mask & IdMat;
  Corresp(id.cam_max,:) = 0;
  [buff, id.cam_to_pair] = max(sum(Corresp')); % find the camera with most correspondences
  idx.corr_to_pair = find(sum(IdMat([id.cam_max,id.cam_to_pair],:))==2);
  % used	   = add(used, id.cam_to_pair, not_used(id.cam_to_pair).pts);
  % not_used = remove(not_used, id.cam_to_pair);
  if size(idx.corr_to_pair,2)<8,
	error('Not enough points to compute epipolar geometry in RANSAC validation')
  end
  Wspair   = [];
  Wspair   = Ws(id.cam_max*3-2:id.cam_max*3, idx.corr_to_pair);
  Wspair   = [Wspair; Ws(id.cam_to_pair*3-2:id.cam_to_pair*3, idx.corr_to_pair)];
  % id
  [F, inls] = rEG(Wspair,tol,tol,0.99);
  IdMatIn(id.cam_max, idx.corr_to_pair(inls)) = 1;
  IdMat(id.cam_max, :)						  = 0;
  IdMat(id.cam_max, idx.corr_to_pair(inls))	  = 1;
end

function list = add(list, id, value)
list(id).pts = value;
return

function list = remove(list, id)
list(id).pts = -1;
return

author: strawlab

repo: MultiCamSelfCal


text:
%% Work out the eye scaling first using the PDM
addpath('../PDM_helpers/');
% also work out a size for each window
% Replace this with the location of in 300 faces in the wild data
if(exist([getenv('USERPROFILE') '/Dropbox/AAM/test data/'], 'file'))
    root_test_data = [getenv('USERPROFILE') '/Dropbox/AAM/test data/'];
else
    root_test_data = 'F:/Dropbox/Dropbox/AAM/test data/';
end

[~, detections, labels] = Collect_wild_imgs(root_test_data, true, true, true, true);

% load('../models/pdm/pdm_68_multi_pie');
load('../models/pdm/pdm_68_aligned_wild');

%%
shapes = zeros(size(labels));
num_points = 68;

for i=1:size(labels,1)

    [ a_orig, R, ~, ~, ~, ~] = fit_PDM_ortho_proj_to_2D(M, E, V, squeeze(labels(i,:,:)));
    view_actual = Rot2Euler(R);

    views = [0,0,0; 0,-30,0; -30,0,0; 0,30,0; 30,0,0];
    views = views * pi/180;

    [~,view_id] = min(sum(abs(bsxfun(@plus, views, -view_actual)), 2));

    rot = Euler2Rot(views(view_id,:));

    rot_m = rot * reshape(M, num_points, 3)';
    width_model = max(rot_m(1,:)) - min(rot_m(1,:));
    height_model = max(rot_m(2,:)) - min(rot_m(2,:));

    bounding_box = detections(i,:);

    a = (((bounding_box(3) - bounding_box(1)) / width_model) + ((bounding_box(4) - bounding_box(2))/ height_model)) / 2;

    tx = (bounding_box(3) + bounding_box(1))/2;
    ty = (bounding_box(4) + bounding_box(2))/2;

    % correct it so that the bounding box is just around the minimum
    % and maximum point in the initialised face
    tx = tx - a*(min(rot_m(1,:)) + max(rot_m(1,:)))/2;
    ty = ty - a*(min(rot_m(2,:)) + max(rot_m(2,:)))/2;

    % visualisation of the initial state
    %hold off;imshow(Image);hold on;plot(a*rot_m(1,:)+tx, a*rot_m(2,:)+ty,'.r');hold on;rectangle('Position', [bounding_box(1), bounding_box(2), bounding_box(3)-bounding_box(1), bounding_box(4)-bounding_box(2)]);
    global_params = [a, 0, 0, 0, tx, ty]';
    global_params(2:4) = views(view_id);

    local_params = zeros(numel(E), 1);

    shape = GetShapeOrtho(M, V, local_params, global_params);
    shape = shape(:,1:2);
    shapes(i,:,:) = shape;

    shapes(i,:,:) = shapes(i,:,:) / a_orig;
    labels(i,:,:) = labels(i,:,:) / a_orig;

end

%%
errs_img = sort(squeeze(mean(mean(abs(labels - shapes), 2), 3)));

errs_lmk = squeeze(mean(mean(abs(labels - shapes), 1), 3));

errs_img_95 = errs_img(round(0.95*end));

scales = [0.25, 0.35, 0.5];

windows = floor(errs_img_95 * 2 * scales + 11) + 1

author: TadasBaltrusaitis

repo: OpenFace


text:
function x = p03_start ( i, x )

%*****************************************************************************80
%
%% P03_START returns a starting point for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the starting point.
%
%    Output, real X, the starting point.
%
  if ( i == 1 )
    x = -1.0;
  elseif ( i == 2 )
    x =  0.5;
  elseif ( i == 3 )
    x =  2.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P03_START - Fatal error!\n' );
    fprintf ( 1, '  Illegal starting point index = %d\n', i );
    error ( 'P03_START - Fatal error!' );
  end

  return
end


text:
function variables( varargin )

%VARIABLES Declares one or more CVX variables.
%   VARIABLES x1 x2 x3 ..., where x1, x2, x3, etc. are valid
%   variable names, declares multiple cvx variables. It is
%   exactly equivalent to issuing a separate VARIABLE command
%   for each x1, x2, x3, ...
%
%   VARIABLES allows the declaration of vector, matrix, and
%   array variables. However, unlike VARIABLE, structure modifiers
%   such as "symmetric", "toeplitz", etc. are NOT permitted. Thus
%   VARIABLES cannot be used with variables with structure.
%
%   Examples:
%      variables x y z;
%      variables x(100) y z(100,10);
%
%   See also VARIABLE, DUAL, DUALS.

if ~iscellstr( varargin ),
    error( 'VARIABLES must be used in command mode.' );
end

for k = 1 : nargin,
    evalin( 'caller', [ 'variable ', varargin{k} ] );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [] = plot_way(ax, parsed_osm, map_img_filename)
%PLOT_WAY   plot parsed OpenStreetMap file
%
% usage
%   PLOT_WAY(ax, parsed_osm)
%
% input
%   ax = axes object handle
%   parsed_osm = parsed OpenStreetMap (.osm) XML file,
%                as returned by function parse_openstreetmap
%   map_img_filename = map image filename to load and plot under the
%                      transportation network
%                    = string (optional)
%
% 2010.11.06 (c) Ioannis Filippidis, jfilippidis@gmail.com
%
% See also PARSE_OPENSTREETMAP, EXTRACT_CONNECTIVITY.

% ToDo
%   add double way roads

if nargin < 3
    map_img_filename = [];
end

[bounds, node, way, ~] = assign_from_parsed(parsed_osm);

disp_info(bounds, size(node.id, 2), size(way.id, 2))
show_ways(ax, bounds, node, way, map_img_filename);

function [] = show_ways(hax, bounds, node, way, map_img_filename)
show_map(hax, bounds, map_img_filename)

%plot(node.xy(1,:), node.xy(2,:), '.')

key_catalog = {};
for i=1:size(way.id, 2)
    [key, val] = get_way_tag_key(way.tag{1,i} );

    % find unique way types
    if isempty(key)
        %
    elseif isempty( find(ismember(key_catalog, key) == 1, 1) )
        key_catalog(1, end+1) = {key};
    end

    % way = highway or amenity ?
    flag = 0;
    switch key
        case 'highway'
            flag = 1;

            % bus stop ?
            if strcmp(val, 'bus_stop')
                disp('Bus stop found')
            end
        case 'amenity'
            % bus station ?
            if strcmp(val, 'bus_station')
                disp('Bus station found')
            end
        otherwise
            disp('way without tag.')
    end

    % plot highway
    way_nd_ids = way.nd{1, i};
    num_nd = size(way_nd_ids, 2);
    nd_coor = zeros(2, num_nd);
    nd_ids = node.id;
    for j=1:num_nd
        cur_nd_id = way_nd_ids(1, j);
        if ~isempty(node.xy(:, cur_nd_id == nd_ids))
             nd_coor(:, j) = node.xy(:, cur_nd_id == nd_ids);
        end
    end

    % remove zeros
    nd_coor(any(nd_coor==0,2),:)=[];

    if ~isempty(nd_coor)
        % plot way (highway = blue, other = green)
        if flag == 1
            plot(hax, nd_coor(1,:), nd_coor(2,:), 'b-')
        else
            plot(hax, nd_coor(1,:), nd_coor(2,:), 'g--')
        end
    end

    %waitforbuttonpress
end
disp(key_catalog.')

function [] = disp_info(bounds, Nnode, Nway)
disp( ['Bounds: xmin = ' num2str(bounds(1,1)),...
    ', xmax = ', num2str(bounds(1,2)),...
    ', ymin = ', num2str(bounds(2,1)),...
    ', ymax = ', num2str(bounds(2,2)) ] )
disp( ['Number of nodes: ' num2str(Nnode)] )
disp( ['Number of ways: ' num2str(Nway)] )

author: johnyf

repo: openstreetmap


text:
% The COBRAToolbox: testWriteSBML.m
%
% Purpose:
%     - read in a model from a .xml file, write it as a .sbml file,
%       read it back in and assess any differences
%
% Authors:
%     - Partial original file: Joseph Kang 04/07/09
%     - CI integration: Laurent Heirendt

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testWriteSBML'));
cd(fileDir);

testModelXML = readCbModel('Ec_iJR904.xml');

% write the model as a .sbml file
writeCbModel(testModelXML, 'sbml', 'testModelSBML.sbml');

% read the model from the newly written .sbml file
testModelSBML = readCbModel('testModelSBML.sbml');

%Due to the original model being SBML 2.1, no proteins/geneNames were
%generated.
%These models are not the same, since the References given by the xml (in
%the notes field), are not PubMed IDs, so they are transferred to the
%rxnNotes field,
testModelSBML = rmfield(testModelSBML,'rxnNotes');
testModelXML = rmfield(testModelXML,'rxnNotes');
testModelXML = rmfield(testModelXML,'rxnReferences');

% check whether both models are the same
 [isSame numDiff fieldNames] = isSameCobraModel(testModelXML, testModelSBML);

% assess any potential differences
assert(~any(numDiff))

% obtain SBML structure
testSBMLstruct = writeSBML(testModelSBML,'testModelSBML.sbml');

% asses if lower/upper bounds have proper SBO terms assigned
for i = 1:length(testSBMLstruct.parameter)
    value = testSBMLstruct.parameter(i).value;
    if abs(value) == 1000 || isinf(value) || value == 0
        assert(testSBMLstruct.parameter(i).sboTerm == 626)
    else
        assert(testSBMLstruct.parameter(i).sboTerm == 625)
    end
end

% remove the written file to clean up
delete 'testModelSBML.sbml.xml';

%read/writeSBML should be the only methods used for cobra sbml conversion.
% change the directory
cd(currentDir)

% try to shut down any left open parpool
try
    poolobj = gcp('nocreate');
    delete(poolobj);
catch
    % pass
end

author: opencobra

repo: cobratoolbox


text:
function cae = caebp(cae, y)

    %%  backprop deltas
    cae.L = 0;
    for i = 1 : numel(cae.o)
        %  error
        cae.e{i} = (cae.o{i} - y{i}) .* cae.edgemask;
        %  loss function
        cae.L = cae.L + 1/2 * sum(cae.e{i}(:) .^2 ) / size(cae.e{i}, 1);
        %  output delta
        cae.od{i} = cae.e{i} .* (cae.o{i} .* (1 - cae.o{i}));

        cae.dc{i} = sum(cae.od{i}(:)) / size(cae.e{i}, 1);
    end

    for j = 1 : numel(cae.a)   %  calc activation deltas
        z = 0;
        for i = 1 : numel(cae.o)
             z = z + convn(cae.od{i}, flipall(cae.ok{i}{j}), 'full');
        end
        cae.ad{j} = cae.a{j} .* (1 - cae.a{j}) .* z;
    end

    %%  calc gradients
    ns = size(cae.e{1}, 1);
    for j = 1 : numel(cae.a)
        cae.db{j} = sum(cae.ad{j}(:)) / ns;
        for i = 1 : numel(cae.o)
            cae.dok{i}{j} = convn(flipall(cae.a{j}), cae.od{i}, 'valid') / ns;
            cae.dik{i}{j} = convn(cae.ad{j}, flipall(cae.i{i}), 'valid') / ns;
        end
    end

end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function obj = setDimension(obj, dim)
    % Sets the dimension of the NLP variable vector
    %
    % Parameters:
    % dim: the dimension of the vector @type double

    assert(isscalar(dim) && dim >=0, 'The dimension must be a scalar positive value.');

    obj.Dimension = dim;

end

author: ayonga

repo: frost-dev

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


text:
function F = in_fread_brainamp(sFile, sfid, SamplesBounds)
% IN_FREAD_BRAINAMP:  Read a block of recordings from BrainVision BrainAmp .eeg file
%
% USAGE:  F = in_fread_brainamp(sFile, sfid, SamplesBounds=[])

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2013

% Parse inputs
if (nargin < 3) || isempty(SamplesBounds)
    SamplesBounds = round(sFile.prop.times .* sFile.prop.sfreq);
end

% BINARY files
if strcmpi(sFile.header.DataFormat, 'BINARY')
    nChan = sFile.header.NumberOfChannels;
    nSamplesToRead = SamplesBounds(2) - SamplesBounds(1) + 1;
    % MULTIPLEXED files
    if strcmpi(sFile.header.DataOrientation, 'MULTIPLEXED')
        % Get start and length of block to read
        offsetData = SamplesBounds(1) * nChan * sFile.header.bytesize;
        % Position file at the beginning of the data block
        fseek(sfid, offsetData, 'bof');
        % Read all values at once
        F = fread(sfid, [nChan, nSamplesToRead], sFile.header.byteformat);
    % VECTORIZED
    elseif strcmpi(sFile.header.DataOrientation, 'VECTORIZED')
        % Get blocks of samples to skip at the beginning and end of each channel
        offsetStart = SamplesBounds(1) * sFile.header.bytesize;
        offsetEnd = (round(sFile.prop.times(2) .* sFile.prop.sfreq) - SamplesBounds(2)) * sFile.header.bytesize;
        % Blocks of samples to skip between two blocks to read
        offsetSkip = offsetStart + offsetEnd;
        % Position file at the beginning of the trial
        fseek(sfid, offsetStart, 'bof');
        % Read the requested samples for all the channels
        % => WARNING: CALL TO FREAD WITH SKIP=0 DOES NOT WORK PROPERLY
        if (offsetSkip == 0)
            F = fread(sfid, [nSamplesToRead, nChan], sFile.header.byteformat)';
        else
            precision = sprintf('%d*%s', nSamplesToRead, sFile.header.byteformat);
            F = fread(sfid, [nSamplesToRead, nChan], precision, offsetSkip)';
        end
    end
% ASCII and VECTORIZED files
elseif (strcmpi(sFile.header.DataFormat, 'ASCII') && strcmpi(sFile.header.DataOrientation, 'VECTORIZED'))
    % Open file
    fid = fopen(sFile.filename, 'r');
    % Initialize data matrix
    F = zeros(sFile.header.NumberOfChannels, sFile.header.DataPoints);
    iChannel = 1;
    % Read the entire file line by line
    while(1)
        % Display message
        % disp(sprintf('BRAINAMP> Reading channel #%d...', iChannel));
        % Reached the end of the file: exit the loop
        if feof(fid)
            break;
        end
        % Read one line
        strChan = strtrim(fgetl(fid));
        if isempty(strChan)
            continue;
        end
        % Find the first separator
        iSep = min([find(strChan == ' ',1), find(strChan == sprintf('\t'),1)]);
        % Replace "," with "." for the numbers
        strChan(strChan == ',') = '.';
        % Read the values
        F(iChannel,:) = sscanf(strChan(iSep+1:end), '%f')';
        iChannel = iChannel+1;
    end
    % Close file
    fclose(fid);
    % Select only the requested time points
    iTime = (SamplesBounds(1):SamplesBounds(2)) - round(sFile.prop.times(1) .* sFile.prop.sfreq) + 1;
    F = F(:,iTime);
end

% Apply gains, if available
if isfield(sFile.header, 'chgain') && (length(sFile.header.chgain) == size(F,1))
    F = bst_bsxfun(@times, F, sFile.header.chgain(:));
% Else: Convert from microVolts to Volts by default
else
    F = F * 1e-6;
end

author: brainstorm-tools