text:
%% Analyzing Investment Strategies with CVaR Portfolio Optimization in MATLAB - Reality
%
% Robert Taylor
% The MathWorks, Inc.

% Copyright (C) 2012 The MathWorks, Inc.

%% Introduction

% This script illustates the simulation of a single realization for both uncovered and covered-call
% positions on a common underlying stock price realization obtained from the file
% |BuyWriteTestData.mat|.

%% Generate a Single Realization for Uncovered and Covered-Call Scenarios

% Example:
%	Assume data are sampled at 30-minute intervals over 22 days. Assume that a "year" is 252 days,
%	that a "day" is from 09:30 to 14:00, and that a "period" is 30 minutes. Consequently, a "day"
%	comprises 13 "periods" with 1 + 22*13 samples and the investment period is 22/252 years. Note
%	that the extra sample is the initial price.
%
%	X(1)				40			initial price of stock
%	numel(X)			1 + 22*13	initial price plus 13 30-minute periods in a day for 22 days
%	T					22/252		investment period has 22 days with 252 days in a year
%	volatility			0.35		35% volatility
%
%	initial_equity		1000000		$1 million initial equity invested in stock (25000 shares)
%	distribution		0.10		10% per year distribution
%
%	option_expiration	45/252		option expiration in 45 days
%	risk_free_rate		0.0015		0.15% risk-free rate
%	strike_cushion		0.05		5% strike price "cushion"
%	exercise_likelihood 0.2			probability of exercise if stock price >= strike price

%	stock_cost			0.06		5 c/share txcost + 1 c/share spread
%	option_cost			0.20		20 c/option contract

%	confirmation_delay	2			1 hour delay to confirm option assignment
%	reinvestment_delay	1*13		1 day to reinvest after confirmation of assignment
%	settlement_delay	3*13		3 days to settle after reinvestment

% fund details

distribution = 0.10;

% initial number of shares invested in stock

initial_holdings = 25000;			% initial number of shares

% investment period
%	periodicity is 30-minute periods, 13 periods in a day, 252 days in a year

T = 22/252;							% duration of investment period is 22 days
N = 22*13;							% 22 days x 13 30-minute periods in a day

% stock information

load BuyWriteTestData X				% realization of asset total return prices X

initial_price = X(1);				% initial price
mu = 0.10;							% stock drift
volatility = 0.35;					% stock volatility

initial_equity = initial_holdings*initial_price;

% option information

contract_expiration = 45/252;		% option contract expiration (45 days)
next_contract_expiration = 112/252;	% next option contract expiration (112 days)
risk_free_rate = 0.0015;			% risk-free rate
strike_cushion = 0.05;				% strike price "cushion"
% exercise_likelihood = NaN;		% period probability of exercise if stock price >= strike price
exercise_likelihood = 0.9;			% period probability of exercise if stock price >= strike price

% costs/frictions

stock_cost = 0.06;					% 5 c/share txcost + 1 c/share spread
contract_cost = 0.20;				% 20 c/option contract + 2 c/option contract spread

% delays

confirmation_delay = 2;				% 1 hour delay to confirm option assignment
reinvestment_delay = 1*13;			% 1 day to reinvest after confirmation of assignment
settlement_delay = 3*13;			% 3 days to settle after reinvestment

no_reinvestment = false;			% true if reinvestment not allowed in an investment period
% no_reinvestment = true;			% true if reinvestment not allowed in an investment period
covered_vs_uncovered = false;

% simulate scenarios

[rU, WU, CU, HU] = uncovered_engine(X, T, ...
	initial_equity, distribution, risk_free_rate, stock_cost);

[rC, WC, CC, HC, KC] = covered_engine(X, T, mu, volatility, ...
	initial_equity, distribution, no_reinvestment, covered_vs_uncovered, ...
	strike_cushion, contract_expiration, next_contract_expiration, risk_free_rate, ...
	stock_cost, contract_cost, exercise_likelihood, ...
	confirmation_delay, reinvestment_delay, settlement_delay);

fprintf('Scenario total returns ...\n');
fprintf('  Uncovered  %10.4f  Total return for uncovered strategy\n',rU);
fprintf('  Covered    %10.4f  Total return for covered strategy\n',rC);

fprintf('Mean of total returns ...\n');
fprintf('  Uncovered  %10.4f\n',N*mean(tick2ret(WU)));
fprintf('  Covered    %10.4f\n',N*mean(tick2ret(WC)));

fprintf('Standard deviation of total returns ...\n');
fprintf('  Uncovered  %10.4f\n',sqrt(N)*std(tick2ret(WU)));
fprintf('  Covered    %10.4f\n',sqrt(N)*std(tick2ret(WC)));

fprintf('Maximum drawdown of total returns ...\n');
fprintf('  Uncovered  %10.4f\n',maxdrawdown(WU));
fprintf('  Covered    %10.4f\n',maxdrawdown(WC));

%% Plot Results

t = linspace(0, T*252, N+1);

subplot(2,2,1);
plot(t, [WC/WC(1), WU/WU(1)]);
title('\bfTotal Wealth');
xlabel('Day');
ylabel('Dollars (Millions)');
h = legend('Covered', 'Uncovered', 'Location', 'SouthEast');
set(h, 'FontSize', 7, 'Box', 'off');

subplot(2,2,2);
plot(t, 1.0e-6*CC);
title('\bfCash');
xlabel('Day');
ylabel('Dollars (Millions)');

subplot(2,2,3);
plot(t, 1.0e-3*HC);
title('\bfShares');
xlabel('Day');
ylabel('Shares (Thousands)');

subplot(2,2,4);
plot(t, KC);
title('\bfStrike Price');
xlabel('Day');
ylabel('Dollars');


text:
function value = hammersley_base_check ( dim_num, base )

