text:
function [p,t] = initmesh(fd,fh,h0,bbox,pfix,varargin)
%% INITMESH initialize a triangulation
%
% Modifed from distmesh. Please refer to distmesh on the usage.
%
%   Example: (Rectangle with circular hole, refined at circle boundary)
%      fd=inline('ddiff(drectangle(p,-1,1,-1,1),dcircle(p,0,0,0.5))','p');
%      fh=inline('min(4*sqrt(sum(p.^2,2))-1,2)','p');
%      [p,t]=initmesh(fd,fh,0.05,[-1,-1;1,1],[-1,-1;-1,1;1,-1;1,1]);
%
%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

% parameters;
geps = 0.001*h0;

% 1. Create initial distribution in bounding box (equilateral triangles)
[x,y] = meshgrid(bbox(1,1):h0:bbox(2,1),bbox(1,2):h0*sqrt(3)/2:bbox(2,2));
x(2:2:end,:) = x(2:2:end,:)+h0/2;                % Shift even rows
p = [x(:),y(:)];                                 % List of node coordinates

% 2. Remove points outside the region, apply the rejection method
p = p(feval(fd,p,varargin{:}) < geps,:);        % Keep only d<0 points
r0 = 1./feval(fh,p,varargin{:}).^2;              % Probability to keep point
p = p(rand(size(p,1),1) < r0./max(r0),:);
p = [pfix; p + (rand(size(p,1),2)-2)*geps];  % Rejection method
% p = unique(p,'rows');
% p(1:size(pfix,1),:) = pfix;

