text:
%--------------------------------------------------------------------------
% HDL Workflow Script
% Generated with MATLAB 9.3 (R2017b) at 20:13:12 on 23/01/2018
% This script was generated using the following parameter values:
%     Filename  : '/tmp/modem-phy/hdlworkflow_hwi.m'
%     Overwrite : true
%     Comments  : true
%     Headers   : true
%     DUT       : 'combinedTxRx_HWInterface/Combined TX and RX'
% To view changes after modifying the workflow, run the following command:
% >> hWC.export('DUT','combinedTxRx_HWInterface/Combined TX and RX');
%--------------------------------------------------------------------------

%% Load the Model
load_system('combinedTxRx_HWInterface');

%% Restore the Model to default HDL parameters
%hdlrestoreparams('combinedTxRx_HWInterface/Combined TX and RX');

%% Model HDL Parameters
%% Set Model 'combinedTxRx_HWInterface' HDL parameters
hdlset_param('combinedTxRx_HWInterface', 'HDLSubsystem', 'combinedTxRx_HWInterface/Combined TX and RX');
hdlset_param('combinedTxRx_HWInterface', 'ReferenceDesign', 'Receive and transmit path');
hdlset_param('combinedTxRx_HWInterface', 'ReferenceDesignParameter', {'ChannelMapping','1','DUTSynthFreqMHz','20'});
hdlset_param('combinedTxRx_HWInterface', 'SynthesisTool', 'Xilinx Vivado');
hdlset_param('combinedTxRx_HWInterface', 'SynthesisToolChipFamily', 'Zynq');
hdlset_param('combinedTxRx_HWInterface', 'SynthesisToolDeviceName', 'xc7z045');
hdlset_param('combinedTxRx_HWInterface', 'SynthesisToolPackageName', 'ffg900');
hdlset_param('combinedTxRx_HWInterface', 'SynthesisToolSpeedValue', '-2');
hdlset_param('combinedTxRx_HWInterface', 'TargetDirectory', 'hdl_prj_hwi/hdlsrc');
hdlset_param('combinedTxRx_HWInterface', 'TargetPlatform', 'ZC706 and FMCOMMS2/3/4');
hdlset_param('combinedTxRx_HWInterface', 'Workflow', 'IP Core Generation');

% Set SubSystem HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX', 'ProcessorFPGASynchronization', 'Free running');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ValidInTx', 'IOInterface', 'Tx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ValidInTx', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/rInTx', 'IOInterface', 'Tx data I1 In [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/rInTx', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/iInTx', 'IOInterface', 'Tx data Q1 In [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/iInTx', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromRadioR', 'IOInterface', 'Rx data I1 In [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromRadioR', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromRadioI', 'IOInterface', 'Rx data Q1 In [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromRadioI', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromRadioValidIn', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromRadioValidIn', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/LoopBW', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/LoopBW', 'IOInterfaceMapping', '');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/EQmu', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/EQmu', 'IOInterfaceMapping', '');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/SelectScope', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/SelectScope', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/IQ Interface Mapper/real', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/IQ Interface Mapper/real', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/IQ Interface Mapper/imag', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/IQ Interface Mapper/imag', 'IOInterfaceMapping', '');

% Set SubSystem HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL', 'ProcessorFPGASynchronization', 'Free running');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/re', 'IOInterface', 'Rx data I1 In [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/re', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/Im', 'IOInterface', 'Rx data Q1 In [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/Im', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/Enable', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/Enable', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/FRLoopBW', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/FRLoopBW', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/EQmu', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/EQmu', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/Scope Select', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/Scope Select', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/DebugSelector', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/DebugSelector', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/bytesOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/bytesOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/validOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/validOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/sync', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/sync', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/payloadLenOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/payloadLenOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/dataRe', 'IOInterface', 'Rx data I1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/dataRe', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/dataIm', 'IOInterface', 'Rx data Q1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/dataIm', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/validIQ', 'IOInterface', 'Rx data Valid Out');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/validIQ', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/crcErr', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/crcErr', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/validCRC', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/validCRC', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/endFlag', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/endFlag', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/DebugSelection', 'IOInterface', 'Rx data Valid Out');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Receiver HDL/DebugSelection', 'IOInterfaceMapping', '[0]');

% Set SubSystem HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL', 'ProcessorFPGASynchronization', 'Free running');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/bytesIn', 'IOInterface', 'Tx data I1 In [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/bytesIn', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/validIn', 'IOInterface', 'Tx data Valid In');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/validIn', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/Packet Generation/trueData', 'IOInterface', 'Tx data Valid Out');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/Packet Generation/trueData', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/real', 'IOInterface', 'Tx data I1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/real', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/imag', 'IOInterface', 'Tx data Q1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/imag', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/validOut', 'IOInterface', 'Tx data Valid Out');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/validOut', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/Need Data', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/Transmitter HDL/Need Data', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ToRadioR', 'IOInterface', 'Tx data I1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ToRadioR', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ToRadioI', 'IOInterface', 'Tx data Q1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ToRadioI', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ToRadioValid', 'IOInterface', 'Tx data Valid Out');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/ToRadioValid', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromIPR', 'IOInterface', 'Rx data I1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromIPR', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromIPI', 'IOInterface', 'Rx data Q1 Out [0:15]');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromIPI', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromIPvalid', 'IOInterface', 'Rx data Valid Out');
hdlset_param('combinedTxRx_HWInterface/Combined TX and RX/FromIPvalid', 'IOInterfaceMapping', '[0]');

%% Workflow Configuration Settings
% Construct the Workflow Configuration Object with default settings
hWC = hdlcoder.WorkflowConfig('SynthesisTool','Xilinx Vivado','TargetWorkflow','IP Core Generation');

% Specify the top level project directory
hWC.ProjectFolder = 'hdl_prj_hwi';
hWC.ReferenceDesignToolVersion = '2017.4';
hWC.IgnoreToolVersionMismatch = false;

% Set Workflow tasks to run
hWC.RunTaskGenerateRTLCodeAndIPCore = true;
hWC.RunTaskCreateProject = true;
hWC.RunTaskGenerateSoftwareInterfaceModel = false;
hWC.RunTaskBuildFPGABitstream = true;
hWC.RunTaskProgramTargetDevice = false;

% Set properties related to 'RunTaskGenerateRTLCodeAndIPCore' Task
hWC.IPCoreRepository = '';
hWC.GenerateIPCoreReport = true;

% Set properties related to 'RunTaskCreateProject' Task
hWC.Objective = hdlcoder.Objective.None;
hWC.AdditionalProjectCreationTclFiles = '';
hWC.EnableIPCaching = false;

% Set properties related to 'RunTaskGenerateSoftwareInterfaceModel' Task
hWC.OperatingSystem = '';

% Set properties related to 'RunTaskBuildFPGABitstream' Task
hWC.RunExternalBuild = true;
hWC.TclFileForSynthesisBuild = hdlcoder.BuildOption.Default;
hWC.CustomBuildTclFile = '';

% Set properties related to 'RunTaskProgramTargetDevice' Task
hWC.ProgrammingMethod = hdlcoder.ProgrammingMethod.Download;

% Validate the Workflow Configuration Object
hWC.validate;

%% Run the workflow
hdlcoder.runWorkflow('combinedTxRx_HWInterface/Combined TX and RX', hWC);

author: analogdevicesinc


text:
function img = imresizecrop(img, M, METHOD)
%
% img = imresizecrop(img, M, METHOD);
%
% Output an image of size M(1) x M(2).