%*****************************************************************************80
%
%% HAMMERSLEY_BASE_CHECK checks BASE for a Hammersley sequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer BASE(1:DIM_NUM), the bases.
%    Bases should not be 0 or 1.
%
%    Output, logical VALUE, is true if BASE is legal.
%
  if ( any ( base(1:dim_num) == 0 ) | any ( base(1:dim_num) == 1 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HAMMERSLEY_BASE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  BASE(I) = 0 or 1 for some I!\n' );
    value = 0;
  else
    value = 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Eugenio Alcala Baselga
% Date: 02/06/2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ u_opt, u_predicted, x_predicted, x_err] = MPC_Kinematic_Computation...
    ( controller, Hp, KC, oldu, u_predicted, x_predicted, x_Real,...
    kin_states, kin_inputs, x_Ref, y_Ref, theta_Ref, vel_Ref, omega_Ref,...
    longitud_vect, car, frozen_based, references_based)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\\ Vehicle kinematic control:
%   controller: object provided by the MPC design
%   Hp: Prediction horizon
%   KC: discrete kinematic time
%   oldu: previous control action
%   u_predicted: vector of predicted control actions in the last iteration
%   x_predicted: vector of predicted states in the last iteration
%   x_Real: vector of measured vehicle states [x y theta v alpha w]
%   kin_states: number of kinematic states: 3
%   kin_inputs: number of kinematic control actions: 2
%   x_Ref: x reference (planner)
%   y_Ref: y reference (planner)
%   theta_Ref: theta reference (planner)
%   vel_Ref: velocity reference (planner)
%   omega_Ref: angular velocity reference (planner)
%   longitud_vect: length of vectors for simulation
%   car: vehicle object
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if KC==1
    % Initializing vectors:
    theta_cos           = ones(1,Hp+1);
    theta_sin           = zeros(1,Hp+1);
    theta_err_cos       = ones(1,Hp);
    U_OPT_VECTOR        = zeros(kin_inputs,Hp,longitud_vect);
    X_OPT_VECTOR        = zeros(kin_states,Hp,longitud_vect);
end

%%%% Constraints %%%%
% Control action constraints:
    % CA
        u_min_ini   = vehicle.u_min_ini;
        u_max_ini   = vehicle.u_max_ini;
        u_min       = vehicle.u_min;
        u_max       = vehicle.u_max;
    % derivative of CA
        du_min_ini  = vehicle.du_min_ini;
        du_max_ini  = vehicle.du_max_ini;
        du_min      = vehicle.du_min;
        du_max      = vehicle.du_max;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Terminal Set LMIs
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %[ PP, WW ] = TerminalSet_Continuous_Knemtic_Cmputtion();
    %[ PP, WW ] = TerminalSet_Discrete_Knemtic_Cmputtion(Ts_Kcontrol);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Errors de verdad (no forecasted from the LPV model):
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x_err(1)    = (x_Ref(1)-x_Real(1))*cos(x_Real(3)) + (y_Ref(1)-x_Real(2))*sin(x_Real(3));
    x_err(2)    = -(x_Ref(1)-x_Real(1))*sin(x_Real(3)) + (y_Ref(1)-x_Real(2))*cos(x_Real(3));
    x_err(3)    = theta_Ref(1) - x_Real(3);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Moving Horizon Constraints:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if KC>1
        for j=1:Hp+1
            theta_cos(j) = cos(theta_Ref(1));
            theta_sin(j) = sin(theta_Ref(1));
            if(j>1)
                theta_err_cos(j-1) = cos(x_predicted(3,1)); % cos (theta_ERROR forecasted)
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Computing PREDICTIONS of KIN-DYN LPV Vehicle Model. Note that
        % this A and B matrices are in continuous time due to the
        % discretization is carried out inside the controller object.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [Ac_pred, Bc_pred, Bref] = KinemError_LPV_Model_predicted(KC,Hp,x_predicted(3,:),...
            omega_Ref, vel_Ref );
        if frozen_based == true
            for i=1:Hp
                Ac_pred(:,:,i) = Ac_pred(:,:,1);
                Bc_pred(:,:,i) = Bc_pred(:,:,1);
            end
        elseif references_based == true
            for i=1:Hp
                Ac_pred(:,:,i) = Ac_pred(:,:,i);
                Bc_pred(:,:,i) = Bc_pred(:,:,i);
            end
        end
    else
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Computing KIN-DYN LPV Vehicle Model for First iteration. Note that
        % this A and B matrices are in continuous time due to the
        % discretization is carried out inside the controller object.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [Ac_pred, Bc_pred, Bref] = KinemError_LPV_Model(KC,Hp,x_err(3),...
            omega_Ref(1), vel_Ref(1) );
    end

    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Optimization stage:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%     if(KC > car.t_ini)
        inputs = {x_err', oldu, vel_Ref, omega_Ref, Ac_pred, Bc_pred,...
            theta_cos, theta_sin, theta_err_cos, u_min, u_max, du_min, du_max};
%     else
%         inputs = {x_err', oldu, vel_Ref, omega_Ref, Ac_pred, Bc_pred,...
%             theta_cos, theta_sin, theta_err_cos, u_min_ini, u_max_ini,...
%             du_min_ini, du_max_ini};
%     end

    [solutions,diagnostics] = controller{inputs};

    if diagnostics == 1
        u_opt = oldu;
        U_OPT_VECTOR(:,:,KC) = U_OPT_VECTOR(:,:,KC-1);
        X_OPT_VECTOR(:,:,KC) = X_OPT_VECTOR(:,:,KC-1);
    else
        U_OPT_VECTOR(:,:,KC) = solutions{1};
        u_opt = U_OPT_VECTOR(:,1,KC);
        X_OPT_VECTOR(:,:,KC) = solutions{2}; % Kinematic errors
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Predictions obtained from the solver:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x_predicted = X_OPT_VECTOR(:,:,KC);
    u_predicted = U_OPT_VECTOR(:,:,KC);

end


text:
function sunDotNormal = computeSunDotNormal(hra)
%     rVectBodyToSun = -1.0 * getPositOfBodyWRTSun(ut, bodyInfo, celBodyData);
%
%     if(norm(rVectBodyToSun) == 0)
%         sunDotNormal = 1.0;
%     else
%         rVectSunECEF = getFixedFrameVectFromInertialVect(ut, rVectBodyToSun, bodyInfo);
%         rECEF = getrVectEcefFromLatLongAlt(zeros(size(long)), long, zeros(size(long)), bodyInfo);
%
%         planarRvectEcef = [rECEF(1,:); rECEF(2,:); zeros(1,size(rECEF,2))];
%         planarRvectSunEcef = [rVectSunECEF(1,:); rVectSunECEF(2,:); zeros(1,size(rVectSunECEF,2))];
%
%         hra = angleNegPiToPi(dang(planarRvectEcef,planarRvectSunEcef));
%
%         sunDotNormal = 0.5 * cos(hra + deg2rad(45)) + 0.5;
%     end
    sunDotNormal = 0.5 * cos(hra + deg2rad(45)) + 0.5;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function b = center_logistic(u,y);
s = mean(y);
b = log( s / (1 - s ) ) - mean(u);
for i=1:20
	sigma = 1./ ( 1 + exp( - u - b ) );
	phider = sum(sigma - y) ;
	phider2 = sum(sigma .* ( 1- sigma));
	bnew = b - phider / phider2;
	if norm(bnew-b) < 1e-16, break; end
	b=bnew;

end

repo: REST


text:
classdef DepthAnalysisWindow < AnalysisWindow
    % DEPTHANALYSISWINDOW shows the number of events, binned per depth
    properties
        binedges    double; % km
    end
    methods
        function obj=DepthAnalysisWindow(ax,binedges)
            obj@AnalysisWindow(ax);
            obj.binedges=binedges;
        end

        function prepare_axes(obj)
            % prepare the depth profile axes
            obj.ax.Tag          = 'dvEventsWidthDepth';
            obj.ax.YDir         = 'reverse';
            obj.ax.Title.String = 'Depth Profile';
            obj.ax.XLabel.String= 'Number of events';
            obj.ax.YLabel.String= 'Depth';
        end

        function [x,y]=calculate(obj,catalog)
            % calculate the number of events (y) for each binned depth(x)
            x = histcounts(catalog.Depth,obj.binedges);
            y=obj.binedges(1:end-1) + diff(obj.binedges)./2;
        end

    end
end

author: CelsoReyes


text:
%compute head speed in the direction perpendicular to the central-head direction
function [data,units]=compute_velheadperch(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velheadperch=cell(1,numlarvae);
centralheadangperp=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    centralheadangperp{1,i}=trx(larva).centralheadang-pi/2;
    velheadperch{1,i}=trx(larva).velmaghead.*(cos(trx(larva).velanghead).*cos(trx(larva).centralheadang(1,1:end-1)+pi/2)+sin(trx(larva).velanghead).*sin(trx(larva).centralheadang(1,1:end-1)+pi/2));
end

units=parseunits('mm/s');
data=velheadperch;

author: kristinbranson

repo: JAABA


text:
function PopDec = EvolEI(Problem,Dec,Obj,model,IFEs)
% Solution update in EGO, where a solution with the best expected
% improvement is re-evaluated

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

    Off   = [OperatorGA(Problem,Dec(TournamentSelection(2,size(Dec,1),Obj),:));OperatorGA(Problem,Dec,{0,0,1,20})];
    N     = size(Off,1);
    EI    = zeros(N,1);
    Gbest = min(Obj);
    E0    = inf;

    while IFEs > 0
        drawnow('limitrate');
        for i = 1 : N
            [y,~,mse] = predictor(Off(i,:),model);
            s         = sqrt(mse);
            EI(i)     = -(Gbest-y)*normcdf((Gbest-y)/s)-s*normpdf((Gbest-y)/s);
        end
        [~,index] = sort(EI);
        if EI(index(1)) < E0
            Best = Off(index(1),:);
            E0   = EI(index(1));
        end
        Parent = Off(index(1:ceil(N/2)),:);
        Off    = [OperatorGA(Problem,Parent(TournamentSelection(2,size(Parent,1),EI(index(1:ceil(N/2)))),:));OperatorGA(Problem,Parent,{0,0,1,20})];
        IFEs   = IFEs - size(Off,1);
    end
    PopDec = Best;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [Mean2, Var2] = online_UpdateWin(Mean, Var, old_x, new_x, n)
%Using this function to update mean and variance
%
%Input
%	Mean	: average of previous data
%	Var     : variance of previous data
%	old_x	: oldest data in sliding window
%	new_x	: new arriving data
%	n		: number of previous data
%
%Output
%	Mean2	: result of average by online updating average with new data
%	Var2	: result of variance by online updating average with new data
%
    dim = size(Var);
	if n <= 0
		Mean = x;
		Var = ones(dim)*0.5;
	else
		Mean2 = ( n * Mean - old_x + new_x ) / n;
		Var2 = max((( n * ( Var + Mean.^2 ) - old_x.^2 + new_x.^2 ) / n - Mean2.^2),0);
	end
end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
filename='Bridge_hexahedra_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef Predict

  methods (Static) % fastPredict

    function fastPredict = fastPredict(nClassifier)
      % fastPredict constructor
      %
      % fastPredict(ibeh).classifier % struct array, reduced/indexed classifier
      % fastPredict(ibeh).windowfeaturescellparams % struct, pffnames -> cell
      % fastPredict(ibeh).wfs % window features
      % fastPredict(ibeh).pffs % all distinct pffs for wfs; fieldnames of windowfeaturescellparams;
      % fastPredict(ibeh).ts % classifier TS for this classifier
      % fastPredict(ibeh).wfidx % indices into 'external'/supplied featurename list corresponding to .wfs
      % fastPredict(ibeh).wfidx_valid = false;

      fastPredict = struct(...
        'classifier',cell(nClassifier,1),...
        'windowfeaturescellparams',[],...
        'wfs',{cell(1,0)},...
        'pffs',{cell(1,0)}, ...
        'ts',[], ...
        'wfidx',[], ...
        'wfidx_valid',false);
    end

    function fastPredict = fastPredictInit(fastPredict,classifiers,classifierTS,feature_names) %#ok<INUSL>
      % classifiers: cell array of classifiers
      % classifierTS: classifier timestamps
      % feature_names: cellstr

        nClassifier = numel(classifiers);
        assert(isequal(nClassifier,numel(classifierTS),numel(feature_names)));

        fastPredict = Predict.fastPredict(nClassifier); % No part of fastPredict is reused

        for iCls = 1:nClassifier

          if isempty(classifiers{iCls})
            continue;
          end

          % which features are actually used
          dims = [classifiers{iCls}.dim];
          feature_names_cls = feature_names{iCls}(dims);

          % put these in with the rest of the classifiers' window features
          wfs = {};
          for j = 1:numel(feature_names_cls),
            wfidxcurr = find(WindowFeatureNameCompare(feature_names_cls{j},wfs),1);
            if isempty(wfidxcurr),
              wfidxcurr = numel(wfs)+1;
              wfs{wfidxcurr} = feature_names_cls{j}; %#ok<AGROW>
            end
          end

          wf2pff = cellfun(@(x)x{1},wfs,'UniformOutput',false);
          [pffs,~,wf2pffidx] = unique(wf2pff);

          windowfeaturescellparams = struct;
          for pfi = 1:numel(pffs),
            pf = pffs{pfi};
            wfidx_cur = wf2pffidx==pfi;
            windowfeaturescellparams.(pf) = WindowFeatureName2Params(wfs(wfidx_cur));
          end

          classifiers_indexed = classifiers{iCls};
          for j = 1:numel(classifiers_indexed),
            classifiers_indexed(j).dim = j;
          end

          fastPredict(iCls).classifier = classifiers_indexed;
          fastPredict(iCls).windowfeaturescellparams = windowfeaturescellparams;
          fastPredict(iCls).wfs = feature_names_cls;
          fastPredict(iCls).pffs = pffs;
          fastPredict(iCls).ts = classifierTS(iCls);
          fastPredict(iCls).wfidx = [];
          fastPredict(iCls).wfidx_valid = false;
        end
    end

    function fastPredict = fastPredictFindWfidx(fastPredict,feature_names)
      feature_names_char = cellfun(@cell2str,feature_names,'UniformOutput',false);
      for i = 1:numel(fastPredict)
        wfs = fastPredict(i).wfs;

        wfs_char = cellfun(@cell2str,wfs,'UniformOutput',false);
        [ism,wfidx] = ismember(wfs_char,feature_names_char);
        if ~all(ism),
          error('Error matching wfs for classifier with window features computed');
        end

  %       wfidx = nan(1,numel(wfs));
  %       for j = 1:numel(wfs),
  %         idxcurr = find(WindowFeatureNameCompare(wfs{j},feature_names));
  %         if numel(idxcurr) ~= 1,
  %           error('Error matching wfs for classifier with window features computed');
  %         end
  %         wfidx(j) = idxcurr;
  %       end

        fastPredict(i).wfidx = wfidx;
        fastPredict(i).wfidx_valid = true;
      end
    end

  end

  methods (Static) % predictblocks

    function pb = predictblocks(nCls)
      pb = struct('t0',cell(nCls,1),'t1',[],'expi',[],'flies',[]);
    end

  end

end


text:
function varargout = process_diff_mean( varargin )
% PROCESS_DIFF_MEAN: Difference of means of sets A and B.

% @=============================================================================
% 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, 2010-2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Difference of means';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Test';
    sProcess.Index       = 100;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Difference';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 1;

    % === Absolue values: legend
    sProcess.options.labelavg.Type = 'label';
    sProcess.options.labelavg.Comment = 'Function to estimate the average across the files:';
    sProcess.options.labelavg.InputTypes = {'data', 'results', 'matrix'};
    % === Absolue values: type
    sProcess.options.avg_func.Comment = {'Arithmetic average
mean(A) - mean(B)', ...
                                         'Absolute value of average
abs(mean(A)) - abs(mean(B))', ...
                                         'Average of absolute values
mean(abs(A)) - mean(abs(B))', ...
                                         'Power
mean(A^2) - mean(B^2)'};
    sProcess.options.avg_func.Type    = 'radio';
    sProcess.options.avg_func.Value   = 2;
    sProcess.options.avg_func.InputTypes = {'data', 'results', 'matrix'};
    % === WEIGHTED AVERAGE
    sProcess.options.weighted.Comment    = 'Weighted average:  mean(x) = sum(Leff_i * x(i)) / sum(Leff_i)';
    sProcess.options.weighted.Type       = 'checkbox';
    sProcess.options.weighted.Value      = 0;
    sProcess.options.weightedlabel.Comment    = '      Leff_i = Effective number of averages for file #i';
    sProcess.options.weightedlabel.Type       = 'label';
    % === MATCH ROWS WITH NAMES
    sProcess.options.matchrows.Comment    = 'Match signals between files using their names';
    sProcess.options.matchrows.Type       = 'checkbox';
    sProcess.options.matchrows.Value      = 1;
    sProcess.options.matchrows.InputTypes = {'timefreq', 'matrix'};
    % === EXCLUDE ZEROS FROM THE AVERAGE
    sProcess.options.iszerobad.Comment    = 'Exclude the flat signals from the average (zero at all times)';
    sProcess.options.iszerobad.Type       = 'checkbox';
    sProcess.options.iszerobad.Value      = 1;
    sProcess.options.iszerobad.InputTypes = {'timefreq', 'matrix'};

    % === UNCONSTRAINED SOURCES
    sProcess.options.label_norm.Comment    = ['Note: For unconstrained sources, "absolute value" refers to the norm
' ...
                                              'of the three orientations: abs(F) = sqrt(Fx^2 + Fy^2 + Fz^2)'];
    sProcess.options.label_norm.Type       = 'label';
    sProcess.options.label_norm.InputTypes = {'results'};
end

%% ===== FORMAT COMMENT =====
function [Comment, strAbs] = FormatComment(sProcess)
    % Weighted
    if isfield(sProcess.options, 'weighted') && isfield(sProcess.options.weighted, 'Value') && ~isempty(sProcess.options.weighted.Value) && sProcess.options.weighted.Value
        Comment = 'Difference of weighted means';
    else
        Comment = 'Difference of means';
    end
    % If sources: averaging option
    if isfield(sProcess.options, 'avg_func')
        switch(sProcess.options.avg_func.Value)
            case 1,  strAbs = '[mean]';
            case 2,  strAbs = '[abs(mean)]';
            case 3,  strAbs = '[mean(abs)]';
            case 4,  Comment = 'Difference of power';    strAbs = '';
        end
    else
        strAbs = '[mean]';
    end
    if ~isempty(strAbs)
        Comment = [Comment, ' ', strAbs];
    end
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputsA, sInputsB) %#ok<DEFNU>
    OutputFiles = {};
    % Make sure that file type is indentical for both sets
    if ~isempty(sInputsA) && ~isempty(sInputsB) && ~strcmpi(sInputsA(1).FileType, sInputsB(1).FileType)
        bst_report('Error', sProcess, sInputsA, 'Cannot process inputs from different types.');
        return;
    end
    % === GET OPTIONS ===
    isResults = strcmpi(sInputsA(1).FileType, 'results');
    if isfield(sProcess.options, 'avg_func')
        switch (sProcess.options.avg_func.Value)
            case 1,  Function = 'mean'; isAbsDiff = 0; strComment = '';
            case 2,  Function = 'mean'; isAbsDiff = 1; strComment = ' [abs(avg)]';
            case 3,  Function = 'norm'; isAbsDiff = 0; strComment = ' [avg(abs)]';
            case 4,  Function = 'rms';  isAbsDiff = 0; strComment = ' [power]';
        end
    else
        strComment = '';
        Function = 'mean';
        isAbsDiff = 0;
    end
    % Weighted
    if isfield(sProcess.options, 'weighted') && isfield(sProcess.options.weighted, 'Value') && ~isempty(sProcess.options.weighted.Value)
        isWeighted = sProcess.options.weighted.Value;
    else
        isWeighted = 0;
    end
    % Match signals between files using their names
    if isfield(sProcess.options, 'matchrows') && isfield(sProcess.options.matchrows, 'Value') && ~isempty(sProcess.options.matchrows.Value)
        isMatchRows = sProcess.options.matchrows.Value;
    else
        isMatchRows = 1;
    end
    % Exclude zero values
    if isfield(sProcess.options, 'iszerobad') && isfield(sProcess.options.iszerobad, 'Value') && ~isempty(sProcess.options.iszerobad.Value)
        isZeroBad = sProcess.options.iszerobad.Value;
    else
        isZeroBad = 1;
    end
    % Read first file of the list
    [sMat, matName] = in_bst(sInputsA(1).FileName);

    % === COMPUTE DIFFERENCE OF AVG ===
    % Compute average of the two sets of files
    isVariance = 0;
    [StatA, MessagesA, iAvgFileA] = bst_avg_files({sInputsA.FileName}, [], Function, isVariance, isWeighted, isMatchRows, isZeroBad);
    [StatB, MessagesB, iAvgFileB] = bst_avg_files({sInputsB.FileName}, [], Function, isVariance, isWeighted, isMatchRows, isZeroBad);
    % Add messages to report
    if ~isempty(MessagesA)
        if isempty(StatA)
            bst_report('Error', sProcess, sInputsA, MessagesA);
            return;
        else
            bst_report('Warning', sProcess, sInputsA, MessagesA);
        end
    end
    if ~isempty(MessagesB)
        if isempty(StatB)
            bst_report('Error', sProcess, sInputsB, MessagesB);
            return;
        else
            bst_report('Warning', sProcess, sInputsB, MessagesB);
        end
    end
    % Absolute values before difference
    if isAbsDiff
        % Unconstrained sources: Convert to flat map (norm of the three orientations)
        if isResults && isfield(sMat, 'nComponents') && (sMat.nComponents ~= 1)
            % Average A
            sMatTmp = sMat;
            sMat.(matName) = StatA.mean;
            sMat = process_source_flat('Compute', sMat, 'rms');
            StatA.mean = sMat.(matName);
            % Average B
            sMat = sMatTmp;
            sMat.(matName) = StatB.mean;
            sMat = process_source_flat('Compute', sMat, 'rms');
            StatB.mean = sMat.(matName);
            % Edit the comment field to indicate we used a norm
            strComment = strrep(strComment, 'abs', 'norm');
        end
        % Enforce absolute values
        StatA.mean = abs(StatA.mean);
        StatB.mean = abs(StatB.mean);
    % Unconstrained or mixed sources: source maps have been flattened by bst_avg_files
    elseif isResults && ismember(Function, {'norm', 'rms'}) && (sMat.nComponents ~= 1)
        sMat = process_source_flat('Compute', sMat, 'rms');
    end
    % Difference of power: square the rms that is returned by the averaging function
    if strcmpi(Function, 'rms')
        StatA.mean = StatA.mean .^ 2;
        StatB.mean = StatB.mean .^ 2;
    end
    % Check timefreq measures
    if ~isempty(StatA.Measure) && ~strcmpi(StatA.Measure, StatB.Measure)
        bst_report('Error', sProcess, [sInputsA(:)',sInputsB(:)'], 'The two sets of files have different measures applied to the time-frequency coefficients.');
        return
    end
    % Compute difference of averages
    StatA.mean = StatA.mean - StatB.mean;

    % === CREATE OUTPUT STRUCTURE ===
    bst_progress('start', 'Difference of means', 'Saving result...');
    [processComment, strAbs] = FormatComment(sProcess);
    % Get output study
    [sStudy, iStudy] = bst_process('GetOutputStudy', sProcess, [sInputsA, sInputsB]);
    % Comment: forced in the options
    if isfield(sProcess.options, 'Comment') && isfield(sProcess.options.Comment, 'Value') && ~isempty(sProcess.options.Comment.Value)
        Comment = sProcess.options.Comment.Value;
    % Comment: process default
    else
        % Get comment for files A and B
        [tmp__, tmp__, CommentA] = bst_process('GetOutputStudy', sProcess, sInputsA, [], 0);
        [tmp__, tmp__, CommentB] = bst_process('GetOutputStudy', sProcess, sInputsB, [], 0);
        if strcmpi(CommentA, CommentB) && ~strcmpi(sInputsA(1).Condition, sInputsB(1).Condition)
            CommentA = sInputsA(1).Condition;
            CommentB = sInputsB(1).Condition;
            % Add file number
            if (length(sInputsA) > 1) || (length(sInputsB) > 1)
                CommentA = [CommentA ' (' num2str(length(sInputsA)) ')'];
                CommentB = [CommentB ' (' num2str(length(sInputsB))  ')'];
            end
        end
        % Comment: difference A-B
        Comment = [CommentA ' - ' CommentB];
        if ~isempty(strAbs)
           Comment = [Comment ' ' strAbs];
        end
    end
    % Copy fields from StatA structure
    sMat.(matName)   = StatA.mean;
    sMat.Comment     = Comment;
    sMat.ChannelFlag = StatA.ChannelFlag;
    sMat.Time        = StatA.Time;
    sMat.nAvg        = StatA.nAvg + StatB.nAvg;
    % Effective number of averages
    % Leff = 1 / sum_i(w_i^2 / Leff_i),  with w1=1 and w2=-1
    %      = 1 / (1/Leff_A + 1/Leff_B))
    sMat.Leff = 1 / (1/StatA.Leff + 1/StatB.Leff);

    % Colormap for recordings: keep the original
    % Colormap for sources, timefreq... : difference (stat2)
    if ~strcmpi(sInputsA(1).FileType, 'data')
        sMat.ColormapType = 'stat2';
    end
    % Time-frequency: Change the measure type
    % if strcmpi(sInputsA(1).FileType, 'timefreq')
    %     sMat.Measure = 'other';
    % end

    % History: Average
    if isfield(sMat, 'History')
        prevHistory = sMat.History;
        sMat = bst_history('reset', sMat);
        sMat = bst_history('add', sMat, 'diff_mean', processComment);
        sMat = bst_history('add', sMat, 'diff_mean', 'History of the first file:');
        sMat = bst_history('add', sMat, prevHistory, ' - ');
    else
        sMat = bst_history('add', sMat, 'diff_mean', Comment);
    end
    % History: Number of files averaged for each channel
    if any(StatA.nGoodSamples(1) ~= StatA.nGoodSamples) && isfield(StatA, 'RowNames') && ~isempty(StatA.RowNames) && iscell(StatA.RowNames) && (length(StatA.RefRowNames) <= 1)
        strInfo = 'Number of files that were averaged, for each signal (FilesA):';
        sMat = bst_history('add', sMat, 'average', strInfo);
        uniqueGood = unique(StatA.nGoodSamples);
        for i = length(uniqueGood):-1:1
            strHistoryGood = [sprintf(' - %d files: ', uniqueGood(i)), sprintf('%s ', StatA.RowNames{StatA.nGoodSamples == uniqueGood(i)})];
            sMat = bst_history('add', sMat, 'average', strHistoryGood);
            strInfo = [strInfo, 10, strHistoryGood];
        end
        % Generate process warning
        bst_report('Warning', sProcess, sInputsA, strInfo);
    end
    % History: List files A
    sMat = bst_history('add', sMat, 'diff_mean', 'List of files in group A:');
    for i = 1:length(iAvgFileA)
        sMat = bst_history('add', sMat, 'diff_mean', [' - ' sInputsA(iAvgFileA(i)).FileName]);
    end
    % History: List files B
    sMat = bst_history('add', sMat, 'diff_mean', 'List of files in group B:');
    for i = 1:length(iAvgFileB)
        sMat = bst_history('add', sMat, 'diff_mean', [' - ' sInputsB(iAvgFileB(i)).FileName]);
    end

    % Averaging results from the different data file: reset the "DataFile" field
    if isfield(sMat, 'DataFile')
        sMat.DataFile = [];
    end
    % Do not keep the events
    if isfield(sMat, 'Events') && ~isempty(sMat.Events)
        sMat.Events = [];
    end
    % Fix surface link for warped brains
    if isfield(sMat, 'SurfaceFile') && ~isempty(sMat.SurfaceFile) && ~isempty(strfind(sMat.SurfaceFile, '_warped'))
        sMat = process_average('FixWarpedSurfaceFile', sMat, sInputsA(1), sStudy);
    end
    % Rownames
    if isfield(StatA, 'RowNames') && ~isempty(StatA.RowNames)
        if strcmpi(matName, 'TF')
            sMat.RowNames = StatA.RowNames;
        elseif strcmpi(matName, 'Value')
            sMat.Description = StatA.RowNames;
        end
    end

    % === SAVE FILE ===
    % Output filename
    fileTag = bst_process('GetFileTag', sInputsA(1).FileName);
    OutputFiles{1} = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), [fileTag, '_diff_mean']);
    % Save on disk
    bst_save(OutputFiles{1}, sMat, 'v6');
    % Register in database
    db_add_data(iStudy, OutputFiles{1}, sMat);