% 3. Retriangulation by the Delaunay algorithm
t = delaunayn(p);                                % List of triangles
pmid = (p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;  % Compute centroids
t = t(feval(fd,pmid,varargin{:})<-geps,:);       % Keep interior triangles

repo: Aero-matlab


text:
function coords = makeSmoothCoords(c)

coords = [];
ii = 1;coordsInd = 1;
buff = 10;
while ii<size(c,2)
    n = c(2,ii);
    if n>100
        x = c(1,ii+1:ii+n);
        y = c(2,ii+1:ii+n);
        x = [x(end-buff:end) x x(1:buff)]; % buffer makes ends meet
        y = [y(end-buff:end) y y(1:buff)];
        x = smooth(x,25,'loess');
        y = smooth(y,25,'loess');
        x = [x(buff+1:end-buff-1);x(buff+1)];
        y = [y(buff+1:end-buff-1);y(buff+1)];
        coords(coordsInd).x = x;
        coords(coordsInd).y = y;
        coordsInd = coordsInd+1;

    end
    ii = ii+n+1;
end

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function metricSelection(keyword, varargin);
%GUI to build a planMetricS object, for passing to graphicalComparison.m
%Uses list of metrics stored in optS.planMetrics, which are the names of
%functions stored in \planMetrics, and designed for this purpose.
%JRA 06.16.03
%
%Bug Note:  If structure name selected for a metric is exactly the same as
%another structure name, the incorrect structure may be used for the
%metric.  [VHC 2/1/06]
%
%Modified:  VHC 2/1/06  Added save metric set feature.
%           VHC 2/6/06  Added criteria feature.
%                       Changed location of Output Range so it doesn't overlap.
%           VHC 2/13/06 Added report feature (not yet complete)
%
%The following function line is for compiler linking. Add functions referenced in callbacks to the list.
%#function meanDose maxDose minDose Vx Dx EUD
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global mSState; %This is for local state info that we dont want clogging up stateS.
global planC;
global stateS;
indexS = planC{end};

uicolor              = [.9 .9 .9];
units = 'normalized';

switch lower(keyword)
    case 'init'

        %Get figure, set its size based on screen size.
        screenSize = get(0,'ScreenSize');
        h = figure('Position', [screenSize(3)/10 screenSize(4)/10 screenSize(3)/10*8 screenSize(4)/10*6]);
        stateS.handle.metricSelectionFig = h;
        set(h, 'NumberTitle', 'off');
        set(h, 'Name', 'Metric Selection');
        set(h, 'MenuBar', 'none');

        mSState = [];
        mSState.handles.figure = h;

        %Init currentMetric vars
        mSState.currentMetric = [];
        mSState.currentMetricType = [];
        mSState.currentMetricIndex = 0;
        mSState.numPlanMetrics = 0;

        %Init saved metric set var
        if ~isfield(planC{end},'metrics')
            savedMetricSetsS = struct([]); %should be empty
            savedMetricSetsNames = {};
        else
            savedMetricSetsS = planC{planC{end}.metrics}.savedMetricSets;
            savedMetricSetsNames = {savedMetricSetsS.name};
        end

        %Init List Boxes
        mSState.handles.metricList = uicontrol('callback', 'metricSelection(''metricList_Callback'');', 'Min', length(stateS.optS.planMetrics), 'Max', length(stateS.optS.planMetrics), 'units',units,'BackgroundColor',uicolor, 'Position',[.05 .5+.45/2 .2 .45/2],'String', stateS.optS.planMetrics, 'Style','listbox','Tag','metricList');
        mSState.handles.myMetricList = uicontrol('callback', 'metricSelection(''myMetricList_Callback'');', 'Min', 1, 'Max', 100, 'units',units,'BackgroundColor',uicolor, 'Position',[.3 .5 .2 .45], 'String', [], 'Style','listbox','Tag','myMetricsList');
        mSState.handles.doseList = uicontrol('Min', 0, 'Max', 100, 'units',units,'BackgroundColor', uicolor, 'Position',[.05 .13 .45 .3],'String', {planC{planC{end}.dose}.fractionGroupID}, 'Style','listbox','Tag','doseList');
        mSState.handles.metricSetList = uicontrol('callback', 'metricSelection(''metricSetList_Callback'');', 'Min', 0, 'Max', 100, 'units',units,'BackgroundColor',uicolor, 'Position',[.05 .5 .2 .45/2-.03],'String', savedMetricSetsNames, 'Style','listbox','Tag','metricSetList');
        set(mSState.handles.myMetricList, 'Value', []);
        set(mSState.handles.metricList, 'Value', 1);
        set(mSState.handles.metricSetList, 'Value', []);

        %Make frame for metric info, and add labels and text inside of it
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.55 .05 .4 .85],'String','frame', 'Style','frame','Tag','metricInfoFrame');
        mSState.handles.metricNameTitle = uicontrol('callback', 'metricSelection(''metricnametitle_callback'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.6 .9 .3 .05],'String', '', 'Style','edit','Tag','metricNameTitle');
        mSState.handles.functionName = uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .83 .35 .04],'String', 'Function callback: ', 'Style','text','Tag','functionNameText', 'horizontalalignment', 'left');
        mSState.handles.functionDesc = uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .1 .35 .08],'String', 'Function description: ', 'Style','text','Tag','functionDescText', 'horizontalalignment', 'left');
        mSState.handles.parameterLabel = uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .77 .1 .04],'String', 'Parameters: ', 'Style','text','Tag','functionNameText', 'horizontalalignment', 'left');
        mSState.handles.criteriaLabel = uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .435 .1 .04],'String', 'Criteria: ', 'Style','text','Tag','functionNameText', 'horizontalalignment', 'left');

        %Make labels for listboxes
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.05 .95 .1 .03],'String', 'Available Metrics: ', 'Style','text','Tag','functionNameText', 'horizontalalignment', 'left');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.3 .95 .1 .03],'String', 'Selected Metrics: ', 'Style','text','Tag','functionNameText', 'horizontalalignment', 'left');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.05 .43 .11 .03],'String', 'Plans to Compare:', 'Style','text','Tag','planSelectionText', 'horizontalalignment', 'left');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.05 .95-.45/2-.03 .13 .03],'String', 'Available Metric Sets:', 'Style','text','Tag','functionNameText', 'horizontalalignment', 'left');
        %uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.18 .95-.45/2-.03 .03 .03],'String', 'Load', 'Style','push','Tag','loadPush', 'horizontalalignment', 'left','callback', 'metricSelection(''metricfileselected'');');

        %make buttons for list boxes, evaluation, canceling, and saving.
        uicontrol('callback', 'metricSelection(''addMetric'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.25 .75 .05 .05],'String','->', 'Style','pushbutton','Tag','metricAdd');
        uicontrol('callback', 'metricSelection(''removeMetric'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.25 .70 .05 .05],'String','<-', 'Style','pushbutton','Tag','metricRemove');
        uicontrol('callback', 'metricSelection(''evaluateMetrics'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.05 .04 .1 .05],'String','Evaluate', 'Style','pushbutton','Tag','evaluateMetrics');
        uicontrol('callback', 'metricSelection(''generateReport'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.17 .04 .1 .05],'String','Report', 'Style','pushbutton','Tag','cancel');
        uicontrol('callback', 'metricSelection(''cancel'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.29 .04 .1 .05],'String','Cancel', 'Style','pushbutton','Tag','cancel');
        uicontrol('callback', 'metricSelection(''saveMetricSet'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.44 .95 .05 .05],'String','Save', 'Style','pushbutton','Tag','saveMetricSet');
        uicontrol('callback', 'metricSelection(''metricfileselected'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.49 .95 .05 .05],'String','Load', 'Style','pushbutton','Tag','loadPush');
        mSState.expandesView = 0;
        return;

    case 'metriclist_callback'
        set(mSState.handles.myMetricList, 'Value', []); %unhighlight the user listbox
        set(mSState.handles.metricSetList, 'Value', []); %unhighlight the builtin set listbox
        mSState.currentMetricType = 'BuiltIn';
        mSState.currentMetricIndex = 0; %built in selected, set index to 0.
        metricSelection('builtInMetricSelected');
        if ~mSState.expandesView && strcmpi(mSState.currentMetric.name,'lkb')
            lkbExpandView('EXPANDEDVIEW')
        elseif mSState.expandesView && ~strcmpi(mSState.currentMetric.name,'lkb')
            lkbExpandView('SIMPLEVIEW')
        end
        return

    case 'mymetriclist_callback'
        set(mSState.handles.metricList, 'Value', 1); %unhighlight the builtin listbox
        set(mSState.handles.metricSetList, 'Value', []); %unhighlight the builtin set listbox
        mSState.currentMetricType = 'User';
        mSState.currentMetricIndex = get(mSState.handles.myMetricList, 'Value');
        metricSelection('usermetricselected');
        if ~mSState.expandesView && strcmpi(mSState.currentMetric.name,'lkb')
            lkbExpandView('EXPANDEDVIEW')
        elseif mSState.expandesView && ~strcmpi(mSState.currentMetric.name,'lkb')
            lkbExpandView('SIMPLEVIEW')
        end
        return

        return

    case 'metricsetlist_callback'
        set(mSState.handles.myMetricList, 'Value', []); %unhighlight the user listbox
        set(mSState.handles.metricList, 'Value', 1); %unhighlight the builtin listbox
        mSState.currentMetricType = 'BuiltInSet';
        mSState.currentMetricIndex = 0; %built in selected, set index to 0.
        metricSelection('builtInMetricSetSelected');
        return

    case 'builtinmetricselected'
        delete(findobj('tag', 'parameter')); %remove all old parameter GUI objects.
        mSState.handles.params = [];
        [planC, metric] = feval(stateS.optS.planMetrics{get(mSState.handles.metricList, 'Value')}, planC, 'getNewMetric'); %special call to get a default metric
        mSState.currentMetric = metric;
        set(mSState.handles.metricNameTitle, 'String', metric.name);
        set(mSState.handles.functionName, 'String', ['Function callback: ' func2str(metric.functionName)]);
        set(mSState.handles.functionDesc, 'String', ['Function description: ' metric.description]);
        set(mSState.handles.parameterLabel, 'String', 'Parameters: ');
        pos = get(mSState.handles.functionName,'Position');
        for i=1:length(metric.params)  %populate parameter window with all existing parameters, using passed default values
            switch metric.params(i).type
                case 'Edit'
                    mSState.handles.params = [mSState.handles.params uicontrol('callback', ['metricSelection(''parameter_callback'',' num2str(i) ');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .70 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).value, 'Style','edit','Tag','parameter')];
                    uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .74 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).name, 'Style','text','Tag','parameter','horizontalalignment', 'left');
                case 'DropDown'
                    mSState.handles.params = [mSState.handles.params uicontrol('callback', ['metricSelection(''parameter_callback'',' num2str(i) ');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .70 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).list, 'Style','popupmenu','Tag','parameter')];
                    uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .74 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).name, 'Style','text','Tag','parameter','horizontalalignment', 'left');
            end
        end

        %criteria
        delete(findobj('tag', 'criteria')); %remove all old criteria GUI objects.
        if ~isfield(mSState.currentMetric,'criteria') %there are no defaults for the criteria, so we'll put them in
            cS = struct();
            cS.passDirectionIndex = 1; cS.passDirectionList = {'above','below'}; cS.passValue = 0;
            cS.marginalDirectionIndex = 1; cS.marginalDirectionList = {'above','below'}; cS.marginalValue = 0;
            cS.priority = 1; cS.priorityTypeIndex = 1; cS.priorityTypeList = {'strict','fuzzy'};
            cS.passStatus = {}; %will be 'passed', 'marginal', or 'failed' for ea dose dist
            mSState.currentMetric.criteria = cS;
        else
            cS = mSState.currentMetric.criteria;
        end
        %if you change something here, then copy it to case 'usermetricselected as well.
        set(mSState.handles.criteriaLabel, 'String', 'Criteria: ');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.04 .12 .04],'String', 'Pass:       Value must be ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.passDirectionIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''passDirectionIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12 .435-.031 .07 .04],'String', cS.passDirectionList, 'Value', cS.passDirectionIndex, 'Style','popupmenu','Tag','criteria');
        mSState.handles.criteria.passValue = uicontrol('callback', ['metricSelection(''criteria_callback'', ''passValue'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12+.07+.005 .435-.031 .1 .04],'String', cS.passValue, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.095 .12 .04],'String', 'Marginal:  Value must be ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.marginalDirectionIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''marginalDirectionIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12 .435-.086 .07 .04],'String', cS.marginalDirectionList, 'Value', cS.marginalDirectionIndex, 'Style','popupmenu','Tag','criteria');
        mSState.handles.criteria.marginalValue = uicontrol('callback', ['metricSelection(''criteria_callback'', ''marginalValue'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12+.07+.005 .435-.086 .1 .04],'String', cS.marginalValue, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.15 .055 .04],'String', 'Priority: ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.priority = uicontrol('callback', ['metricSelection(''criteria_callback'', ''priority'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.060 .435-.141 .05 .04],'String', cS.priority, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.145 .435-.15 .15 .04],'String', '1 = highest, 2 = second-highest, etc.', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.205 .1 .04],'String', 'Priority Type: ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.priorityTypeIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''priorityTypeIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.07 .435-.196 .07 .04],'String', cS.priorityTypeList, 'Value', cS.priorityTypeIndex, 'Style','popupmenu','Tag','criteria');

        %dont forget range dialogue, defaults again
        %         mSState.handles.range = [uicontrol('callback', ['metricSelection(''range_callback'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[.57 .1 .1 .04],'String', num2str(metric.range), 'Style','edit','Tag','parameter')];
        %         uicontrol('units',units,'BackgroundColor',uicolor, 'Position', [.57 .06 .1 .04],'String', 'Output Range', 'Style','text','Tag','parameter','horizontalalignment', 'center');
        %Double clicks automatically add metric to userlist
        if strcmpi(get(mSState.handles.figure, 'SelectionType'), 'open')
            set(mSState.handles.figure, 'SelectionType', 'normal');
            metricSelection('addmetric');
        end

    case 'metricfileselected'
        delete(findobj('tag', 'parameter')); %remove all old parameter GUI objects.
        mSState.handles.params = [];
        [fname, pathname] = uigetfile({'*.xls'}, 'Select Excel file to load stored Metric');
        if isnumeric(fname) && fname == 0
            return;
        end
        file = fullfile(pathname, fname);
        [numeric,txt,rawData] = xlsread(file);
        [planC, metric] = feval(rawData{2,2}, planC, 'getNewMetric'); %special call to get a default metric
        delete(findobj('tag', 'criteria')); %remove all old criteria GUI objects.
        cS = struct();
        cS.passDirectionIndex = 1; cS.passDirectionList = {'above','below'}; cS.passValue = 0;
        cS.marginalDirectionIndex = 1; cS.marginalDirectionList = {'above','below'}; cS.marginalValue = 0;
        cS.priority = 1; cS.priorityTypeIndex = 1; cS.priorityTypeList = {'strict','fuzzy'};
        cS.passStatus = {}; %will be 'passed', 'marginal', or 'failed' for ea dose dist

        %Fill-in user defined params and criteria
        structIndex = strmatch(lower(rawData{7,2}),lower({planC{indexS.structures}.structureName}),'exact');
        switch lower(metric.name)
            case {'meandose','mindose','maxdose','stdDevDose'}
                %Set params
                metric.params(1).value = structIndex;
                if strcmpi(lower(rawData{8,2}),'no')
                    metric.params(2).value = 2;
                else
                    metric.params(2).value = 1;
                end
                critInd = 11;

            case 'vx'
                %Set params
                metric.params(1).value = structIndex;
                metric.params(2).value = rawData{8,2};
                if strcmpi(lower(rawData{9,2}),'absolute')
                    metric.params(3).value = 2;
                else
                    metric.params(3).value = 1;
                end
                if strcmpi(lower(rawData{10,2}),'no')
                    metric.params(4).value = 2;
                else
                    metric.params(4).value = 1;
                end
                critInd = 13;

            case 'dx'
                %Set params
                metric.params(1).value = structIndex;
                metric.params(2).value = rawData{8,2};
                if strcmpi(lower(rawData{9,2}),'no')
                    metric.params(3).value = 2;
                else
                    metric.params(3).value = 1;
                end
                critInd = 12;

            case 'eud'
                %Set params
                metric.params(1).value = structIndex;
                if strcmpi(lower(rawData{8,2}),'no')
                    metric.params(2).value = 2;
                else
                    metric.params(2).value = 1;
                end
                metric.params(3).value = rawData{9,2};
                critInd = 12;

            case 'erp'
                %Set params
                metric.params(1).value = structIndex;
                if strcmpi(lower(rawData{8,2}),'no')
                    metric.params(2).value = 2;
                else
                    metric.params(2).value = 1;
                end
                metric.params(3).value = rawData{9,2};
                metric.params(4).value = rawData{10,2};
                metric.params(5).value = rawData{11,2};
                critInd = 13;

            case 'lkb'
                %Set params
                metric.params(1).value = structIndex;
                if strcmpi(lower(rawData{8,2}),'no')
                    metric.params(2).value = 2;
                else
                    metric.params(2).value = 1;
                end
                metric.params(3).value = rawData{9,2};
                metric.params(4).value = rawData{10,2};
                metric.params(5).value = rawData{11,2};
                critInd = 14;

        end

        %Set criteria
        if strcmpi(rawData{critInd,2},'above')
            cS.passDirectionIndex = 1;
        else
            cS.passDirectionIndex = 2;
        end
        cS.passValue = rawData{critInd+1,2};
        if strcmpi(rawData{critInd+2,2},'above')
            cS.marginalDirectionIndex = 1;
        else
            cS.marginalDirectionIndex = 2;
        end
        if strcmpi(rawData{critInd+4,2},'high')
            cS.priority = 1;
        else
            cS.priority = 2;
        end
        if strcmpi(rawData{critInd+5,2},'strict')
            cS.priorityTypeIndex = 1;
        else
            cS.priorityTypeIndex = 2;
        end
        metric.criteria = cS;

        mSState.currentMetric = metric;
        set(mSState.handles.metricNameTitle, 'String', metric.name);
        set(mSState.handles.functionName, 'String', ['Function callback: ' func2str(metric.functionName)]);
        set(mSState.handles.functionDesc, 'String', ['Function description: ' metric.description]);
        set(mSState.handles.parameterLabel, 'String', 'Parameters: ');
        for i=1:length(metric.params)  %populate parameter window with all existing parameters, using passed default values
            switch metric.params(i).type
                case 'Edit'
                    mSState.handles.params = [mSState.handles.params uicontrol('callback', ['metricSelection(''parameter_callback'',' num2str(i) ');'], 'units',units,'BackgroundColor',uicolor, 'Position',[.57 + mod(i-1,3)*.11 .70 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).value, 'Style','edit','Tag','parameter')];
                    uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 + mod(i-1,3)*.11 .74 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).name, 'Style','text','Tag','parameter','horizontalalignment', 'left');
                case 'DropDown'
                    mSState.handles.params = [mSState.handles.params uicontrol('callback', ['metricSelection(''parameter_callback'',' num2str(i) ');'], 'units',units,'BackgroundColor',uicolor, 'Position',[.57 + mod(i-1,3)*.11 .70 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).list, 'Style','popupmenu','Tag','parameter')];
                    uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 + mod(i-1,3)*.11 .74 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).name, 'Style','text','Tag','parameter','horizontalalignment', 'left');
            end
        end

        %if you change something here, then copy it to case 'usermetricselected as well.
        set(mSState.handles.criteriaLabel, 'String', 'Criteria: ');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .435-.04 .15 .04],'String', 'Pass:       Value must be ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.passDirectionIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''passDirectionIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.57+.15 .435-.031 .07 .04],'String', cS.passDirectionList, 'Value', cS.passDirectionIndex, 'Style','popupmenu','Tag','criteria');
        mSState.handles.criteria.passValue = uicontrol('callback', ['metricSelection(''criteria_callback'', ''passValue'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[.57+.15+.07+.005 .435-.031 .1 .04],'String', cS.passValue, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .435-.095 .15 .04],'String', 'Marginal:  Value must be ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.marginalDirectionIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''marginalDirectionIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.57+.15 .435-.086 .07 .04],'String', cS.marginalDirectionList, 'Value', cS.marginalDirectionIndex, 'Style','popupmenu','Tag','criteria');
        mSState.handles.criteria.marginalValue = uicontrol('callback', ['metricSelection(''criteria_callback'', ''marginalValue'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[.57+.15+.07+.005 .435-.086 .1 .04],'String', cS.marginalValue, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .435-.15 .055 .04],'String', 'Priority: ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.priority = uicontrol('callback', ['metricSelection(''criteria_callback'', ''priority'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[.57+.060 .435-.141 .1 .04],'String', cS.priority, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57+.165 .435-.15 .21 .04],'String', '1 = highest, 2 = second-highest, etc.', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[.57 .435-.205 .1 .04],'String', 'Priority Type: ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.priorityTypeIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''priorityTypeIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.57+.10 .435-.196 .07 .04],'String', cS.priorityTypeList, 'Value', cS.priorityTypeIndex, 'Style','popupmenu','Tag','criteria');

        mSState.currentMetricType = 'BuiltIn';
        metricSelection('addmetric');

    case 'builtinmetricsetselected'  %need to remove more things than just this
        delete(findobj('tag', 'parameter')); %remove all old parameter GUI objects.
        delete(findobj('tag', 'criteria')); %remove all old criteria GUI objects.
        mSState.handles.params = [];
        mSState.currentMetric = [];
        set(mSState.handles.metricNameTitle, 'String', planC{planC{end}.metrics}.savedMetricSets(get(mSState.handles.metricSetList, 'Value')).name);
        set(mSState.handles.functionName, 'String', '');
        set(mSState.handles.functionDesc, 'String', '');
        set(mSState.handles.parameterLabel, 'String', '');
        set(mSState.handles.criteriaLabel, 'String', '');
        %Double clicks automatically add metric set to userlist
        if strcmpi(get(mSState.handles.figure, 'SelectionType'), 'open')
            set(mSState.handles.figure, 'SelectionType', 'normal');
            metricSelection('addmetric');
        end

    case 'usermetricselected'
        delete(findobj('tag', 'parameter')); %remove all old parameter GUI objects.
        mSState.handles.params = [];
        metric = mSState.planMetricsS(get(mSState.handles.myMetricList, 'Value'));
        if length(metric)>=1
            metric = metric(1);
        else
            return;
        end
        set(mSState.handles.metricNameTitle, 'String', metric.name);
        set(mSState.handles.functionName, 'String', ['Function callback: ' func2str(metric.functionName)]);
        set(mSState.handles.functionDesc, 'String', ['Function description: ' metric.description]);
        set(mSState.handles.parameterLabel, 'String', 'Parameters: ');
        mSState.currentMetric = metric;
        pos = get(mSState.handles.functionName,'Position');
        for i=1:length(metric.params)  %populate parameter window with all existing parameters, using passed default values
            switch metric.params(i).type
                case 'Edit'
                    mSState.handles.params = [mSState.handles.params uicontrol('callback', ['metricSelection(''parameter_callback'',' num2str(i) ');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .70 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).value, 'Style','edit','Tag','parameter')];
                    uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .74 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).name, 'Style','text','Tag','parameter','horizontalalignment', 'left');
                case 'DropDown'
                    mSState.handles.params = [mSState.handles.params uicontrol('callback', ['metricSelection(''parameter_callback'',' num2str(i) ');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .70 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).list, 'value', metric.params(i).value, 'Style','popupmenu','Tag','parameter')];
                    uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) + mod(i-1,3)*.11 .74 - mod(floor((i-1)/3), 6)*.1 .1 .04],'String', metric.params(i).name, 'Style','text','Tag','parameter','horizontalalignment', 'left');
            end
        end

        %criteria
        delete(findobj('tag', 'criteria')); %remove all old criteria GUI objects.
        if ~isfield(mSState.currentMetric,'criteria') %there are no defaults for the criteria, so we'll put them in
            cS = struct();
            cS.passDirectionIndex = 1; cS.passDirectionList = {'above','below'}; cS.passValue = 0;
            cS.marginalDirectionIndex = 1; cS.marginalDirectionList = {'above','below'}; cS.marginalValue = 0;
            cS.priority = 1; cS.priorityTypeIndex = 1; cS.priorityTypeList = {'strict','fuzzy'};
            cS.passStatus = {}; %will be 'passed', 'marginal', or 'failed' for ea dose dist
            mSState.currentMetric.criteria = cS;
        else
            cS = mSState.currentMetric.criteria;
        end
        %if you change something here, then copy it to case 'usermetricselected as well.
        set(mSState.handles.criteriaLabel, 'String', 'Criteria: ');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.04 .12 .04],'String', 'Pass:       Value must be ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.passDirectionIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''passDirectionIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12 .435-.031 .07 .04],'String', cS.passDirectionList, 'Value', cS.passDirectionIndex, 'Style','popupmenu','Tag','criteria');
        mSState.handles.criteria.passValue = uicontrol('callback', ['metricSelection(''criteria_callback'', ''passValue'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12+.07+.005 .435-.031 .1 .04],'String', cS.passValue, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.095 .12 .04],'String', 'Marginal:  Value must be ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.marginalDirectionIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''marginalDirectionIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12 .435-.086 .07 .04],'String', cS.marginalDirectionList, 'Value', cS.marginalDirectionIndex, 'Style','popupmenu','Tag','criteria');
        mSState.handles.criteria.marginalValue = uicontrol('callback', ['metricSelection(''criteria_callback'', ''marginalValue'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.12+.07+.005 .435-.086 .1 .04],'String', cS.marginalValue, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.15 .055 .04],'String', 'Priority: ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.priority = uicontrol('callback', ['metricSelection(''criteria_callback'', ''priority'');'], 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.060 .435-.141 .05 .04],'String', cS.priority, 'Style','edit','Tag','criteria');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.145 .435-.15 .15 .04],'String', '1 = highest, 2 = second-highest, etc.', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position',[pos(1) .435-.205 .1 .04],'String', 'Priority Type: ', 'Style','text','Tag','criteria', 'horizontalalignment', 'left');
        mSState.handles.criteria.priorityTypeIndex = uicontrol('callback', 'metricSelection(''criteria_callback'', ''priorityTypeIndex'');', 'units',units,'BackgroundColor',uicolor, 'Position',[pos(1)+.07 .435-.196 .07 .04],'String', cS.priorityTypeList, 'Value', cS.priorityTypeIndex, 'Style','popupmenu','Tag','criteria');

        %dont forget range dialogue, stored values
        mSState.handles.range = uicontrol('callback', 'metricSelection(''range_callback'');', 'units',units,'BackgroundColor',uicolor, 'Position',[.57 .1-.015 .1 .04],'String', num2str(metric.range), 'Style','edit','Tag','parameter');
        uicontrol('units',units,'BackgroundColor',uicolor, 'Position', [pos(1) .06-.015+.01 .1 .04-.01],'String', 'Output Range', 'Style','text','Tag','parameter','horizontalalignment', 'center');
        %Double clicks automatically remove metric from userlist
        if strcmpi(get(mSState.handles.figure, 'SelectionType'), 'open')
            set(mSState.handles.figure, 'SelectionType', 'normal');
            metricSelection('removemetric');
        end

    case 'addmetric' %add current metric to the user metrics list, and store--select the stored user metric as the current metric.
        if(strcmp(mSState.currentMetricType, 'BuiltIn'))
            mSState.planMetricsS(mSState.numPlanMetrics + 1) = mSState.currentMetric;
            mSState.numPlanMetrics = mSState.numPlanMetrics + 1;
            set(findobj('tag', 'myMetricsList'), 'String', {mSState.planMetricsS.name});
            set(findobj('tag', 'myMetricsList'), 'Value', mSState.numPlanMetrics);
            metricSelection('mymetriclist_callback');
        elseif (strcmp(mSState.currentMetricType, 'BuiltInSet'))
            numMetricSet = get(mSState.handles.metricSetList, 'Value');
            for i=1:size(planC{planC{end}.metrics}.savedMetricSets(numMetricSet).metricSet,2)
                mSState.currentMetric = planC{planC{end}.metrics}.savedMetricSets(numMetricSet).metricSet(i);
                mSState.currentMetricType = 'BuiltIn';
                metricSelection('addmetric');
            end
            mSState.currentMetricType = 'BuiltInSet';
        end

        return;

    case 'parameter_callback' %parameter has changed, update the current and stored version if it exists
        switch get(mSState.handles.params(varargin{1}), 'Style');
            case 'edit'
                mSState.currentMetric.params(varargin{1}).value = get(mSState.handles.params(varargin{1}), 'String');
            case 'popupmenu'
                mSState.currentMetric.params(varargin{1}).value = get(mSState.handles.params(varargin{1}), 'Value');
        end
        if(strcmp(mSState.currentMetricType, 'User'))
            mSState.planMetricsS(mSState.currentMetricIndex) = mSState.currentMetric;
        end
        if strcmpi(mSState.currentMetric.name,'lkb')
            lkbExpandView('UPDATECURVE')
        end

    case 'criteria_callback' %criteria input has changed, update the current and stored version if it exists
        tempHandle = eval(['mSState.handles.criteria.' varargin{1}]);
        switch get(tempHandle, 'Style');
            case 'edit'
                eval(['mSState.currentMetric.criteria.' varargin{1} ' = str2num(get(tempHandle, ''String''));']);
            case 'popupmenu'
                eval(['mSState.currentMetric.criteria.' varargin{1} ' = get(tempHandle, ''Value'');']);
        end
        if(strcmp(mSState.currentMetricType, 'User'))
            mSState.planMetricsS(mSState.currentMetricIndex) = mSState.currentMetric;
        end

    case 'range_callback' %range has changed, update current and stored version if it exists
        mSState.currentMetric.range = str2num(get(mSState.handles.range, 'String'));
        if(strcmp(mSState.currentMetricType, 'User'))
            mSState.planMetricsS(mSState.currentMetricIndex) = mSState.currentMetric;
        end

    case 'metricnametitle_callback' %metric name has changed, update current and stored version if it exists
        mSState.currentMetric.name = get(findobj('tag', 'metricNameTitle'), 'String');
        if(strcmp(mSState.currentMetricType, 'User'))
            mSState.planMetricsS(mSState.currentMetricIndex) = mSState.currentMetric;
            set(findobj('tag', 'myMetricsList'), 'String', {mSState.planMetricsS.name});
        end

    case 'removemetric' %current metric, if it is in the user list, is removed from list & planMetricsS
        if(strcmp(mSState.currentMetricType, 'User'))
            mSState.planMetricsS(mSState.currentMetricIndex) = [];
            if(mSState.currentMetricIndex == mSState.numPlanMetrics)
                set(mSState.handles.myMetricList, 'Value', [mSState.currentMetricIndex - 1]);
                mSState.currentMetricIndex = mSState.currentMetricIndex - 1;
            end
            mSState.numPlanMetrics = mSState.numPlanMetrics - length(mSState.currentMetricIndex);
            set(findobj('tag', 'myMetricsList'), 'String', {mSState.planMetricsS.name});
        end

    case 'evaluatemetrics' %init graphical comparison, and give it planMetricsS and a list of doses/names to perform calcs on
        try
            graphicalComparison('init', planC);
            dosesToCompare = get(mSState.handles.doseList, 'Value');

            badDosesV = [];
            for i=1:length(dosesToCompare)
                if isCompressed(planC{indexS.dose}(dosesToCompare(i)).doseArray)
                    badDosesV = [badDosesV i];
                end
            end
            if ~isempty(badDosesV)
                warning(['Cannot calculate metrics on compressed dosearrays.  Skipping doses [' num2str(badDosesV) '].']);
            end
            dosesToCompare(badDosesV) = [];

            for i=1:length(mSState.planMetricsS) %set doses to operate on in metrics.
                mSState.planMetricsS(i).doseSets = dosesToCompare;
            end
            planC = graphicalComparison('add', planC, stateS.optS, mSState.planMetricsS, {planC{indexS.dose}(mSState.planMetricsS(1).doseSets).fractionGroupID});
        catch
            warning(['Error in graphicalComparison: ' lasterr]);
        end

    case 'generatereport'
        try

            %get name of html report file to generate
            [fname, pathname] = uiputfile({'*.html', 'HTML files (*.html)';'*.*', 'All Files (*.*)'}, 'Select a location and filename to save the generated html report file'); %,'Location',[100,100]);

            if (fname == 0) %then user canceled
                return;
            end

            %get doses to compare, and put them into appropriate place in planMetricsS
            dosesToCompare = get(mSState.handles.doseList, 'Value');

            badDosesV = [];
            for i=1:length(dosesToCompare)
                if isCompressed(planC{indexS.dose}(dosesToCompare(i)).doseArray)
                    badDosesV = [badDosesV i];
                end
            end
            if ~isempty(badDosesV)
                warning(['Cannot calculate metrics on compressed dosearrays.  Skipping doses [' num2str(badDosesV) '].']);
            end
            dosesToCompare(badDosesV) = [];

            for i=1:length(mSState.planMetricsS) %set doses to operate on in metrics.
                mSState.planMetricsS(i).doseSets = dosesToCompare;
            end

            %call function to generate report
            [planC, mSState.planMetricsS] = generateReport(planC, mSState.planMetricsS, fname, pathname, {planC{indexS.dose}(mSState.planMetricsS(1).doseSets).fractionGroupID}); % stateS.optS

            %display html file in help browser
            open([pathname fname]);

        catch
            warning(['Error in generateReport: ' lasterr]);
        end

    case 'savemetricset' %save a set of metrics (in planC) for later reference.
        if ~isfield(mSState, 'planMetricsS') | isempty(mSState.planMetricsS)
            errordlg('Cannot save a metric set with no metrics in it.');
            return;
        end
        if ~isfield(planC{end},'metrics')
            metricSelection('createMetricsStructInPlanC');
        end
        metricSetNameC = inputdlg('Name to save metric set as:');
        if isempty(metricSetNameC) %they cancelled
            return;
        end
        planC{planC{end}.metrics}.savedMetricSets(end+1).name = metricSetNameC{1};
        planC{planC{end}.metrics}.savedMetricSets(end).metricSet = mSState.planMetricsS;
        set(mSState.handles.metricSetList, 'String', {planC{planC{end}.metrics}.savedMetricSets.name});

    case 'createmetricsstructinplanc'
        planC{end}.metrics = size(planC,2);
        planC{end}.indexS = size(planC,2)+1;
        planC{planC{end}.indexS} = planC{end};
        planC{planC{end}.metrics} = struct('savedMetricSets',struct([]));

    case 'cancel'
        close;

    otherwise
        metricSelectionError = 1
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function q3 = quat_mul ( q1, q2 )