if nargin < 3
    METHOD = 'bilinear';
end

if length(M) == 1
    M = [M(1) M(1)];
end

scaling = max([M(1)/size(img,1) M(2)/size(img,2)]);

%scaling = M/min([size(img,1) size(img,2)]);

newsize = round([size(img,1) size(img,2)]*scaling);
img = imresize(img, newsize, METHOD);

[nr nc cc] = size(img);

sr = floor((nr-M(1))/2);
sc = floor((nc-M(2))/2);

img = img(sr+1:sr+M(1), sc+1:sc+M(2),:);

author: CSAILVision

repo: LabelMeToolbox


text:
function str = gencode_substructcode(subs, name)

% GENCODE_SUBSTRUCTCODE  Create code for a subscript structure
% Generate MATLAB code (using SUBSTRUCT) to create subscript structure
% subs. See help on SUBSTRUCT, SUBSASGN and SUBSREF for details how
% subscript structures are used.
%
% str = gencode_substructcode(subs, name)
% Input arguments:
%  subs - a subscript structure
%  name - optional: name of variable
% Output arguments:
%  str  - a one-line cellstr containing a call to SUBSTRUCT that returns
%         an substruct equivalent to subs.
% If name is supplied as input argument, the generated code will assign
% the output of SUBSTRUCT to the variable 'name'.
% then only the rhs of the expression will be returned.
% For '()' and '{}' also pseudo subscripts are allowed: if subs.subs{...}
% is a string, it will be printed literally, even if it is not equal to
% ':'. This way, one can create code snippets that contain e.g. references
% to a loop variable by name.
%
% See also GENCODE, GENCODE_RVALUE, GENCODE_SUBSTRUCT.
%
% This code has been developed as part of a batch job configuration
% system for MATLAB. See
%      http://sourceforge.net/projects/matlabbatch
% for details about the original project.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: gencode_substructcode.m 3355 2009-09-04 09:37:35Z volkmar $

rev = '$Rev: 3355 $'; %#ok

ind = 1;
if nargin < 2
    name = '';
end

if ~isstruct(subs) || ~all(isfield(subs, {'type','subs'}))
    if any(exist('cfg_message') == 2:6)
        cfg_message('matlabbatch:usage', 'Item is not a substruct.');
    else
        warning('gencode_substructcode:usage', ...
                'Item is not a substruct.');
    end
else
    if isempty(subs)
        str = {'struct(''type'',{},''subs'',{})'};
    else
        str = {'substruct('};
        for k = 1:numel(subs)
            str{1} = sprintf('%s''%s'',', str{1}, subs(k).type);
            switch subs(k).type
                case '.',
                    substr = sprintf('''%s''', subs(k).subs);
                case {'()','{}'},
                    substr = '{';
                    for l = 1:numel(subs(k).subs)
                        if ischar(subs(k).subs{l})
                            if strcmp(subs(k).subs{l},':')
                                substr = sprintf('%s'':'', ', substr);
                            else
                                substr = sprintf('%s%s, ', substr, subs(k).subs{l});
                            end
                        else
                            substr1 = sprintf('%d ', subs(k).subs{l});
                            if numel(subs(k).subs{l}) > 1
                                substr1 = sprintf('[%s]', substr1(1:end-1));
                            else
                                substr1 = substr1(1:end-1);
                            end
                            substr = sprintf('%s%s, ', substr, substr1);
                        end
                    end
                    substr = sprintf('%s}', substr(1:end-2));
            end
            str{1} = sprintf('%s%s, ', str{1}, substr);
        end
        str{1} = sprintf('%s)', str{1}(1:end-2));
    end
    if ~isempty(name)
        str{1} = sprintf('%s = %s;', name, str{1});
    end
end

author: spm


text:
function A = spm_mesh_adjacency(F)
% Compute the adjacency matrix of a triangle mesh
% FORMAT A = spm_mesh_adjacency(F)
% F        - a [fx3] faces array or a patch structure
%
% A        - adjacency matrix as a sparse [vxv] array
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_adjacency.m 4035 2010-08-05 18:54:32Z guillaume $

if ~isnumeric(F) && isfield(F,'vertices')
    N = size(F.vertices,1);
    F = double(F.faces);
    A = sparse([F(:,1); F(:,1); F(:,2); F(:,2); F(:,3); F(:,3)], ...
           [F(:,2); F(:,3); F(:,1); F(:,3); F(:,1); F(:,2)], 1, N, N);
else
    if isstruct(F), F = F.faces; end
    F = double(F);
    A = sparse([F(:,1); F(:,1); F(:,2); F(:,2); F(:,3); F(:,3)], ...
           [F(:,2); F(:,3); F(:,1); F(:,3); F(:,1); F(:,2)], 1);
end

A = double(A > 0);


text:
function r = polya_sample(a,n)
% POLYA_SAMPLE      Sample from Dirichlet-multinomial (Polya) distribution.
% POLYA_SAMPLE(a,n) returns a matrix of histograms.
% If A is a row, the histograms are the rows, others they are the columns.
% N is a vector whose length is the number of histograms.
% N(i) will be the total count in histogram i.

row = (rows(a) == 1);

a = a(:);
p = dirichlet_sample(a,length(n));
r = zeros(size(p));
for i = 1:length(n)
  r(:,i) = sample_hist(p(:,i),n(i));
end
if row
  r = r';
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
classdef MaterialDesignExperiment < handle

   properties (Access = public)

    end

    properties (Access = private)

    end

    properties (Access = private)
        fileNames
        topOptSet
        topOptProblem
    end

    methods (Access = public)

        function obj = MaterialDesignExperiment()
            obj.init();
            for icases = 1:numel(obj.fileNames)
                obj.createSettings(icases);
                obj.solveProblem();
                close all
            end

        end

    end

    methods (Access = private)

        function init(obj)
            obj.fileNames = {...
                             'HorizontalMaterialDesign';
                          %   'CompositeMaterialDesignTriDensityP1';
                           %  'CompositeMaterialDesignTriDensityPDE';
                           %  'CompositeMaterialDesignQuadDensityP1';
                            %  'CompositeMaterialDesignQuadDensityPDE';
                            % 'CompositeMaterialDesignTriLevelSetP1';
                            % 'CompositeMaterialDesignTriLevelSetPDE';
                           % 'CompositeMaterialDesignQuadLevelSetP1';
                           %  'CompositeMaterialDesignQuadLevelSetPDE';
                             };
        end

        function createSettings(obj,icases)
            s = SettingsTopOptProblem(obj.fileNames{icases});
            obj.topOptSet = s;
        end

        function solveProblem(obj)
            obj.topOptProblem = TopOpt_Problem(obj.topOptSet);
            obj.topOptProblem.computeVariables();
            obj.topOptProblem.postProcess();
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: get_rate_diff_ci1.m
% one-sided confidence intervals for the difference of rates

% 030113 tdr created

function ci = get_rate_diff_ci1(x1,A1,x2,A2,alpha,method,verbose)

r1_hat = x1/A1; r2_hat = x2/A2; delta_r_hat = r1_hat - r2_hat;
a = find_lower_limit(x1,A1,x2,A2, alpha);
b = find_upper_limit(x1,A1,x2,A2, alpha);
ci = [delta_r_hat a b];

% ------------------------------------------------------------
% find lower limit

function limit = find_lower_limit(x1,A1,x2,A2,alpha)
% alpha is 1-Pr{delta_r >= a}
r1_hat = x1/A1; r2_hat = x2/A2; delta_r_hat = r1_hat - r2_hat;