end

author: brainstorm-tools


text:
function F = constraint(X,quantifier,Y)
% Internal class for constraint list

superiorto('sdpvar');
superiorto('double');
superiorto('logical');

% Evaluate
switch quantifier
    case '>'
        Z = X - Y;
    case '>='
        Z = X - Y;
    case '<'
        Z = Y - X;
    case '<='
        Z = Y - X;
    case '=='
        Z = Y - X;
    otherwise
        error('Quantifier not supported')
end

if isa(Z,'double')

    if  size(Z,1)==size(Z,2) &&  norm(Z-Z',inf)<1e-12 && ~isequal(quantifier,'==')
        checkSDP = 1;
    else
        checkSDP = 0;
    end

    if checkSDP
        if min(eig(Z))>=0
            warning('Inequality constraint evaluated to trivial true.')
            F = [];
            return
        else
            error('Inequality constraint evaluated to trivial false (no decision variable in constraint)')
        end
    else
        Z = Z(:);
        switch quantifier
            case '=='
                if all(Z)==0
                    warning('Equality constraint evaluated to trivial true.')
                    F = [];
                    return
                else
                    error('Equality constraint evaluated to trivial false (no decision variable in constraint)')
                end
            case {'<=','>='}
                if all(Z>=0)
                    warning('Inequality constraint evaluated to trivial true.')
                    F = [];
                    return
                else
                    error('Inequality constraint evaluated to trivial false (no decision variable in constraint)')
                end
            case {'<','>'}
                if all(Z>0)
                    warning('Inequality constraint evaluated to trivial true.')
                    F = [];
                    return
                else
                    error('Inequality constraint evaluated to trivial false (no decision variable in constraint)')
                end
        end
    end
end

switch quantifier
    case {'>','<'}
        F.strict(1) = 1;
    case {'>=','<=','=='}
        F.strict(1) = 0;
    otherwise
        error('Quantifier not supported')
end
if isa(Z,'sdpvar') && ~isequal(quantifier,'==')
    if issquare(Z)
        if ~ishermitian(Z)
            warning('YALMIP:SuspectNonSymmetry','Suspect non-symmetry in square constaint (Learn more) ')
        end
    end
end
F.List={X,quantifier,Y};
F.Evaluated{1} = Z;
F.ConstraintID = yalmip('ConstraintID');
F.tag{1} = '';
F = class(F,'constraint');

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function candoit = canGetHessian(problem)
% Checks whether the Hessian can be computed for a problem structure.
%
% function candoit = canGetHessian(problem)
%
% Returns true if the Hessian of the cost function can be computed given
% the problem description, false otherwise.
%
% See also: canGetCost canGetDirectionalDerivative canGetGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    candoit = isfield(problem, 'hess') || ...
              (isfield(problem, 'ehess') && canGetEuclideanGradient(problem));

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function z = times( x, y, oper )

%   Disciplined convex programming information for TIMES:
%      In general, disciplined convex programs must obey the "no-product
%      rule" which states that two non-constant expressions cannot be
%      multiplied together. Under this rule, the following combinations
%      are allowed:
%         {constant} .* {non-constant}  {non-constant} .* {constant}
%      Furthermore, if the non-constant expression is nonlinear, then
%      the constant term must be real.
%
%      A lone exception to the no-product rule is made for quadratic
%      forms: two affine expressions may be multiplied together if the
%      result is convex or concave. For example, the construction
%         variable x(n)
%         x.*x  <= 1;
%      would be permitted because each element of x.*x is convex.
%
%   Disciplined geometric programming information for TIMES:
%      Both terms in a multiplication must have the same log-curvature,
%      so the following products are permitted:
%         {log-convex} .* {log-convex}  {log-concave} .* {log-concave}
%         {log-affine} .* {log-affine}
%      Note that log-affine expressions are both log-convex and
%      log-concave.
%
%   For vectors, matrices, and arrays, these rules are verified
%   indepdently for each element.

error( nargchk( 2, 3, nargin ) );
if nargin < 3, oper = '.*'; end

%
% Check sizes
%

sx = size( x );
sy = size( y );
xs = all( sx == 1 );
ys = all( sy == 1 );
if xs,
    sz = sy;
elseif ys,
    sz = sx;
elseif ~isequal( sx, sy ),
    error( 'Matrix dimensions must agree.' );
else
    sz = sx;
end
nn = prod( sz );
if nn == 0,
    z = zeros( sz );
    return
end

%
% Determine the computation methods
%

persistent remap_m remap_l remap_r
if isempty( remap_r ),
    % zero .* valid, valid .* zero, constant .* constant
    temp_1    = cvx_remap( 'constant' );
    temp_2    = cvx_remap( 'zero' );
    temp_3    = temp_2' * cvx_remap( 'valid' );
    remap_1   = ( temp_1' * temp_1 ) | temp_3 | temp_3';
    remap_1n  = ~remap_1;
    % constant / nonzero, zero / log-concave, zero / log-convex
    temp_4    = temp_1 & ~temp_2;
    remap_1r  = ( temp_1' * temp_4 ) | ( temp_2' * cvx_remap( 'log-valid' ) );
    remap_1rn = ~remap_1r;

    % constant * affine, real * convex/concave/log-convex, positive * log-concave
    temp_5   = cvx_remap( 'real' );
    temp_6   = cvx_remap( 'positive' );
    temp_7   = cvx_remap( 'affine' );
    temp_1n   = ~temp_1;
    temp_8   = cvx_remap( 'log-concave' ) & temp_1n;
    remap_2  = ( temp_1' * temp_7 ) | ...
               ( temp_5' * cvx_remap( 'convex', 'concave', 'log-convex' ) ) | ...
               ( temp_6' * temp_8 );
    remap_2  = remap_2 & remap_1n;
    % real / log-concave, positive / log-convex
    temp_9   = cvx_remap( 'log-convex' ) & temp_1n;
    remap_2r = ( temp_5' * temp_8 ) | ...
               ( temp_6' * temp_9 );
    remap_2r = remap_2r & remap_1rn;

    % Affine * constant, convex/concave/log-convex * real, log-concave * positive
    remap_3  = remap_2';
    % Affine / nonzero, convex/concave/log-convex / nzreal, log-concave / positive
    remap_3r = remap_3;
    remap_3r(:,2) = 0;

    % Affine * affine
    remap_4  = temp_7 & ~temp_1;
    remap_4  = remap_4' * +remap_4;

    % log-concave * log-concave, log-convex * log-convex
    remap_5  = ( temp_8' * +temp_8 ) | ( temp_9' * +temp_9 );
    % log-concave / log-convex, log-convex / log-concave
    remap_5r = temp_8' * +temp_9;
    remap_5r = remap_5r | remap_5r';

    remap_m = remap_1 + 2 * remap_2  + 3 * remap_3 + 4 * remap_4 + 5 * remap_5;
    remap_r = remap_1r + 2 * remap_2r + 3 * remap_3r + 5 * remap_5r;
    remap_r(:,2) = 0;
    remap_l = remap_r';
end
switch oper,
    case '.*',
        remap = remap_m;
        r_recip = 0;
        l_recip = 0;
    case './',
        remap = remap_r;
        r_recip = 1;
        l_recip = 0;
    case '.\',
        remap = remap_l;
        r_recip = 0;
        l_recip = 1;
end
vx = cvx_classify( x );
vy = cvx_classify( y );
vr = remap( vx + size( remap, 1 ) * ( vy - 1 ) );
vu = sort( vr(:) );
vu = vu([true;diff(vu)~=0]);
nv = length( vu );
if vu(1) == 1 && nv > 1,
    vr(vr==1) == vu(2);
    nv = nv - 1;
    vu(1) = [];
end

%
% Process each computation type separately
%

x   = cvx( x );
y   = cvx( y );
xt  = x;
yt  = y;
if nv ~= 1,
    z = cvx( sz, [] );
end
for k = 1 : nv,

    %
    % Select the category of expression to compute
    %

    if nv ~= 1,
        t = vr == vu( k );
        if ~xs,
            xt = cvx_subsref( x, t );
            sz = size( xt );
        end
        if ~ys,
            yt = cvx_subsref( y, t );
            sz = size( yt );
        end
    end

    %
    % Apply the appropriate computation
    %

    switch vu( k ),
    case 0,

        % Invalid
        error( 'Disciplined convex programming error:\n    Cannot perform the operation: {%s} %s {%s}', cvx_class( xt, true, true, true ), oper, cvx_class( yt, true, true, true ) );

    case 1,

        % constant .* constant
        xb = cvx_constant( xt );
        yb = cvx_constant( yt );
        if l_recip,
            cvx_optval = xb .\ yb;
        elseif r_recip,
            cvx_optval = xb ./ yb;
        else
            cvx_optval = xb .* yb;
        end
        if nnz( isnan( cvx_optval ) ),
            error( 'Disciplined convex programming error:\n    This expression produced one or more invalid numeric values (NaNs).', 1 ); %#ok
        end
        cvx_optval = cvx( cvx_optval );

    case 2,

        % constant .* something
        xb = cvx_constant( xt );
        if l_recip, xb = 1.0 ./ xb; end
        yb = yt;
        if r_recip && nnz( xb ), yb = exp( - log( yb ) ); end
        yb = yb.basis_;
        if ~xs,
            nn = numel(  xb  );
            if ys,
                xb = cvx_reshape( xb, [ 1, nn ] );
                if issparse( yb ) && ~issparse( xb ),
                    xb = sparse( xb );
                end
            else
                n1 = 1 : nn;
                xb = sparse( n1, n1, xb( : ), nn, nn );
            end
        end
        cvx_optval = cvx( sz, yb * xb );

    case 3,

        % something .* constant
        yb = cvx_constant( yt );
        if r_recip, yb = 1.0 ./ yb; end
        xb = xt;
        if l_recip && any( xb ), xb = exp( - log( xb ) ); end
        xb = xb.basis_;
        if ~ys,
            nn = numel(  yb  );
            if xs,
                yb = cvx_reshape( yb, [ 1, nn ] );
                if issparse( xb ) && ~issparse( yb ),
                    yb = sparse( yb );
                end
            else
                n1 = 1 : nn;
                yb = sparse( n1, n1, yb( : ), nn, nn );
            end
        end
        cvx_optval = cvx( sz, xb * yb );

    case 4,

        % affine .* affine
        nn = prod( sz );
        xA = xt.basis_; yA = yt.basis_;
        if xs && ~ys, xA = xA( :, ones( 1, nn ) ); end
        if ys && ~xs, yA = yA( :, ones( 1, nn ) ); end
        mm = max( size( xA, 1 ), size( yA, 1 ) );
        if size( xA, 1 ) < mm, xA( mm, end ) = 0; end
        if size( yA, 1 ) < mm, yA( mm, end ) = 0; end
        xB = xA( 1, : ); xA( 1, : ) = 0;
        yB = yA( 1, : ); yA( 1, : ) = 0;
        cyA   = conj( yA );
        alpha = sum( real( xA .* yA ), 1 ) ./ max( sum( cyA .* yA, 1 ), realmin );
        adiag = sparse( 1 : nn, 1 : nn, alpha, nn, nn );
        if all( sum( abs( xA - cyA * adiag ), 2 ) <= 2 * eps * sum( abs( xA ), 2 ) ),
            beta  = xB - alpha .* conj( yB );
            alpha = reshape( alpha, sz );
            if isreal( y ),
                cvx_optval = alpha .* square( y ) + reshape( beta, sz ) .* y;
            elseif all( abs( beta ) <= 2 * eps * abs( xB ) ),
                cvx_optval = alpha .* square_abs( y );
            else
                error( 'Disciplined convex programming error:\n    Invalid quadratic form(s): product is not real.\n', 1 ); %#ok
            end
        else
            error( 'Disciplined convex programming error:\n    Invalid quadratic form(s): not a square.\n', 1 ); %#ok
        end

    case 5,

        % posynomial .* posynomial
        xb = log( xt );
        if l_recip, xb = - xb; end
        yb = log( yt );
        if r_recip, yb = - yb; end
        cvx_optval = exp( xb + yb );

    otherwise,

        error( 'Shouldn''t be here.' );

    end

    %
    % Store the results
    %

    if nv == 1,
        z = cvx_optval;
    else
        z = cvx_subsasgn( z, t, cvx_optval );
    end

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.


text:
% VL_LBP  Local Binary Patterns
%   F = VL_LBP(IM, CELLSIZE) computes the Local Binary Pattern (LBP)
%   features for image I.
%
%   IM is divided in cells of size CELLSIZE. F is a three-dimensional
%   array containing one histograms of quantized LBP features per
%   cell. The witdh of F is FLOOR(WIDTH/CELLSIZE), where WIDTH is the
%   width of the image. The same for the height. The third dimension
%   is 58.
%
%   See also: LBP, VL_LBPFLIPLR(),
%   VL_HELP().

% Copyright (C) 2013 Andrea Vedaldi.
% Copyright (C) 2010-11 Andrea Vedaldi.
% 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).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function lbl = loadLbl(lbl_file)

tname = tempname;
mkdir(tname);
try
  fprintf('Untarring project into %s\n',tname);
  untar(lbl_file,tname);
  fprintf('... done with untar.\n');
  rawLblFile = fullfile(tname,'label_file.lbl');
catch ME
  if strcmp(ME.identifier,'MATLAB:untar:invalidTarFile') || ...
      strcmp(ME.identifier,'MATLAB:io:archive:untar:invalidTarFile')
    warningNoTrace('Label file %s is not bundled. Using it in raw (mat) format.',lbl_file);
    rawLblFile = lbl_file;
  else
    ME.rethrow();
  end
end

warning('off','MATLAB:load:classNotFound');
warning('off','MATLAB:load:cannotInstantiateLoadedVariable');
lbl = load(rawLblFile,'-mat');
warning('on','MATLAB:load:classNotFound');
warning('on','MATLAB:load:cannotInstantiateLoadedVariable');

[success, message, ~] = rmdir(tname,'s');
if ~success
  error('Could not clear the temp directory %s\n',message);
else
  fprintf('Cleared out temp directory %s\n',tname);
end


text:
 function [sino, Hk, hn, nn] = fbp2_sino_filter(type, sino, varargin)
%function [sino, Hk, hn, nn] = fbp2_sino_filter(type, sino, [options])
%|
%| Apply ramp-like filters to sinogram(s) for 2D FBP image reconstruction.
%| Both parallel-beam and fan-beam tomographic geometries are supported.
%| This approach of sampling the band-limited ramp avoids the aliasing that
%| would be caused by sampling the ramp directly in the frequency domain.
%|
%| in
%|	type		'arc' (3rd generation CT) or 'flat' (for parallel too)
%|	sino	[nb (L)] sinogram(s)
%| options
%|	dr | ds	(real)	sample spacing (in distance units, e.g., cm) (default 1)
%|	dsd	(real)	source-to-detector distance, for 'arc' case only.
%|	extra		# of extra sinogram radial samples to keep (default: 0)
%|	npad		# of padded samples. (default: 0, means next power of 2)
%|	decon1		deconvolve effect of linear interpolator? (default: 1)
%|	window	[npad]	samples of apodization window function
%|			for [-np/2,...,np/2-1]. (default: '' = plain ramp)
%|			or a string like 'hann' for some predefined windows
%| out
%|	sino	[nb (L)] filtered sinogram rows
%|	Hk	[npad]	apodized ramp filter frequency response
%|	hn	[npad]	samples of band-limited ramp filter
%|	nn	[npad]	[-np/2,...,np/2-1] vector for convenience
%|
%| Copyright 2005-12-19, Jeff Fessler, University of Michigan

if nargin == 1 && streq(type, 'test'), fbp2_sino_filter_test, return, end
if nargin < 2, ir_usage, end

arg.ds = 1;
arg.dsd = [];
arg.window = '';
arg.extra = 0;
arg.npad = 0;
arg.decon1 = false;
arg = vararg_pair(arg, varargin, 'subs', {'dr', 'ds'; 'Dsd', 'dsd'});

dims = size(sino);
sino = reshape(sino, dims(1), []);
[sino Hk hn nn] = fbp2_sino_filter_do(type, sino, ...
	arg.ds, arg.dsd, arg.window, arg.extra, arg.npad, arg.decon1);
sino = reshape(sino, [size(sino, 1) dims(2:end)]);

%
% fbp2_sino_filter_do()
%
function [sino, Hk, hn, nn] = fbp2_sino_filter_do(type, sino, ...
	ds, dsd, window, extra, npad, decon1);

[nb na] = size(sino);
if ~npad
	npad = 2^ceil(log2(2*nb-1)); % padded size
%	printm('nb=%d npad=%d', nb, npad)
end
sino = [sino; zeros(npad-nb,na)]; % padded sinogram

[hn nn] = fbp_ramp(type, npad, ds, dsd);

Hk = reale(fft(fftshift(hn)));

Hk = Hk .* fbp2_window(npad, window);

Hk = ds * Hk; % differential for discrete-space convolution vs integral

% linear interpolation is like blur with a triangular response,
% so we can compensate for this approximately in frequency domain
if decon1
	Hk = Hk ./ fftshift(nufft_sinc(nn / npad).^2);
end

sino = ifft_sym( fft(sino, [], 1) .* repmat(Hk, [1 na]), [], 1); % apply filter

% trick: possibly keep extra column(s) for zeros!
sino = sino([1:(nb+extra)],:);
sino([(nb+1):(nb+extra)],:) = 0;

%
% test
%
function fbp2_sino_filter_test
nb = 2^5;
sino = zeros(nb,2); sino(nb/2+1,1) = 1;
ds = 1;
dsd = 30;
[sino1 H1 h1 nn] = fbp2_sino_filter('arc', sino, 'ds', ds, 'dsd', dsd);
[sino2 H2 h2 nn] = fbp2_sino_filter('flat', sino, 'ds', ds);
[sino2 H3 h2 nn] = fbp2_sino_filter('flat', sino, 'ds', ds, 'decon1', 1);
%max_percent_diff(h1, h2) % small!
%max_percent_diff(H1, H2) % small!

if im
	clf, subplot(121)
	plot(nn, h1, 'o', nn, h2, '+')
	xlabel 'n', ylabel 'h[n]'
	axis tight
	legend('arc', 'flat')

	subplot(122)
	plot(nn, fftshift(H1), 'o', nn, fftshift(H2), '+', nn, fftshift(H3), '.')
	xlabel 'k', ylabel 'H[k]', axis tight, axisy(0, 0.5)
	legend('arc', 'flat', 'decon', 'location', 'north')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function eddyscan_compiled_script(file_name, file_path, save_path, varargin)%#codegen
if ~strcmp(file_path(end), '/')
    file_path = strcat(file_path, '/');
end
if ~strcmp(save_path(end), '/')
    save_path = strcat(save_path, '/');
end
vars = load('area_map.mat');
area_map = vars.area_map;
[dir, rem] = strtok(file_name, '/');
filename = strtok(rem, '/');
eddy_dir = [save_path, dir];
indices = regexp(filename, '[0-9]');
nums = filename(indices);
date = nums(1:8);
eddy_file = ['anticyc_', date, '.mat'];
if exist(eddy_dir, 'dir')
    cd(eddy_dir);
    if exist(eddy_file, 'file')
        disp('Eddy file detected, quitting.');
        %quit;
    end
end
cd([file_path, dir]);
if exist([file_path, dir, '/', filename], 'file') && ~exist([save_path, dir, '/', eddy_file], 'file')
    ssh = ncread(filename, 'sla')';
    lat = double(ncread(filename, 'lat'));
    lon = double(ncread(filename, 'lon'));
    cd(save_path);
    ant_eddies = scan_single(ssh, lat, lon, date, 'anticyc', 'v2', area_map, varargin{:});%#ok
    cyc_eddies = scan_single(ssh, lat, lon, date, 'cyclonic', 'v2', area_map, varargin{:});%#ok
    if ~exist([save_path, dir], 'dir')
        mkdir([save_path, dir]);
    end
    cd([save_path, dir]);
    save(['anticyc_', date, '.mat'], 'ant_eddies');
    save(['cyclonic_', date, '.mat'], 'cyc_eddies');
else

end
%quit;
end

function [ eddies ] = scan_single( ssh, lat, lon, date, cyc, scan_type, areamap, varargin )
%SCAN_SINGLE Wrapper function to do scanning
% ssh: ssh slice with nans for land, size should be [length(lat) length(lon)]
% lat: 1D array of the latitudes of ssh grid
% lon: 1D array of the longitudes of ssh grid
% cyc: 'anticyc' or 'cyclonic'
% scan_type: 'v1', 'v2', 'hybrid'
%         v1: Will run top-down scanning (only works with full data of 0.25 x 0.25 ssh grid)
%         v2: Will run bottom-up scanning from the minima of the field
%     hybrid: Will run v2 and v1 scanning and will take the union of the
%             two sets where, for common features, v2 bodies will be used
% areamap: A 2D array that refer to the area of each pixel in SSH data (should have same size as ssh), or 1D array
%   that refer to area of each pixel for a specific lat in a regular grid (pixeld have same area for the same
%   latitude)
% Optional parameters (only applicable for v2 eddyscan):
%   'minimumArea': minimum number of pixels for an eddy, used for validating eddies, default value is 9
%   'thresholdStep': the minimum step for thresholding, the unit is SSH's unit, default value is 0.05
%   'isPadding': whether or not to pad SSH data, should be true when scanning SSH data with the longitudes expanding the
%   whole world dmap. Set to false if only partial SSH data is used. Default value is true
    if ~any(isnan(ssh(:)))
        error('Invalid ssh data, must contain NaNs for land values');
    end

    if ~all(size(ssh) == [length(lat) length(lon)])
        error('Invalid ssh data size, should be [length(lat) length(lon]');
    end

    if ~all(size(areamap) == size(ssh))
        % Not a 2d array with same size as ssh
        if ~any(size(areamap) == [1 1]) || length(areamap) ~= length(lat)
            disp('Invalid areamap, using NaN for eddy surface area');
        end
    end

    stype = get_stype(scan_type);
    ctype = get_ctype(cyc);
    disp('About to start scanning eddies');

    %oldpath = addpath('lib');
    switch stype
        case 1
            eddies = top_down_single(ssh, lat, lon, areamap, ctype, varargin{:});
        case 2
            eddies = bottom_up_single(ssh, lat, lon, areamap, ctype, varargin{:});
%        case 0
%            scanners = {@top_down_single, @bottom_up_single};
%            eddies_out = {[], []};
%            parfor i = 1:2
%                eddies_out{i} = scanners{i}(ssh, lat, lon, areamap, ctype, varargin{:});
%            end
%            eddies = get_combined_eddy_frames(eddies_out{2}, eddies_out{1}, ssh);
    end
    [eddies.Date] = deal(date);

    %path(oldpath);

end
%end

function cyc_t = get_ctype(cyc)
    switch cyc
        case 'anticyc'
            cyc_t = 1;
        case 'cyclonic'
            cyc_t = -1;
        otherwise
            error('cyc must be anticyc or cyclonic');
    end
end

function scan_t = get_stype(scan_type)
    switch scan_type
        case 'v1'
            scan_t = 1;
        case 'v2'
            scan_t = 2;
        case 'hybrid'
            scan_t = 0;
        otherwise
            error('scan_type must be v1, v2, or hybrid');
    end
end

function [ eddies ] = bottom_up_single(ssh_data, lat, lon, areamap, cyc, varargin)
%BOTTOM_UP_SINGLE Finds eddies using the Bottom Up method
%   Will return an array of struct's that contain the eddy data.
%   ssh_data: A 2D array of double's that contain the sea surface heights (latsxlons)
%   lat: A 1D array of double's that gives the latitude for a given index (dimension should match
%         that of ssh_data)
%   lon: A 1D array of double's that gives the longitude for a given index (dimension should match
%         that of ssh_data)
%   areamap: A 2D array that refer to the area of each pixel in SSH data (should have same size as ssh), or 1D array
%   that refer to area of each pixel for a specific lat in a regular grid (pixeld have same area for the same
%   latitude)
%   cyc: Pass 1 to output anticyclonic eddies or -1 to output cyclonic eddies
%   Optional parameters:
%   'minimumArea': minimum number of pixels for an eddy, used for validating eddies
%   'thresholdStep': the minimum step for thresholding, the unit is SSH's unit
%   'isPadding': whether or not to pad SSH data, should be true when scanning SSH data of the whole map. Set to false if
%   only partial SSH data is used.
%   'sshUnits': The units the SSH data is in. bottom_up_single is built to work natively on centimeter SSH data.
%   Valid parameters are 'meters' and 'centimeters'. If the paramater passed in is 'meters', the SSH data will
%   be multiplied by 100. No changes will be made if the paramater passed in is 'centimeters'.
%   The default value of 'sshUnits' is centimeters.
    p = inputParser;
    defaultMinPixelSize = 9;
    defaultThresholdStep = 0.05;
    defaultSSHUnits = 'centimeters';
    defaultPaddingFlag = true;
    addRequired(p, 'ssh_data');
    addRequired(p, 'lat');
    addRequired(p, 'lon');
    addRequired(p, 'areamap');
    addRequired(p, 'cyc');
    addParameter(p, 'minimumArea', defaultMinPixelSize);%, @isnumeric);
    addParameter(p, 'thresholdStep', defaultThresholdStep);%, @isnumeric);
    addParameter(p, 'isPadding', defaultPaddingFlag);
    addParameter(p, 'sshUnits', defaultSSHUnits);
    parse(p, ssh_data, lat, lon, areamap, cyc, varargin{:});
    minimumArea = p.Results.minimumArea;
    thresholdStep = p.Results.thresholdStep;
    isPadding = p.Results.isPadding;
    SSH_Units = p.Results.sshUnits;
    disp(minimumArea);
    if isa(minimumArea, 'char')
        disp('Minimum area was a string, converting to double.');
        minimumArea = str2double(minimumArea);
        disp(minimumArea);
    end
    if isa(thresholdStep, 'char')
        disp('Threshold step was a string, converting to double.');
        thresholdStep = str2double(thresholdStep);
        disp(thresholdStep);
    end
    if strcmp(SSH_Units, 'meters')
        ssh_data = ssh_data * 100;
    elseif strcmp(SSH_Units, 'centimeters')
        max_val = max(ssh_data(:));
        min_val = max(ssh_data(:));
        if max_val < 1 && min_val > -1
            ssh_data = ssh_data * 100;
        elseif max_val < 100 && min_val > -100

        else
            error('Could not figure out what units the SSH data provided is in. Please specify it as an additional parameter: sshUnits');
        end
    end

    %Check if the grid is regular (differences between lats and lons are equal)
    lat_diffs = lat(2:end) - lat(1:end-1);
    lat_diffs2 = lat_diffs(2:end) - lat_diffs(1:end-1);
    lon_diffs = lon(2:end) - lon(1:end-1);
    lon_diffs(lon_diffs <= -180) = lon_diffs(lon_diffs <= -180) + 360;
    lon_diffs(lon_diffs >= 180) = lon_diffs(lon_diffs >= 180) - 360;
    lon_diffs = abs(lon_diffs);
    lon_diffs2 = lon_diffs(2:end) - lon_diffs(1:end-1);
    if all(lat_diffs2 == 0) && all(lon_diffs2 == 0)
        % Regular grid, create a georasterref object to get eddy's centroid
        geo_raster_lat_limit = [lat(1) lat(end)];
        if lon(1) > lon(end)
            geo_raster_lon_limit = [lon(1) (360 + lon(end))];
        else
            geo_raster_lon_limit = [lon(1) lon(end)];
        end

        R = georasterref('LatLim', geo_raster_lat_limit, 'LonLim', geo_raster_lon_limit, 'RasterSize', ...
         size(ssh_data), 'ColumnsStartFrom', 'south', 'RowsStartFrom', 'west');
    else
        % Use normal indexing to get eddy's centroid
        R = [];
    end
    disp('About to get extrema');

    extrema = get_extrema(ssh_data, cyc);
    disp('Got extrema');
    if isPadding
        origExtrema = extrema;
        extrema = [zeros(size(extrema, 1), 200), extrema, zeros(size(extrema, 1), 200)];
        sshExtended = [ssh_data(:, end-199:end), ssh_data(:, :), ssh_data(:, 1:200)];
        [extrema_lat_indexes, extrema_lon_indexes] = ind2sub(size(extrema), find(extrema == 1));

        extrema(:, 1:200) = origExtrema(:, end-199:end);
        extrema(:, end-199:end) = origExtrema(:, 1:200);
    else
        [extrema_lat_indexes, extrema_lon_indexes] = ind2sub(size(extrema), find(extrema == 1));
        sshExtended = ssh_data;
    end

    disp('Scanning');
    eddies = new_eddy();
    eddies(length(extrema_lat_indexes)).Date = NaN;
    cyc_sshExtended = sshExtended * cyc;
    parfor i = 1:length(extrema_lat_indexes) % Normally a parfor. Modified to a for loop solely for itasca testing
        curr_lat_index = extrema_lat_indexes(i); curr_lon_index = extrema_lon_indexes(i);
        e = thresholdBU(cyc, curr_lat_index-5, curr_lat_index+5, curr_lon_index-5, curr_lon_index+5, ...
            sshExtended, extrema, curr_lat_index, curr_lon_index, sshExtended(curr_lat_index, curr_lon_index), ...
            thresholdStep, NaN, ...
            zeros(size(sshExtended)), lat, lon, R, areamap, minimumArea, isPadding, cyc_sshExtended);
        %disp(e);
        if ~isempty(e.Stats)
            %disp(['Adding eddy at index: ', num2str(i)]);
            %disp(e);
            eddies(i) = e;
        end
    end
    %lat_array = [eddies.Lat];
    %disp(lat_array);
    mask = false(1, length(eddies));
    for i = 1:length(eddies)
        if isempty(eddies(i).Lat)
    %        disp(['Empty lat array index at ', num2str(i)]);
            mask(i) = true;
        end
    end
    disp('Scanned');
    %mask = cellfun('isempty', {eddies.Lat});
    eddies = eddies(~mask);
    %lat_array = [eddies.Lat];
    %disp(lat_array);
end

function [eddy] = thresholdBU(cyc, block_bottom_index, block_top_index, block_left_index, block_right_index, ...
        ssh, extrema, extrema_lat_index, extrema_lon_index, thresh, threshold_step, last_step, previous, ...
        lat, lon, R, areamap, min_pixel_size, is_padding, cyc_ssh)
% THRESHOLDBU Get an eddy by bottom up method
%   cyc: 1 for anticyclonic and -1 for cyclonic
%   block_bottom(top/left/right)_index: index of the bottom/top/left/right of the block that will be used for
%   thresholding
%   ssh: ssh data(extended if is_padding is true)
%   extrema: logical index of ssh extrema
%   extrema_lat/lon_index: lat/lon index of the extremum that is being used to find an eddy
%   thresh: current threshold that is being used to find an eddy
%   threshold_step: the step to increase/decrease threshold value, based on eddy type
%   last_step: the last step was used for thresholding
%   previous: 2d logical array of the connected component that contains the extremum in the last thresholdBU call
%   lat: 1d array of latitudes of the ssh grid
%   lon: 1d array of longitudes of the ssh grid
%   R: the georasterref object for SSH grid
%   areamap: 1D or 2D array of reference to area of each pixel in SSH grid
%   min_pixel_size: minimum number of pixels for an eddy
%   is_padding: whether or not the SSH data is padded

    switch cyc
        case 1
            intensity = 'MaxIntensity';
        case -1
            intensity = 'MinIntensity';
        otherwise
            error('Invalid cyc');
    end

    if block_bottom_index < 1 || block_top_index > size(ssh, 1) || ...
            block_left_index < 1 || block_right_index > size(ssh, 2)
        edgeOfWorld = true;

        while block_bottom_index < 1 || block_top_index > size(ssh, 1) || ...
            block_left_index < 1 || block_right_index > size(ssh, 2)
            % Make sure that the block is inside the grid

            block_bottom_index = block_bottom_index + 1;
            block_top_index = block_top_index - 1;
            block_left_index = block_left_index + 1;
            block_right_index = block_right_index - 1;
            if block_top_index <= block_bottom_index + 2 || block_right_index <= block_left_index + 2
                % If the block is too small, just return an empty eddy
                eddy = new_eddy();
                return;
            end
        end

        block = ssh(block_bottom_index:block_top_index, block_left_index:block_right_index);

        extremaBlock = extrema(block_bottom_index:block_top_index, block_left_index:block_right_index);

    else
        edgeOfWorld = false;
        block = ssh(block_bottom_index:block_top_index, block_left_index:block_right_index);
        extremaBlock = extrema(block_bottom_index:block_top_index, block_left_index:block_right_index);
    end

    if isnan(last_step)
        step = threshold_step;
    else
        step = last_step;
    end

    iter = 1;
    while true
        iter = iter+1;
        if iter > 5000

            perim = imdilate(logical(current), ones(3)) & ~logical(current);
            if all(isnan(block(perim)))
                eddy = new_eddy();
                return;
            end
            disp('potential infinite loop')
        end

        bw = cyc .* block >= cyc .* thresh;
        labels = bwlabel(bw);

        extrema_label = labels(extrema_lat_index - block_bottom_index + 1, extrema_lon_index - block_left_index + 1);
        current = labels == extrema_label;
        currentExtrema = extremaBlock(current);

        existing_pixel_at_box_edge = outterRing(labels, extrema_label);

        if sum(currentExtrema) > 1 || ( edgeOfWorld && existing_pixel_at_box_edge)

            if step ~= threshold_step
                % Go back to last threshold
                thresh = thresh + cyc*step;
                step = threshold_step;
                thresh = thresh - cyc*step;
                continue;
            end

            if size(block, 1) ~= size(previous(block_bottom_index:block_top_index, block_left_index:block_right_index), 1)
                prevBlock = block(2:end-1, 2:end-1);
            else
                prevBlock = block;
            end

            perim = imdilate(logical(previous(block_bottom_index:block_top_index, block_left_index:block_right_index)), ones(3)) ...
                & ~logical(previous(block_bottom_index:block_top_index, block_left_index:block_right_index));
            nan = isnan(prevBlock(perim));
            if sum(nan) / length(nan) > .3
                %if more than half of your perimeter is land, then throw it out.
                eddy = new_eddy();
                return;
            end

            if sum(previous(:)) < min_pixel_size
                    eddy = new_eddy();
                    return;
            end

            perim = bwperim(previous);
            meanPerim = mean(ssh(logical(perim)));
            amp = cyc * (ssh(extrema_lat_index, extrema_lon_index)-meanPerim);

            stats = regionprops(previous, ssh, 'Area', 'Extrema',...
                'PixelIdxList', intensity, 'ConvexImage', 'PixelList', ...
                'Solidity', 'Extent', 'Orientation', 'MajorAxisLength', ...
                'MinorAxisLength');

            stats.Intensity = stats.(intensity);
            stats = rmfield(stats, intensity);

            if is_padding
                [idx, r, c] = extidx2original(stats.PixelIdxList, [length(lat) length(lon)], size(ssh));
                stats.PixelIdxList = idx;
            else
                [r, c] = ind2sub(size(ssh), stats.PixelIdxList);
            end

            stats.PixelList = [r, c];

            % Getting geodesic speed
            if is_padding
                geoSpeed = mean_geo_speed(ssh(:, 201:end-200), stats.PixelIdxList, lat, lon);
                if ~isempty(R)
                    [elat, elon] = weighted_centroid(cyc_ssh(:, 201:end-200), stats.PixelList, stats.PixelIdxList, R);
                else
                    [elat, elon] = weighted_centroid_irregular_grid(cyc_ssh(:, 201:end-200), stats.PixelList, stats.PixelIdxList, lat, lon);
                end
            else
                geoSpeed = mean_geo_speed(ssh, stats.PixelIdxList, lat, lon);
                if ~isempty(R)
                    [elat, elon] = weighted_centroid(cyc_ssh, stats.PixelList, stats.PixelIdxList, R);
                else
                    [elat, elon] = weighted_centroid_irregular_grid(cyc_ssh, stats.PixelList, stats.PixelIdxList, lat, lon);
                end
            end

            % weighted_centroid returns lon from 0-360, fix this
            % TODO: should we also fix lat lon -270 to 80?
            elon = (elon > 180).*(elon - 360) + (elon <= 180).*elon;

            % Getting surface area of the eddy
            if all(size(areamap) == [length(lat) length(lon)])
                % area is 2D array for areas of pixels at [lat, lon]
                sarea = sum(areamap(stats.PixelIdxList));
            elseif any(size(areamap) == [1 1]) && length(areamap) == length(lat)
                % Area is 1D array for areas of pixels at a specific latitude
                sarea = sum(areamap(stats.PixelList(:, 1)));
            else
                % Invalid areamap
                sarea = NaN;
            end

            eddy = new_eddy(rmfield(stats, 'PixelList'), amp, elat, elon, thresh, sarea, cyc, geoSpeed, 'ESv2');

            return
        end

        if existing_pixel_at_box_edge
            %disp('expanding size');
            eddy = thresholdBU(cyc, block_bottom_index-1, block_top_index+1, block_left_index-1,...
                block_right_index+1, ssh, extrema, extrema_lat_index, extrema_lon_index, thresh, threshold_step, ...
                step, previous, lat, lon, R, areamap, min_pixel_size, ...
                is_padding, cyc_ssh);
            return
        end

        previous(block_bottom_index:block_top_index, block_left_index:block_right_index) = current;

        step = step * 2; % double the step for less number of iterations
        thresh = thresh - cyc*step;
    end

end

function [ extrema ] = get_extrema( ssh, cyc )
%GET_EXTREMA Returns a matrix containing all of the minima or maxima
%(depending on the value of cyc) in a 5x5 matrix within the 2D ssh field.
% ssh: ssh slice containing NaNs for land
% cyc: 1 for anticyclonic, -1 for cyclonic

    padded = [ssh(:,end-1:end) ssh(:,:) ssh(:,1:2)];
    padded(isnan(padded)) = cyc*-Inf;
    padded = padarray(padded, [1, 1], cyc*-Inf);
    n = ones(5); n(3, 3) = 0;
    padded = cyc .* padded; % Want to find right extrema for cyclonic and anticyc eddies

    extrema = padded > imdilate(padded, n);
    extrema = extrema(2:end-1, 4:end-3);
end

function [res] = outterRing(box, val)
    ring = [box(1, 1:end)'; box(end, 1:end)'; box(1:end, 1); box(1:end, end)];
    res = any(ring == val);
end

function [idx, row, col] = extidx2original(idx, original_size, extended_size)
%EXTIDX2ORIGINAL Convert from extended indexes to original indexes
    [row, col] = ind2sub(extended_size,idx);

    offright = col > (extended_size(2) + original_size(2)) / 2;
    offleft = col < (extended_size(2) - original_size(2)) / 2 + 1;
    notoff = ~(offleft | offright);
    col(offright)=col(offright) - (extended_size(2) + original_size(2)) / 2;
    col(offleft) = col(offleft) + original_size(2) - (extended_size(2) - original_size(2)) / 2;
    col(notoff)=col(notoff) - (extended_size(2) - original_size(2)) / 2;

    idx = sub2ind(original_size,row,col);
end

function mean_speed = mean_geo_speed(ssh, pixels, lat, lon)
%MEAN_GEO_SPEED Returns the mean geostrophic speed for pixels.
% lat and lon should yield the correct values for indices of ssh.

    g = 980.665; % cm/s
    omega = 7.2921e-5;
    [x, y] = ind2sub(size(ssh), pixels);
    lats = lat(x);
    f = 2*omega*sin(lats);
    f(f == 0) = 2*omega*sin(0.25); % f is coriolis frequency

    dSSH_y = ssh(sub2ind(size(ssh), x, mod(y, size(ssh, 2))+1)) - ...
        ssh(sub2ind(size(ssh), x, mod(y-2, size(ssh, 2))+1));
    dSSH_x = ssh(sub2ind(size(ssh), min(x+1, zeros(size(x)) + size(ssh, 1)), y))...
        - ssh(sub2ind(size(ssh), max(x-1, ones(size(x))), y));

    dy = deg2km(distance(lat(x), lon(mod(y, size(ssh, 1))+1), lat(x), lon(mod(y-2, size(ssh, 1))+1))) * 100000;
    dx = (min(x+1, size(ssh, 1)) - max(1, x-1)) .* 111.12 .* 100000 .* 180 ./ (length(lats) - 1);

    vs = -g .* (dSSH_y) ./ (2.*f .* dy);
    us = g .* (dSSH_x) ./ (2 .* f .* dx );
    speeds = sqrt(us .^2 + vs .^2);
    mean_speed = nanmean(speeds);
end

function [lat, lon] = weighted_centroid(cyc_ssh, pixellist, pixelidxlist, R)
%WEIGHTED_CENTROID Returns the location of the weighted centroid for the
%pixels provided
    %ssh = cyc * ssh;
    shift = min(cyc_ssh(pixelidxlist));

    x = pixellist(:, 1);
    y = pixellist(:, 2);

    if min(y) == 1 && max(y) == size(cyc_ssh,2)
        y(y > size(cyc_ssh,2)/2) = y(y > size(cyc_ssh,2)/2) - size(cyc_ssh,2);
    end

    mask = ~isnan(cyc_ssh(pixelidxlist));
    intensities = cyc_ssh(pixelidxlist)+shift;
    intensities = intensities - min(intensities); % should start from 0

    xbar = sum(x(mask) .* intensities(mask).^2) / sum(intensities(mask).^2);
    ybar = sum(y(mask) .* intensities(mask).^2) / sum(intensities(mask).^2);

    if ybar <= 0
        ybar = ybar + size(cyc_ssh,2);
    end

    [lat, lon] = pix2latlon(R, xbar, ybar);

end

function [lat, lon] = weighted_centroid_irregular_grid(cyc_ssh, pixellist, pixelidxlist, lats, lons)
%WEIGHTED_CENTROID_IRREGULAR_GRID Returns the location of the weighted centroid for the
%pixels provided
    %ssh = cyc * ssh;
    shift = min(cyc_ssh(pixelidxlist));

    x = pixellist(:, 1);
    y = pixellist(:, 2);

    if min(y) == 1 && max(y) == size(cyc_ssh,2)
        y(y > size(cyc_ssh,2)/2) = y(y > size(cyc_ssh,2)/2) - size(cyc_ssh,2);
    end

    mask = ~isnan(cyc_ssh(pixelidxlist));
    intensities = cyc_ssh(pixelidxlist)+shift;
    intensities = intensities - min(intensities); % should start from 0
    xbar = sum(x(mask) .* intensities(mask).^2) / sum(intensities(mask).^2);
    ybar = sum(y(mask) .* intensities(mask).^2) / sum(intensities(mask).^2);

    if ybar <= 0
        ybar = ybar + size(cyc_ssh,2);
    end

    x_lower = floor(xbar);
    x_upper = ceil(xbar);
    if x_lower == 0
        lat = lats(1);
    elseif x_upper == length(lats) + 1
        lat = lats(end);
    else
        lat_lower = lats(x_lower);
        lat_upper = lats(x_upper);
        lat = lat_lower + (lat_upper - lat_lower) * (xbar - x_lower) / (x_upper - x_lower);
        if isnan(lat)
            lat = 90;
        end
    end
    y_lower = floor(ybar);
    y_upper = ceil(ybar);
    if y_lower == 0
        lon = lons(1);
    elseif y_upper == length(lons) + 1
        lon = lons(end);
    else
        lon_lower = lons(y_lower);
        lon_upper = lons(y_upper);
        if lon_upper - lon_lower > 180
            lon_upper = lon_upper - 360;
        elseif lon_lower - lon_upper > 180
            lon_lower = lon_lower - 360;
        end
        lon = lon_lower + (lon_upper - lon_lower) * (ybar - y_lower) / (y_upper - y_lower);
    end

end

function eddy = new_eddy(STATS, amplitude, lat, lon, thresh, sa, cyc, geospeed, detect)
%NEW_EDDY Initializes new eddy objects, run with no arguments to create an
%empty matrix
    if nargin
        eddy = struct('Stats', STATS, ...
            'Lat', lat, ...
            'Lon', lon, ...
            'Amplitude', amplitude, ...
            'ThreshFound', thresh, ...
            'SurfaceArea', sa, ...
            'Date', NaN, ...
            'Cyc', cyc, ...
            'MeanGeoSpeed', geospeed, ...
            'DetectedBy', detect);
    else
        eddy = struct('Stats', [], ...
            'Lat', [], ...
            'Lon', [], ...
            'Amplitude', [], ...
            'ThreshFound', [], ...
            'SurfaceArea', [], ...
            'Date', [], ...
            'Cyc', [], ...
            'MeanGeoSpeed', [], ...
            'DetectedBy', []);
    end
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
%% Anisotropic Diffusion
%
% This sample demonstrates Perona-Malik anisotropic diffusion.
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.3.0/modules/ximgproc/samples/filterdemo.cpp>
%

function varargout = anisodiff_demo_gui(im)
    % load source image
    if nargin < 1
        src = cv.imread(fullfile(mexopencv.root(),'test','fruits.jpg'), 'Color',true);
    elseif ischar(im)
        src = cv.imread(im, 'Color',true);
    else
        src = im;
    end

    % not too big
    if size(src,1) > 480
        src = cv.resize(src, [round(480*size(src,2)/size(src,1)), 480]);
    end

    % create the UI
    h = buildGUI(src);
    if nargout > 0, varargout{1} = h; end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    niters = round(get(h.slid, 'Value'));
    set(h.txt, 'String',sprintf('No. of time steps: %2d',niters));

    % apply filtering
    dst = cv.anisotropicDiffusion(h.src, ...
        'Alpha',1.0, 'K',0.02, 'Iterations',niters);

    % show result
    set(h.img, 'CData',dst);
    drawnow;
end

function h = buildGUI(img)
    %BUILDGUI  Creates the UI

    % parameters
    niters = 10;
    sz = size(img);
    sz(2) = max(sz(2), 250);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Anisotropic Diffusion', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+30-1]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 30 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(img, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(img);
    end
    h.txt = uicontrol('Parent',h.fig, 'Style','text', ...
        'Position',[5 5 150 20], 'FontSize',11, ...
        'String',sprintf('No. of time steps: %2d',niters));
    h.slid = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Position',[155 5 sz(2)-155-5 20], 'Value',niters, ...
        'Min',0, 'Max',30, 'SliderStep',[1 5]./(30-0));

    % hook event handlers, and trigger default start
    set(h.slid, 'Callback',{@onChange,h}, ...
        'Interruptible','off', 'BusyAction','cancel');
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function varargout = prtPlotUtilMatrixTable(X,cLim,cMap,num2strFormat,textCMap)
% prtPlotUtilMatrixTable displays a matrix as a table. The table is actually an
%   image with lines between each pixel.
%
% Syntax: [axesHandle, textHandles, verticleLineHandles, horizontalLineHandles]
%           = prtPlotUtilMatrixTable(X,cLim,cMap,num2strFormat,textCMap)
%
% Inputs:
%   X - 2-D Matrix to plot
%   cLim - The limits for image scaling default = [min(X(:)), max(X(:))]
%   cMap - The colormap for image default = gray;
%   num2strFormat - The string format for each matrix element
%       default = '%.2f'
%   textCMap - The colorMap to use for the text Color.
%       default - [cMap(end,:); cMap(1,:)];
%           If cMap has only 1 color the default = 1-cMap.
%
% Outputs: These names are all fairly descriptive
%   axesHandle
%   textHandles
%   verticleLineHandles
%   horizontalLineHandles
%
% Example:
%   prtPlotUtilMatrixTable(rand(5),[0 1],summer(128),'%.2f',[0 0 0])
%   prtPlotUtilMatrixTable(rand(5),[0 1],summer(128),'%.2f')
%   prtPlotUtilMatrixTable(rand(5),[0 1],[1 1 1],[],[0 0 0])
%
% Other m-files required: none
% Subfunctions: getBestFontSize
% MAT-files required: none
%
% See also: plotConfusionMatrix.m

if nargin < 2 || isempty(cLim)
    cLim = [min(X(:)) max(X(:))];
end

if nargin < 3 || isempty(cMap)
    cMap = gray;
end

if nargin < 4 || isempty(num2strFormat)
    num2strFormat = [];
end

if nargin < 5 || isempty(textCMap)
    if size(cMap,1) == 1;
        textCMap = 1-cMap;
    else
        textCMap = [cMap(end,:); cMap(1,:)];
    end
end

[nRows, nCols] = size(X);

imageAxes = gca;
imagesc(X,cLim);
colormap(cMap);
fontSize = max(getBestFontSize(imageAxes),1);
textHandles = prtUtilPreAllocateHandles(size(X));
verticleLineHandles = prtUtilPreAllocateHandles(nRows,1);
horizontalLineHandles = prtUtilPreAllocateHandles(nCols,1);
hold on;

[dontNeed, textCMapInds] = histc( (X-cLim(1))./(cLim(2)-cLim(1)) , linspace(0,1+eps,size(textCMap,1)+1)); %#ok<ASGLU>
textCMapInds(textCMapInds==0 | ~isfinite(textCMapInds)) = size(textCMap,1);
textCMapInds(textCMapInds==0 | ~isfinite(textCMapInds)) = 1;
textCMapInds(textCMapInds > size(textCMap,1)) = size(textCMap,1);

for iRow = 1:nRows
    for jCol = 1:nCols
        cNum = X(iRow,jCol);

        cTextColor = textCMap(textCMapInds(iRow,jCol),:);

        cTextString = num2str(cNum,num2strFormat);

        % Some decimal place pruning
        done = false;
        while ~done
            if length(cTextString) > 1 && strcmpi(cTextString(end),'0') && any(cTextString == '.')
                cTextString(end) = [];
            else
                done = true;
            end
        end
        % Remove last decimal if necessary
        if strcmpi(cTextString(end),'.')
            cTextString(end) = [];
        end

        textHandles(iRow,jCol) = text(jCol,iRow,cTextString,...
            'color',cTextColor,'horizontalAlignment','center',...
            'fontsize',fontSize,'clipping','on','visible','on');

        if iRow == 1
            horizontalLineHandles(jCol) = plot([jCol jCol]+0.5,[0.5 0.5+nRows],'k','linewidth',1);
        end
    end
    verticleLineHandles(iRow) = plot([0.5 0.5+nCols],[iRow, iRow]+0.5,'k','linewidth',1);
end
hold off;

set(imageAxes,'YTick',1:nRows,'XTick',1:nCols,...
    'XTickLabel',num2str((1:nCols)'),'YTickLabel',num2str((1:nRows)'),...
    'Xlim',[0.5 nCols+0.5],'Ylim',[0.5 nRows+0.5],...
    'TickLength',[0 0]);

varargout = {};
if nargout > 0
    varargout = {imageAxes, textHandles, verticleLineHandles, horizontalLineHandles};
end

f = ancestor(imageAxes,'figure');
set(f,'ResizeFcn',@(src,evt)setBestFontSize(imageAxes,textHandles));
setBestFontSize(imageAxes,textHandles);

function setBestFontSize(imAxes,textHandles)

if ~ishandle(imAxes)
    return
end
try
    fs = getBestFontSize(imAxes);
catch %#ok<CTCH>
    return
end

for iHandle = 1:numel(textHandles)
    if ishandle(textHandles(iHandle))
        if fs > 0
            set(textHandles(iHandle),'FontSize',fs);
            set(textHandles(iHandle),'visible','on');
        else
            set(textHandles(iHandle),'visible','off');
        end
    end
end

function fs = getBestFontSize(imAxes)
% I found this little gem in a MATLAB central file heatmaptext.m
% This should solve alot of our problems with this function
% It adjusts the font size relative to the number of col and rows.
%
% 31-Jan-2008 - I have modified this so it checks the extent of the axes
% instead of the figure. This allows the fontsize to change relative to a
% subplot. - KDM
%
% 17-May-2011 - I have modified this to check the minimumum of both ratios
% so that resized plots look good. Also added a listener above so that the
% resize will call this function again.
%
%
% Apparently this is copyright the MathWorks as is stated in the funciton.
%
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=15877&objectType=file#
%

% Try to keep font size reasonable for text
nrows = diff(get(imAxes,'YLim'));
ncols = diff(get(imAxes,'XLim'));

set(imAxes,'units','pixels');
extent = get(imAxes,'Position');
set(imAxes,'units','normalized');
ratioNum = extent(3:4);

magicNumber = 80;
if ncols < magicNumber && nrows < magicNumber
    ratio = min(ratioNum./[nrows,ncols]);
elseif ncols < magicNumber
    ratio = ratioNum(2)/ncols;
elseif nrows < magicNumber
    ratio = ratioNum(1)/nrows;
else
    ratio = 1;
end

%fs = min(maxFontSize,ceil(ratio/4));    % the gold formula
fs = ceil(ratio/4);
if fs < 4 % Font sizes less than 4 still look like crap
    fs = 0;
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function y = cvx_isvalid( x, full )
y = isnan( x ) | isinf( x );
if nargin < 2 || ~full,
	y = ~any( y(:) );
else
	y = ~y;
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function cdf = quasigeometric_cdf ( x, a, b )

%*****************************************************************************80
%
%% QUASIGEOMETRIC_CDF evaluates the Quasigeometric CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the number of trials.
%
%    Input, real A, the probability of 0 successes.
%    0.0 <= A <= 1.0.
%
%    Input, real B, the depreciation constant.
%    0.0 <= B < 1.0.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x < 0 )
    cdf = 0.0;
  elseif ( x == 0 )
    cdf = a;
  elseif ( b == 0.0 )
    cdf = 1.0;
  else
    cdf = a + ( 1.0 - a ) * ( 1.0 - b^x );
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rls.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [x,act]=rls(A,b,fac);
% function [x,act]=rls(A,b,x0,r);
% minimizes ||Ax-b||_2^2  s.t. |x-x0|<=r
% act=(is some constraint active?)
%
% if nargin<4, x0=0 and r(k)=fac*||b||/||A(:,k)||
% (default: fac=1000);
%
function [x,act]=rls(A,b,x0,r);