%*****************************************************************************80
%
%% QUAT_MUL multiplies two quaternions.
%
%  Discussion:
%
%    A quaternion is a quadruplet (A,B,C,D) of real numbers, which
%    may be written as
%
%      Q = A + Bi + Cj + Dk.
%
%    To multiply two quaternions, use the relationships:
%
%      i * j = -j * i = k
%      j * k = -k * j = i
%      k * i = -i * k = j
%      i * i =  j * j = k * k = -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real Q1(4), Q2(4), the two quaternions to be multiplied.
%
%    Output, real Q3(4), the product of the two quaternions.
%
  q3(1) = q1(1) * q2(1) - q1(2) * q2(2) - q1(3) * q2(3) - q1(4) * q2(4);
  q3(2) = q1(1) * q2(2) + q1(2) * q2(1) + q1(3) * q2(4) - q1(4) * q2(3);
  q3(3) = q1(1) * q2(3) - q1(2) * q2(4) + q1(3) * q2(1) + q1(4) * q2(2);
  q3(4) = q1(1) * q2(4) + q1(2) * q2(3) - q1(3) * q2(2) + q1(4) * q2(1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function region = p12_region ( )

%*****************************************************************************80
%
%% P12_REGION returns the name of the integration region for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end


text:
function figure_num = elman_contour ( figure_num )

%*****************************************************************************80
%
%% ELMAN_CONTOUR displays a contour plot of a 2D stochastic diffusivity function.
%
%  Discussion:
%
%    The diffusivity function is compute by DIFFUSIVITY_2D_ELMAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Howard Elman, Darran Furnaval,
%    Solving the stochastic steady-state diffusion problem using multigrid,
%    IMA Journal on Numerical Analysis,
%    Volume 27, Number 4, 2007, pages 675-688.
%
%  Parameters:
%
%    Input/output, integer FIGURE_NUM, the current number of figures.
%
  if ( nargin < 1 )
    figure_num = 0;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELMAN_CONTOUR\n' );
  fprintf ( 1, '  Display contour or surface plots of the stochastic\n' );
  fprintf ( 1, '  diffusivity function defined by DIFFUSIVITY_2D_ELMAN.\n' );
%
%  Initialize the normal random number generator.
%
  seed = -1;
  [ ~, ~ ] = r8vec_normal_01 ( -1, seed );
%
%  Specify the X and Y evaluation points.
%
  nx = 51;
  a = 1.0;
  xvec = linspace ( -a, +a, nx );
  yvec = linspace ( -a, +a, nx );

  [ xmat, ymat ] = meshgrid ( xvec, yvec );
%
%  Sample OMEGA.
%
  m_1d = 5;
  msq = m_1d * m_1d;
  seed = 123456789;
  if ( seed == 0 )
    omega = randn ( msq, 1 );
  else
    [ omega, seed ] = r8vec_normal_01 ( msq, seed );
  end
%
%  Compute the diffusivity field.
%
  cl = 0.1;
  dc0 = 10.0;
  dc = diffusivity_2d_elman ( a, cl, dc0, m_1d, omega, nx, nx, xmat, ymat );
%
%  Make a surface plot.
%
  figure_num = figure_num + 1;
  figure ( figure_num );
  surf ( ymat, xmat, dc, 'EdgeColor', 'interp' )
  xlabel ( 'Y' )
  ylabel ( 'X' )
  zlabel ( 'DC(X,Y)' )
  title ( 'ELMAN Stochastic diffusivity function' )

  filename = 'elman_contour.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file stored as "%s".\n', filename );

  return