% set limits of search
too_big = 1e12; too_small = -1e12;

%check that too_big is indeed too big etc.
alpha_too_big = 1 - rate_diff(x1,A1,x2,A2,too_big);
alpha_too_small = 1 - rate_diff(x1,A1,x2,A2,too_small);
if (alpha_too_big < alpha | alpha_too_small > alpha), error('Failure to set bounds for search'); end;

tolerance = 1e-6;
max_count = 50;
count = 0;

% initial guess for search
guess = (too_big + too_small)/2;
alpha_guess = 1- rate_diff(x1,A1,x2,A2,guess);

% binary search
while (abs(alpha_guess - alpha) > tolerance & (count < max_count))
    if (alpha_guess > alpha)
        too_big = guess;
        guess = (too_small + guess)/2;
    else
        too_small=guess;
        guess = (too_big + guess)/2;
    end;
    count = count+1;
    alpha_guess = 1 - rate_diff(x1,A1,x2,A2,guess);
end;

limit = guess;

% ------------------------------------------------------------
% find upper limit

function limit = find_upper_limit(x1,A1,x2,A2,alpha, bounds)
% alpha is Pr{delta_r >= b}
r1_hat = x1/A1; r2_hat = x2/A2; delta_r_hat = r1_hat - r2_hat;

% set limits for search
too_big = 1e12;  too_small = -1e12;

%check that too_big is indeed too big etc.
alpha_too_big = rate_diff(x1,A1,x2,A2,too_big);
alpha_too_small = rate_diff(x1,A1,x2,A2,too_small);
if (alpha_too_big > alpha | alpha_too_small < alpha), error('Failure to set bounds for search'); end;

tolerance = 1e-6;
max_count = 50;
count = 0;

% initial guess for search
guess = (too_big + too_small)/2;
alpha_guess = rate_diff(x1,A1,x2,A2,guess);

% binary search
while (abs(alpha_guess - alpha) > tolerance & (count < max_count))
    if (alpha_guess < alpha)
        too_big = guess;
        guess = (too_small + guess)/2;
    else
        too_small=guess;
        guess = (too_big + guess)/2;
    end;
    count = count+1;
    alpha_guess = rate_diff(x1,A1,x2,A2,guess);
end;

limit = guess;

author: Sable

repo: mcbench-benchmarks


text:
im = imgRead(101085,'gray');
fb = fbCreate(6,1,1,3);
ntex = 32;
[tmap,tex] = computeTextons(fbRun(fb,im),ntex);
[tim,tperm] = visTextons(tex,fb);
wt = zeros(ntex,1);
for i = 1:ntex,
  wt(i) = sum(abs(tim{i}(:))); % L1 norm of texton
end
wt = wt / max(wt(:));
tsim = zeros(ntex);
for i = 1:ntex,
  for j = 1:ntex,
    tsim(i,j) = sum(sum(abs(tim{i}-tim{j})));
  end
end
r = 10;
norient = 6;
tic; [tg,theta] = tgmo(tmap,ntex,r,norient,tsim); toc;
aa = cell(size(tg));
bb = cell(size(tg));
cc = cell(size(tg));
for i = 1:numel(tg),
  tic; [c,b,a] = fitparab(tg{i},r,theta(i)); toc;
  aa{i}=a; bb{i}=b; cc{i}=c;
end
tgs = cell(size(tg));
pb = zeros(size(tmap));
for i = 1:numel(tgs),
  tgs{i} = max(0,cc{i}) .* (aa{i}<0) .* exp(-abs(bb{i})/0.1);
  pb = max(pb,tgs{i});
end
pb2 = zeros(size(tmap));
for i = 1:numel(tgs),
  pb2 = max(pb2,(tgs{i}==pb).*nonmax(tgs{i},theta(i)));
end

figure(1); clf;
imshow(im);

figure(2); clf;
imagesc(mymontage({tim{tperm}}));
axis image; colorbar;

figure(3); clf;
imagesc(tmap);
truesize;

figure(4); clf;
imagesc(pb);
truesize;

figure(5); clf;
imagesc(pb2)
truesize;

author: zouchuhang

repo: LayoutNet


text:
% OUTER_HULL Compute the "outer hull" of a potentially non-manifold mesh (V,F)
% whose intersections have been "resolved" (e.g. using `cork` or
% `igl::selfintersect`). The outer hull is defined to be all facets (regardless
% of orientation) for which there exists some path from infinity to the face
% without intersecting any other facets. For solids, this is the surface of the
% solid. In general this includes any thin "wings" or "flaps".  This
% implementation largely follows Section 3.6 of "Direct repair of
% self-intersecting meshes" [Attene 2014].
%
% [G,J,flip] = outer_hull(V,F);
%
% Inputs:
%   V  #V by 3 list of vertex positions
%   F  #F by 3 list of triangle indices into V
% Outputs:
%   G  #G by 3 list of output triangle indices into V
%   J  #G list of indices into F
%   flip  #F list of whether facet was added to G **and** flipped orientation
%     (false for faces not added to G)

author: alecjacobson

repo: gptoolbox


text:
function tauc = taucv(alfa,f,n)
%TAUCV  Upper percentage point of the tau distribution.
%       tauc=taucv(alfa,f,n) computes the critical value of tau distribution
%       for the Type I error--significance level (alfa), degree of freedom (f)
%       and the number of observations (n). It is used in Pope test for outliers.
%
%       References:
%       [1] Koch KR,1999, Parameter estimation and hypothesis testing in linear models,
%       Springer Verlag, Berlin
%       [2] Pope AJ,1976, The  statistics  of residuals and the detection of outliers.
%       NOAA Technical Report NOS65 NGS1, US Department of Commerce, National Geodetic
%       Survey,Rockville, Maryland
%       -----------------------------------------------------------------------------
%       Written by Cuneyt AYDIN, 2003, Yildiz Technical University, Geodesy Division
%       e-mail:caydin@yildiz.edu.tr
%       -----------------------------------------------------------------------------

if nargin < 3,
    error('requires at least three input arguments; alfa,f and n');
 end

if f==1
   error('degree of freedom must be larger than 1 (f>1)')

end

fdist=finv((1-alfa/n),1,(f-1));
tauc=sqrt((f*fdist)/(f-1+fdist));


text:
function g = fix(f)
%FIX   SINGFUN does not support FIX.

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

error('CHEBFUN:SINGFUN:fix:notSupported', 'SINGFUN does not support FIX.');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function bound = compute_shape_boundary(M)

% compute_shape_boundary - extract boundary points
%
% bound = compute_shape_boundary(M);
%
%   bound is oriented counter clockwise.
%   This is the 8 connectivity boundary.
%
%   Copyright (c) 2007 Gabriel Peyre

n = size(M,1);

% enforce only 1 connected component
epsilon = 1e-10;
options.nb_iter_max = Inf;
options.Tmax = 2*n;
W = epsilon+M;
I = find(M(:)==1); x = linspace(-1,1,n);
[Y,X] = meshgrid(x,x);
[tmp,c] = min( X(I).^2 + Y(I).^2 );
[a b] = ind2sub([n n],I(c));
options.constraint_map = [];
[d,tmp] = perform_fast_marching(W, [a;b], options);
I = find(d>1e5 | d<0); M(I) = 0;

% points on the boundary
if 0
    % outside the shape
    h = [0 1 0; 1 1 1; 0 1 0];
    Mh = perform_convolution(M,h);
    A = Mh>0 & M==0;
