text:
function varargout = panel_femcond(varargin)
% PANEL_FEMCOND Edit FEM conductivity for a list of named layers (isotropic/anisotropic).
%
% USAGE:  bstPanel = panel_femcond('CreatePanel', OPTIONS)           : Call from the interactive interface
%         bstPanel = panel_femcond('CreatePanel', sProcess, sFiles)  : Call from the process editor
%                s = panel_femcond('GetPanelContents')

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

eval(macro_method);
end

%% ===== CREATE PANEL =====
function [bstPanelNew, panelName] = CreatePanel(sProcess, sFiles) %#ok<DEFNU>
    panelName = 'FemCondOptions';
    % Java initializations
    import java.awt.*;
    import javax.swing.*;

    % GUI CALL:  panel_femcond('CreatePanel', OPTIONS)
    if (nargin == 1)
        OPTIONS = sProcess;
    % PROCESS CALL:  panel_femcond('CreatePanel', sProcess, sFiles)
    else
        OPTIONS = sProcess.options.femcond.Value;
        % Get FEM files
        sSubject = bst_get('Subject', sProcess.options.subjectname.Value);
        if isempty(sSubject.iFEM)
            error('No available FEM mesh file for this subject.');
        end
        OPTIONS.FemFile = file_fullpath(sSubject.Surface(sSubject.iFEM).FileName);
    end

    % ==== GET MESH INFO ====
    % Load tissue labels
    FemMat = load(OPTIONS.FemFile, 'TissueLabels');
    % Get default conductivities
    OPTIONS.FemNames = FemMat.TissueLabels;
    OPTIONS.FemCond = panel_duneuro('GetDefaultCondutivity', OPTIONS.FemNames);

    % ==== FRAME STRUCTURE ====
    jPanelNew = java_create('javax.swing.JPanel');
    jPanelNew.setLayout(BoxLayout(jPanelNew, BoxLayout.PAGE_AXIS));
    jPanelNew.setBorder(BorderFactory.createEmptyBorder(12,12,12,12));

    % ===== FEM LAYERS =====
    jPanelLayers = gui_river([6,6], [-5,6,15,6], 'FEM conductivities');
        nLayers = length(OPTIONS.FemNames);
        jRadioLayerIso = javaArray('javax.swing.JRadioButton', nLayers);
        jRadioLayerAniso = javaArray('javax.swing.JRadioButton', nLayers);
        jTextCond = javaArray('javax.swing.JTextField', nLayers);
        % Loop on each layer
        for i = 1:nLayers
            gui_component('label', jPanelLayers, 'br', [OPTIONS.FemNames{i} ':'], [], [], [], []);
            jTextCond(i) = gui_component('texttime', jPanelLayers, 'tab', '', [], [], [], []);
            gui_validate_text(jTextCond(i), [], [], {0.0001,1000,10000}, 'S/m', [], OPTIONS.FemCond(i), []);
            gui_component('label', jPanelLayers, 'tab', 'S/m     ', [], [], [], []);
            jGroupRadio = ButtonGroup();
            jRadioLayerIso(i) = gui_component('radio', jPanelLayers, 'tab', 'Isotropic', jGroupRadio, [], @(h,ev)UpdatePanel(), []);
            jRadioLayerAniso(i) = gui_component('radio', jPanelLayers, 'tab', 'Anisotropic', jGroupRadio, [], @(h,ev)UpdatePanel(), []);
            jRadioLayerIso(i).setSelected(1);
        end
    jPanelNew.add(jPanelLayers);

    % ===== ANISOTROPY OPTIONS =====
    jPanelAniso = gui_river([2,2], [6,6,6,6], 'Anisotropy method');
        jGroupAniso = ButtonGroup();
        jRadioMethodEma = gui_component('radio', jPanelAniso, '', 'EMA: Effective Medium Approach (k=0.736) (DTI)
[Rullmann et al., 2009, Tuch et al., 1999, Haueisen et al., 2002]', jGroupAniso, [], @(h,ev)UpdatePanel(), []);
        jRadioMethodEmaVc = gui_component('radio', jPanelAniso, 'br', 'EMA+VC: EMA with volume constraint (DTI)
[Rullmann et al., 2009, Vorwerk et al., 2014]', jGroupAniso, [], @(h,ev)UpdatePanel(), []);
        jRadioMethodSim = gui_component('radio', jPanelAniso, 'br', 'Simulated: Artificial anisotropy (no DTI)', jGroupAniso, [], @(h,ev)UpdatePanel(), []);
        jRadioMethodEma.setSelected(1);
    jPanelNew.add(jPanelAniso);

    % ===== SIMULATED OPTIONS =====
    jPanelSim = gui_river([2,2], [6,6,6,6], 'Simulated anisotropy');
        % Ratio
        jLabelRatio = gui_component('label', jPanelSim, '', 'Ratio longitudinal/transversal: ', [], [], [], []);
        SimRatio = 10;
        jTextSimRatio = gui_component('texttime', jPanelSim, '', '', [], [], [], []);
        gui_validate_text(jTextSimRatio, [], [], {0.0001,100,100}, '', 2, SimRatio, []);
        % Constraint
        jGroupConstr = ButtonGroup();
        jRadioConstrWang = gui_component('radio', jPanelSim, 'br', 'Wang''s constraint (Wolters, 2003):
[sig_r*sig_t = sig_iso^2]', jGroupConstr, [], [], []);
        jRadioConstrWolters = gui_component('radio', jPanelSim, 'br', 'Volume''s constraint (Wolters, 2003):
[(4/3)*pi*(sig_r*sig_t^2) = (4/3)*pi*(sig_iso^3)]', jGroupConstr, [], [], []);
        jRadioConstrWolters.setSelected(1);
    jPanelNew.add(jPanelSim);

    % ===== VALIDATION BUTTONS =====
    jPanelValidation = gui_river([10 0], [6 10 0 10]);
        gui_component('Button', jPanelValidation, 'br right', 'Cancel', [], [], @ButtonCancel_Callback, []);
        gui_component('Button', jPanelValidation, [], 'OK', [], [], @ButtonOk_Callback, []);
    jPanelNew.add(jPanelValidation);

    % ===== PANEL CREATION =====
    % Return a mutex to wait for panel close
    bst_mutex('create', panelName);
    % Create the BstPanel object that is returned by the function
    ctrl = struct('jTextCond',           jTextCond, ...
                  'jRadioLayerIso',      jRadioLayerIso, ...
                  'jRadioLayerAniso',    jRadioLayerAniso, ...
                  'jPanelAniso',         jPanelAniso, ...
                  'jRadioMethodEma',     jRadioMethodEma, ...
                  'jRadioMethodEmaVc',   jRadioMethodEmaVc, ...
                  'jRadioMethodSim',     jRadioMethodSim, ...
                  'jPanelSim',           jPanelSim, ...
                  'jRadioConstrWang',    jRadioConstrWang, ...
                  'jRadioConstrWolters', jRadioConstrWolters, ...
                  'jLabelRatio',         jLabelRatio, ...
                  'jTextSimRatio',       jTextSimRatio);
    ctrl.FemNames = OPTIONS.FemNames;
    % Create the BstPanel object that is returned by the function
    bstPanelNew = BstPanel(panelName, jPanelNew, ctrl);
    % Redraw panel
    UpdatePanel();

%% =================================================================================
%  === LOCAL CALLBACKS  ============================================================
%  =================================================================================
    %% ===== BUTTON: CANCEL =====
    function ButtonCancel_Callback(varargin)
        % Close panel
        gui_hide(panelName);
    end

    %% ===== BUTTON: OK =====
    function ButtonOk_Callback(varargin)
        % Release mutex and keep the panel opened
        bst_mutex('release', panelName);
    end

    %% ===== UPDATE PANEL ======
    function UpdatePanel()
        % FEM Layers
        isIsotropic = false(1, nLayers);
        for j = 1:nLayers
            isIsotropic(j) = jRadioLayerIso(j).isSelected();
            % jTextCond(j).setEnabled(isIsotropic(j));
        end
        % All isotropic: disable anisotropy options
        isAllIso = all(isIsotropic);
        jPanelAniso.setVisible(~isAllIso);
        % Simulated
        isSim = ~isAllIso && jRadioMethodSim.isSelected();
        jPanelSim.setVisible(isSim);

        % Get panel
        [bstPanel iPanel] = bst_get('Panel', 'FemCondOptions');
        container = get(bstPanel, 'container');
        % Re-pack frame
        if ~isempty(container)
            jFrame = container.handle{1};
            if ~isempty(jFrame)
                jFrame.pack();
            end
        end
    end
end

%% =================================================================================
%  === EXTERNAL CALLBACKS  =========================================================
%  =================================================================================
%% ===== GET PANEL CONTENTS =====
function s = GetPanelContents() %#ok<DEFNU>
    % Get panel controls handles
    ctrl = bst_get('PanelControls', 'FemCondOptions');
    if isempty(ctrl)
        s = [];
        return;
    end
    % FEM layers
    for i = 1:length(ctrl.jTextCond)
        s.FemCond(i) = str2double(char(ctrl.jTextCond(i).getText()));
        s.isIsotropic(i) = ctrl.jRadioLayerIso(i).isSelected();
    end
    % Anisotropy options
    if ctrl.jRadioMethodEma.isSelected()
        s.AnisoMethod = 'ema';
    elseif ctrl.jRadioMethodEmaVc.isSelected()
        s.AnisoMethod = 'ema+vc';
    elseif ctrl.jRadioMethodSim.isSelected()
        s.AnisoMethod = 'simulated';
    end
    % Simulated: Ratio
    s.SimRatio = str2double(char(ctrl.jTextSimRatio.getText()));
    % Simulated: Constraint method
    if ctrl.jRadioConstrWang.isSelected()
        s.SimConstrMethod = 'wang';
    elseif ctrl.jRadioConstrWolters.isSelected()
        s.SimConstrMethod = 'wolters';
    end
end

author: brainstorm-tools


text:
function [d] = spm_kl_eig_normal (m_q,c_q,c_p)
% KL divergence between normal densities using eigendecomposition
% function [d] = spm_kl_eig_normal (m_q,c_q,c_p)
%
% Calculate the KL distance
%
% KL (Q||P) = <log Q/P> where avg is wrt Q
%
% between two Normal densities Q and P where P is
% zero mean and has a diagonal covariance.
%
% m_q, c_q    Mean and covariance of first Normal density
% c_p         Covariance of second (zero-mean) Normal density
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_kl_eig_normal.m 1143 2008-02-07 19:33:33Z spm $

d=length(m_q);
m_q=m_q(:);

ldcp=0.5*sum(log(diag(c_p)));
[v,eigvals]=eig(c_q);
ldcq=-0.5*sum(log(diag(eigvals)));
Term1=ldcp+ldcq;

inv_c_p=diag(1./diag(c_p));
Term2=0.5*trace(inv_c_p*c_q)+0.5*m_q'*inv_c_p*m_q;
d=Term1+Term2-0.5*d;

author: spm


text:
function p04_story ( )

%*****************************************************************************80
%
%% P04_STORY prints the "story" for problem p04.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Larry Irvine, Samuel Marin, Philip Smith,
%    Constrained Interpolation and Smoothing,
%    Constructive Approximation,
%    Volume 2, Number 1, December 1986, pages 129-151.
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This example is due to Larry Irvine, Samuel Marin and Philip Smith.\n' );
  fprintf ( 1, '  This data can cause problems for interpolation methods.\n' );
  fprintf ( 1, '  There are sudden changes in direction, and at the same time,\n' );
  fprintf ( 1, '  sparsely-placed data.  This can cause an interpolant to overshoot\n' );
  fprintf ( 1, '  the data in a way that seems implausible.\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function r = YZRectangle(x,y1,y2,z1,z2)

%YZRECTANGLE Rectangle in the YZ plane.

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

s1 = ySegment(x,y1,y2,z1);
s2 = zSegment(x,y2,z1,z2);
s3 = ySegment(x,y2,y1,z2);
s4 = zSegment(x,y1,z2,z1);

r = [s1 s2 s3 s4];

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% This script depends on
% - ex_subspace_angles.m
% - ex_omp_rep_subspace_angles.m
% Please execute the dependencies before
% running this script.

close all; clear all; clc;

export = true;

orig = load('bin/principal_angles');
mf = spx.graphics.Figures();
mf.new_figure('Original Principal angles');
hist(orig.off_diag_angles, 200);
xlabel('Principal angle (degrees)');
ylabel('Number of subspace pairs');
title('Distribution of principal angles over subspace pairs in signal space');
grid on;

if export
export_fig bin/images/signal_space_principal_angles.png -r120 -nocrop;
export_fig bin/images/signal_space_principal_angles.pdf;
end

rep = load('bin/omp_rep_principal_angles');
mf.new_figure('Representation Principal angles');
hist(rep.off_diag_angles, 200);
xlabel('Principal angle (degrees)');
ylabel('Number of subspace pairs');
title('Distribution of principal angles over subspace pairs in representations');
grid on;

if export
export_fig bin/images/representation_space_principal_angles.png -r120 -nocrop;
export_fig bin/images/representation_space_principal_angles.pdf;
end

author: indigits

repo: sparse-plex


text:
function child = mutateParent(parent)
% Oren Rosen
% The MathWorks
% 8/29/2007
%
% This custom mutation function is written to work on a population of
% vectors of zeros and ones with the same amount of ones in each vector.
% The mutated child is formed by randomly permuting the elements of the
% parent.
% Note: Performance-wise this hasn't worked out to be that efficient. A
% better implementation may swap only two of the elements.

% *** Calculate dimensions of parent ***
numBits = length(parent);

% *** Randomly permute bits
child = parent( randperm(numBits) );

% *** Display results ***
disp(' ');
disp(['Given: parent = [ ',num2str(parent),' ]']);
disp(' ');
disp(['        child = [ ',num2str(child),' ]']);

author: Sable

repo: mcbench-benchmarks


text:
function [EC,ec,degij] = edge_nei_overlap_bd(CIJ)
%EDGE_NEI_OVERLAP_BD     overlap amongst neighbors of two adjacent nodes
%
%   [EC,ec,degij] = edge_nei_bd(CIJ);
%
%   This function determines the neighbors of two nodes that are linked by
%   an edge, and then computes their overlap.  Connection matrix must be
%   binary and directed.  Entries of 'EC' that are 'inf' indicate that no
%   edge is present.  Entries of 'EC' that are 0 denote "local bridges",
%   i.e. edges that link completely non-overlapping neighborhoods.  Low
%   values of EC indicate edges that are "weak ties".
%
%   If CIJ is weighted, the weights are ignored. Neighbors of a node can be
%   linked by incoming, outgoing, or reciprocal connections.
%
%   Inputs:     CIJ,      directed (binary/weighted) connection matrix
%
%   Outputs:    EC,     edge neighborhood overlap matrix
%               ec,     edge neighborhood overlap per edge, in vector format
%               degij,  degrees of node pairs connected by each edge
%
%   Reference:
%
%       Easley and Kleinberg (2010) Networks, Crowds, and Markets.
%           Cambridge University Press, Chapter 3
%
%   Olaf Sporns, Indiana University, 2012

[ik,jk,ck] = find(CIJ);
lel = length(ck);
N = size(CIJ,1);

[~,~,deg] = degrees_dir(CIJ);

ec = zeros(1,lel);
degij = zeros(2,lel);
for e=1:lel
    neiik = setdiff(union(find(CIJ(ik(e),:)),find(CIJ(:,ik(e))')),[ik(e) jk(e)]);
    neijk = setdiff(union(find(CIJ(jk(e),:)),find(CIJ(:,jk(e))')),[ik(e) jk(e)]);
    ec(e) = length(intersect(neiik,neijk))/length(union(neiik,neijk));
    degij(:,e) = [deg(ik(e)) deg(jk(e))];
end;

ff = find(CIJ);
EC = 1./zeros(N);
EC(ff) = ec;                            %#ok<FNDSB>

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function idx = getclassindex(cnames,g)
%GETCLASSINDEX Find indices for class names in another list of names
%   IDX = GETCLASSINDEX(CNAMES,G) takes a list CNAMES of class names
%   (such as the grouping variable values in the treefit or classify
%   function) and another list G of group names (as might be supplied
%   in the "prior" argument to those functions), and finds the indices
%   of the CNAMES names in the G list.  CNAMES should be a cell array
%   of strings.  G can be numbers, a string array, or a cell array of
%   strings

%   Copyright 1993-2002 The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 2002/03/22 22:03:02 $

% Convert to common string form, whether input is char, cell, or numeric
if isnumeric(g)
   g = cellstr(strjust(num2str(g(:)), 'left'));
elseif ~iscell(g)
   g = cellstr(g);
end

nclasses = length(cnames);
idx = zeros(1,nclasses);

% Look up each class in the grouping variable.
for i = 1:nclasses
   j = strmatch(cnames{i}, g, 'exact');
   if ~isempty(j)
      idx(i) = j(1);
   end
end

author: zouchuhang

repo: LayoutNet


text:
function [TableHHStats] = getBasicHarveyStats(male,female)
% This function compiles basic statistics on the male and female whole-body
% metabolic model
%
% function [TableHHStats] = getBasicHarveyStats(male,female)
%
% INPUT
% male          Model structure containing the male model
% female        Model structure containing the female model
%
% OUTPUT
% TableHHStats  Table containing the statistics
%
% Ines Thiele 2016- 2019

% build table
clear TableHHStats
c = 1;

SL = (find(~cellfun(@isempty,strfind(female.mets,'slack_'))));
SL_female = SL;
SL = (find(~cellfun(@isempty,strfind(male.mets,'slack_'))));
SL_male = SL;

TableHHStats{c , 1} = '';
TableHHStats{c , 2} = 'Harvetta';
TableHHStats{c , 3} = 'Harvey'; c = c+1;

TableHHStats{c , 1} = 'Number of Reactions';

TableHHStats{c , 2} = num2str(length(female.rxns));
TableHHStats{c , 3} = num2str(length(male.rxns)); c = c+1;

TableHHStats{c , 1} = 'Number of Metabolites';

TableHHStats{c , 2} = num2str(length(female.mets)-length(SL_female));
TableHHStats{c , 3} = num2str(length(male.mets)-length(SL_male)); c = c+1;

TableHHStats{c , 1} = 'Number of Genes (transcripts)';
TableHHStats{c , 2} = num2str(length(female.genes));
TableHHStats{c , 3} = num2str(length(male.genes)); c = c+1;

TableHHStats{c , 1} = 'Number of Genes (unique)';
[a,b]=strtok(female.genes,'.');

TableHHStats{c , 2} = num2str(length(unique(a)));
[a,b]=strtok(male.genes,'.');
TableHHStats{c , 3} = num2str(length(unique(a))); c = c+1;

TableHHStats{c , 1} = 'Number of Subsystems';

% Feb 2020, moved this to loadPSCMfile
% Oct 2017 % clean up final remaining subsystem names
% female.subSystems(strmatch('Transport, endoplasmic reticular',female.subSystems,'exact'))={'Transport, endoplasmic reticulum'};
% female.subSystems(strmatch('Arginine and Proline Metabolism',female.subSystems,'exact'))={'Arginine and proline Metabolism'};
% female.subSystems(strmatch(' ',female.subSystems,'exact'))={'Miscellaneous'};
%
% male.subSystems(strmatch('Transport, endoplasmic reticular',male.subSystems,'exact'))={'Transport, endoplasmic reticulum'};
% male.subSystems(strmatch('Arginine and Proline Metabolism',male.subSystems,'exact'))={'Arginine and proline Metabolism'};
% male.subSystems(strmatch(' ',male.subSystems,'exact'))={'Miscellaneous'};

TableHHStats{c , 2} = num2str(length(unique(female.subSystems))-1);%do not count dummy reactions
TableHHStats{c , 3} = num2str(length(unique(male.subSystems))-1); c = c+1;%do not count dummy reactions

TableHHStats{c , 1} = 'Blood compartment metabolites';
BC_mets = (find(~cellfun(@isempty,strfind(female.mets,'[bc]'))));
SL_female = (find(~cellfun(@isempty,strfind(female.mets,'slack_'))));
BC_mets_female = setdiff(BC_mets,SL_female);

BC_mets = (find(~cellfun(@isempty,strfind(male.mets,'[bc]'))));
SL_male = (find(~cellfun(@isempty,strfind(male.mets,'slack_'))));
BC_mets_male = setdiff(BC_mets,SL_male);

TableHHStats{c , 2} = num2str(length(BC_mets_female));
TableHHStats{c , 3} = num2str(length(BC_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Urine metabolites';
U_mets = (find(~cellfun(@isempty,strfind(female.mets,'[u]'))));
U_mets_female = setdiff(U_mets,SL_female);

U_mets = (find(~cellfun(@isempty,strfind(male.mets,'[u]'))));
U_mets_male = setdiff(U_mets,SL_male);

TableHHStats{c , 2} = num2str(length(U_mets_female));
TableHHStats{c , 3} = num2str(length(U_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Portal vein metabolites';
BP_mets = (find(~cellfun(@isempty,strfind(female.mets,'[bp]'))));
BP_mets_female = setdiff(BP_mets,SL_female);

BP_mets = (find(~cellfun(@isempty,strfind(male.mets,'[bp]'))));
BP_mets_male = setdiff(BP_mets,SL_male);

TableHHStats{c , 2} = num2str(length(BP_mets_female));
TableHHStats{c , 3} = num2str(length(BP_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Bile duct metabolites';

BD_mets = (find(~cellfun(@isempty,strfind(female.mets,'[bd]'))));
BD_mets_female = setdiff(BD_mets,SL_female);

BD_mets = (find(~cellfun(@isempty,strfind(male.mets,'[bd]'))));
BD_mets_male = setdiff(BD_mets,SL_male);

TableHHStats{c , 2} = num2str(length(BD_mets_female));
TableHHStats{c , 3} = num2str(length(BD_mets_male)); c = c+1;

TableHHStats{c , 1} = 'CSF metabolites';

CSF_mets = (find(~cellfun(@isempty,strfind(female.mets,'[csf]'))));
CSF_mets_female = setdiff(CSF_mets,SL_female);

CSF_mets = (find(~cellfun(@isempty,strfind(male.mets,'[csf]'))));
CSF_mets_male = setdiff(CSF_mets,SL_male);

TableHHStats{c , 2} = num2str(length(CSF_mets_female));
TableHHStats{c , 3} = num2str(length(CSF_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Diet metabolites';

D_mets = (find(~cellfun(@isempty,strfind(female.mets,'[d]'))));
D_mets_female = setdiff(D_mets,SL_female);
D_mets = (find(~cellfun(@isempty,strfind(male.mets,'[d]'))));
D_mets_male = setdiff(D_mets,SL_male);

TableHHStats{c , 2} = num2str(length(D_mets_female));
TableHHStats{c , 3} = num2str(length(D_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Fecal metabolites';

Fe_mets = (find(~cellfun(@isempty,strfind(female.mets,'[fe]'))));
Fe_mets_female = setdiff(Fe_mets,SL_female);

Fe_mets = (find(~cellfun(@isempty,strfind(male.mets,'[fe]'))));
Fe_mets_male = setdiff(Fe_mets,SL_male);

TableHHStats{c , 2} = num2str(length(Fe_mets_female));
TableHHStats{c , 3} = num2str(length(Fe_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Sweat metabolites';

Fe_mets = (find(~cellfun(@isempty,strfind(female.mets,'[sw]'))));
Fe_mets_female = setdiff(Fe_mets,SL_female);

Fe_mets = (find(~cellfun(@isempty,strfind(male.mets,'[sw]'))));
Fe_mets_male = setdiff(Fe_mets,SL_male);

TableHHStats{c , 2} = num2str(length(Fe_mets_female));
TableHHStats{c , 3} = num2str(length(Fe_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Air metabolites';

A_mets = (find(~cellfun(@isempty,strfind(female.mets,'[a]'))));
A_mets_female = setdiff(A_mets,SL_female);

A_mets = (find(~cellfun(@isempty,strfind(male.mets,'[a]'))));
A_mets_male = setdiff(A_mets,SL_male);

TableHHStats{c , 2} = num2str(length(A_mets_female));
TableHHStats{c , 3} = num2str(length(A_mets_male)); c = c+1;

TableHHStats{c , 1} = 'Milk metabolites';

Fe_mets = (find(~cellfun(@isempty,strfind(female.mets,'[mi]'))));
Fe_mets_female = setdiff(Fe_mets,SL_female);

Fe_mets = (find(~cellfun(@isempty,strfind(male.mets,'[mi]'))));
Fe_mets_male = setdiff(Fe_mets,SL_male);

TableHHStats{c , 2} = num2str(length(Fe_mets_female));
TableHHStats{c , 3} = num2str(length(Fe_mets_male)); c = c+1;

author: opencobra

repo: cobratoolbox


text:
function alpha = p00_alpha ( problem )

%*****************************************************************************80
%
%% P00_ALPHA returns the value of ALPHA for any problem.
%
%  Discussion:
%
%    ALPHA is the lower, finite limit of integration in the integral.
%
%    The typical or default value is 0.0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the index of the problem.
%
%    Output, real ALPHA, the value of ALPHA.
%
  if ( problem == 1 )
    alpha = p01_alpha ( );
  elseif ( problem == 2 )
    alpha = p02_alpha ( );
  elseif ( problem == 3 )
    alpha = p03_alpha ( );
  elseif ( problem == 4 )
    alpha = p04_alpha ( );
  elseif ( problem == 5 )
    alpha = p05_alpha ( );
  elseif ( problem == 6 )
    alpha = p06_alpha ( );
  elseif ( problem == 7 )
    alpha = p07_alpha ( );
  elseif ( problem == 8 )
    alpha = p08_alpha ( );
  elseif ( problem == 9 )
    alpha = p09_alpha ( );
  elseif ( problem == 10 )
    alpha = p10_alpha ( );
  elseif ( problem == 11 )
    alpha = p11_alpha ( );
  elseif ( problem == 12 )
    alpha = p12_alpha ( );
  elseif ( problem == 13 )
    alpha = p13_alpha ( );
  elseif ( problem == 14 )
    alpha = p14_alpha ( );
  elseif ( problem == 15 )
    alpha = p15_alpha ( );
  elseif ( problem == 16 )
    alpha = p16_alpha ( );
  elseif ( problem == 17 )
    alpha = p17_alpha ( );
  elseif ( problem == 18 )
    alpha = p18_alpha ( );
  elseif ( problem == 19 )
    alpha = p19_alpha ( );
  elseif ( problem == 20 )
    alpha = p20_alpha ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_ALPHA - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', problem );
    error ( 'P00_ALPHA - Fatal error!' );
  end

  return
end


text:
function val = m00Q0011(F00, F11)
%------------------------------------------------------------------------------
%
% Integrates interpolated gridfunction F on quincunx grid (F00 united with F11)
% where interpolation is assumed piecewise constant.
%
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: February 1, 2001.
%  2001 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
[hx, hy] = Q0011gridfdims(F00, F11);
Qhxy = 2 * hx * hy;
val = ( sum(sum(F00)) + sum(sum(F11)) ) * Qhxy;
%------------------------------------------------------------------------------


text:
function r4mat_uniform_ab_test ( )

%*****************************************************************************80
%
%% R4MAT_UNIFORM_AB_TEST tests R4MAT_UNIFORM_AB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 Deember 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;
  a = -1.0;
  b = +5.0;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4MAT_UNIFORM_AB_TEST\n' );
  fprintf ( 1, '  R4MAT_UNIFORM_AB computes a random R4MAT.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %g <= x <= %g\n', a, b );
  fprintf ( 1, '  Initial seed is %d\n', seed );

  [ v, seed ] = r4mat_uniform_ab ( m, n, a, b, seed );

  r4mat_print ( m, n, v, '  Uniform R4MAT:' );

  return
end


text:
function b = rhs(disc, f)
%RHS   Discretize the right-hand side of a linear system for ULTRAS.
%   B = RHS(DISC, F) returns a discrete version, B,  of the function (or
%   chebmatrix) F, as defined by the discretization DISC.
%
% See also MATRIX, REDUCE.

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

% Developers note: This method works as follows.
%   The original function(s) are discretized, then reduced in dimension the same
%   way as the operator. The constraints are prepended to the top of the vector.

% Create an ULTRAS object from the CHEBMATRIX F:
fDisc = ultraS(f, disc.dimension, disc.domain);
fDisc.outputSpace = disc.outputSpace;

% Instantiate (discretize) the ULTRAS discretisation.
b = cell2mat(instantiate(fDisc));

% Developer note:
%   The continuity conditions go above the constraints. See getConstraints().

% Prepend the values of the constraints and continuity conditions.
L = disc.source;
if ( ~isempty(L.constraint) )
    b = [ L.constraint.values ; b ];
end
if ( ~isempty(L.continuity) )
    b = [ L.continuity.values ; b ];
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y = getFieldSafe(data,fld)
% y = getFieldSafe(data,fld)
% mode: optional, either 'safe' or 'safest' (default).
% y: col cell

if isempty(data)
    y = cell(0,1);
else
    y = {data.(fld)}';
    assert(isequal(size(y),[numel(data) 1]));
end

author: kristinbranson

repo: JAABA


text:
function [solution, nIterations, bestApprox] = sparseLP(model, approximation, params)
% DC programming for solving the sparse LP
% :math:`min ||x||_0` subject to linear constraints
% See `Le Thi et al., DC approximation approaches for sparse optimization,
% European Journal of Operational Research, 2014`;
% http://dx.doi.org/10.1016/j.ejor.2014.11.031
%
% USAGE:
%
%    [solution, nIterations, bestApprox] = sparseLP(model, approximation, params);
%
% INPUTS:
%    model:       Structure containing the following fields describing the linear constraints:
%
%                        * .A - `m x n` LHS matrix
%                        * .b - `m x 1` RHS vector
%                        * .lb - `n x 1` Lower bound vector
%                        * .ub - `n x 1` Upper bound vector
%                        * .csense - `m x 1` Constraint senses, a string containting the model sense for
%                          each row in `A` ('E', equality, 'G' greater than, 'L' less than).
%
% OPTIONAL INPUTS
%    approximation:    appoximation type of zero-norm. Available approximations:
%
%                        * 'cappedL1' : Capped-L1 norm
%                        * 'exp'      : Exponential function
%                        * 'log'      : Logarithmic function
%                        * 'SCAD'     : SCAD function
%                        * 'lp-'      : `L_p` norm with `p < 0`
%                        * 'lp+'      : `L_p` norm with `0 < p < 1`
%                        * 'l1'       : L1 norm
%                        * 'all'      : try all approximations and return the best result
%
% OPTIONAL INPUTS:
%    params:           Parameters structure:
%
%                        * .nbMaxIteration - stopping criteria - number maximal of iteration (Defaut value = 1000)
%                        * .epsilon - stopping criteria - (Defaut value = 10e-6)
%                        * .theta - parameter of the approximation (Defaut value = 0.5)
%
% OUTPUT:
%    solution:         Structure containing the following fields:
%
%                        * .x - `n x 1` solution vector
%                        * .stat - status:
%
%                          * 1 =  Solution found
%                          * 2 =  Unbounded
%                          * 0 =  Infeasible
%                          * -1=  Invalid input
%
%   nIterations:       Number of iterations
%   bestApprox:        Best approximation
%
% .. Author: - Hoai Minh Le,	20/10/2015
%              Ronan Fleming,    2017

availableApprox = {'exp','log','SCAD','lp-','lp+','l1','cappedL1','all'};

if ~exist('approximation','var')
    approximation='cappedL1';
end

% Check inputs
if nargin < 3
    params.nbMaxIteration = 1000;
    params.epsilon = 1e-6;
    params.theta   = 0.5;
    params.pNeg = -1;
    params.pPos = 0.5;
else
    if ~isfield(params,'nbMaxIteration')
        params.nbMaxIteration = 1000;
    end

    if ~isfield(params,'epsilon')
        params.epsilon = 1e-6;
    end

    if ~isfield(params,'theta')
        params.theta   = 0.5;
    end

    if ~isfield(params,'pNeg')
        params.pNeg = -1;
    end

    if ~isfield(params,'pPos')
        params.pPos = 0.5;
    end

end

if ~isfield(model,'A')
    error('Error:LHS matrix is not defined');
    solution.stat = -1;
    return;
end
if ~isfield(model,'b')
    error('RHS vector is not defined');
    solution.stat = -1;
    return;
end
if ~isfield(model,'lb')
    error('Lower bound vector is not defined');
    solution.stat = -1;
    return;
end
if ~isfield(model,'ub')
    error('Upper bound vector is not defined');
    solution.stat = -1;
    return;
end
if ~isfield(model,'csense')
    error('Constraint sense vector is not defined');
    solution.stat = -1;
    return;
end

if ~ismember(approximation,availableApprox)
    error('Approximation is not valid');
    solution.stat = -1;
    return;
end

bestApprox = '';

switch approximation
    case 'all'
        approximations = setdiff(availableApprox,'all','stable');
        bestResult = size(model.A,2);
        bestSolution.x = [];
        bestSolution.stat = 0;
        bestIterations = 0;
        feasTol = getCobraSolverParams('LP','feasTol');
        for i=1:length(approximations)
            %disp(approximations(i))
            %try
            [candSolution,candIterations] = sparseLP(model,approximations{i},params);
            %catch
            %fail gracefully
            %solutionL0.stat = 0;
            %end
            if candSolution.stat == 1
                candResult = nnz(abs(candSolution.x) > feasTol);
                if bestResult >= candResult
                    bestResult = candResult;
                    bestApprox = approximations{i};
                    bestSolution = candSolution;
                    bestIterations = candIterations;
                end
            end
        end
        solution = bestSolution;
        nIterations = bestIterations;

    otherwise
        [nbMaxIteration,epsilon,theta,pNeg,pPos] = deal(params.nbMaxIteration,params.epsilon,params.theta,params.pNeg,params.pPos);
        [A,b,lb,ub,csense] = deal(model.A,model.b,model.lb,model.ub,model.csense);

        %Parameters
        nbIteration = 0;
        epsilonP = 10e-2;
        alpha = 3;
        [m,n] = size(A);

        %Create the linear sub-programme that one needs to solve at each iteration, only its
        %objective function changes, the constraints set remains.

        % Constraints
        % Ax <=b
        % t >= x
        % t >= -x
        A2 = [A         sparse(m,n);
            speye(n) -speye(n);
            -speye(n) -speye(n)];
        b2 = [b; zeros(2*n,1)];
        csense2 = [csense;repmat('L',2*n, 1)];

        % Bound;
        % lb <= x <= ub
        % 0  <= t <= max(|lb|,|ub|)
        lb2 = [lb;zeros(n,1)];
        ub2 = [ub;max(abs(lb),abs(ub))];

        %Define the linear sub-problem
        subLPproblem = struct('osense',1,'A',A2,'csense',csense2,'b',b2,'lb',lb2,'ub',ub2);

        %Initialisation
        x = zeros(n,1);
        obj_old = evalObj(x,theta,pNeg,pPos,epsilonP,alpha,approximation);

        %DCA
        tic
        while nbIteration < nbMaxIteration

            x_old = x;

            %Solve the linear problem
            subLPproblem.c = updateObj(x,theta,pNeg,pPos,epsilonP,alpha,approximation);
            solution = solveCobraLP(subLPproblem);

            if solution.stat == 1
                x = solution.full(1:n);

                %Check stopping criterion
                error_x = norm(x - x_old);
                obj_new = evalObj(x,theta,pNeg,pPos,epsilonP,alpha,approximation);
                error_obj = abs(obj_new - obj_old);
                if (error_x < epsilon) || (error_obj < epsilon)
                    break;
                else
                    obj_old = obj_new;
                end
                % Automatically update the approximation parameter theta
                if theta < 1000
                    theta = theta * 1.5;
                end
                nbIteration = nbIteration + 1;
                %             disp(strcat('DCA - Iteration: ',num2str(nbIteration)));
                %             disp(strcat('Obj:',num2str(obj_new)));
                %             disp(strcat('Stopping criteria error: ',num2str(min(error_x,error_obj))));
                %             disp('=================================');
            else
                x = [];
                break;
            end
        end

        time = toc;
        solution.x = x;
        solution.time = time;
        nIterations = nbIteration;

end


text:
function b = legcoeffs(f, n)
%LEGCOEFFS   Compute Legendre series coefficients of a CHEBTECH object.
%   B = LEGCOEFFS(F) returns the Legendre series coefficients of CHEBTECH F, so
%   that F = B(1)*P_0 + ... + B(N)*P_(N-1), where P_k is the kth Legendre
%   polynomial. B is a vector of the same length as that of F.
%
%   B = LEGCOEFFS(F, N) returns the first N coefficients. If length(F) < N
%   then the additional entries of B are padded with zeros.
%
%   If F is an array-valued CHEBTECH, then a matrix of coefficients is returned
%   so that F(:,k) = B(1,k)*P_0 + ... + B(N,k)*P_(N-1).
%
% See also CHEBCOEFFS.

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

b = cheb2leg(f.coeffs);

if ( nargin > 1 )
    s = size(b);
    if ( s(1) > n )
        b = b(1:n, :);
    else
        b = [b; zeros(n-s(1), s(2))];
    end
end

end


text:
function lidar(fileName, sensorType, varName)

    %% Check if input is valid
    if ~isstring(fileName) && ~ischar(fileName)
        error("File name must a string or char array\n")
    end

    if ~strcmp(sensorType, "lidar")
       error("Wrong sensor selected\n")
    end

    if ~isvarname(varName)
        error("Invalid variable name\n")
    end

    %% Generate the python file
    file = fopen(strcat(fileName, '.py'), 'w');

    % Automatically genrates a python file containing the sensor call back
    % bindings
    fprintf(file, 'import numpy as np\n');
    fprintf(file, '\n');
    fprintf(file, 'def bindSensor(sensor):\n');
    fprintf(file, '    sensor.listen(lambda _image: do_something(_image))\n');
    fprintf(file, '\n');
    fprintf(file, 'def do_something(_image):\n');
    fprintf(file, '    global %s\n', varName);
    fprintf(file, '    data = np.frombuffer(_image.raw_data, dtype="float32")\n');
    fprintf(file, '\n');
    fprintf(file, '    # Pair up in [x,y,z] format\n');
    fprintf(file, '    data = np.reshape(data, (-1,4))\n');
    fprintf(file, '\n');
    fprintf(file, '    # Convert the data into MATLAB cast compatible type\n');
    fprintf(file, '    %s = np.ascontiguousarray(data)\n', varName);

    fclose(file);
end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
%--- help for sym/jacobian ---
%
% JACOBIAN Jacobian matrix.
%    JACOBIAN(f,x) computes the Jacobian of the scalar or vector f
%    with respect to the vector x. The (i,j)-th entry of the result
%    is df(i)/dx(j). Note that when f is scalar, the Jacobian of f
%    is the gradient of f. Also, note that scalar x is allowed,
%    although this is just DIFF(f,x).
%
%    Example:
%        syms x y z u v; jacobian([x*y*z; y; x+z],[x y z])
%        returns  [y*z, x*z, x*y; 0, 1, 0; 1, 0, 1]
%
%        jacobian(u*exp(v),[u;v])
%        returns  [exp(v), u*exp(v)]
%
%    See also SYM/CURL, SYM/DIFF, SYM/DIVERGENCE, SYM/GRADIENT, SYM/HESSIAN,
%    SYM/POTENTIAL, CURL, DIVERGENCE, HESSIAN, LAPLACIAN, VECTORPOTENTIAL,
%    SUBS.
%

author: jmaih


text:
% closest fly, based on dell2nose
function [data,units,mind] = compute_closestfly_ell2nose(trx,n,dosave_d)

if nargin < 3,
  dosave_d = true;
end

flies = trx.exp2flies{n};
nflies = numel(flies);
closestfly = cell(1,nflies);
mind = cell(1,nflies);

parfor i1 = 1:nflies,
  fly1 = flies(i1);
  fprintf('fly1 = %d\n',fly1);
  flies2 = flies(trx.roi(fly1)==trx.roi(flies));
  d = nan(numel(flies2),trx(fly1).nframes);

  % use dcenter2nose and major axis length to compute upper and lower bounds on
  % dell2nose
  [mindupper,dlower] = dell2nose_bounds(trx,fly1,flies2);

  for i2 = 1:numel(flies2),
    fly2 = flies2(i2);
    if fly1 == fly2,
      continue;
    end
    % only try for frames where lower bound is smaller than min upper bound
    idx1try = find(mindupper >= dlower(i2,:));
    d(i2,:) = dell2nose_pair(trx,fly1,fly2,idx1try);
  end
  [mind{i1},closesti] = min(d,[],1);
  closestfly{i1} = flies2(closesti);
  closestfly{i1}(isnan(mind{i1})) = nan;
end

% so that we don't compute dcenter twice
if dosave_d,
  data = mind; %#ok<NASGU>
  units = parseunits('mm'); %#ok<NASGU>
  filename = trx.GetPerFrameFile('dell2nose',n);
  save(filename,'data','units');
end

data = closestfly;
units = parseunits('unit');


text:
function newROI = dtiMergeROIs(roi1,roi2)
%
%   newROI = dtiMergeROIs(roi1,roi2)
%
% Author: RFD, BW
% Purpose:
%  Merge ROIs
%
%Example:
% handles = guidata(gcf); roi1 = handles.rois(1); roi2 = handles.rois(2);
% newROI = dtiMergeROIs(roi1,roi2);
% handles = dtiAddROI(newROI,handles,1);
% guidata(gcf,handles);
% dtiRefresh ... or something

newROI = roi1;
newROI.name = sprintf('%s_%s',roi1.name,roi2.name);

newROI.coords = unique(vertcat(roi1.coords,roi2.coords),'rows');

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [dyn] = N2dyn(N)
% Convert force from newtons to dyne.
% Chad A. Greene 2012
dyn = N*100000;


text:
%NUMROWS Number of rows in matrix
%
% NR = NUMROWS(M) is the number of rows in the matrix M.
%
% Notes::
% - Readable shorthand for SIZE(M,1);
%
% See also NUMCOLS, SIZE.

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

function r = numrows(m)

	r = size(m, 1);

author: Allopart

repo: rbpf-gmapping


text:
function view = editMenu(view)
%
% view = editMenu(view)
%
% Set up the callbacks of the EDIT menu.
%
% djh, 2/28/2001
editMenu = uimenu('Label','Edit','Separator','on');

% edit mrSESSION:
% [mSESSION,ok] = EditSession(mrSESSION);
% if ok, saveSession(1), end;
cb=['[mrSESSION,ok] = EditSession(mrSESSION,0); ',...
    'if ok, saveSession(1); end;'];
uimenu(editMenu,'Label','Edit/view mrSESSION','Separator','off',...
   'CallBack',cb);

% Edit alignment(using mrRx):
uimenu(editMenu,'Label','Edit/View Alignment...','Separator','off',...
    'Callback','rxAlign;', 'Accelerator', '7');

datatypeMenu = uimenu(editMenu,'Label','Data Type','Separator','off');

% Edit data type:
% c = curDataType(view)
% [dataTYPES(c),ok] = EditDataType(dataTYPES(c))
% if ok, saveSession(1), end;
cb=['c = ',view.name,'.curDataType; ',...
	'[dataTYPES(c), ok] = EditDataType(dataTYPES(c)); ',...
    'if ok, saveSession(1), end;'];
uimenu(datatypeMenu,'Label','Edit data type','Separator','off',...
   'CallBack',cb);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assign Visual Field Map Parameters submenu:
retinoMapMenu = uimenu(datatypeMenu, 'Label', 'Set Retinotopy Parameters...', ...
                        'Separator', 'off');

% submenus for assigning visual field map params
% retinoSetParams(view);
cb = sprintf('retinoSetParams(%s); ', view.name);
uimenu(retinoMapMenu, 'Label', 'Current Scan', 'Callback', cb);

% scans = er_selectScans(view); retinoSetParams(view, [], scans);
cb = [sprintf('scans = er_selectScans(%s); ', view.name) ...
      sprintf('retinoSetParams(%s, [], scans); ', view.name)];
uimenu(retinoMapMenu, 'Label', 'Select Scans', 'Callback', cb);

% retinoSetParams(view, [], 1:numScans(view));
cb = sprintf('retinoSetParams(%s, [], 1:numScans(%s)); ', ...
               view.name, view.name);
uimenu(retinoMapMenu, 'Label', 'All Scans', 'Callback', cb);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Edit Event-Related Parameters callback:
%   params = er_getParams(view);
%   params = er_editParams(params);
%   er_setParams(view, params);
cb = sprintf('params = er_getParams(%s); ', view.name);
cb = [cb 'params = er_editParams(params); '];
cb = [cb sprintf('er_setParams(%s, params); ', view.name)];
uimenu(datatypeMenu, 'Label', 'Edit Event-Related Parameters (cur scan)', ...
    'Separator', 'off', 'CallBack', cb);

% Duplicate data type:
%   duplicateDataType(view);
cb = ['duplicateDataType(',view.name,');'];
uimenu(datatypeMenu, 'Label', 'Duplicate current data type', ...
        'Separator', 'on', 'CallBack', cb);

% Group scans into data type:
%   groupScans(view));
cb = ['groupScans(',view.name,');'];
uimenu(datatypeMenu, 'Label', 'Group scans into data type', ...
    'Separator', 'off', 'CallBack',cb);

% Remove data type:
%   removeDataType(getDataTypeName(view));
cb = ['removeDataType(getDataTypeName(',view.name,'));'];
uimenu(datatypeMenu, 'Label', 'Remove data type', 'Separator', 'on',...
   'CallBack', cb);

% Remove last scan from data type:
%   removeScan(view, numScans(view), '', 1);
cb = sprintf('removeScan(%s, numScans(%s), '''', 1); ', view.name, view.name);
uimenu(datatypeMenu, 'Label', 'Remove last scan in data type', 'Separator', 'off',...
   'CallBack', cb);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Edit Preferences callback:
% editPreferences;
uimenu(editMenu, 'Label', 'Edit Preferences', 'Accelerator', ';', ...
    'Separator', 'off', 'CallBack','editPreferences;');

return


text:
function sftpack_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests R8VEC_SFTB and R8VEC_SFTF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
  n = 36;
  alo = 0.0;
  ahi = 5.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  For real slow Fourier transforms,\n' );
  fprintf ( 1, '  R8VEC_SFTF computes the forward transform.\n' );
  fprintf ( 1, '  R8VEC_SFTB computes the backward transform.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of data values, N = %d\n', n );

  seed = 123456789;

  [ x, seed ] = r8vec_uniform ( n, alo, ahi, seed );

  r8vec_print_part ( n, x, 10, '  The original data:' );
%
%  Compute the slow Fourier transform of the data.
%
  [ azero, a, b ] = r8vec_sftf ( n, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A (cosine) coefficients:\n' );
  fprintf ( 1, '\n' );

  fprintf ( 1, '  %4d  %f\n', 0, azero );

  for i = 1 : floor ( n / 2 )
   fprintf ( 1, '  %4d  %f\n', i, a(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  B (sine) coefficients:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : floor ( n / 2 )
    fprintf ( 1, '  %4d  %f\n', i, b(i) );
  end
%
%  Now try to retrieve the data from the coefficients.
%
  x = r8vec_sftb ( n, azero, a, b );

  r8vec_print_part ( n, x, 10, '  The retrieved data:' );

  return
end


text:
function y = abs( x )

%Disciplined convex/geometric programming information for ABS:
%   ABS(X) is convex and nonmonotonic in X. Therefore, when used in
%   DCPs, X must be affine. ABS(X) is not useful in DGPs, since all
%   log-convex and log-concave expressions are already positive.

%
% Determine the expression types
%

% 0 : convex, concave, invalid
% 1 : constant
% 2 : real affine
% 3 : complex affine
persistent remap
if isempty( remap ),
    remap_1 = cvx_remap( 'constant' );
    remap_2 = cvx_remap( 'real-affine' );
    remap_3 = cvx_remap( 'complex-affine' );
    remap_4 = cvx_remap( 'log-valid' );
    remap = remap_1 + ( 2 * remap_2 + 3 * remap_3 + 4 * remap_4 ) .* ~remap_1;
end
v = remap( cvx_classify( x ) );

%
% Process each type of expression one piece at a time
%

vu = sort( v(:) );
vu = vu([true;diff(vu)~=0]);
nv = length( vu );
sx = x.size_;
if nv ~= 1,
    y = cvx( sx, [] );
end
for k = 1 : nv,

    %
    % Select the category of expression to compute
    %

    vk = vu( k );
    if nv == 1,
        xt = x;
    else
        t = v == vk;
        xt = cvx_subsref( x, t );
    end

    %
    % Perform the computations
    %

    switch vk,
        case 0,
            % Invalid
            error( 'Disciplined convex programming error:\n    Illegal operation: abs( {%s} ).', cvx_class( xt ) );
        case 1,
            % Constant
            cvx_optval = cvx( builtin( 'abs', cvx_constant( xt ) ) );
        case 2,
            % Real affine
            st = size( xt );
            cvx_begin
                epigraph variable w( st )
                { xt, w } == lorentz( st, 0 );
            cvx_end
        case 3,
            % Complex affine
            st = size( xt );
            cvx_begin
                epigraph variable w( st )
                { xt, w } == complex_lorentz( st, 0 );
            cvx_end
        case 4,
            % log-affine, log-convex
            cvx_optval = xt;
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    %
    % Store the results
    %

    if nv == 1,
        y = cvx_optval;
    else
        y = cvx_subsasgn( y, 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.

repo: REST


text:
function updatelarvaspecies(hfly,hfly_extra,pos)

set(hfly,'XData',pos.xspine([1,6,11]),'YData',pos.yspine([1,6,11]));
%set(hfly,'XData',[pos.xcontour;nan;pos.xspine],'YData',[pos.ycontour;nan;pos.yspine]);
% xhead = pos.x + 2*pos.a*cos(pos.theta);
% yhead = pos.y + 2*pos.a*sin(pos.theta);
xhead = pos.xspine(1);
yhead = pos.yspine(1);
set(hfly_extra,'XData',xhead,'YData',yhead);


text:
function [dat] = read_nmc_archive_k_data(datafile, hdr, begsample, endsample, channelsel)

% READ_NMC_ARCHIVE_K_DATA reads data from nmc_archive_k datasets
%
% Used in read_data as
%   dat = read_nmc_archive_k_data(datafile, hdr, begsample, endsample, channelsel);
%
%
% This function specifically only reads data from one of the archived
% datasets of the Neurophysiological Mechanisms of Cognition group of
% Eric Maris, at the Donders Centre for Cognition, Radboud University,
% Nijmegen, the Netherlands. It should not be used for any other data
% format.
%
%
%

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

% Sanity checks (to be added)

% Getting data-directory out of data-filename and check
datadir = datafile(1:(strfind(datafile, 'eeg.noreref')+11));
if exist(datadir,'dir') ~= 7
    ft_error('no proper data-directory provided, please check your paths');
end

% Getting session name + path out of data-filename
sessionpath = datafile(1:(end-4));

% Build array containing channel file extensions to be used in reading the data
channelext = [];
% Remove the 'CH' from channel labels
chanlabnum = hdr.label;
for ichan = 1:length(chanlabnum)
    chanlabnum{ichan} = chanlabnum(3:end);
end
for ichan = 1:length(channelsel)
    if length(hdr.label{channelsel(ichan)}) == 3
        channelext{ichan} = ['.00' hdr.label{channelsel(ichan)}(3:end)];
    elseif length(hdr.label{channelsel(ichan)}) == 4
        channelext{ichan} = ['.0' hdr.label{channelsel(ichan)}(3:end)];
    elseif  length(hdr.label{channelsel(ichan)}) == 5
        channelext{ichan} = ['.' hdr.label{channelsel(ichan)}(3:end)];
    end
end

% Loop over channels while reading in data from file, 1 file per channel
dat = zeros(length(channelsel),(endsample-begsample+1));
for ichan = 1:length(channelsel)
    channelfile = [sessionpath channelext{ichan}];
    datafid = fopen_or_error(channelfile,'r','l');
    fseek(datafid,(hdr.nBytes*(begsample-1)),'bof');
    dat(ichan,:) = fread(datafid,(endsample-begsample+1),hdr.dataformat)';
    fclose(datafid);
end


text:
function MRI = out_mri_bst( MRI, MriFile )
% OUT_MRI_BST: Save a Brainstorm MRI structure.
%
% USAGE:  MRI = out_mri_bst( MRI, MriFile )
%
% INPUT:
%     - MRI     : Brainstorm MRI structure
%     - MriFile : full path to file where to save the MRI in brainstorm format
% OUTPUT:
%     - MRI : Modificed MRI structure
%
% NOTES:
%     - MRI structure:
%         |- Voxsize:   [x y z], size of each MRI voxel, in millimeters
%         |- Cube:      MRI volume
%         |- SCS:       Subject Coordinate System definition
%         |    |- NAS:    [x y z] coordinates of the nasion, in voxels
%         |    |- LPA:    [x y z] coordinates of the left pre-auricular point, in voxels
%         |    |- RPA:    [x y z] coordinates of the right pre-auricular point, in voxels
%         |    |- R:      Rotation to convert MRI coordinates -> SCS coordinatesd
%         |    |- T:      Translation to convert MRI coordinates -> SCS coordinatesd
%         |- Landmarks[]: Additional user-defined landmarks
%         |- NCS:         Normalized Coordinates System (Talairach, MNI, ...)
%         |    |- AC:             [x y z] coordinates of the Anterior Commissure, in voxels
%         |    |- PC:             [x y z] coordinates of the Posterior Commissure, in voxels
%         |    |- IH:             [x y z] coordinates of any Inter-Hemispheric point, in voxels
%         |- Comment:     MRI description

% @=============================================================================
% 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, 2008-2012

% ===== Clean-up MRI structure =====
% Remove (useless or old fieldnames)
Fields2BDeleted     = {'Origin','sag','ax','cor','hFiducials','header','filename'};
SCSFields2BDeleted  = {'Origin','Comment'};

for k = 1:length(Fields2BDeleted)
    if isfield(MRI,Fields2BDeleted{k})
        MRI = rmfield(MRI,Fields2BDeleted{k});
    end
end

if isfield(MRI,'SCS')
    for k = 1:length(SCSFields2BDeleted)
        if isfield(MRI.SCS,SCSFields2BDeleted{k})
            MRI.SCS = rmfield(MRI.SCS,SCSFields2BDeleted{k});
        end
    end
end

if isfield(MRI,'Landmarks') && ~isempty(MRI.Landmarks)
    tmpLandmarks = MRI.Landmarks;
    nLandmarks2Remove = 0;
    if isfield(MRI,'SCS')
        nLandmarks2Remove = 4; %Remove SCS fiducials from Landmark list
    end
    if isfield(MRI,'talCS')
        nLandmarks2Remove = nLandmarks2Remove + 3; % Remove TAL/MNI fiducials
    end
    if isfield(MRI.Landmarks, 'Names') && ~isempty(MRI.Landmarks.Names)
        MRI.Landmarks.Names    = MRI.Landmarks.Names(nLandmarks2Remove+1:end);
    end
    if isfield(MRI.Landmarks, 'MRImmXYZ') && ~isempty(MRI.Landmarks.MRImmXYZ)
        MRI.Landmarks.MRImmXYZ = MRI.Landmarks.MRImmXYZ(:,nLandmarks2Remove+1:end);
    end
    if isfield(MRI.Landmarks, 'Handles')
        MRI.Landmarks = rmfield(MRI.Landmarks,'Handles');
    end
end

if isfield(MRI,'SCS2Landmarks') % don't need to be stored in file
    tmpSCS2Landmarks = MRI.SCS2Landmarks;
    MRI = rmfield(MRI,'SCS2Landmarks');
end

% Remove FileName field
if isfield(MRI,'FileName')
    MRI = rmfield(MRI, 'FileName');
end

% SAVE .mat file
try
    bst_save(MriFile, MRI, 'v7');
catch

end


text:
function [grade,X,SR,improve,AppSet] = LocalSearch3(Problem,X,SR,improve,AppSet)
% Local Search 3 of MTS

%------------------------------- 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".
%--------------------------------------------------------------------------

    grade   = 0;
    SearchL = Problem.lower;
    SearchU = Problem.upper;
    Disp    = (SearchU-SearchL)/10;
    Best    = X;
    while any(Disp<1e-2)
        for i = randperm(length(SR))
            value = SearchL(i) : Disp(i) : SearchU(i);
            Decs  = repmat(Best.dec,length(value),1);
            Decs(:,i) = value;
            Y     = Problem.Evaluation(Decs);
            for j = 1 : length(Y)
                % 'grade', 'SR' and 'improve' will not be updated in local
                % search 3
                [~,~,AppSet] = Grading(Y(j),Best,grade,improve,AppSet,Problem.N);
                if all(Y(j).obj<=Best.obj)
                    Best = Y(j);
                end
            end
            SearchL(i) = max(Best.dec(i)-2*Disp(i),Problem.lower);
            SearchU(i) = min(Best.dec(i)+2*Disp(i),Problem.upper);
            Disp(i)    = (SearchU(i)-SearchL(i))/10;
        end
    end
    X = Best;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% CS294A/CS294W Linear Decoder Exercise

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear decoder exericse. For this exercise, you will only need to modify
%  the code in sparseAutoencoderLinearCost.m. You will not need to modify
%  any code in this file.

addpath '../library/'
addpath '../library/minFunc/'

%%======================================================================
%% STEP 0: Initialization
%  Here we initialize some parameters used for the exercise.

imageChannels = 3;     % number of channels (rgb, so 3)

patchDim = 8;          % patch dimension
numPatches = 100000;   % number of patches

visibleSize = patchDim * patchDim * imageChannels;  % number of input units
outputSize = visibleSize;   % number of output units
hiddenSize = 400;           % number of hidden units

sparsityParam = 0.035; % desired average activation of the hidden units.
lambda = 3e-3;         % weight decay parameter
beta = 5;              % weight of sparsity penalty term

epsilon = 0.1;	       % epsilon for ZCA whitening

% %%======================================================================
% %% STEP 1: Create and modify sparseAutoencoderLinearCost.m to use a linear decoder,
% %          and check gradients
% %  You should copy sparseAutoencoderCost.m from your earlier exercise
% %  and rename it to sparseAutoencoderLinearCost.m.
% %  Then you need to rename the function from sparseAutoencoderCost to
% %  sparseAutoencoderLinearCost, and modify it so that the sparse autoencoder
% %  uses a linear decoder instead. Once that is done, you should check
% % your gradients to verify that they are correct.
%
% % NOTE: Modify sparseAutoencoderCost first!
%
% % To speed up gradient checking, we will use a reduced network and some
% % dummy patches
% hiddenSize = 5;
% visibleSize = 8;
% patches = rand([8 10]);
% theta = initializeParameters(hiddenSize, visibleSize);
%
% [cost, grad] = sparseAutoencoderLinearCost(theta, visibleSize, hiddenSize, ...
%                                            lambda, sparsityParam, beta, ...
%                                            patches);
%
% % Check gradients
% numGrad = computeNumericalGradient( @(x) sparseAutoencoderLinearCost(x, visibleSize, hiddenSize, ...
%                                                   lambda, sparsityParam, beta, ...
%                                                   patches), theta);
%
% % Use this to visually compare the gradients side by side
% disp([numGrad grad]);
%
% diff = norm(numGrad-grad)/norm(numGrad+grad);
% % Should be small. In our implementation, these values are usually less than 1e-9.
% disp(diff);
%
% assert(diff < 1e-9, 'Difference too large. Check your gradient computation again');
%
% % NOTE: Once your gradients check out, you should run step 0 again to
% %       reinitialize the parameters
% %}

%%======================================================================
%% STEP 2: Learn features on small patches
%  In this step, you will use your sparse autoencoder (which now uses a
%  linear decoder) to learn features on small patches sampled from related
%  images.

%% STEP 2a: Load patches
%  In this step, we load 100k patches sampled from the STL10 dataset and
%  visualize them.

patches = load('../data/STL10Patches100k.mat');
patches = patches.patches;

displayColorNetwork(patches(:, 1:100));

%% STEP 2b: Apply preprocessing
%  In this sub-step, we preprocess the sampled patches, in particular,
%  ZCA whitening them.
%
%  In a later exercise on convolution and pooling, you will need to replicate
%  exactly the preprocessing steps you apply to these patches before
%  using the autoencoder to learn features on them. Hence, we will save the
%  ZCA whitening and mean image matrices together with the learned features
%  later on.

% Scale data to range [0, 1]
patches = patches / 255;

% Subtract mean patch (hence zeroing the mean of the patches)
meanPatch = mean(patches, 2);
patches = bsxfun(@minus, patches, meanPatch);

% Apply ZCA whitening
sigma = patches * patches' / numPatches;
[u, s, v] = svd(sigma);
ZCAWhite = u * diag(1 ./ (diag(s) + epsilon)) * u';
patches = ZCAWhite * patches;

displayColorNetwork(patches(:, 1:100));

%% STEP 2c: Learn features
%  You will now use your sparse autoencoder (with linear decoder) to learn
%  features on the preprocessed patches. This should take around 45 minutes.

theta = initializeParameters(hiddenSize, visibleSize);

% Use minFunc to minimize the function

clear('options');
options.Method = 'lbfgs';
options.maxIter = 400;
options.display = 'on';

[optTheta, cost] = minFunc( @(p) sparseAutoencoderLinearCost(p, ...
                                   visibleSize, hiddenSize, ...
                                   lambda, sparsityParam, ...
                                   beta, patches), ...
                              theta, options);

% Save the learned features and the preprocessing matrices for use in
% the later exercise on convolution and pooling
fprintf('Saving learned features and preprocessing matrices...\n');
save('STL10Features.mat', 'optTheta', 'ZCAWhite', 'meanPatch');
fprintf('Saved\n');

%% STEP 2d: Visualize learned features

W = reshape(optTheta(1:visibleSize * hiddenSize), hiddenSize, visibleSize);
b = optTheta(2*hiddenSize*visibleSize+1:2*hiddenSize*visibleSize+hiddenSize);
displayColorNetwork(W');

author: zellyn


text:
function VAR = RenameClusGroup(VAR,i_fish,i_ClusGroup,name)
VAR(i_fish).ClusGroupName{i_ClusGroup} = name;
end

author: xiuyechen

repo: FishExplorer


text:
function b=bitstream(l,p)

% b=bitstream(l,p)
%
% Return a stream of l bits where the probability of a '1' is p.
% p defaults to 0.5.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

error(nargchk(1,2,nargin));
error(chk_param(l,'l','scalar','integer','real','>=',0));
if (nargin<2)
  p=0.5;
end
error(chk_param(p,'p','scalar','real','>=',0,'<=',1));

b=double(rand(1,l)<p);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%

classdef UMAP_extra_results < handle
    properties
        supervisorMatchedLabels;
                        %vector of labels matched from a supervised template
                        % 1 value per row of input data matrix
        fig;            % main figure with UMAP output plot
        qft;            % instance of qf tree object for input data.
                        % The fig property contains the view of the tree.
        qftSupervisors; % instance of qf tree object for supervisor data.
                        % The fig property
                        % contains the view of the tree.
        qfd={};         % array qf QFTable instances.
                        % The function getAverages summarizes the result.
                        % Relevant figure properties are fig for table view,
                        % qHistFig for view of dismilarity score histogram
                        % fHistFig for view of f measure histogram
       matchHtmlHead1;
       matchHtmlHead2;
       matchHtmlHead3;
       matchHtmlBody;
       matchCsvHead;
       matchCsvBody;
       falsePosNegFile;
       falsePosNegHead;
       falsePosNegBody='';
       app;
       htmlFile;
       csvFile;
       args;
    end

    methods

        function closeMatchFigs(this)
            N=length(this.qfd);
            for i=1:N
                this.qfd{i}.closeFigs;
            end
        end

        function closeTreeFigs(this)
            if ~isempty(this.qft)
                Gui.CloseFig(this.qft.fig);
            end
            if ~isempty(this.qftSupervisors)
                Gui.CloseFig(this.qftSupervisors.fig);
            end
        end

        function avgs=getMatchAverages(this)
            N=length(this.qfd);
            avgs=[];
            for i=1:N
                qf=this.qfd{i};
                [medianDissimilarity, meanDissimilarity, ...
                    medianOverlap, meanOverlap]=qf.getAverages;
                avgs(end+1,:)=[medianDissimilarity meanDissimilarity...
                    medianOverlap meanOverlap];
            end
        end

        function fig=seeFalsePosNeg(this)
            if ~isempty(this.falsePosNegFile) && ~isempty(this.falsePosNegHead)
                File.SaveTextFile(this.falsePosNegFile, ...
                    [this.falsePosNegHead this.falsePosNegBody])
                [fig, that]=FalsePositiveNegative.Plot([0 1], ...
                    this.falsePosNegFile);
                that.fcnMoreHtml=@()getFalsePosNegMatrixHtml(this);
            else
                fig=[];
            end
        end

        function saveMatchFiles(this, h1)
            if nargin<2
                h1='';
            end
            if ~isempty(this.csvFile) && ~isempty(this.matchCsvHead)
                File.SaveTextFile(this.csvFile, [this.matchCsvHead ...
                    this.matchCsvBody])
            end
            if ~isempty(this.falsePosNegFile) && ~isempty(this.falsePosNegHead)
                File.SaveTextFile(this.falsePosNegFile, ...
                    [this.falsePosNegHead this.falsePosNegBody])
            end
            if ~isempty(this.htmlFile) && ~isempty(this.matchHtmlHead1)
                File.SaveTextFile(this.htmlFile, this.getMatchHtml(h1));
            end
        end

        function seeMatches(this, how, h1)
            if nargin<3
                h1='';
            end
            if how==-1
                msg(this.getMatchHtml(h1), 0, 'east++', ['Matches: '...
                    this.describeRedution ' reduction'], 'none')
            elseif how==1
                Html.BrowseString(this.getMatchHtml(h1));
            elseif how==2
                Html.BrowseFile(this.htmlFile);
            end
        end

        function this=UMAP_extra_results()
            this.app=BasicMap.Global;
        end

        function doMatchOutput(this, D)
            if ~isempty(this.qfd)
                this.doMatchHtmlHead(D);
                this.doMatchHtmlBody;
                this.doMatchCsv;
                this.doFalsePosNeg;
            end
        end

        function htmls=getFalsePosNegMatrixHtml(this)
            args_=this.args;
            qfds=this.qfd;
            N=length(qfds);
            if N<1 || isempty(find(args_.match_scenarios==4, 1))
                return;
            end
            reduction=args_.reduction;
            if ~isempty(args_.training_set)
                trainingSet=args_.training_set;
            elseif ischar(args_.template_file)
                [~,fn]=fileparts(args_.template_file);
                trainingSet=fn;
            else
                trainingSet='';
            end
            if ~isempty(args_.test_set)
                testSet=args_.test_set;
            elseif ischar(args_.csv_file_or_data)
                [sampleFldr,fn]=fileparts(args_.csv_file_or_data);
                testSet=fn;
            else
                testSet='';
            end
            if ~isempty(args_.sample_set)
                sampleSet=args_.sample_set;
            elseif exist('sampleFldr', 'var')
                [~, sampleSet]=fileparts(sampleFldr);
            else
                sampleSet='';
            end
            htmls=cell(1,N);
            for i=1:N
                qt=qfds{i};
                N2=length(qt.qf.falsePosNegs);
                if N2>0
                    mt=qt.context.matchType;
                    cd=qt.context.clusterDetail;
                    cd=StringArray.IndexOf(Density.DETAILS, cd);
                    htmls{i}=FalsePositiveNegative.MatrixHtml(qt.qf);
                end
            end
        end

        function doFalsePosNeg(this, doPlotNow)
            args_=this.args;
            qfds=this.qfd;
            N=length(qfds);
            if N<1 || isempty(find(args_.match_scenarios==4, 1))
                return;
            end
            fullBody=this.falsePosNegBody;
            this.falsePosNegFile=[UMAP_extra_results.FileName(...
                args_, 'falsePosNeg') '.txt'];
            reduction=args_.reduction;
            if ~isempty(args_.training_set)
                trainingSet=args_.training_set;
            elseif ischar(args_.template_file)
                [~,fn]=fileparts(args_.template_file);
                trainingSet=fn;
            else
                trainingSet='';
            end
            if ~isempty(args_.test_set)
                testSet=args_.test_set;
            elseif ischar(args_.csv_file_or_data)
                [sampleFldr,fn]=fileparts(args_.csv_file_or_data);
                testSet=fn;
            else
                testSet='';
            end
            if ~isempty(args_.sample_set)
                sampleSet=args_.sample_set;
            elseif exist('sampleFldr', 'var')
                [~, sampleSet]=fileparts(sampleFldr);
            else
                sampleSet='';
            end
            for i=1:N
                qt=qfds{i};
                N2=length(qt.qf.falsePosNegs);
                if N2>0
                    mt=qt.context.matchType;
                    cd=qt.context.clusterDetail;
                    cd=StringArray.IndexOf(Density.DETAILS, cd);
                    [body, notFound]=FalsePositiveNegative.TabRows(...
                        qt.qf.falsePosNegs, reduction, sampleSet, ...
                        trainingSet, testSet, args_.n_neighbors, ...
                        args_.hiD, args_.n_components, mt, cd);
                    if ~isempty(notFound)
                        body=[body newline notFound];
                    end
                    fullBody=[fullBody body newline];
                end
            end
            this.falsePosNegHead=FalsePositiveNegative.TabHead;
            this.falsePosNegBody=fullBody;
            if nargin>1 && doPlotNow
                this.seeFalsePosNeg;
            end
        end

        function html=getMatchHtml(this, h1)
            if nargin<2
                h1='';
            end
            html=['' h1 '

' ...
                '
' this.matchHtmlHead1 '
' ...
                '
' this.matchHtmlHead2 '
' ...
                '
' this.matchHtmlHead3 '

'...
                '
' this.matchHtmlBody '
'...
                '
'];
        end

        function s=describeRedution(this)
            s=UmapUtil.GetReductionLongText(this.args.reductionType);
        end

        function records=getFalsePosNeg(this, verbose)
            N=length(this.qfd);
            records={};
            for i=1:N
                qf=this.qfd{i};
                N2=length(qf.qf.falsePosNegs);
                if N2>0
                    if nargin>1 && verbose
                        for j=1:N2
                            record=qf.qf.falsePosNegs{j};
                            fprintf(['#%d "%s"; false pos=%s '...
                                'false neg=%s; match 1 x %d\n'], ...
                                j, record.trainingClass, ...
                                String.encodePercent(record.falsePos, ...
                                record.testSize, 2), ...
                                String.encodePercent(record.falseNeg, ...
                                record.trainingSize, 2), ...
                                length(record.testIds));
                        end
                    end
                    records=[records qf.qf.falsePosNegs];
                end
            end
        end
    end

    methods(Access=private)
        function doMatchHtmlHead(this, D)
            args_=this.args;
            nMatchTypes=length(args_.match_supervisors);
            nMatchScenarios=length(args_.match_scenarios);
            if any(args_.match_scenarios==1)
                nUstMatches=nMatchScenarios-1;
            else
                nUstMatches=nMatchScenarios;
            end
            this.htmlFile=[UMAP_extra_results.FileName(args_) '.html'] ;
            colspan=['colspan="' num2str(nUstMatches*3) '"'];
            if isempty(args_.description)
                htmlDsc='	';
            else
                htmlDsc='Context	';
            end
            sm1=this.app.smallStart;
            sm2=this.app.smallEnd;
            html1='';
            html2='';
            html3='';
            has1=false;
            cluDtls=args_.cluster_detail;
            nCluDtls=length(cluDtls);
            if length(args_.ust_test_components)>1
                loD=[];
            else
                loD=args_.n_components;
            end
            for c=1:nCluDtls
                for i=1:nMatchTypes
                    matchType=args_.match_supervisors(i);
                    if c>1
                        if matchType>=3 % nearest neighbor no clustering
                            continue;
                        end
                    end
                    mt=UmapUtil.GetMatchTypeLongText(matchType, ...
                        args_.reductionType, loD, D);
                    if matchType<3
                        mt=[mt ' ' cluDtls{c}];
                    end
                    html1=[html1 '' mt '	'];
                    for j=1:nMatchScenarios
                        scenario=args_.match_scenarios(j);
                        sc=UmapUtil.GetMatchScenarioText(scenario, args_.reductionType) ;
                        if scenario==1
                            if ~has1
                                has1=true;
                                html1_1=['Prior
'...
                                    sm1 'classifications' sm2 '	' ];
                                html2_1=['' sc '
' ...
                                    sm1 'Dissimilarity' sm2 '	'];
                                html3_1='Unmatched	Median	Mean	';
                            end
                        else
                            if scenario==4
                                htmlStat=['
' sm1 'Overlap' sm2];
                            else
                                htmlStat=['
' sm1 'Dissimilarity' sm2];
                            end
                            html2=[html2 '' sc htmlStat '	' ];
                            html3=[html3 'Unmatched	Median	Mean	'];
                        end
                    end
                end
            end
            if has1
                if nUstMatches>0
                    this.matchHtmlHead1=[htmlDsc html1_1  html1];
                    this.matchHtmlHead2=['	' html2_1 html2];
                    this.matchHtmlHead3=['	' html3_1 html3];
                else
                    this.matchHtmlHead1=[htmlDsc html1_1];
                    this.matchHtmlHead2=['	' html2_1];
                    this.matchHtmlHead3=['	' html3_1];
                end
            else
                this.matchHtmlHead1=[htmlDsc html1];
                this.matchHtmlHead2=['	'  html2];
                this.matchHtmlHead3=['	' html3];
            end
        end

        function doMatchHtmlBody(this)
            args_=this.args;
            qfds=this.qfd;
            dsc=args_.description;
            td='';
            td_='	';
            N_=length(qfds);
            htmlDsc=['' dsc td_];
            html='';
            mdn=zeros(1, N_);
            mn=zeros(1,N_);
            adjMdn=zeros(1, N_);
            adjMn=zeros(1,N_);
            matchRate=zeros(1,N_);
            matches=cell(1, N_);
            scenarios=zeros(1,N_);
            isDis=true(1, N_);
            for i=1:N_
                qf=qfds{i};
                matchStrategy=qf.context.matchStrategy;
                [tUnmatched, tN, sUnmatched, sN]=qf.getMatchCounts;
                matches{i}=sprintf('%d/%d, %d/%d', tUnmatched, tN,...
                    sUnmatched, sN);
                if matchStrategy==1
                    scenarios(i)=qf.context.matchScenario;
                    [~,mdn(i), mn(i)]=qf.getData(true);
                else
                    scenarios(i)=4;
                    [~,mdn(i), mn(i)]=qf.getData(false);
                    isDis(i)=false;
                end
                tM=tN-tUnmatched;
                sM=sN-sUnmatched;
                matchedSubsetCnt=tM+sM;
                subsetCnt=tN+sN;

                if tUnmatched>0 || sUnmatched>0
                    if matchStrategy==1
                        addDs=(tUnmatched+sUnmatched)*100;
                        adjMdn(i)=(matchedSubsetCnt*mdn(i)+addDs)/subsetCnt;
                        adjMn(i)=(matchedSubsetCnt*mn(i)+addDs)/subsetCnt;
                    else
                        adjMdn(i)=matchedSubsetCnt*mdn(i)/subsetCnt;
                        adjMn(i)=matchedSubsetCnt*mn(i)/subsetCnt;
                    end
                else
                    adjMdn(i)=mdn(i);
                    adjMn(i)=mn(i);
                end
                matchRate(i)=matchedSubsetCnt/subsetCnt;
            end
            usc=unique(scenarios);
            nSc=length(usc);
            bestMedian=zeros(1, 4);
            bestMean=zeros(1, 4);
            bestMatches=zeros(1,4);
            worstMatches=zeros(1,4);
            for i=1:nSc
                u=usc(i);
                if u==4
                    bestMedian(u)=max(adjMdn(scenarios==u));
                    bestMean(u)=max(adjMn(scenarios==u));
                else
                    bestMedian(u)=min(adjMdn(scenarios==u));
                    bestMean(u)=min(adjMn(scenarios==u));
                end
                bestMatches(u)=max(matchRate(scenarios==u));
                worstMatches(u)=min(matchRate(scenarios==u));
            end
            has1=false;
            for i=1:N_
                scenario=scenarios(i);
                sMdn=String.encodeRounded(mdn(i), 1, true);
                sMn=String.encodeRounded(mn(i), 1, true);
                match=matches{i};
                ms1='';
                ms2='';

                if scenario>1
                    if matchRate(i)==bestMatches(scenario)
                        ms1='';
                        ms2='';
                    elseif matchRate(i)==worstMatches(scenario)
                        ms1='';
                        ms2='';
                    end
                end
                if scenario>1 && adjMdn(i)==bestMedian(scenario)
                    s1='';
                    s2='';
                else
                    s1='';
                    s2='';
                end
                if scenario==1
                    html1=[td ms1 match ms2 td_ td s1 sMdn  s2 td_];
                else
                    html=[html td ms1 match ms2 td_ td s1 sMdn  s2 td_];
                end
                if scenario>1 && adjMn(i)==bestMean(scenario)
                    s1='';
                    s2='';
                else
                    s1='';
                    s2='';
                end
                if scenario==1
                    has1=true;
                    html_1=[html1 td s1 sMn s2 td_ ];
                else
                    html=[html td s1 sMn s2 td_ ];
                end
            end
            if has1
                html=[htmlDsc html_1 html];
            else
                html=[htmlDsc html];
            end
            this.matchHtmlBody=html;
        end

        function doMatchCsv(this)
            args_=this.args;
            qfds=this.qfd;
            N=length(qfds);
            if N>0
                this.csvFile=[UMAP_extra_results.FileName(args_) '.csv'] ;
                lf=newline;
                this.matchCsvHead=[ StringArray.toString(fieldnames(...
                    qfds{1}.context), ',') 'trainingUnmatched,'...
                    'trainingSubsets,testUnMatched,testSubsets,'...
                    'median,mean' lf];
                csv='';
                for i=1:N
                    qf=qfds{i};
                    [tUnmatched, tN, sUnmatched, sN]=qf.getMatchCounts;

                    if qf.context.matchStrategy==1
                        [~,mdn, mn]=qf.getData(true);
                    else
                        [~,mdn, mn]=qf.getData(false);
                    end
                    csv=[csv String.toString(qf.context) ...
                        num2str(tUnmatched) ',' num2str(tN) ',' ...
                        num2str(sUnmatched) ',' num2str(sN) ',' ...
                        num2str(mdn) ',' num2str(mn) ',' lf];
                end
                this.matchCsvBody=csv;
            end
        end
    end

    methods(Static)
        function file=FileName(args, prefix)
            if nargin<2
                prefix='match';
            end
            if isfield(args, [prefix '_file'])
                fn=getfield(args, [prefix '_file']);
            else
                fn='';
            end
            if isempty(fn)
                strMatchTypes=strrep(strrep(num2str(...
                    args.match_supervisors), ' ', '_'), '__', '_');
                strMatchScenarios=strrep(strrep(num2str(...
                    args.match_scenarios), ' ', '_'), '__', '_');
                fn=[prefix '_' strMatchTypes...
                    '_for_' strMatchScenarios ];
            end
            file=fullfile(args.result_folder, fn);
            File.mkDir(args.result_folder);
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Test file for singfun/feval.m

function pass = test_feval(pref)

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

% Generate a few random points to use as test values.
seedRNG(786);
x = -1 + 2*rand(100, 1);

%%
% Check feval on empty set of points
data.exponents = [0, 0];
data.singType = {'none', 'none'};
f = singfun(@(x) x, data, pref);
pass(1) = isempty(feval(f, []));
%%
% Check feval on a SINGFUN without exponents
fh = @(x) sin(cos(10*x.^2));
data.exponents = [0, 0];
data.singType = {'none', 'none'};
f = singfun(fh, data, pref);
pass(2) = norm(feval(f,x) - feval(fh,x), inf) < 1e1*eps;

%%
% Check feval on a SINGFUN with negative exponents
a = 1 + rand();
b = 1 + rand();
fh = @(x) sin(cos(10*x.^2))./((1+x).^a.*(1-x).^b);
data.exponents = [-a, -b];
data.singType = {'sing', 'sing'};
f = singfun(fh, data, pref);
err = norm(feval(f,x) - feval(fh,x), inf);
tol = 2e3*eps;
pass(3) = err < tol;

%%
% Check feval on a SINGFUN with positive exponents
a = rand();
b = rand();
fh = @(x) sin(cos(10*x.^2)).*(1+x).^a.*(1-x).^b;
data.exponents = [a, b];
data.singType = {'root', 'root'};
f = singfun(fh, data, pref);
pass(4) = norm(feval(f,x) - feval(fh,x), inf) < 1e1*eps;

end


text:
function linpack_d_test10 ( )

%*****************************************************************************80
%
%% TEST10 tests DGEFA and DGESL.
%
%  Discussion:
%
%    Solve A*x = b where A is a given matrix, and B a right hand side.
%
%    We will also assume that A is stored in the simplest
%    possible way.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST10\n' );
  fprintf ( 1, '  For a general matrix,\n' );
  fprintf ( 1, '  DGEFA computes the LU factors;\n' );
  fprintf ( 1, '  DGESL solves a factored linear system;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of equations is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  a = [ 1.0, 2.0, 3.0;
        4.0, 5.0, 6.0;
        7.0, 8.0, 0.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %12f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Set the values of the right hand side vector B.
%
  b(1:3) = [ 6.0, 15.0, 15.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The right hand side B is\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12f\n', b(i) );
  end
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix\n' );

  [ a, ipvt, info ] = dgefa ( a, lda, n );

  if ( info ~= 0 )
    fprintf ( 1, '  DGEFA returned an error flag INFO = %d\n', info );
    return
  end
%
%  If no error occurred, now use DGESL to solve the system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  job = 0;
  b = dgesl ( a, lda, n, ipvt, b, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  DGESL returns the solution:\n' );
  fprintf ( 1, '  (Should be (1,1,1))\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %12f\n', b(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F=FPolarCoordTurn2D(T,x,turnType,discPoint,tauTurn,tauLinAccel)
%%FPOLARCOORDTURN2D The state transition matrix for a two-dimensional
%              coordinated turn model where the velocity is specified in
%              terms of a heading and a speed. The turn rate can be
%              specified in terms of a turn rate in radians per second, or
%              in terms of a transversal acceleration. Additionally, a
%              linear acceleration can be given. The turn rate and linear
%              acceleration can optionally have time constants associated
%              with them, like in the Singer model, modelling a tendancy to
%              eventually want to return to non-accelerating, straight-line
%              motion.
%
%INPUTS: T The time-duration of the propagation interval in seconds.
%        x The target state for 2D motion where the velocity is given in
%          terms of heading and speed components. If there is no linear
%          acceleration (acceleration along the direction of motion), then
%          x can either be x=[x;y;h;v;omega], where h is the heading in
%          terms of radians counterclockwise from the x-axis, v is the
%          speed, and omega is the turn rate (the derivative of h with
%          respect to time) or  x=[x;y;h;v;at] where at is the transversal
%          acceleration, which is orthogonal to the velocity and is defined
%          such that positive values of at map to positive values of omega.
%          If there is a linear acceleration, then the target state is
%          either x=[x;y;h;v;omega;al] where omega is the turn rate and al
%          is the linear acceleration or the target state is
%          x=[x;y;h;v;at;al] if the turn is expressed in terms of a
%          transversal acceleration. The dimensionality of the state is
%          used to determine whether a linear acceleration component is
%          present. The linear acceleration component changes the speed.
%          That means that it is the derivative of the speed.
% turnType A string specifying whether the turn is given in terms of a
%          turn rate in radians per second or a transversal acceleration in
%          m/s^2. Possible values are
%          'TurnRate'   The turn is specified in terms of a turn rate
%                       (The default if this parameter is omitted).
%          'TransAccel' The turn is specified in terms of a transversal
%                       acceleration.
% discPoint This optional parameter specified what value of the turn rate
%          is used for the discretized state prediction. The three possible
%          values were suggested in Li's paper, [3]. Possible values are:
%          0 (The default if omitted) use omega=x(5), or the equivalent
%            value when specifying the turn rate using a transverse
%            acceleration, from the non-predicted target state for
%            building the state transition matrix. \omega_k
%          1 Use the average value of the predicted omega (or the average
%            value of the transverse acceleration) over the interval T
%            for building the state transition matrix. \bar{\omega}
%          2 Use the approximate average value of the predicted omega
%            over the interval T for building the state transition
%            matrix. \bar{\omega} This is the suggestion of using half
%            the prior and prediction, as was given in Li's paper. When
%            given a transverse acceleration instead of a turn rate, half
%            of the prior and predicted accelerations is used.
%          3 Use the forward-predicted omega (or transverse acceleration)
%            for building the state transition matrix. \omega_{k+1}
%  tauTurn The correlation time constant for the turn rate in seconds.
%          tau must be positive but does not have to be finite. If this
%          parameter is omitted, then tauTurn is set to infinity.
% tauLinAccel The correlation time constant for the linear acceleration (if
%          present) in seconds. This parameter is not needed if there is
%          no linear acceleration. If a linear acceleration is present
%          and this parameter is omitted, then tauLinAccel is set to
%          infinity.
%
%OUTPUT: F The state transition matrix under a possibly linearlly
%          accelerating coordinated turn model where the velocity is given
%          by a heading and a speed.
%
%The state transition matrix is, in part, a realization of the unnumered
%transition Equation after equation 1b in [1], which comes from the initial
%derivation in [2], where a different definition of heading direction is
%used. Note that the order of the components in the state here is different
%than the ordering of the components in Blackman's paper. Also, Blackman's
%paper does not consider linear acceleration terms.
%
%Despite the form of the transition matrix without a linear acceleration
%term, the result is mathematically equivalent to equation 75 in [3]. After
%being modified to account for a decaying turn rate. It is assumed that the
%continuous-time turn rate model is
%omegaDot=-(1/tauTurn)*omega+noise, which discretizes to
%omega[k+1]=exp(-T/tauTurn)*omega[k]+noise.
%Analogously, the optional linear acceleration discretizes to
%al[k+1]=exp(-T/tauLinAccel)*al[k]+noise
%
%More information on discrete-time turning models is given in the comments
%to the function FCoordTurn2D.
%
%This state transition matrix goes with the process noise covariance matrix
%given by QPolarCoordTurn2D.The corresponding continuous-time drift
%function are aPolarCoordTurn2DOmega and aPolarCoordTurn2Dtrans with the
%diffusion matrix DPolarCoordTurn2D. Note that this model is a direct-
%discrete-time model and is not just a discretization of the continuous-
%time model.
%
%REFERENCES:
%[1] M. Busch and S. Blackman, "Evaluation of IMM filtering for an air
%    defense system application," in Proceedings of SPIE: Signal and Data
%    Processing of Small Targets, vol. 2561, 9 Jul. 1995, pp. 435-447.
%[2] J. L. Gertz, "Multisensor surveillance for improved aircraft
%    tracking," The Lincoln Laboratory Journal, vol. 2, no. 3, pp. 381-396,
%    1989.
%[3] X. R. Li and V. P. Jilkov, "Survey of maneuvering target tracking.
%    Part I: Dynamic models," IEEE Transactions on Aerospace and Electronic
%    Systems, vol. 39, no. 4, pp. 1333-1364, Oct. 2003.
%
%July 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The time constant for the linear acceleration.
if(nargin<6)
    tauLinAccel=Inf;
end

%The time constant for the turn.
if(nargin<5)
   tauTurn=Inf;
end

if(nargin<4)
    discPoint=0;
end

if(nargin<3)
    turnType='TurnRate';
end

beta=exp(-T/tauTurn);
switch(discPoint)
    case 0%Use the prior value of at/omega for the linearization.
        turnVal=x(5);
    case 1
        %Use the average value of at.omega over the prediction
        %interval for the linearization.
        tau=param3;
        turnVal0=x(5);
        turnVal=(tau-beta*tau+T*turnVal0)/T;

        %Assume tau=Inf or T=0 was used, in which case the
        %asymptotic average value of at/omega should be used.
        if(~isfinite(tauTurn)||~isfinite(turnVal))
            turnVal=turnVal0;
        end
    case 2%Use the simple mean of at/omega for the linearization.
        turnVal=(beta+1)*x(5)/T;
    case 3%Use the predicted value of at/omega for the linearization.
        turnVal=beta*x(5);
    otherwise
        error('Invalid value entered for discPoint');
end

v=x(4);%The speed

switch(turnType)
    case 'TransAccel'%The turn is expressed in terms of a transverse
                     %acceleration.

        omega=turnVal/v;
        if(~isfinite(omega))%Deal with zero velocity.
            omega=0;
        end
    case 'TurnRate'%The turn is expressed in terms of a turn rate.
        omega=turnVal;
    otherwise
        error('Unknown turn type specified.');
end
%We now have the omega term for the turn.

theta=x(3);%The heading (counterclockwise from the x axis).

sinHead=sin(theta);
cosHead=cos(theta);

sinVal=sin(omega*T);
cosVal=cos(omega*T);

sinRat=sinVal/omega;
if(~isfinite(sinRat))
    sinRat=T;%The limit as omega goes to zero.
end
cosRat=(1-cosVal)/omega;
if(~isfinite(cosRat))
    cosRat=0;%The limit as omega goes to zero.
end

switch(turnType)
    case 'TransAccel'%The turn is expressed in terms of a transverse
        F=[1,0,0,cosHead*sinRat-sinHead*cosRat,0;%The x-component.
            0,1,0,sinHead*sinRat+cosHead*cosRat,0;%The y-component.
            0,0,1,0,                            T/v;%The heading component.
            0,0,0,1,                            0;%The speed component
            0,0,0,0,                            beta];%The transverse accel component.
    case 'TurnRate'
        F=[1,0,0,cosHead*sinRat-sinHead*cosRat,0;%The x-component.
            0,1,0,sinHead*sinRat+cosHead*cosRat,0;%The y-component.
            0,0,1,0,                            T;%The heading component.
            0,0,0,1,                            0;%The speed component
            0,0,0,0,                            beta];%The turn rate component.
end
if length(x)==6
    betaLinAccel=exp(-T/tauLinAccel);
    F(4,6)=T;
    F(6,6)=betaLinAccel;
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