end


text:
function varargout=rdmimic2wave(varargin)
%
% [tm,signal,Fs,recList,sigInfo]=rdmimic2wave(subjectID,clinicalTimeStamp,dataType,beginMinute,endMinute)
%
%
%
% Output Parameters:
%
% tm
%       Nx1 vector of doubles representing time in minutes if a match is
%       found, otherwise tm is empty. The time of the first index is an offset
%       from beginMinute. So that the total duration in minutes is:
%
%       tm(end) - tm(1) = endMinute - beginMinute
%
%       And the time of the clinical event (clinicalTimeStamp) in the waveform signal is given by
%       (to the closest minute):
%
%       timeStamp = tm(1) + beginMinute
%
%
% signal
%       NxM vector of doubles representing the signals from dataType that
%       match the timeStamp. If no match is found an empty matrix is
%       returned.
%
% Optional Ouput:
%
% Fs
%       A 1x1 double representing the sampling frequency (in Hz).
%
%recList
%       A Lx1 double specifying a list of valid subjectIDs  or matched record. You can use this
%       field to find which IDs are in the matched waveform list (see below). If the
%       subjectID is defined in the function signature, recList will be
%       equal to the first found matched record if it exists.
%
%
%sigInfo
%       A Mx1 structure containing meta information about the waverforms in the 'signal' output.
%
%
% Input Parameters:
%
% subjectID
%       A 1x1 Double specifying a valid MIMIC II subject ID. For a list
%       of valid subjectID with matched waveform use this to query:
%
%       [~,~,~,recList]=rdmimic2wave([],[],dataType);
%
%       Once you have a valid subjectID and pass it to RDMIMI2WAVE, recList
%       will return the string name of the first matched record if any (empty otherwise).
%
%
% clinicalTimeStamp
%      String specifying the clinical  time of the event. This string
%      should have the following format (as described in
%       http://www.physionet.org/physiobank/database/mimic2wdb/matched/) :
%
%      'YYYY-MM-DD-hh-mm'
%
%      Where:
%           YYYY = surrogate year
%           MM = month (01-12)
%           DD = day (01-31)
%           hh = real hour (00-23)
%           mm = minute (00-59)
%
% dataType (Optional)
%     String specifying what time of high resolution waveform to fetch.
%     Options are: 'numerics' and 'waveform'. Default is 'numerics'.
%
%
% beginMinute (Optional)
%     1x1 Double specifying time in minutes of
%     how much of the signal to get before clinicalTimeStamp occured.
%     Default is 60 minutes.
%
% endMinute (Optional)
%     1x1 Double specifying time in minutes of
%     how much of the signal to get after clinicalTimeStamp occured.
%     Default is 60 minutes.
%
%
%
% % Example:
%[tm,signal,Fs,recList,sigInfo]=rdmimic2wave(32805,'2986-12-15-10-00',[],0,2);
%plot(tm,signal(:,2))
%title(['Found data in record: ' recList])
%legend(sigInfo(2).Description)
%
% Written by Ikaro Silva, 2013
% Last Modified: December 1, 2014
%
% Version 1.1
%
% Since 0.9.0
%
%
%
% See also rdsamp, wfdbdesc

%endOfHelp

inputs={'subjectID','clinicalTimeStamp','dataType','beginWindow','endWindow'};
outputs={'tm','signal','Fs','recList','sigInfo'};
subjectID=[];
beginWindow=60; %beginWindow in minutes!
endWindow=60; %endWindow in minutes!
dataType='numerics';
dBName='mimic2wdb/matched/';
tm=[];
signal=[];
Fs=[];
sigInfo=[];
recList=[];
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

persistent cachedDataType matched_id matched

if(~strcmp(cachedDataType,dataType))
    switch (dataType)
        case 'numerics'
            matched=urlread('http://www.physionet.org/physiobank/database/mimic2wdb/matched/RECORDS-numerics');
            %case 'all'
            % Dont want to mess around with mixed sample rates yet...
            %%matched=urlread('http://www.physionet.org/physiobank/database/mimic2wdb/matched/RECORDS');
        case 'waveform'
            matched=urlread('http://www.physionet.org/physiobank/database/mimic2wdb/matched/RECORDS-waveforms');
        otherwise
            error(['Unknow dataType. Options are: ''numerics'', ''all'', or ''waveform''']);
    end
    cachedDataType=dataType;

    %Get patient ID list
    matched_id=regexprep(regexp(matched,'(\<s\d\d\d\d\d)','match'),'s','');
    for n=1:length(matched_id)
        matched_id{n}=str2num(matched_id{n});
    end
    matched_id=unique(cell2mat(matched_id));

   recList=matched_id;

end

if(isempty(subjectID))
    %In this case the user is just querying for a list of matched
    %records
    recList=matched_id;
    for n=1:nargout
        eval(['varargout{n}=' outputs{n} '']);
    end
    return

end

%Convert timestam to serial data (in days)
dateFormat='yyyy-mm-dd-HH-MM';
clinicalDateNum=datenum(clinicalTimeStamp,dateFormat);

%Window to include before and after the measurement (in days)
beginTime=clinicalDateNum - (beginWindow/(60*24));
endTime=clinicalDateNum + (endWindow/(60*24));

%If id exists loop through the files to see if any file is within the
%specific time range of the clinial event
matched_pid=find(matched_id==subjectID,1);
if(~isempty(matched_pid))

    %Get all records that match the patient and search for a record
    %that includes the clinical time
    strSubjectID=sprintf('%05d',subjectID);
    eval(['recs=regexp(regexp(matched,''(\<s' strSubjectID '.*)'',''match''),''\n'',''split'');']);
    recs=recs{:};
    N=length(recs);
    for n=1:N

        recName=recs(n); %File name
        thisTimeStamp=regexp(recName,'/s\d\d\d\d\d-','end');
        if(isempty(recName{1}))
            %Empty name, keep moving...
           continue;
        end

        thisTimeStamp=recName{1,1}(1,thisTimeStamp{1}+1:end);
        if(strcmp(thisTimeStamp(end),'n'))
            thisTimeStamp(end)=[];
        end

        %Convert to datenum
        thisDateNum=datenum(thisTimeStamp,dateFormat);

        %Continue if this record started at or before beginTime
        if(thisDateNum <= beginTime)

            %Get record length and check its duration
            recList=[dBName recName{:}];
            sigInfo=wfdbdesc(recList);
            Fs=sigInfo(1).SamplingFrequency;
            waveLengthSamples=sigInfo(1).LengthSamples;

            %Get the starting/ending offsets in samples wrt the beginning of
            %the waveform
            start_offset=round(Fs*(beginTime-thisDateNum)*24*60*60) + 1;
            end_offset=round(Fs*(endTime-thisDateNum)*24*60*60) + 1;

            if(end_offset <= waveLengthSamples)
                %Found a match. Get the waveform and exit the the search
                %Get signal and exit
                [signal,Fs,tm]=rdsamp(recList,[],end_offset,start_offset);
                break;
            end
        end %of thisDateNum< beginTime
    end %of matched record sublit
end

for n=1:nargout
    eval(['varargout{n}=' outputs{n} ';']);
end

author: ikarosilva

repo: wfdb-app-toolbox


text:
function f = comp(varargin)
% FUNC = COMP(FUNCTIONS+)
%   Form the composition of FUNCTIONS, resulting in FUNC.
%

  functions = varargin;
  function varargout = func(varargin)
    in = varargin;
    for i = 1:numel(functions)
      func = functions{i};
      if iscell(func)
        [func, nout] = deal(func{:});
      else
        nout = max(nargout(func), 1);
      end
      out = cell(1, nout);
      [out{:}] = func(in{:});
      in = out;
    end

    [varargout{1:max(nargout, 1)}] = deal(in{:});
  end

  f = @func;
end

author: Sable

repo: mcbench-benchmarks


text:
function signal = flt_project(varargin)
% Spatially project the given data set, e.g. to apply an IC decomposition
% Signal = flt_project(Signal, ProjectionMatrix, ChannelNames)
%
% In:
%   Signal           : epoched or continuous EEGLAB data set
%
%   ProjectionMatrix : projection matrix to apply; can be one of the following:
%                      * '.icaweights*.icasphere' : apply ICA forward projection
%                      * '.icawinv' : apply ICA back-projection
%                      * string : use the named variable or expression in the base workspace
%                      * matrix : use the given matrix as-is
%                      * anonymous function of 1 argument as string: compute spatial filter with N channels
%
%   ChannelNames     : optional cell array of new channel names (default: {'1','2','3',...})
%                      If this is set to false, the old channel labels will be retained.
%
%   ComponentSubset  : List of component indicies to which the result shall be restricted; can also
%                      be expressed as fractional intervals as in [0.1 0.3; 0.7 0.9], denoting components
%                      from 10% to 30%, and 70% to 90% of the number of components. (default: [] = retain all);
%
%   ChannelSubset    : List of channel indices to which the data shall be restricted prior to
%                      application of the matrix (default: [] = retain all);
%
% Out:
%   Signal : projected EEGLAB data set
%
% Examples:
%   % project onto a 10-dimensional random subspace, assuming that there are 32 channels in the data
%   eeg = flt_project(eeg,randn(10,32))
%
%   % project onto a 10-dimensional random subspace, and pass some alphabetic channel names
%   eeg = flt_project(eeg,randn(10,32),{'A','B','C','D','E','F','G','H','I','J'})
%
%   % project using some independent component decomposition
%   eeg = flt_project(eeg,eeg.icaweights*eeg.icasphere)
%
%   % project using some independent component decomposition, passing arguments by name
%   eeg = flt_project('Signal',eeg,'ProjectionMatix',eeg.icaweights*eeg.icasphere)
%
% See also:
%   flt_ica, flt_stationary, flt_laplace, flt_selchans
%
% TODO:
%   Try to clean up somewhat.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-03-28

% flt_project_version<1.0> -- for the cache

if ~exp_beginfun('filter') return; end;

% would be reverted by an ICA
declare_properties('name','Projection', 'follows',{'flt_ica','flt_selvolume'}, 'independent_channels',false, 'independent_trials',true);

arg_define(varargin,...
    arg_norep({'signal','Signal'}), ...
    arg({'projmat','ProjectionMatrix'}, '.icaweights*.icasphere', {'.icaweights*.icasphere','.icawinv'}, 'Projection matrix. The data is multiplied by this matrix, which can therefore implement any linear spatial filter. If left empty, flt_project will try to apply the ICA projection matrix, if present.','type','expression'),...
    arg({'newchans','ChannelNames'}, [], [], 'New channel names. Cell array of new channel names, if known. If empty, channels will be named 1:n; if false, the old channel labels will be retained.','type','cellstr','shape','row'), ...
    arg({'subcomps','ComponentSubset'}, [], [], 'Component subset. List of component indices to which the result shall be restricted (or []).','shape','row'), ...
    arg({'subchans','ChannelSubset'}, [], [], 'Channel subset. List of channel indices (or names) to which the data shall be restricted prior to application of the matrix.','type','expression','shape','row'));

append_online = {}; % this is a set of arguments that shall be appended during online use
if strcmp(projmat,'.icaweights*.icasphere') || isempty(projmat)
    if isfield(signal.etc,'amica')
        % project using AMICA weights
        subchans = signal.icachansind;
        for m=1:size(signal.etc.amica.W,3)
            tmp{m} = signal.etc.amica.W(:,:,m)*signal.etc.amica.S; end %#ok<AGROW>
        projmat = cat(1,tmp{:});
    elseif isfield(signal,'icaweights') && ~isempty(signal.icaweights)
        % project using ICA weights
        subchans = signal.icachansind;
        projmat = signal.icaweights*signal.icasphere;
    else
        error('No icaweights field is present.')
    end
    % make sure that we know what the used projection matrix & channel subset was during online use
    append_online = {'projmat',projmat,'subchans',subchans};
elseif strcmp(projmat,'.icawinv')
    if isfield(signal,'icaweights') && ~isempty(signal.icaweights) && isfield(signal,'icawinv') && ~isempty(signal.icawinv)
        subchans = signal.icachansind;
        projmat = signal.icawinv*signal.icaweights*signal.icasphere;
        newchans = 1:signal.nbchan;
    else
        error('No icaweights or icawinv fields are present.')
    end
elseif isa(projmat,'function_handle')
    randn('seed',23456);
    projmat = projmat(signal.nbchan);
elseif ischar(projmat)
    projmat = evalin('base',projmat);
end

if ~isempty(subchans)
    subset = set_chanid(signal,subchans);
    if ~isequal(subset,1:signal.nbchan)
        signal.data = signal.data(subset,:,:,:,:,:,:,:);
        signal.chanlocs = signal.chanlocs(subset);
        signal.nbchan = size(signal.data,1);
    end
end

% project data
[C,S,T] = size(signal.data); %#ok<*NODEF>
if size(projmat,2) ~= C
    error('The given projection matrix needs to have the same number of rows as the data has channels (%i), but had: %i',C,size(projmat,2)); end
signal.data = reshape(projmat*reshape(signal.data,C,[]),[],S,T);
signal.nbchan = size(signal.data,1);

% rewrite chanlocs
if isempty(signal.urchanlocs)
    signal.urchanlocs = signal.chanlocs; end
if isempty(newchans)
    signal.chanlocs = struct('labels',cellfun(@num2str,num2cell(1:signal.nbchan,1),'UniformOutput',false));
elseif length(newchans) == signal.nbchan
    if isfield(newchans,'labels')
        signal.chanlocs = newchans;
    elseif iscellstr(newchans)
        signal.chanlocs = struct('labels',newchans);
    elseif isnumeric(newchans)
        signal.chanlocs = signal.chanlocs(newchans);
    else
        error('The chanlocs format is unsupported.');
    end
elseif isequal(newchans,false) && size(projmat,1)==size(projmat,2)
    % retain current labels
else
    error('The number of provided channel labels does not match the data dimension.');
end

if ~isempty(subcomps)
    if isnumeric(subcomps) && any((subcomps(:)-floor(subcomps(:)) ~= 0)) && size(subcomps,2) == 2
        % components are given as fractional intervals
        subcomps = 1 + floor(subcomps*(signal.nbchan-1));
        tmp = [];
        for k=1:size(subcomps,1)
            tmp = [tmp subcomps(k,1):subcomps(k,2)]; end %#ok<AGROW>
        subset = set_chanid(signal,tmp);
    else
        % components are given as indices
        subset = set_chanid(signal,subcomps);
    end
    if ~isequal(subset,1:signal.nbchan)
        signal.data = signal.data(subset,:,:,:,:,:,:,:);
        signal.chanlocs = signal.chanlocs(subset);
        signal.nbchan = size(signal.data,1);
    end
end

exp_endfun('append_online',append_online);

repo: REST


text:
function value = p37_autonomous ( )

%*****************************************************************************80
%
%% P37_AUTONOMOUS reports whether problem P37 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 0;

  return
end


text:
clear all, close all, clc
%% Generate Data
Beta = [10; 28; 8/3]; % Lorenz's parameters (chaotic)
n = 3;
x0=[-8; 8; 27];  % Initial condition
tspan=[.01:.01:50];
options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,3));
[t,x]=ode45(@(t,x) lorenz(t,x,Beta),tspan,x0,options);