% set defaults
if     nargin<3, fac=1000;x0=zeros(size(A,2),1);
elseif nargin<4, fac=x0;x0=zeros(size(A,2),1);
else             b=b-A*x0;
end;

x=x0;
aa=diag(A'*A);

% handle zero columns directly
ind=find(aa>0);
n=length(ind);
if n==0, act=1; return; end;

A=A(:,ind);
gamma=b'*b;
if gamma==0, act=1; return; end;
if nargin<4, r=fac*sqrt(gamma./aa(ind));
else         r=r(ind);
end;

% now we need to minimize ||Az-b||^2  s.t. |z|<=r
prt=0;
[z,fct,ier]=minq(0,-A'*b,A'*A,-r,r,prt);
act=max(abs(z)./r);
x(ind)=x(ind)+z;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function varargout = minandmax2est(varargin)
%MINANDMAX2EST   Estimates the minimum and maximum of a DISKFUN.
%   mM = MINANDMAX2EST(F) returns estimates for the minimum and maximum of the
%   DISKFUN F over its domain.  mM is a vector of length 2 such that mM(1) is
%   the estimated minimum and mM(2) is the estimated maximum.
%
%   mM = MINANDMAX2EST(F, N) returns estimates for the minimum and maximum of
%   the DISKFUN F over its domain, based on samples on an N by N grid
%   (N = 33 by default).
%
% See also DISKFUN/MINANDMAX2.

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

[varargout{1:nargout}] = minandmax2est@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function C = ecov_music_1d(design, wavelength, doas, P, noise_var, snapshot_count)
%ECOV_MUSIC_1D Asymptotic covariance matrix of the estimation errors
%of the classical MUSIC algorithm.
%Syntax:
%   C = ECOV_MUSIC_1D(design, wavelength, doas, P, noise_var[, snapshot_count])
%Inputs:
%   design - Array design.
%   wavelength - Wavelength.
%   doas - DOA vector in radians.
%   P - Source covariance matrix. If all sources are uncorrelated and
%       shares the same power, you can just pass in a scalar. If all
%       sources are uncorrelated but have different powers, you can just
%       pass in a vector.
%   noise_var - Noise power.
%   snapshot_count - (Optional) number of snapshots. Default is one.
%Output:
%   C - Asymptotic error covariance matrix.
%Reference:
%   [1] P. Stoica and A. Nehorai, "MUSIC, maximum likelihood, and
%       Cramer-Rao bound: further results and comparisons," IEEE
%       Transactions on Acoustics, Speech and Signal Processing, vol. 38,
%       no. 12, pp. 2140-2150, Dec. 1990.
%   [2] P. Stoica and A. Nehorai, "MUSIC, maximum likelihood, and
%       Cramer-Rao bound," IEEE Transactions on Acoustics, Speech and
%       Signal Processing, vol. 37, no. 5, pp. 720-741, May 1989.
if design.dim ~= 1
    error('1D array expected.');
end
if nargin <= 5
    snapshot_count = 1;
end
m = design.element_count;
k = length(doas);
P = unify_source_power_matrix(P, k);
[A, D] = steering_matrix(design, wavelength, doas);
H = D'*(eye(m) - A/(A'*A)*A')*D;
B = P\(P + noise_var*eye(k)/(A'*A))/P;
h = real(1 ./ diag(H));
C = (noise_var/2/snapshot_count) * (real(H .* B.') .* (h*h'));
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function pass = test_undampedNewton(pref)
%TEST_UNDAMPEDNEWTON    Test what happens when we turn off damping in Newton

% This test tries to solve steady state Allen-Cahn using undamped Newton
% iteration.

if ( nargin == 0 )
    pref = cheboppref();
end

% Turn off damping
pref.damping = 0;

%% Setup and solve
dom = [0, 10];
x = chebfun(@(x) x, dom);
f = sin(x);
N = chebop(@(u) diff(u,2) + u - u.^3, dom, 1, -1);
u = solvebvp(N, f, pref);

% Happy?
err = norm(N(u) - f);
pass = ( err < 1e-8 );

% end


text:
function idx=diagElIdx(N,minorDiag)
%%DIAGELIDX For an NXN matrix, determine the linear indices of the elements
%           on the main or minor diagonal of the matrix.
%
%INPUTS:  N The size of the NXN matrix.
% minorDiag Optionally, this can specify whether the main or minor diagonal
%           of the matrix is desired. minorDiag=false (the default if
%           omitted or an empty matrix is passed) chooses the main
%           diagonal, and minorDiag=true chooses the minor diagonal of the
%           matrix.
%
%OUTPUTS: idx A 1XN vector of the indices on the main (or minor) diagonal
%             of the matrix. Minor diagonal elements go from the lower-left
%             corner of the matrix to the upper-right corner.
%
%The major diagonal is also known as the principal diagonal, the primary
%diagonal, the major diagonal and the leading diagonal. The minor diagonal
%of a matrix is also known as the antidiagonal, the counterdiagonal, the
%trailing diagonal, and the secondary diagonal.
%
%EXAMPLE:
% N=8;
% M=magic(N)
% mainDiagEls=M(diagElIdx(N))
% minorDiagEls=M(diagElIdx(N,true))
%One will see that mainDiagEls contains the main diagonal elements of M,
%and minorDiagEls contains the minor diagonal elements, which increase
%sequentialy.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(minorDiag))
    minorDiag=false;
end

if(minorDiag==false)
    idx=1:(N+1):(N^2);
else
    idx=N:(N-1):(N^2-N+1);
end
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:
% Test file for trigtech/poly.m

function pass = test_poly(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

%%
% Check a few simple examples.

f = testclass.make(@(x) zeros(size(x)), [], pref);
p = poly(f);
pass(1) = (norm(p, inf) <= 10*vscale(f)*eps);

f = testclass.make(@(x) 3*ones(size(x)), [], pref);
p = poly(f);
pass(2) = (norm(p - 3, inf) < 10*vscale(f)*eps);

f = testclass.make(@(x) 1+cos(pi*x), [], pref);
p = poly(f);
pass(3) = (norm(p - [0.5 1 0.5], inf) < 10*vscale(f)*eps);

f = testclass.make(@(x) 1 + exp(2*1i*pi*x) + exp(-1i*pi*x), [], pref);
p = poly(f);
pass(4) = (norm(p - [0 1 1 0 1], inf) ...
    < 10*vscale(f)*eps);

%%
% Verify operation for array-valued chebtech objects.

f = testclass.make(@(x) [3*ones(size(x)), 1+cos(pi*x), ...
    1 + exp(2*1i*pi*x) + exp(-1i*pi*x)], [], pref);
p = poly(f);
p_exact = [0 0   0;...
           0 0.5 1;...
           3 1   1;...
           0 0.5 0;...
           0 0   1].';
pass(5) = (norm(p(:) - p_exact(:), inf) < 10*max(vscale(f)*eps));

end


text:
%Script cmssnn;
%Center of mass calculation
%xc, yc; xcc, ycc - coordinates of center
%of mass of curve and spot
%(interior of the curve) correspondingly
%mss - mass (square) of a spot
mss=sum(sum(~J3));
xcc=0;ycc=0;
  for i=1:si;
    for j=1:sj;
      if ~J3(i,j)
        xcc=xcc+i;
        ycc=ycc+j;
      end
    end
  end
    xcc=xcc/mss;
    ycc=ycc/mss;

xc=mean(X0);
yc=mean(Y0);

    kz=[];
for iz=1:numel(X0);
  kz=[kz max(sqrt((X0(iz)-X0(iz+1:end)).^2+(Y0(iz)-Y0(iz+1:end)).^2))];
end
[uz1,uz]=max(kz);
[kz1,kz2]=max(sqrt((X0(uz)-X0(1:end)).^2+(Y0(uz)-Y0(1:end)).^2));

author: Sable

repo: mcbench-benchmarks


text:
%ADDCIRCLE  Add a circle to the current plot
%
%   addcircle(center, radius)
%   addcircle(center, radius, linestyle)
%
%  Returns the graphics handle for the circle.

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

function h = addcircle(center, radius, varargin)

	n = 100;

	th = [0:n]'/n*2*pi;
	cth = cos(th);
	sth = sin(th);

    ih = ishold;
    hold on
    handles = [];   % list of handles for circles

    for i=1:numrows(center)
        handles = [handles; patch(radius(i)*cos(th)+center(i,1), radius(i)*sin(th)+center(i,2), varargin{:})];
    end
    if ih == 0,
        hold off
    end

    if nargout > 0,
        h = handles;
    end

author: petercorke

repo: machinevision-toolbox-matlab