else
    % inside the shape
    h = ones(3);
    Mh = perform_convolution(M,h);
    A = Mh>0 & Mh<9 & M==1;
end

% 8 connectivity
wx = [1 -1 0  0 -1 1 -1  1];
wy = [0  0 1 -1 -1 1  1 -1];

% order the points
a = find(A(:)==1); a = a(1);
[x,y] = ind2sub(size(A),a);

boundx = x;
boundy = y;
A(x,y) = 0;

while sum(A(:))>0
    % find closest point
    for k=1:8
        if x+wx(k)>0 && x+wx(k)<=n && y+wy(k)>0 && y+wy(k)<=n && A(x+wx(k),y+wy(k))>0
            x = x+wx(k);
            y = y+wy(k);
            break;
        end
    end
    if A(x,y)==0
        % try to find the existing point the closest to the last one
        I = find(A>0);
        [xlist,ylist] = ind2sub(size(A),I);
        d = (xlist-x).^2 + (ylist-y).^2;
        [tmp,I] = min(d(:));
        x = xlist(I(1)); y = ylist(I(1));
    end
    A(x,y) = A(x,y)-1;
    boundx(end+1) = x;
    boundy(end+1) = y;
end
bound = cat(2,boundx(:), boundy(:));

% reorient the curve so that it is clockwise oriented
p = size(bound,1);
c = mean(bound,1);
v = bound - repmat(c, [p 1]);
% cross product
z = v(1:end-1,1).*v(2:end,2) - v(1:end-1,2).*v(2:end,1);
if sum(z<0)>sum(z<0)
    % should flip the curve
    bound = bound(end:-1:1,:);
end


text:
function [tracks, revived, dropped_criteria, dropped_vicinity] = tolerance_track_lnn( eddies_path, type, time_frequency, tolerance, minimumArea ) % old args: eddies, time_frequency, tolerance
%TOLERANCE_TRACK_LNN Tracking eddies by lnn method with loss tolerance

%   ---- NOTE ----: This is not the regular version of track_lnn. This
%   version of track_lnn will interpolate fake eddies in-between two close
%   tracks. Our theory is that these nearby tracks are the same eddy, but
%   the eddy on the missing day just wasn't strong enough to be picked up
%   by the original track_lnn. We created this function to try to recover
%   some of the dropped eddies so that we could connect the broken tracks.

%   For each eddy in a timestep, all eddies in next timestep within a boundary are checked to see if there's an eddy
%   that qualifies some conditions to be stitched to current eddy. The north, south and east bounds are the gate
%   distance while the west bound is computed based on rossby-phase_speed.
% Input:
%   eddies_path: The path to the directory where eddies that were detected
%                by eddyscan were saved
%   type: Cyclone type of eddies ('anticyc' or 'cyclonic')
%   time_frequency: Number of days between timesteps. For weekly data, put
%                   7, for daily data, put 1
%   tolerance: Number of timesteps an eddy can be lost for
%   minimumArea: Minimum area of the eddies we want to include in the
%                tracking. I.E. eddy size is 4+, but we only want to
%                include eddies of size 9+ in the tracks. Set minimumArea
%                to 9
%
% Understanding the columns of the tracks
% Column 1 is the latitude of the centroid of the eddy
% Column 2 is the longitude of the centroid of the eddy
% Column 3 is the time slice (day) that the eddy is from. To use this data to
% get the actual day that the eddy is from, use this value as an index to a
% dates array, and the return value of your dates array will be the day the
% eddy is from.
% Column 4 is the index in the struct array of eddies for that particular
% day. If you get the correct eddies struct that corresponds to the date in
% column 3, the value in this column is the index in that struct where this
% actual eddy resides. It is through getting the appropriate struct and
% accessing the appropriate entry that we relate this track data to actual
% eddies in the struct.
% Column 5 is whether the eddy in the tracking data is a "fake" eddy or
% not. Normal eddies will be denoted by a '0' in the fifth column. Fake
% eddies that are generated by this tolerance_track_lnn function will have
% a '1' in the fifth column.

%% Setting constants
temp = load('rossby_daily_phase_speeds.mat');
rossby_phase_speed = temp.rossby_daily_phase_speeds;
rossby_phase_speed(:, 3) = rossby_phase_speed(:, 3) * time_frequency;

gate_distance = 150/7 * time_frequency; % 150km is default maximum travel distance for 1 week
gate_distance_in_degree = km2deg(gate_distance);

% Make lat and lon to be in [-180 180]
rossby_phase_speed(rossby_phase_speed(:, 2) > 180, 2) = rossby_phase_speed(rossby_phase_speed(:, 2) > 180, 2) - 360;
rossby_phase_speed(rossby_phase_speed(:, 2) < -180, 2) = rossby_phase_speed(rossby_phase_speed(:, 2) < -180, 2) + 360;
max_rossby_lat = max(rossby_phase_speed(:, 1));
min_rossby_lat = min(rossby_phase_speed(:, 1));

rangesearch_radius = 5; % euclidean radius in degree for rangesearch

% All next eddies with longitude in range [min_lon (min_lon + pad_range)] and [(max_lon - pad_range) max_lon] will be
% duplicated with longitude lon+360 or lon-360 to make sure Euclidean distances can still get the eddies on the other
% side of the map
pad_range = 5;
revived = 0;
dropped_criteria = 0;
dropped_vicinity = 0;

eddies_names = get_eddies_names(eddies_path, type);
eddy_length = length(eddies_names);
fake_eddy_count = ones(eddy_length, 1);
eddies = cell(eddy_length, 1);
indices = cell(eddy_length, 1);
orig_count = zeros(eddy_length, 1);
[eddies, indices, orig_count] = load_eddies_cell(eddies, indices, orig_count, eddies_names, minimumArea);

%% Begin tracking
tracks = {};

curr_track_indexes = nan(length(eddies{1}), 1); % track index of first eddies

for curr_time_step = 1:(length(eddies) - 1)
    curr_eddies_coordinates = [[eddies{curr_time_step}.Lat]' [eddies{curr_time_step}.Lon]'];
    next_eddies_coordinates = [[eddies{curr_time_step + 1}.Lat]' [eddies{curr_time_step + 1}.Lon]'];
    orig_next_eddies_coordinates_length = size(next_eddies_coordinates, 1);

    curr_real_indices = indices{curr_time_step};
    next_real_indices = indices{curr_time_step + 1};
    % eddy entries for connecting tracks.

    disp(curr_time_step)

    % The availability of next eddies, will be updated when an in-ranged eddy matches current eddy and the tracks are updated
    next_eddies_availabilities = true(size(next_eddies_coordinates, 1), 1);
    % Index of the tracks of the tracked next eddies, NaN for untracked eddies
    next_eddies_track_indexes = NaN(size(next_eddies_coordinates, 1), 1);

    % First, pad the next eddies so that euclidean distances don't mess up eddy distances at the edge of world map
    [padding_eddy_index, padding_eddies] = pad_eddies(next_eddies_coordinates, pad_range);
    padded_next_eddies = [next_eddies_coordinates; padding_eddies];

    padding_eddies_length = size(padding_eddies, 1);

    % Get in-ranged next eddies
    [next_padded_eddy_in_range_indexes, ~] = ...
        rangesearch(padded_next_eddies(:, 1:2), curr_eddies_coordinates(:, 1:2), rangesearch_radius);

    %next_padded_eddy_in_range_indexes =

    % Get distances from current eddies to in-ranged next eddies
    [dists, first_curr_eddy_index_in_distance_data] = get_eddy_distances(curr_eddies_coordinates, padded_next_eddies, ...
        next_padded_eddy_in_range_indexes);

    % Get westbound of curr eddies
    westbounds = get_west_bounds(curr_eddies_coordinates(:, 1:2), gate_distance, rossby_phase_speed, ...
        max_rossby_lat, min_rossby_lat);
    % For each current eddy, check if any in-ranged next eddy matches and update tracks
    for curr_eddy_index = 1:size(curr_eddies_coordinates, 1)
        % Index of the distance between current eddy and next eddy in the dists array
        curr_distance_index = first_curr_eddy_index_in_distance_data(curr_eddy_index);
        eddies_within_range = 0;
        matched = false;

        if curr_eddy_index <= length(curr_real_indices)
            curr_real_index = curr_real_indices(curr_eddy_index);
        else
            curr_real_index = orig_count(curr_time_step) + (curr_eddy_index - length(curr_real_indices));
        end

        for next_padded_eddy_index = next_padded_eddy_in_range_indexes{curr_eddy_index}
            % Get real next eddy index
            if next_padded_eddy_index > orig_next_eddies_coordinates_length && next_padded_eddy_index < orig_next_eddies_coordinates_length + ...
                    padding_eddies_length + 1
                % This is index of a padded eddy
                next_eddy_index = padding_eddy_index(next_padded_eddy_index - orig_next_eddies_coordinates_length);
            else
                next_eddy_index = next_padded_eddy_index;
            end

            % Check if the next eddy matches the current eddy
            if next_eddies_availabilities(next_eddy_index)
                % Only check available next eddies
                curr_eddy_lat = curr_eddies_coordinates(curr_eddy_index, 1);
                curr_eddy_lon = curr_eddies_coordinates(curr_eddy_index, 2);
                next_eddy_lat = next_eddies_coordinates(next_eddy_index, 1);
                next_eddy_lon = next_eddies_coordinates(next_eddy_index, 2);

                next_real_index = next_real_indices(next_eddy_index);

                if ~is_within_range(curr_eddy_lat, curr_eddy_lon, next_eddy_lat, next_eddy_lon, ...
                        dists(curr_distance_index), westbounds(curr_eddy_index), gate_distance, gate_distance_in_degree)
                    continue;
                else
                    eddies_within_range = eddies_within_range + 1;
                end

                if is_matched(eddies{curr_time_step}(curr_eddy_index), eddies{curr_time_step + 1}(next_eddy_index));

                    % Update tracks and get out of the loop
                    if isnan(curr_track_indexes(curr_eddy_index))
                        % new track
                        track_id = length(tracks) + 1;
                        tracks{track_id}(1, :) = ...
                            [curr_eddy_lat curr_eddy_lon curr_time_step curr_real_index];
                        tracks{track_id}(2, :) = ...
                            [next_eddy_lat next_eddy_lon (curr_time_step + 1) next_real_index];
                    else
                        % old track
                        track_id = curr_track_indexes(curr_eddy_index);
                        %check if track is flagged
                        if length(tracks{track_id}(end,:)) == 5
                            if tracks{track_id}(end,5) ~= 0
                                revived = revived + 1;
                                %need to interpolate new locations for flagged
                                %track entries
                                for i = 1:size((tracks{track_id}), 1)
                                    if tracks{track_id}(i,5) ~= 0
                                        flag = tracks{track_id}(i,5);
                                        break
                                    end
                                end

                                number_of_flags_at_end = 0;
                                for j = size(tracks{track_id},1):1
                                    if tracks{track_id}(j,5) ~= 0
                                        number_of_flags_at_end = number_of_flags_at_end + 1;
                                    else
                                        break;
                                    end
                                end

                                last_real_track_index = size(tracks{track_id},1) - number_of_flags_at_end;

                                lat_increment = (next_eddy_lat - tracks{track_id}(last_real_track_index, 1))/(number_of_flags_at_end + 1);
                                lon_increment = (next_eddy_lon - tracks{track_id}(last_real_track_index, 2))/(number_of_flags_at_end + 1);
                                %now just increment the flagged tracks' lat and lon by these amounts
                                for k = last_real_track_index + 1: size(tracks{track_id},1)
                                    tracks{track_id}(k,1) = tracks{track_id}(last_real_track_index,1) + (k - last_real_track_index)*lat_increment;
                                    tracks{track_id}(k,2) = tracks{track_id}(last_real_track_index,2) + (k - last_real_track_index)*lon_increment;
                                end
                            end
                            tracks{track_id}(end + 1, :) = ...
                                [next_eddy_lat next_eddy_lon (curr_time_step + 1) next_real_index 0 ];
                        else
                            tracks{track_id}(end + 1, :) = ...
                                [next_eddy_lat next_eddy_lon (curr_time_step + 1) next_real_index];
                        end
                    end
                    next_eddies_track_indexes(next_eddy_index) = track_id;
                    next_eddies_availabilities(next_eddy_index) = false;
                    matched = true;
                    break;
                end
                % getting here means that it had at least one eddy within
                % range but it wasn't a match for criteria

            end

            curr_distance_index = curr_distance_index + 1;

        end

        if ~matched
            if ~isnan(curr_track_indexes(curr_eddy_index))
                % getting here means there were no matches
                if ~eddies_within_range
                    dropped_vicinity = dropped_vicinity + 1;
                else
                    %there was at least one within range but none were matches
                    dropped_criteria = dropped_criteria + 1;
                end
            end

            if tolerance > 0
                %if a track already exists, add flag for potential
                %missing track.
                %This section checks if it has the max tolerated
                %number of flagged tracks already
                if ~isnan(curr_track_indexes(curr_eddy_index))
                    track_id = curr_track_indexes(curr_eddy_index);

                    if length(tracks{track_id}(end,:)) == 5  % protect against out of bounds indexing (too large)
                        if size(tracks{track_id}, 1) > tolerance %protect against negative indexing
                            if tracks{track_id}(end - tolerance + 1, 5) ~= 0
                                %remove flagged tracks

                                tracks{track_id} = tracks{track_id}(1: (end - tolerance),:);
                                continue;
                            end
                        end
                    else
                        %track isn't flagged at all, so add 5th column in order to add flags
                        tracks{track_id}(:,5) = 0;
                    end
                    %won't get here if it got trimmed due to the
                    %continue statement above
                    %now make entry to next coordinates and eddies arrays so we can
                    %search from this fake eddy in next timestep
                    lat_increment = tracks{track_id}(end, 1) - tracks{track_id}(end - 1, 1);
                    lon_increment = tracks{track_id}(end, 2) - tracks{track_id}(end - 1, 2);
                    next_eddies_coordinates(end+1,:) = [curr_eddies_coordinates(curr_eddy_index, 1) + lat_increment ...
                        curr_eddies_coordinates(curr_eddy_index, 2) + lon_increment];
                    %now add the stats of the current eddy to the
                    %end of the next time steps eddies, but alter
                    %the coordinates
                    eddies{curr_time_step + 1}(end+1) = eddies{curr_time_step}(curr_eddy_index);
                    eddies{curr_time_step + 1}(end).Lat = curr_eddies_coordinates(curr_eddy_index,1) + lat_increment;%next_eddies_coordinates(end,1);
                    eddies{curr_time_step + 1}(end).Lon = curr_eddies_coordinates(curr_eddy_index,2) + lon_increment;%next_eddies_coordinates(end,2);
                    %make a flagged track entry with the above interpolated stats
                    if ~eddies_within_range
                        flag = -1; %dropped for vicinity
                    else
                        flag = 1; % dropped for criteria
                    end
                    next_size = orig_count(curr_time_step + 1) + fake_eddy_count(curr_time_step + 1);
                    fake_eddy_count(curr_time_step + 1) = fake_eddy_count(curr_time_step + 1) + 1;
                    tracks{track_id}(end + 1, :) = ...
                        [next_eddies_coordinates(end, 1) next_eddies_coordinates(end, 2)...
                        (curr_time_step + 1) next_size flag];
                    next_eddies_track_indexes(length(eddies{curr_time_step + 1})) = track_id;
                    next_eddies_availabilities(size(next_eddies_coordinates), 1) = false;
                end
            end

        end
    end
    curr_track_indexes = next_eddies_track_indexes;
    eddies{curr_time_step} = [];
    [eddies, indices, orig_count] = load_eddies_cell(eddies, indices, orig_count, eddies_names, minimumArea);
end

for i = 1:length(tracks)
    track = tracks{i};
    [~,b] = size(track);
    if b < 5
        track(:, 5) = 0;
        tracks{i} = track;
    end
end

end

function within_range = is_within_range(curr_lat, curr_lon, next_lat, next_lon, km_distance, westbound, ...
    gate_distance, gate_distance_in_degree)
% An approximate way of checking if the eddy in next timestep is within edddy in current timestep's range,
% following the description in D.B. Chelton et al. / Progress in Oceanography 91 (2011) page 208
% The north, east and south bounds are defined by the gate distance while the westbound is precomputed by rossby
% phase speed.
% Input:
%   curr_lat: latitude of the eddy in current timestep
%   curr_lon: longitude of the eddy in current timestep
%   next_lat: latitude of the eddy in next timestep
%   next_lon: longitude of the eddy in next timestep
%   km_distance: distance in km between the eddies in current and next timestep
%   gate_distance: default value of how far an eddy can travel in one timestep in km
%   gate_distance: default value of how far an eddy can travel in one timestep in degree

% Check north and south bound
if next_lat > curr_lat + gate_distance_in_degree || next_lat < curr_lat - gate_distance_in_degree
    within_range = false; % Avoid using km2deg or deg2km because
else
    % Check east and west bound
    if abs(next_lon - curr_lon) > 180
        if next_lon > curr_lon
            curr_lon = curr_lon + 360;
        else
            next_lon = next_lon + 360;
        end
    end
    if next_lon >= curr_lon
        is_on_east_side = true;
    else
        is_on_east_side = false;
    end

    if is_on_east_side
        within_range = km_distance <= gate_distance;
    else
        within_range = km_distance <= westbound;
    end

end

end

function westbounds = get_west_bounds(eddy_coordinates, gate_distance, phase_speed_data, max_rossby_lat, min_rossby_lat)
% Getting westbound of eddies by eddy coordinates and rossby phasespeed data
% West bound of an eddy is defined as the maximum value between 1.75 * rossby_phase_speed and default gate_distance
%   eddy_coordinates: coordinates of eddies in format [lat lon]
%   gate_distance: the default maximum distance an eddy can travel in a single timestep
%   phase_speed_data: rossby phase speed data in format [lat lon speed]
%   max_rossby_lat: to make sure eddies with latitude out of range get the default gate distance
%   min_rossby_lat: to make sure eddies with latitude out of range get the default gate distance

curr_eddy_lats = eddy_coordinates(:, 1);
curr_eddy_lons = eddy_coordinates(:, 2);
% Convert eddy longtidue format to -180 to 180 to make sure knn search work
curr_eddy_lons(curr_eddy_lons > 180) = curr_eddy_lons(curr_eddy_lons > 180) - 360;
curr_eddy_lons(curr_eddy_lons < -180) = curr_eddy_lons(curr_eddy_lons < -180) + 360;

nearest_ind = knnsearch(phase_speed_data(:, 1:2), [curr_eddy_lats curr_eddy_lons]);
westbounds = 1.75 * phase_speed_data(nearest_ind, 3);
westbounds(westbounds < gate_distance) = gate_distance;
% Make sure out of range eddies get the default value
westbounds(curr_eddy_lats > max_rossby_lat | curr_eddy_lats < min_rossby_lat) = gate_distance;

end

function [dists, first_curr_eddy_index_in_distance_data] = get_eddy_distances(curr_eddies, next_eddies, ...
    next_eddies_indexes)
% Get distances from current eddies to next eddies with indexes in the cell array next_eddies_indexes. This function
% computes the distance in one distance and deg2km to increase the tracking function performance
% Also return an array of indexes of the first appearance of a curr eddy in the distance array
% Input:
%   curr_eddies: eddies in current timestep with format [lat lon ...]
%   next_eddies: eddies in next timestep with format [lat lon ...]
%   next_eddies_indexes: a cell array of indexes of next eddies that will be computed the distances to current eddies

next_eddy_indexes_lengths = cellfun('length', next_eddies_indexes);
vectorized_next_eddy_indexes = [next_eddies_indexes{:}]';
vectorized_curr_eddy_indexes = zeros(sum(next_eddy_indexes_lengths), 1);
curr_index = 1;
for i = 1:length(next_eddy_indexes_lengths)
    vectorized_curr_eddy_indexes(curr_index:(curr_index + next_eddy_indexes_lengths(i) - 1)) = i;
    curr_index = curr_index + next_eddy_indexes_lengths(i);
end

dists = deg2km(distance(curr_eddies(vectorized_curr_eddy_indexes, 1:2), ...
    next_eddies(vectorized_next_eddy_indexes, 1:2)));

first_curr_eddy_index_in_distance_data = ones(length(next_eddy_indexes_lengths), 1);
for i = 2:length(next_eddy_indexes_lengths)
    first_curr_eddy_index_in_distance_data(i) = first_curr_eddy_index_in_distance_data(i - 1) + next_eddy_indexes_lengths(i - 1);
end

end

function [padding_eddy_index, padding_eddies] = pad_eddies(eddies, pad_range)
% Duplicate eddies in range [min_lon (min_lon + pad_range)] and [(max_lon - pad_range) max_lon] to make sure knnsearch
% work correctly to get in-range eddies in next timestep
%   eddies: eddies in format [lat lon ...]
%   pad_range: how far we want to pad (in degree)

min_lon = min(eddies(:, 2));
max_lon = max(eddies(:, 2));
padding_eddy_index = find( (eddies(:, 2) >= min_lon & eddies(:, 2) <= (min_lon + pad_range)) ...
    | (eddies(:, 2) >= (max_lon - pad_range) & eddies(:, 2) <= max_lon ) );

padding_eddies = eddies(padding_eddy_index, :);
padding_eddies(padding_eddies(:, 2) <= (min_lon + pad_range), 2) = ...
    padding_eddies(padding_eddies(:, 2) <= (min_lon + pad_range), 2) + 360;
padding_eddies(padding_eddies(:, 2) >= (max_lon - pad_range), 2) = ...
    padding_eddies(padding_eddies(:, 2) >= (max_lon - pad_range), 2) - 360;

end

function matched = is_matched(curr_eddy, next_eddy)
% An eddy is consider matched with current eddy if the amplitude and surface area is in range [0.25 2.75] *
% curr_eddy_amp/surface_area
%   curr_eddy: an eddy in current timestep in format [lat lon amp surface_area]
%   next_eddy: an eddy in next timestep in format [lat lon amp surface_area]

curr_amplitude = curr_eddy.Amplitude;
next_amplitude = next_eddy.Amplitude;
curr_surface_area = curr_eddy.SurfaceArea;
next_surface_area = next_eddy.SurfaceArea;
%if isempty(curr_amplitude) | isempty(next_amplitude) | isempty
%end
matched_amp = curr_amplitude > .25* next_amplitude && curr_amplitude < 2.75 * next_amplitude;
matched_area = curr_surface_area > .25 * next_surface_area && curr_surface_area < 2.75 * next_surface_area;

% matched = curr_amplitude > 0.25 * next_amplitude && curr_amplitude < 2.75 * next_amplitude ...
%     && curr_surface_area > 0.25 * next_surface_area && curr_surface_area < 2.75 * next_surface_area;
matched = matched_amp && matched_area;
end

function [eddies_names] = get_eddies_names(path, type)
% path is the path to the eddies directory
% type is anticyclonic or cyclonic
if ~strcmp(path(end), '/')
    path = strcat(path, '/');
end
files = dir(path);
x = 0;
for i = 1:length(files)
    if ~isempty(strfind(files(i).name, [type, '_']))
        x = x + 1;
    end
end
eddies_names = cell(x, 1);
x = 1;
for i = 1:length(files)
    if files(i).isdir && ~isequal(files(i).name, '.') && ~isequal(files(i).name, '..')
        rec_names = get_eddies_names([path, files(i).name, '/'], type);
        for j = 1:length(rec_names)
            eddies_names{x} = rec_names{j};
            x = x + 1;
        end
        continue;
    end
    file = files(i).name;
    [~, name, ext] = fileparts([path, file]);
    if ~isempty(strfind(name, [type, '_'])) && strcmp(ext, '.mat')
        eddies_names{x} = [path, file];
        x = x + 1;
    end
end
end

function [modified_eddies_cell, modified_indices_cell, modified_orig_count_array] = load_eddies_cell(eddies_cell, real_indices_cell, orig_count_array, eddies_names, minimumArea)
count_to_load = 10;
names_length = length(eddies_names);
x = 0;
pos_1 = 1;
for i = 1:length(eddies_cell)
    if ~isempty(eddies_cell{i})
        if x == 0
            pos_1 = i;
        end
        x = x + 1;
    end
end
iterations = count_to_load - x;
first_empty_pos = pos_1 + x;
last_empty_pos = first_empty_pos + iterations - 1;
if last_empty_pos > names_length
    last_empty_pos = names_length;
end
for i = first_empty_pos:last_empty_pos
    eddy_name = eddies_names{i};
    vars = load(eddy_name);
    names = fieldnames(vars);
    eddies = vars.(names{1});
    [f_eddies, real_indices, orig_count] = filter_eddies(eddies, minimumArea);
    eddies_cell{i} = f_eddies;
    real_indices_cell{i} = real_indices;
    orig_count_array(i) = orig_count;
end
modified_eddies_cell = eddies_cell;
modified_indices_cell = real_indices_cell;
modified_orig_count_array = orig_count_array;
end

function [filtered_eddies, real_indices, original_eddy_count] = filter_eddies(eddies, minimumArea)
stats = [eddies.Stats];
areas = [stats.Area];
x = 1;
original_eddy_count = length(eddies);
for i = 1:length(eddies)
    if areas(i) >= minimumArea
        filtered_eddies(x) = eddies(i);
        real_indices(x) = i;
        x = x + 1;
    end
end
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
% This script file load a data set using fscanf
% The default reads Northern California Hypoellipse Format
%

report_this_filefun(mfilename('fullpath'));
disp('Please make sure the has 88 characters for each line')
disp('and all blanks have been substituted by zeros')

% Lets make sure the file is closed...
safe_fclose(fid);

% reset paramteres
a = []; b = []; n = 0;

if inda == 1
    % initial selection option
    tmin   = 10.0001;
    tmax   = 98.000;
    lonmin = -180.0;
    lonmax =  180.0;
    latmin =  -90.00;
    latmax =  90.0 ;
    Mmin   = -4.0;
    Mmax   = 10.;
    mindep = -10;
    maxdep = 700;

    % call the pre-selection window
    presel
    return
end

% open the file and read 10000 lines at a time
[file1,path1] = uigetfile([ '*.dat'],' Earthquake Datafile');
if length(file1) >1
    fid = fopen([path1 file1],'r') ;;
else
    disp('Data import canceled'); return
end

while  ferror(fid) == ''
    n = n+1;
    % vari name   yr mo da hr mi se lat   la  lon    lo de ma1     ma he   hz
    % variabl #   1  2  3  4  5  6  7      8  9      10 11 12      13 14   15
    % position    2  4  6  8  10 14 16 17  21 24 25  29 34 36 67   69 84   88
    l = fscanf(fid,'%2d%2d%2d%2d%2d%4d%2d%*1c%4d%3d%*1c%4d%5f%2d%*31c%2d%*15c%4d',...
        [14 10000]) ;
    %if ferror(fid) ~= '' ; break; end

    b = [ -l(9,:)-l(10,:)/6000 ; l(7,:)+l(8,:)/6000 ; l(1,:);l(2,:);l(3,:);
        l(13,:)/10;l(11,:)/100;l(4,:);l(5,:); l(14,:)/100;l(12,:)/10];
    b = b';
    l =  b(:,6) >= Mmin & b(:,1) >= lonmin & b(:,1) <= lonmax & ...
        b(:,2) >= latmin & b(:,2) <= latmax & b(:,3) <= tmax  & ...
        b(:,3) >= tmin  & b(:,7) <= maxdep & b(:,7) >= mindep;
    a = [a ; b(l,:)];

    disp([ num2str(n*10000) ' earthquakes scanned, ' num2str(length(a)) ' EQ found'])
    if max(b(:,3)) >  tmax ; break; end

end
ferror(fid)
fclose(fid);

% Convert the third column into time in decimals
if length(a(1,:))== 7
    a.Date = decyear(a(:,3:5));
elseif length(a(1,:))>=9       %if catalog includes hr and minutes
    a.Date = decyear(a(:,[3:5 8 9]));
end

% save the data
[file1,path1] = uiputfile(fullfile(hodi, 'eq_data', '*.mat'), 'Save Earthquake Datafile');
sapa2 = ['save ' path1 file1 ' a'];
if length(file1) > 1; eval(sapa2);end

% call the map window
update(mainmap())

author: CelsoReyes


text:
% demo for binary data classification

close all
clear variables

dim.p = 32;
dim.n_t = 1;
dim.n_phi = 16;
dim.n_theta = 0;
dim.n = 0;

g_fname = @g_classif;
options.sources = struct('type',1 ,'out', 1); % one binomial observation;
options.priors.muPhi = zeros(dim.n_phi,1);
options.priors.SigmaPhi = 1e0*eye(dim.n_phi);
options.isYout = zeros(dim.p,1);
options.DisplayWin = 0;
options.verbose = 0;

Nmcmc = 64;
p = cell(2,2,Nmcmc);
o = cell(2,2,Nmcmc);
F = zeros(2,2,Nmcmc);
ner = zeros(2,2,Nmcmc); % proportion of correct predictions
mner = zeros(2,Nmcmc); % maximum performance rate

fprintf(1,['MCMC simulations... ']);
fprintf(1,'%6.2f %%',0)

for ii=1:Nmcmc
    options.inG.X = randn(dim.n_phi-1,dim.p);
    for i=1:2
        % simulate data with and without real mapping
        phi = (2-i)*randn(dim.n_phi,1);
        [y,x,x0,eta,e] = VBA_simulate (dim.n_t,[],g_fname,[],phi,[],[],[],options);
        g = y-e;
        g = g>0.5; % denoised data
        mner(i,ii) = sum(g.*y + (1-g).*(1-y))./dim.p; % max performance rate
        for j=1:2
            % invert model with and without the 2nd half of the data
            options.isYout(dim.p/2:dim.p) = 2-j;
            [p{j,i,ii},o{j,i,ii}] = VBA_NLStateSpaceModel(y,[],[],g_fname,dim,options);
            F(j,i,ii) = o{j,i,ii}.F - VBA_LMEH0(y,o{j,i,ii}.options);
            % proportion of correct predictions on 2nd half of the data
            gx = o{j,i,ii}.suffStat.gx(dim.p/2:dim.p);
            g0 = y(dim.p/2:dim.p);
            ner(j,i,ii) = sum(gx.*g0 + (1-gx).*(1-g0))./(dim.p/2);
        end

    end

    fprintf(1,repmat('\b',1,8))
    fprintf(1,'%6.2f %%',100*ii/Nmcmc)

end

fprintf(1,repmat('\b',1,8))
fprintf(' OK.')
fprintf('\n')

hf = figure('color',[1 1 1]);
pos = get(hf,'position');
set(hf,'position',pos.*[1 1 1.5 1]);
test = {'test','train'};
data = {'half dataset','full dataset'};
model = {'H1','H0'};
ylim = [0;0];

for j=1:2

    ha(j,1) = subplot(2,2,(j-1)*2+1,'parent',hf);
    mr = squeeze(mean(ner(j,:,:),3));
    vr = squeeze(var(ner(j,:,:),[],2)./Nmcmc);
    plotUncertainTimeSeries(mr',vr',[],ha(j,1));
    set(ha(j,1),'xlim',[0,3],'xtick',[1,2],'xticklabels',model)
    xlabel(ha(j,1),'type of simulated data')
    ylabel(ha(j,1),['P[correct prediction]'])
    title(ha(j,1),test{j})
    box(ha(j,1),'off')
    hold(ha(j,1),'on')
    plot(ha(j,1),[0,3],[0.5,0.5],'r--')
    mmr = mean(mner,2);
    plot(ha(j,1),[1,2],mmr,'go')

    ha(j,2) = subplot(2,2,(j-1)*2+2,'parent',hf);
    dF = F(j,:,:);
    mdF = squeeze(mean(dF,3));
    vdF = squeeze(var(dF,[],3)./Nmcmc);
    plotUncertainTimeSeries(mdF',vdF',[],ha(j,2));
    set(ha(j,2),'xlim',[0,3],'xtick',[1,2],'xticklabels',model)
    xlabel(ha(j,2),'type of simulated data')
    ylabel(ha(j,2),['log p(y|H1) - log p(y|H0)'])
    title(ha(j,2),'evidence for a mapping')
    box(ha(j,2),'off')
    title(ha(j,2),data{j})

    ylim(1) = min([ylim(1),get(ha(j,2),'ylim')]);
    ylim(2) = max([ylim(2),get(ha(j,2),'ylim')]);
end

for j=1:2
    set(ha(j,2),'ylim',ylim)
end

VBA_getSubplots ();

author: MBB-team

repo: VBA-toolbox


text:
function svgOffset(s, source, offset, result)
% Adds a feOffset SVG filter
% PRELIMINARY IMPLEMENTATION (Parameters may change)
%
% svgOffset(s, source, offset, result)
% Parameters:
%   s : Array of plot object handles
%   source : Any previous defined filter result string, 'SourceGraphic',
%            or 'SourceAlpha'.
%   offset : Offset value [x y]
%   result : String that identifies the filter result for following filter
%            stages.
for i = 1:length(s)
    userdata = get(s(i),'UserData');
    if isfield(userdata, 'svg') && isfield(userdata.svg, 'Filter')
        next = length(userdata.svg.Filter) + 1;
    else
        next = 1;
    end
    userdata.svg.Filter(next).Subfilter.Type = 'feOffset';
    userdata.svg.Filter(next).Subfilter.Source = source;
    userdata.svg.Filter(next).Subfilter.Offset = offset;
    userdata.svg.Filter(next).Subfilter.Result = result;
    set(s(i),'UserData', userdata);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function kernel = calibrate2D(obj, AtA, nCha, coil)
% retrieve calibration kernel from 2D calibration region
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

sampling = ones([obj.kernelSize,nCha],obj.measPara.precision);

dummyK = zeros(obj.kernelSize(1),obj.kernelSize(2),nCha,obj.measPara.precision);
dummyK((end+1)/2,(end+1)/2,coil) = 1;
idxY = find(dummyK);
sampling(idxY) = 0;
idxA = find(sampling);

Aty = AtA(:,idxY); Aty = Aty(idxA); % correlation values to target point, take complete neighbourhood and not just aquired ones
AtA = AtA(idxA,:); AtA =  AtA(:,idxA); % kick out the searched point

kernel = sampling*0;

lambda = norm(AtA,'fro')/size(AtA,1)*obj.calibTyk;

cnd = cond(AtA + eye(size(AtA))*lambda);
if(cnd == inf)
    rawkernel = pinv(AtA + eye(size(AtA))*lambda)*Aty; % grappa weighting values
else
    rawkernel = (AtA + eye(size(AtA))*lambda)\Aty; % grappa weighting values
end
kernel(idxA) = rawkernel;

end

author: thomaskuestner


text:
report_this_filefun(mfilename('fullpath'));

figure

l  = get(h1,'XLim');
s1_east = l(2); s2_west = l(1);
l  = get(h1,'YLim');
s3_north = l(2); s4_south = l(1);

cl =  get(h1,'Clim');

m_proj('lambert','long',[s2_west s1_east],'lat',[s4_south s3_north]);

l = re4 >= cl(2);
re4(l) = re4(l)*0+cl(2)-0.01;
l = re4 <= cl(1);
re4(l) = re4(l)*0+cl(1)+0.01;

cstep = abs(cl(2) - cl(1))/5 ;
[m,c] = m_contourf(gx,gy,re4,(cl(1):cstep:cl(2)));
%set(c,'LineColor','r');

hold on

if isempty(coastline) ==  0
    lico = m_line(coastline(:,1),coastline(:,2),'color','k','LineWidth',0.5);
end

if isempty(faults) == 0
    %lifa = m_line(faults(:,1),faults(:,2),'color','k');
end

m_grid('box','on','linestyle','none','tickdir','out','color','k','linewidth',2,...
    'fontsize',10,'fontname','Helveticabold');
hold on
shading flat

caxis([cl(1) cl(2)-cstep])

hold on

%li = m_line(a.Longitude,a.Latitude);
%set(li,'Linestyle','none','Marker',ty1,'MarkerSize',ms6,'color',co)
set(gcf,'Color','w')

vx =  (cl(1):cstep:cl(2));
v = [vx ; vx]; v = v';
rect = [0.93 0.4 0.015 0.25];
axes('position',rect)
pcolor((1:2),vx,v)
shading flat
set(gca,'XTickLabels',[])
set(gca,'FontSize',10,'FontWeight','normal',...
    'LineWidth',1.0,'YAxisLocation','right',...
    'Box','on','SortMethod','childorder','TickDir','out')

g = gray;
g = g(60:-1:1,:);

colormap(g)
brighten(0.2)
caxis([cl(1) cl(2)-cstep])