%% Compute Derivative
for i=1:length(x)
    dx(i,:) = lorenz(0,x(i,:),Beta);
end

%% Pool Data  (i.e., build library of nonlinear time series)
polyorder = 3;
Theta = poolData(x,n,polyorder);
m = size(Theta,2);

%% Compute Sparse regression: sequential least squares
lambda = 0.025;      % lambda is our sparsification knob.
Xi = sparsifyDynamics(Theta,dx,lambda,n)
poolDataLIST({'x','y','z'},Xi,n,polyorder);

author: dynamicslab


text:
%
% Usage:  [W [W2]]=mexStochasticProx(y,X,W0,param);
%
% Name: mexStochasticProx
%
% Description: mexStochasticProx implements a proximal MM stochastic algorithm
% for composite optimization in a large scale setting.
%        X is a design matrix of size p x n
%        y is a vector of size n
% WARNING, X is transposed compared to the functions mexFista*, and y is a vector
%        param.lambda is a vector that contains nlambda different values of the
%        regularization parameter (it can be a scalar, in that case nlambda=1)
%        W0: is a dense matrix of size p x nlambda   It is in fact ineffective
%        in the current release
%        W: is the output, dense matrix of size p x nlambda
%
%         - if param.loss='square' and param.regul corresponds to a regularization
%           function (currently 'l1' or 'l2'), the following problem is solved
%           w = argmin (1/n)sum_{i=1}^n 0.5(y_i- x_i^T w)^2 + lambda psi(w)
%         - if param.loss='logistic' and param.regul corresponds to a regularization
%           function (currently 'l1' or 'l2'), the following problem is solved
%           w = argmin (1/n)sum_{i=1}^n log(1+ exp(-y_ix_i^T w)) + lambda psi(w)
%           Note that here, the y_i's should be -1 or +1
%
%         The current release does not handle intercepts
%
% Inputs: y: double dense vector of size n
%         X: dense or sparse matrix of size p x n
%         W0: dense matrix of size p x nlambda
%         param: struct
%           param.loss (choice of loss)
%           param.regul (choice of regularization function)
%           param.lambda : vector of size nlambda
%           param.iters : number of iterations (n corresponds to one pass over the data)
%           param.minibatches: size of the mini-batches: recommended value is 1
%           param.normalized : (optional, can be set to true if the x_i's have
%              unit l2-norm, false by default)
%           param.weighting_mode : (optional, 1 by default),
%                0:  w_t = (t_0+1)/(t+t_0)
%                1:  w_t = ((t_0+1)/(t+t_0))^(0.75)
%                2:  w_t = ((t_0+1)/(t+t_0))^(5)
%           param.averaging_mode: (optional, false by default)
%                0: no averaging
%                1: first averaging mode for W2
%                2: second averaging mode for W2
%                WARNING: averaging can be very slow for sparse solutions
%           param.determineEta (optional, automatically choose the parameters of the
%             learning weights w_t, true by default)
%           param.t0 (optional, set up t0 for weights w_t = ((1+t0)/(t+t0))^(alpha)
%           param.numThreads (optional, number of threads)
%           param.verbose (optional)
%           param.seed (optional, choice of the random seed)
%
% Output:  W:  double dense p x nlambda matrix (contains the solution without averaging)
%          W2:  double dense p x nlambda matrix (contains the solution with averaging)
%
% Author: Julien Mairal, 2013

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [Mean2, Var2] = online_Update(Mean, Var, x, n)
%Using this function to update mean and variance
%
%Input
%	Mean	: average of previous data
%	Var     : variance of previous data
%	x		: new arriving data
%	n		: number of previous data
%
%Output
%	Mean2	: result of average by online updating average with new data
%	Var2	: result of variance by online updating average with new data
%

	dim = size(Var);
	if n <= 0
		Mean = x;
		Var = ones(dim)*0.5;
	elseif n > 0
		Mean2 = ( n * Mean + x ) / ( n + 1);
		Var2 = max((( n * ( Var + Mean.^2 ) + x.^2 ) / (n + 1) - Mean2.^2),0);
	end
end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function mask = compute_inpainting_mask(type,n,options)

% compute_inpainting_mask - compute an inpainting mask.
%
%   mask = compute_inpainting_mask(type,n,options);
%
%   mask=1 : data removed
%
%   Copyright (c) 2008 Gabriel Peyre

options.null = 0;
rho = getoptions(options, 'rho',.8);
qsub = getoptions(options, 'qsub',2);
name = getoptions(options, 'name','');

rand('state',123456);
switch type
    case 'rand'
        sel = randperm(n^2);
        mask = zeros(n);
        mask(sel(1:round(rho*end)))=1;
    case 'mask'
        n0 = getoptions(options, 'n0',1,1);
        mask = load_image([name '-mask'], n0);
        mask = mean(mask,3);
        mask = rescale( crop(mask,n), 0,1 );
        mask = rescale(mask)<.5;
    case 'checkboard'
        x = floor(0:1/(qsub-1):n); x= x(1:n);
        [Y,X] = meshgrid(x,x);
        mask = mod( X+Y,2 )==0;
    case 'superresol'
        [Y,X] = meshgrid(0:n-1,0:n-1);
        mask = (mod(X,qsub)==0) & (mod(Y,qsub)==0);
        mask = 1-mask;
    case {'grids' 'grids-inv'}
        [Y,X] = meshgrid(0:n-1,0:n-1);
        ngrid = getoptions(options, 'ngrid', 6);
        width = getoptions(options, 'width', 8);
        delta = n/ngrid;
        s = round(delta/2:delta:n);
        mask = zeros(n);
        for i=1:length(s);
            mask( abs(X-s(i))<=width/2 ) = 1;
        end
        mask = (mask + mask')>0;
        if strcmp(type, 'grids-inv')
            mask = 1-mask;
        end
    otherwise
        error('Unknown mask');
end

author: gpeyre

repo: matlab-toolboxes


text:
% TP model transformation based controller design for the parallel-type double inverted pendulum
% Szabolcs Nagy, Zoltan Petres, Peter Baranyi
% FUZZ-IEEE 2008 p1374-1380

m_k = 1.0; % kg
m_1 = 0.3; % kg
m_2 = 0.1; % kg
g   = 9.8; % m/s^2
l_1 = 0.6; % m
l_2 = 0.2; % m

F = @(p)(1-3/4*cos(p(1))^2)*m_1;
G = @(p)(1-3/4*cos(p(2))^2)*m_2;
H = @(p)4/3*(m_k+F(p)+G(p));

% state vector: pos, a, b, dpos, da, db
% parameter vector: a, b, da, db
LPV = {...
  @(p)0 @(p)0 @(p)0 @(p)1 @(p)0 @(p)0 @(p)0;
  @(p)0 @(p)0 @(p)0 @(p)0 @(p)1 @(p)0 @(p)0;
  @(p)0 @(p)0 @(p)0 @(p)0 @(p)0 @(p)1 @(p)0;
  @(p)0 @(p)-m_1*g*sinc(p(1)/pi)*cos(p(1))/H(p)                  @(p)-m_2*g*sinc(p(2)/pi)*cos(p(2))/H(p)                  @(p)0 @(p)4/3*m_1*l_1*p(3)*sin(p(1))/H(p)            @(p)4/3*m_2*l_2*p(4)*sin(p(2))/H(p)            @(p)4/3/H(p);
  @(p)0 @(p)(m_k+m_1+G(p))*g*sinc(p(1)/pi)/l_1/H(p)              @(p)3/4*m_2*g*sinc(p(2)/pi)*cos(p(2))*cos(p(1))/l_1/H(p) @(p)0 @(p)-m_1*p(3)*sin(p(1))*cos(p(1))/H(p)         @(p)-m_2*l_2*p(4)*sin(p(2))*cos(p(1))/l_1/H(p) @(p)-cos(p(1))/l_1/H(p);
  @(p)0 @(p)3/4*m_1*g*sinc(p(1)/pi)*cos(p(1))*cos(p(2))/l_2/H(p) @(p)(m_k+m_2+F(p))*g*sinc(p(2)/pi)/l_2/H(p)              @(p)0 @(p)-m_1*l_1*p(3)*sin(p(1))*cos(p(2))/l_2/H(p) @(p)-m_2*p(4)*sin(p(2))*cos(p(2))/H(p)         @(p)-cos(p(2))/l_2/H(p);
};

%Parameter relation tensor
dep = zeros([size(LPV) 4]);
dep(4,2,:) = [1 1 0 0];
dep(4,3,:) = [1 1 0 0];
dep(4,5,:) = [1 1 1 0];
dep(4,6,:) = [1 1 0 1];
dep(4,7,:) = [1 1 0 0];
dep(5,2,:) = [1 1 0 0];
dep(5,3,:) = [1 1 0 0];
dep(5,5,:) = [1 1 1 0];
dep(5,6,:) = [1 1 0 1];
dep(5,7,:) = [1 1 0 0];
dep(6,2,:) = [1 1 0 0];
dep(6,3,:) = [1 1 0 0];
dep(6,5,:) = [1 1 1 0];
dep(6,6,:) = [1 1 0 1];
dep(6,7,:) = [1 1 0 0];

n = 6;


text:
function l = legendre_symbol ( q, p )

%*****************************************************************************80
%
%% LEGENDRE_SYMBOL evaluates the Legendre symbol (Q/P).
%
%  Discussion:
%
%    Let P be an odd prime.  Q is a QUADRATIC RESIDUE modulo P
%    if there is an integer R such that R**2 = Q ( mod P ).
%    The Legendre symbol ( Q / P ) is defined to be:
%
%      + 1 if Q ( mod P ) /= 0 and Q is a quadratic residue modulo P,
%      - 1 if Q ( mod P ) /= 0 and Q is not a quadratic residue modulo P,
%        0 if Q ( mod P ) == 0.
%
%    We can also define ( Q / P ) for P = 2 by:
%
%      + 1 if Q ( mod P ) /= 0
%        0 if Q ( mod P ) == 0
%
%  Example:
%
%    (0/7) =   0
%    (1/7) = + 1  ( 1^2 = 1 mod 7 )
%    (2/7) = + 1  ( 3^2 = 2 mod 7 )
%    (3/7) = - 1
%    (4/7) = + 1  ( 2^2 = 4 mod 7 )
%    (5/7) = - 1
%    (6/7) = - 1
%
%  Note:
%
%    For any prime P, exactly half of the integers from 1 to P-1
%    are quadratic residues.
%
%    ( 0 / P ) = 0.
%
%    ( Q / P ) = ( mod ( Q, P ) / P ).
%
%    ( Q / P ) = ( Q1 / P ) * ( Q2 / P ) if Q = Q1 * Q2.
%
%    If Q is prime, and P is prime and greater than 2, then:
%
%      if ( Q == 1 ) then
%
%        ( Q / P ) = 1
%
%      else if ( Q == 2 ) then
%
%        ( Q / P ) = + 1 if mod ( P, 8 ) = 1 or mod ( P, 8 ) = 7,
%        ( Q / P ) = - 1 if mod ( P, 8 ) = 3 or mod ( P, 8 ) = 5.
%
%      else
%
%        ( Q / P ) = - ( P / Q ) if Q = 3 ( mod 4 ) and P = 3 ( mod 4 ),
%                  =   ( P / Q ) otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Pinter,
%    A Book of Abstract Algebra,
%    McGraw Hill, 1982, pages 236-237.
%
%    Daniel Zwillinger,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996, pages 86-87.
%
%  Parameters:
%
%    Input, integer Q, an integer whose Legendre symbol with
%    respect to P is desired.
%
%    Input, integer P, a prime number, greater than 1, with respect
%    to which the Legendre symbol of Q is desired.
%
%    Output, integer L, the Legendre symbol (Q/P).
%    Ordinarily, L will be -1, 0 or 1.
%    L = -2, P is less than or equal to 1.
%    L = -3, P is not prime.
%    L = -4, the internal stack of factors overflowed.
%    L = -5, not enough factorization space.
%

%
%  P must be greater than 1.
%
  if ( p <= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_SYMBOL - Fatal error!\n' );
    fprintf ( 1, '  P must be greater than 1.\n' );
    l = -2;
    error ( 'LEGENDRE_SYMBOL - Fatal error!' );
  end
%
%  P must be prime.
%
  if ( ~i4_is_prime ( p ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_SYMBOL - Fatal error!\n' );
    fprintf ( 1, '  P is not prime.\n' );
    l = -3;
    error ( 'LEGENDRE_SYMBOL - Fatal error!' );
  end
%
%  ( k*P / P ) = 0.
%
  if ( mod ( q, p ) == 0 )
    l = 0;
    return
  end
%
%  For the special case P = 2, (Q/P) = 1 for all odd numbers.
%
  if ( p == 2 )
    l = 1;
    return
  end
%
%  Make a copy of Q, and force it to be nonnegative.
%
  qq = q;

  while ( qq < 0 )
    qq = qq + p;
  end

  nstack = 0;
  pp = p;
  l = 1;

  while ( 1 )

    qq = mod ( qq, pp );
%
%  Decompose QQ into factors of prime powers.
%
    [ nfactor, factor, power, nleft ] = i4_factor ( qq );

    if ( nleft ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'LEGENDRE_SYMBOL - Fatal error!\n' );
      fprintf ( 1, '  Not enough factorization space.\n' );
      l = -5;
      error ( 'LEGENDRE_SYMBOL - Fatal error!' );
    end
%
%  Each factor which is an odd power is added to the stack.
%
    nmore = 0;

    for i = 1 : nfactor

      if ( mod ( power(i), 2 ) == 1 )

        nmore = nmore + 1;
        nstack = nstack + 1;

        pstack(nstack) = pp;
        qstack(nstack) = factor(i);

      end

    end

    if ( nmore ~= 0 )

      qq = qstack(nstack);
      nstack = nstack - 1;
%
%  Check for a QQ of 1 or 2.
%
      if ( qq == 1 )

        l = + 1 * l;

      elseif ( qq == 2 & ( mod ( pp, 8 ) == 1 | mod ( pp, 8 ) == 7 ) )

        l = + 1 * l;

      elseif ( qq == 2 & ( mod ( pp, 8 ) == 3 | mod ( pp, 8 ) == 5 ) )

        l = - 1 * l;

      else

        if ( mod ( pp, 4 ) == 3 & mod ( qq, 4 ) == 3 )
          l = - 1 * l;
        end

        [ pp, qq ] = i4_swap ( pp, qq );

        continue

      end

    end
%
%  If the stack is empty, we're done.
%
    if ( nstack == 0 )
      break
    end
%
%  Otherwise, get the last P and Q from the stack, and process them.
%
    pp = pstack(nstack);
    qq = qstack(nstack);
    nstack = nstack - 1;

  end

  return
end


text:
function [C_pos,C_neg,Ctot_pos,Ctot_neg] = clustering_coef_wu_sign(W,coef_type)
%CLUSTERING_COEF_WU_SIGN     Multiple generalizations of the clustering coefficient
%
%   [C_pos,C_neg,Ctot_pos,Ctot_neg] = clustering_coef_wu_sign(W,coef_type);
%
%   The weighted clustering coefficient is the average weight or intensity
%   of all triangles associated with each node.
%
%   Inputs:
%       W,
%           Weighted undirected connection matrix
%
%       corr_type,
%           Desired type of clustering coefficient.
%           Options:
%           1,  (default) Onnela et al. formula, used in original
%               clustering_coef_wu.m. Computed separately for positive &
%               negative weights.
%           2,  Zhang & Horvath formula, similar to Onnela formula except
%               denominator of Onnela formula relies on binarizing the
%               network whereas this denominator is based on weight value,
%               which reduces the sensitivity of this measure to the
%               weights directly connected to the node of interest.
%               Computed separately for positive & negative weights.
%           3,  Constantini & Perugini's generalization of the Zhang &
%               Horvath formula. This formula takes both positive &
%               negative weights into account simultaneously, & is
%               particularly sensitive to non-redundnacy in path
%               information based on sign (i.e., when two weights are
%               positive & one negative, or all three are negative, both of
%               which indicate that the weight of the third path is not
%               redundant information). Produces only one value.
%
%
%   Outputs:
%       C_pos/C_neg,
%           Clustering coefficient vector for positive/negative weights.
%           For the third option, only one vector is outputted (as C_pos).
%       Ctot_pos/Ctot_neg,
%           Mean clustering coefficient for positive and negative weights.
%
%   References:
%       Onnela et al. (2005) Phys Rev E 71:065103
%       Zhang & Horvath (2005) Stat Appl Genet Mol Biol 41:1544-6115
%       Costantini & Perugini (2014) PLOS ONE 9:e88669
%
%
%   Contributor: Jeff Spielberg, Boston University, 2014-2015
%                (script based on clustering_coef_wu.m)

%
%   Modification History:
%   May 2014: Added computation of pos & neg weights separately &
%             computation of mean coefficient (Jeff Spielberg)
%   May 2015: Added computation of Zhang & Horvath and Constantini &
%             Perugini formulas (Jeff Spielberg)

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

n            = length(W);                   %number of nodes
W(1:n+1:end) = 0;

switch coef_type
    case 1
        W_pos                = W.*(W>0);
        K_pos                = sum(W_pos~=0,2);
        cyc3_pos             = diag((W_pos.^(1/3))^3);
        K_pos(cyc3_pos == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_pos                = cyc3_pos./(K_pos.*(K_pos-1));         %clustering coefficient
        Ctot_pos             = mean(C_pos);

        W_neg                = -W.*(W<0);
        K_neg                = sum(W_neg~=0,2);
        cyc3_neg             = diag((W_neg.^(1/3))^3);
        K_neg(cyc3_neg == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_neg                = cyc3_neg./(K_neg.*(K_neg-1));         %clustering coefficient
        Ctot_neg             = mean(C_neg);
    case 2
        W_pos    = W.*(W>0);
        cyc3_pos = zeros(n,1);
        cyc2_pos = zeros(n,1);
        for i = 1:n
            for j = 1:n
                for q = 1:n
                    cyc3_pos(i) = cyc3_pos(i)+(W_pos(j,i)*W_pos(i,q)*W_pos(j,q));
                    if j~=q
                        cyc2_pos(i) = cyc2_pos(i)+(W_pos(j,i)*W_pos(i,q));
                    end
                end
            end
        end
        cyc2_pos(cyc3_pos == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_pos                   = cyc3_pos./cyc2_pos;         %clustering coefficient
        Ctot_pos                = mean(C_pos);

        W_neg    = -W.*(W<0);
        cyc3_neg = zeros(n,1);
        cyc2_neg = zeros(n,1);
        for i = 1:n
            for j = 1:n
                for q = 1:n
                    cyc3_neg(i) = cyc3_neg(i)+(W_neg(j,i)*W_neg(i,q)*W_neg(j,q));
                    if j~=q
                        cyc2_neg(i) = cyc2_neg(i)+(W_neg(j,i)*W_neg(i,q));
                    end
                end
            end
        end
        cyc2_neg(cyc3_neg == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_neg                   = cyc3_neg./cyc2_neg;         %clustering coefficient
        Ctot_neg                = mean(C_neg);
    case 3
        cyc3         = zeros(n,1);
        cyc2         = zeros(n,1);

        for i = 1:n
            for j = 1:n
                for q = 1:n
                    cyc3(i) = cyc3(i)+(W(j,i)*W(i,q)*W(j,q));
                    if j~=q
                        cyc2(i) = cyc2(i)+(W(j,i)*W(i,q));
                    end
                end
            end
        end

        cyc2(cyc3 == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_pos           = cyc3./cyc2;         %clustering coefficient
        Ctot_pos        = mean(C_pos);
        C_neg           = nan(size(C_pos));
        Ctot_neg        = nan(size(Ctot_pos));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function out = getqr(a)
% =======================================================================
% QR factorization of a matrix a.
% =======================================================================
% out = getqr(a)
% -----------------------------------------------------------------------
% INPUT
%	- a: Matrix to factorize
% -----------------------------------------------------------------------
% OUTPUT
%   - out: factorized matrix such that out*out'=I
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% December 2021
% -----------------------------------------------------------------------

[q,r]=qr(a);
for i=1:size(q,1)
    if r(i,i)<0
        q(:,i)=-q(:,i);
    end
end
out=q;
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function CPD = tabular_CPD(bnet, self, varargin)
% TABULAR_CPD Make a multinomial conditional prob. distrib. (CPT)
%
% CPD = tabular_CPD(bnet, node) creates a random CPT.
%
% The following arguments can be specified [default in brackets]
%
% CPT - specifies the params ['rnd']
%   - T means use table T; it will be reshaped to the size of node's family.
%   - 'rnd' creates rnd params (drawn from uniform)
%   - 'unif' creates a uniform distribution
% adjustable - 0 means don't adjust the parameters during learning [1]
% prior_type - defines type of prior ['none']
%  - 'none' means do ML estimation
%  - 'dirichlet' means add pseudo-counts to every cell
%  - 'entropic' means use a prior P(theta) propto exp(-H(theta)) (see Brand)
% dirichlet_weight - equivalent sample size (ess) of the dirichlet prior [1]
% dirichlet_type - defines the type of Dirichlet prior ['BDeu']
%  - 'unif' means put dirichlet_weight in every cell
%  - 'BDeu' means we put 'dirichlet_weight/(r q)' in every cell
%    where r = self_sz and q = prod(parent_sz) (see Heckerman)
% trim - 1 means trim redundant params (rows in CPT) when using entropic prior [0]
% entropic_pcases - list of assignments to the parents nodes when we should use
%      the entropic prior; all other cases will be estimated using ML [1:psz]
% sparse - 1 means use 1D sparse array to represent CPT [0]
%
% e.g., tabular_CPD(bnet, i, 'CPT', T)
% e.g., tabular_CPD(bnet, i, 'CPT', 'unif', 'dirichlet_weight', 2, 'dirichlet_type', 'unif')
%
% REFERENCES
% M. Brand - "Structure learning in conditional probability models via an entropic  prior
%   and parameter extinction", Neural Computation 11 (1999): 1155--1182
% M. Brand - "Pattern discovery via entropy minimization" [covers annealing]
%   AI & Statistics 1999. Equation numbers refer to this paper, which is available from
%   www.merl.com/reports/docs/TR98-21.pdf
% D. Heckerman, D. Geiger and M. Chickering,
%   "Learning Bayesian networks: the combination of knowledge and statistical data",
%   Microsoft Research Tech Report, 1994

if nargin==0
  % This occurs if we are trying to load an object from a file.
  CPD = init_fields;
  CPD = class(CPD, 'tabular_CPD', discrete_CPD(0, []));
  return;
elseif isa(bnet, 'tabular_CPD')
  % This might occur if we are copying an object.
  CPD = bnet;
  return;
end
CPD = init_fields;

ns = bnet.node_sizes;
ps = parents(bnet.dag, self);
fam_sz = ns([ps self]);
psz = prod(ns(ps));
CPD.sizes = fam_sz;
CPD.leftright = 0;
CPD.sparse = 0;

% set defaults
CPD.CPT = mk_stochastic(myrand(fam_sz));
CPD.adjustable = 1;
CPD.prior_type = 'none';
dirichlet_type = 'BDeu';
dirichlet_weight = 1;
CPD.trim = 0;
selfprob = 0.1;
CPD.entropic_pcases = 1:psz;

% extract optional args
args = varargin;
% check for old syntax CPD(bnet, i, CPT) as opposed to CPD(bnet, i, 'CPT', CPT)
if ~isempty(args) & ~isstr(args{1})
  CPD.CPT = myreshape(args{1}, fam_sz);
  args = [];
end

for i=1:2:length(args)
  switch args{i},
   case 'CPT',
    T = args{i+1};
    if ischar(T)
      switch T
       case 'unif', CPD.CPT = mk_stochastic(myones(fam_sz));
       case 'rnd',  CPD.CPT = mk_stochastic(myrand(fam_sz));
       otherwise,   error(['invalid CPT ' T]);
      end
    else
      CPD.CPT = myreshape(T, fam_sz);
    end
   case 'prior_type', CPD.prior_type = args{i+1};
   case 'dirichlet_type', dirichlet_type = args{i+1};
   case 'dirichlet_weight', dirichlet_weight = args{i+1};
   case 'adjustable', CPD.adjustable = args{i+1};
   case 'clamped', CPD.adjustable = ~args{i+1};
   case 'trim', CPD.trim = args{i+1};
   case 'entropic_pcases', CPD.entropic_pcases = args{i+1};
   case 'sparse', CPD.sparse = args{i+1};
   otherwise, error(['invalid argument name: ' args{i}]);
  end
end

switch CPD.prior_type
 case 'dirichlet',
  switch dirichlet_type
   case 'unif', CPD.dirichlet = dirichlet_weight * myones(fam_sz);
   case 'BDeu',  CPD.dirichlet = (dirichlet_weight/psz) * mk_stochastic(myones(fam_sz));
   otherwise, error(['invalid dirichlet_type ' dirichlet_type])
  end
 case {'entropic', 'none'}
  CPD.dirichlet = [];
 otherwise, error(['invalid prior_type ' prior_type])
end

% fields to do with learning
if ~CPD.adjustable
  CPD.counts = [];
  CPD.nparams = 0;
  CPD.nsamples = [];
else
  %CPD.counts = zeros(size(CPD.CPT));
  CPD.counts = zeros(prod(size(CPD.CPT)), 1);
  psz = fam_sz(1:end-1);
  ss = fam_sz(end);
  if CPD.leftright
    % For each of the Qps contexts, we specify Q elements on the diagoanl
    CPD.nparams = Qps * Q;
  else
    % sum-to-1 constraint reduces the effective arity of the node by 1
    CPD.nparams = prod([psz ss-1]);
  end
  CPD.nsamples = 0;
end

CPD.trimmed_trans = [];
fam_sz = CPD.sizes;

%psz = prod(fam_sz(1:end-1));
%ssz = fam_sz(end);
%CPD.trimmed_trans = zeros(psz, ssz); % must declare before reading

%sparse CPT
if CPD.sparse
   CPD.CPT = sparse(CPD.CPT(:));
end

CPD = class(CPD, 'tabular_CPD', discrete_CPD(~CPD.adjustable, fam_sz));

%%%%%%%%%%%

function CPD = init_fields()
% This ensures we define the fields in the same order
% no matter whether we load an object from a file,
% or create it from scratch. (Matlab requires this.)

CPD.CPT = [];
CPD.sizes = [];
CPD.prior_type = [];
CPD.dirichlet = [];
CPD.adjustable = [];
CPD.counts = [];
CPD.nparams = [];
CPD.nsamples = [];
CPD.trim = [];
CPD.trimmed_trans = [];
CPD.leftright = [];
CPD.entropic_pcases = [];
CPD.sparse = [];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function plot(a,d)
%
% Usage :
%           plot(a)
%
% a -- The algorithm already trained
%      The colormap can be changed afterwards.
%
% Example:
%   d=gen(spiral({'n=1','m=50'}));
%   [r s0]=train(svm(kernel('rbf',1)),d)
%   plot(s0);

%clf
%d=a.keep;
x=get_x(d);
y=get_y(d);
ax1=min(x); ax2=max(x);

granul=100;
minX=floor(ax1(1));
maxX=ceil(ax2(1));
minY=floor(ax1(2));
maxY=ceil(ax2(2));
axis_sz=[minX maxX minY maxY];
minX=minX-2; maxX=maxX+2; minY=minY-2; maxY=maxY+2;
gridx=[]; gridy=[];

mx=zeros(granul*granul,2);
for i=1:granul
 for j=1:granul
   mx((i-1)*granul+j,:)= [minX+(i-1)*(maxX-minX)/granul minY+(j-1)*(maxY-minY)/granul ] ;
   gridx(i)=minX+(i-1)*(maxX-minX)/granul;
   gridy(j)=minY+(j-1)*(maxY-minY)/granul;
end
end

temp=zeros(granul,granul);
dx=data(mx);
N=get_dim(dx);
blocks=N/50;
rr=[];
last=1;

% a.algorithm.use_signed_output=0;
%a.return_indices=1;
allx=[];
for i=blocks:blocks:N
       resX=test(a,get(dx,[last:i]));
       rr=[rr;resX.X'];allx=[allx;resX.X];
       last=i+1;
end

for i=1:granul
 for j=1:granul
  temp(i,j)= rr( (i-1)*granul+j);
  end
end
hold on;
%surf(1:granul,1:granul,temp'-1000);

% Ipos=find(allx(:,1)==+1);
% Ineg=find(allx(:,1)==-1);

% plot(mx(Ipos,1),mx(Ipos,2),'r.');
% plot(mx(Ineg,1),mx(Ineg,2),'b.');

%clf

% colormap('cool')
 pcolor(gridx, gridy, temp) ;
 shading interp;

%surf(gridx,gridy,temp');
%view(2);
%shading interp;
%[c,h]=contour(gridx,gridy,temp',[0 0],'k');
%   clabel(c,h);
colorbar

if(1)

    Ipos=find(y==+1);
    Ineg=find(y==-1);
    h=plot(x(Ipos,1),x(Ipos,2),'rx'); hold on;
    set(h,'LineWidth',2,'MarkerSize',5);

    h=plot(x(Ineg,1),x(Ineg,2),'ko'); hold on;
    set(h,'LineWidth',2,'MarkerSize',5);
end

axis(axis_sz);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [alpha, beta, gamma, stability] = evalABGParam(process, noisy, dt)
% evalABGParam - evaluates alpha, beta and gamma parameters
% With this parameters, alpha-beta filter becomes a steady-state Kalman filter
%
% Syntax:  [alpha,beta,gamma,stability] = evalABGParam(process,noisy,dt)
%          [alpha,beta,gamma,stability] = evalABGParam([alpha,beta,gamma])
%
% Inputs:
%   process - real system state
%     noisy - measured system state
%        dt - delta time (sample rate)
%
% Outputs:
%     alpha - alpha parameter
%      beta - beta parameter
%     gamma - gamma parameter
% stability - Juri's Stability Test
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: abgFilter;

% Author: Marco Borges, Ph.D. Student, Computer/Biomedical Engineer
% UFMG, PPGEE, Neurodinamica Lab, Brazil
% email address: marcoafborges@gmail.com
% Website: http://www.cpdee.ufmg.br/
% References:
%   NEAL, S. R., Parametric relations for a-b-g filter predictor, IEEE
%      Trans. on Automatic Control, AC-12, June 1967
%   TENNE, D. and Singh, T., Characterizing Performance a-b-g Filters, IEEE
%      Transactions on Aerospace and Electronic Systems, 38, 2002
% September 2013; Version: v1; Last revision: 2013-09-18
% Changelog:
%
%------------------------------- BEGIN CODE -------------------------------

if nargin == 3
    varProcess = var(process);
    varNoise = var(noisy);
    l = varProcess * dt / varNoise; % lambda
    r = (4+l-sqrt(8*l+l^2))/4;
    alpha = 1 - r^2;
    beta = 2*(2-alpha)-4*sqrt(1-alpha);
    gamma = beta^2/(2*alpha);
elseif nargin == 1 && length(process) == 3
    alpha = process(1);
    beta = process(2);
    gamma = process(3);
else
    error('evalABGParam : Incorrect Parameters!');
end

if ( alpha > 0 && alpha < 2 && beta > 0 && beta < (4-2*alpha) && ...
        gamma > 0 && gamma < (4*alpha*beta)/(2*alpha) )
    stability = 'Stable';
else
    stability = 'Unstable';
end
%-------------------------------- END CODE --------------------------------


text:
function w = wf_fft(w);

% WAVEFORM = WF_FFT(WAVEFORM)  calculate frequency spectrum of a waveform.
% The results of the fast fourier transform is added as new fields
% in the waveform:
%   FFT_FREQ is a vector of frequencies with N samples
%   FFT_AMP is a vector of spetral amplitudes
%   FFT_PHASE is a vector of phases
%   FFT_DOM is the scalar frequency of the maximum amplitude
%           peak (or dominant frequency. This could change)

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

admin.deprecated(mfilename,'wf_fft.compute');

% CHECK ARGUMENTS
if ~strcmpi(class(w),'waveform')
    error('First input must be a waveform object');
end

% STEP THROUGH WAVFORMS ADDING NEW FIELDS
[N,M] = size(w);
for i = 1:N*M
    Fn = get(w(i),'NYQ');
    x = get(w(i),'DATA');
    NFFT=2.^(ceil(log(length(x))/log(2)));  % Next highest power of 2
    FFTX=fft(x,NFFT);                       % Take fft, padding with zeros.
    NumUniquePts = ceil((NFFT+1)/2);
    FFTX=FFTX(1:NumUniquePts);              % throw out neg frequencies
    MX=abs(FFTX);                           % Take magnitude of X
    MX=MX*2;                                % Multiply by 2
    MX=MX/length(x);
    PX=phase(FFTX);                           % Take magnitude of X
    f=(0:NumUniquePts-1)*2/NFFT;
    f=f*Fn;
    w(i) = addfield(w(i),'FFT_FREQ',f');
    w(i) = addfield(w(i),'FFT_AMP',MX);
    w(i) = addfield(w(i),'FFT_PHASE',PX);
    a = find(MX == max(MX));
    w(i) = addfield(w(i),'FFT_DOM',f(a));
end;

author: geoscience-community-codes

repo: GISMO


function xyz = readVelodyne(filename)

fid = fopen(filename,'rb');
if fid == -1
  error('%s does not exist', filename);
end
velo = uint8(fread(fid,[32 inf],'uint8'));  % 32 /18
xyz = velo(1:12,:);
xyz = xyz(:);
xyz = typecast(xyz, 'single');
xyz = reshape(xyz, 3, []);
xyz(4,:) = 1;

fclose(fid);


text:
function toms097 ( )

%*****************************************************************************80
%
%% MAIN is the main program for TOMS097_PRB.
%
%  Discussion:
%
%    TOMS097_PRB tests the TOMS097 library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS097_PRB\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TOMS097 library.\n' );

  toms097_test01 ( );
  toms097_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS097_PRB\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% Compute hub disruption indices (HDIs) using BCT functions. See Achard et
% al 2012; Mansour et al 2016.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2020 Yoni Ashar
%
%     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/>.
% ..
%
% :Inputs:
%
%   **bs:**
%        brainpathway_multisubject object. Will use connectivity.regions.r
%        for computing the HDIs.
%
%   **refgroup:**
%        Two options: 1) a brainpathway_multisubject object to use
%        as a reference group. Will compute BCT measures on refgroup.connectivity.regions.r
%        and take the mean value as a the reference distribution.
%        2) a vector on integers to use for cross-validated estimation of
%        HDIs, where each integer represents one fold. The matlab functions
%        crossvalind and cvpartition can be helpful in generating these
%        indices. The hold-out set is used as the reference distribution
%
%   **thresh:**
%        Followed by a link density, ranging from .01 - .99. Will use .1 by default.
%
% :Optional Inputs:
%    **'doweighted'**
%       Compute weighted metrics
%
% :Outputs:
%
%   **bs:**
%        Will save the HDI values in bs.graph_prop, overwriting any
%        existing data in that field
%
%   **refgroup:**
%        With the BCT measures saved in the graph_properties.region field
%
% :Examples:
% ::
%
%    load('brainpathway_data_gsr_censoring_pipeline_nosmoothing.mat');
%    load('paingen_reference_group.mat');
%
%    [bs, refgroup] = compute_HDIs(bs, bs_paingen)
%
% :References:
%   Achard et al 2012, Mansour et al 2016
%
%
% ..
%    Programmers' notes:
%    Initial commit -- Yoni Ashar, April 2020
% ..
function [bs, refgroup] = compute_HDIs(bs, refgroup, thresh, varargin)

doweighted = 0;
if any(strcmp(varargin, 'doweighted')), doweighted = true; end

%% compute BCT measures on my subjects
print_header(sprintf('Computing BCT measures on %d subjects in test group', size(bs.connectivity.regions.r,3)));
bs = bct_toolbox_undirected_graph_metrics(bs, thresh, varargin{:});

%% Compute BCT measures on ref group and regress

% what kind of reference group are we dealing with?

% CV
if isnumeric(refgroup)
    error('CV not implemented yet')

% External reference group
elseif isa(refgroup, 'brainpathway_multisubject')

    print_header(sprintf('Computing BCT measures on %d subjects in external reference group', size(refgroup.connectivity.regions.r,3)));

    %% compute BCT measures on ref group
    refgroup = bct_toolbox_undirected_graph_metrics(refgroup, thresh, varargin{:});

    %% for each subject, regress on mean referece group to get an HDI (i.e., the slope)

    for i=1:size(bs.connectivity.regions.r,3)

        % dont compute HDI on global metrics, only nodal
        % TODO: put in subfunction
        mymetrics = bs.graph_properties.regions.Properties.VariableNames;
        metrics = {};
        for m=1:length(mymetrics)
            if numel(bs.graph_properties.regions.(mymetrics{m})(1,:)) > 1
                metrics{end+1} = mymetrics{m};
            end
        end

        % compute the HDIs
        for m=1:length(metrics) % for each metric

            betas = regress(mean(refgroup.graph_properties.regions.(metrics{m}))', [ones(489, 1) bs.graph_properties.regions.(metrics{m})(i,:)' - mean(refgroup.graph_properties.regions.(metrics{m}))']);

            %figure; scatter(mean(graph_prop_reference_group.(metrics{m})), graph_prop.(metrics{m})(i,:) - mean(graph_prop_reference_group.(metrics{m}))), lsline, title(metrics{m})
            bs.HDIs.regions.(metrics{m})(i) = betas(2);
        end
    end

% Unknown
else
    error('Unknown type of reference group')
end


text:
function res = mtimes(a,x)
% performs the SPIRIT operator
% see mtimes

res = mtimes(a,x);

author: thomaskuestner


text:
function y = mrdivide(X,Y)

if (isa(Y,'sdpvar'))
    if Y.dim(1)*Y.dim(2) == 1
        y = X*Y^-1;
        return
    else
        error('Division of matrix variables not possible (maybe you meant ./).')
    end
end

try

    % Quick exit on scalar division
    if length(Y)==1
        if isinf(Y)
            y = zeros(X.dim);
        else
            y = X;
            y.basis = y.basis/Y;
        end
        y.conicinfo = [0 0];
        return
    end

    lmi_variables = getvariables(X);
    nv = length(lmi_variables);
    y  = X;
    n = X.dim(1);
    m = X.dim(2);
    if (n==1) & (m==1) % SPECIAL CODE FOR FREAKY MATLAB BUG
        y.basis = sparse(reshape(reshape(full(X.basis(:,1)),n,m)/Y,n*m,1));
        for i = 1:nv
            temp = reshape(full(X.basis(:,i+1)),n,m)/Y;
            y.basis(:,i+1) = sparse(temp(:));
        end;
    else
        y.basis = reshape(reshape(X.basis(:,1),n,m)/Y,n*m,1);
        for i = 1:nv
            temp = reshape(X.basis(:,i+1),n,m)/Y;
            y.basis(:,i+1) = temp(:);
        end;
    end

    y.dim(1) = size(temp,1);
    y.dim(2) = size(temp,2);
catch
    error(lasterr);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP