text:
function varargout = doseManagementGui(command, varargin)
%"doseManagementGui" GUI
%   Create a GUI to manage doses.
%
%   JRA 2/27/04
%
%Usage:
%   doseManagementGui()
%
% 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/>.

%The cursed globals. Need em.
global planC stateS
indexS = planC{end};

%Use a static window size, by pixels.  Do not allow resizing.
screenSize = get(0,'ScreenSize');
y = 380;
x = 640;
units = 'normalized';

%Height of a single row of text.
rowHeight = .06;

%If no command given, default to init.
if ~exist('command','var') || isempty(command)
    command = 'init';
end

%Find handle of the gui figure.
%h = findobj('tag', 'doseManagementGui');
h = [];
if isfield(stateS.handle,'doseManagementFig') && ishandle(stateS.handle.doseManagementFig)
    h = stateS.handle.doseManagementFig;
end

%Set framecolor for uicontrols and pseudoframes.
frameColor = [0.8314 0.8157 0.7843];

switch upper(command)
    case 'INIT'
        %If gui doesnt exist, create it, else refresh it.
        if isempty(h)
            %Set up a new GUI window.
            h = figure('units', 'pixels', 'position',[(screenSize(3)-x)/2 (screenSize(4)-y)/2 x y], 'MenuBar', 'none', 'NumberTitle', 'off', 'resize', 'off', 'Tag', 'doseManagementGui', 'Color', [.75 .75 .75], 'WindowButtonUpFcn', 'doseManagementGui(''FIGUREBUTTONUP'')');
            stateS.handle.doseManagementFig = h;
            set(h, 'Name','Dose Management');

            %Create pseudo frames.
            axes('Position', [.02 .05 + (rowHeight + .02) .96 .87 - (rowHeight + .02)/2], 'color', frameColor, 'ytick',[],'xtick',[], 'box', 'on');
            line([.5 .5], [0 1], 'color', 'black');
            axes('Position', [.02 .03 .96 rowHeight+.02], 'color', frameColor, 'ytick',[],'xtick',[], 'box', 'on');

            %Create controls for displaying the GUI status.
            ud.handles.status = uicontrol(h, 'units', units, 'string', 'Status:', 'Position', [.04 .03 .05 rowHeight], 'Style', 'text', 'BackgroundColor', frameColor);
            ud.handles.status = uicontrol(h, 'units', units, 'string', '', 'Position', [.1 .03 .8 rowHeight], 'Style', 'text', 'HorizontalAlignment', 'left', 'foregroundcolor', 'red', 'BackgroundColor', frameColor);
            ud.currentDose = 1;
            ud.doseBlock = 1;
            set(h, 'userdata', ud);
        end
        doseManagementGui('refresh', h);
        figure(h);

    case 'REFRESHFIELDS'
        %Refresh the field list for the current dose.
        ud = get(h, 'userdata');
        doseNum = ud.currentDose;
        set(ud.handles.name, 'string', planC{indexS.dose}(doseNum).fractionGroupID);
        set(ud.handles.units, 'string', planC{indexS.dose}(doseNum).doseUnits);
        [dA, isCompress, isRemote] = getDoseArray(doseNum, planC);
        if isCompress
            set(ud.handles.compbutton, 'string', 'Decompress');
            set(ud.handles.storageMethod, 'string', 'Memory, Compressed');
        else
            set(ud.handles.compbutton, 'string', 'Compress');
        end

        if isRemote
            set(ud.handles.remotebutton, 'string', 'Use Memory');
            set(ud.handles.storageMethod, 'string', 'On Disk');
        else
            set(ud.handles.remotebutton, 'string', 'Use Disk');
        end

        if ~isCompress & ~isRemote
            set(ud.handles.storageMethod, 'string', 'Memory');
        end
        set(ud.handles.maxdose, 'string', ud.maxDoses{doseNum});
        doseSize = getByteSize(planC{indexS.dose}(doseNum));
        set(ud.handles.dosesize, 'string', [num2str(doseSize/(1024*1024), '%6.2f') 'MB']);

    case 'REFRESH'
        %Recreate and redraw the entire doseManagementGui.
        if isempty(h)
            return;
        end

        %Save the current figure so focus can be returned.
        hOld = gcf;

        %Focus on doseManagementGui for the moment.
        set(0, 'CurrentFigure', h);
        ud = get(h, 'userdata');
        try
            %Scrap all old buttons and axes.
            delete(ud.handles.thumbaxis);
            delete(ud.handles.name);
            delete(ud.handles.dosesize);
            delete(ud.handles.storageMethod);
            delete(ud.handles.maxdose);
            delete(ud.handles.units);
            delete(ud.handles.compbutton);
            delete(ud.handles.copybutton);
            delete(ud.handles.delbutton);
            ud.handles.thumbaxis = [];
            ud.handles.name = [];
            ud.handles.dosesize = [];
            ud.handles.storageMethod = [];
            ud.handles.maxdose = [];
            ud.handles.units = [];
            ud.handles.compbutton = [];
            ud.handles.copybutton = [];
            ud.handles.delbutton = [];
        end

        i = 1;

        txtLeft = .05;
        textWidth = .13;
        fieldLeft = .35;
        fieldWidth = .13;

        %Make text to describe uicontrols.
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.15 textWidth rowHeight],'String', 'Edit Dose name:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor,'TooltipString','Change/Edit dose name string');
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.21 textWidth rowHeight],'String', 'RAM Used:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor,'TooltipString','Memory used');
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.27 textWidth rowHeight],'String', 'Storage Method:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor,'TooltipString','Tells if dose is stored on disk OR RAM');
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.33 textWidth rowHeight],'String', 'Max Dose:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor,'TooltipString','Maximum Dose Value');
        uicontrol(h, 'units',units,'Position',[txtLeft 1-.39 textWidth rowHeight],'String', 'Units:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor,'TooltipString','Dose Units');
        uicontrol(h, 'units',units,'Position',[0.24 1-.39 textWidth rowHeight],'String', 'Scale by:', 'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor, 'TooltipString', 'Scale Factor');

        %Make uicontrols for managing the doses, and displaying info.
        ud.handles.name          = uicontrol(h, 'units',units,'Position',[fieldLeft-fieldWidth 1-.15+.02 fieldWidth*2 rowHeight-.01],'String','', 'Style', 'edit', 'callback', 'doseManagementGui(''NAMEFIELD'');', 'userdata', i, 'horizontalAlignment', 'right','TooltipString','Change/Edit dose name string');
        ud.handles.dosesize      = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.21 fieldWidth rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor,'TooltipString','Memory used');
        ud.handles.storageMethod = uicontrol(h, 'units',units,'Position',[fieldLeft-.05 1-.27 fieldWidth+.05 rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor,'TooltipString','Tells if dose is stored on disk OR RAM');
        ud.handles.maxdose       = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.33+.02 fieldWidth rowHeight-.01],'String', '', 'Style', 'edit', 'callback', 'doseManagementGui(''DOSESCALE'');', 'userdata', i, 'horizontalAlignment', 'right','TooltipString','Maximum Dose Value');
        %ud.handles.units         = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.39 fieldWidth rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'right', 'BackgroundColor', frameColor);
        ud.handles.units         = uicontrol(h, 'units',units,'Position',[0.1 1-.39+0.02 fieldWidth-0.02 rowHeight-0.01],'String', '',  'Style', 'edit', 'callBack','doseManagementGui(''UNITFIELD'');','horizontalAlignment', 'left','TooltipString','Dose Units');
        ud.handles.scale         = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.39+0.02 fieldWidth rowHeight-0.01],'String', '',  'Style', 'edit', 'callBack','doseManagementGui(''SCALEBYFACTOR'');', 'horizontalAlignment', 'left','TooltipString','Scale Factor');

        ud.handles.remotebutton  = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.59 fieldWidth rowHeight],'String', 'Make Remote', 'Style', 'pushbutton', 'callback', 'doseManagementGui(''REMOTE'');', 'userdata', i,'TooltipString','Store Dose Remotely, save memory used');
        ud.handles.compbutton    = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.67 fieldWidth rowHeight],'String', 'Compress', 'Style', 'pushbutton', 'callback', 'doseManagementGui(''COMPRESS'');', 'userdata', i,'TooltipString','Compress Remote Dose');
%         ud.handles.copybutton    = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.75 fieldWidth rowHeight],'String','Save', 'Style', 'pushbutton', 'callback', 'doseManagementGui(''SAVE'');', 'userdata', i);
        ud.handles.delbutton     = uicontrol(h, 'units',units,'Position',[fieldLeft 1-.83 fieldWidth rowHeight],'String','Delete', 'Style', 'pushbutton', 'callback', 'doseManagementGui(''DELETE'');', 'userdata', i,'TooltipString','Delete Dose Permantely/ Save Plan to take effect');

        ud.handles.previewAxis   = axes('position', [txtLeft 1-.83 .25 .25*x/y], 'parent', h, 'box', 'on', 'ytick', [], 'xtick', [], 'buttondownfcn', 'doseManagementGui(''PREVIEWBUTTONDOWN'')');
        ud.handles.previewSliceNum = uicontrol(h, 'units',units,'Position',[fieldLeft-0.03 1-.47 fieldWidth rowHeight],'String', '',  'Style', 'text', 'horizontalAlignment', 'left', 'BackgroundColor', frameColor);

        text(.5, .6, 'After selecting a dose, ', 'fontsize', 8, 'hittest', 'off', 'parent', ud.handles.previewAxis, 'horizontalAlignment', 'center');
        text(.5, .5, 'click and drag mouse up', 'fontsize', 8, 'hittest', 'off', 'parent', ud.handles.previewAxis, 'horizontalAlignment', 'center');
        text(.5, .4, 'and down here to preview.', 'fontsize', 8, 'hittest', 'off', 'parent', ud.handles.previewAxis, 'horizontalAlignment', 'center');

        % Create Prev and Next Buttons
        nDoses = length(planC{indexS.dose});
        if nDoses > 9
            ud.handles.prev = uicontrol(h, 'units',units,'Position',[fieldLeft+0.18 1-.86 0.1 rowHeight],'String','<--', 'Style', 'pushbutton', 'callback', 'doseManagementGui(''PREV_BLOCK'');', 'TooltipString','Show Previous Set of Doses','fontWeight','bold');
            ud.handles.next = uicontrol(h, 'units',units,'Position',[fieldLeft+0.5 1-.86 0.1 rowHeight],'String','-->', 'Style', 'pushbutton', 'callback', 'doseManagementGui(''NEXT_BLOCK'');', 'TooltipString','Show Next Set of Doses','fontWeight','bold');
            uicontrol(h, 'units',units,'Position',[fieldLeft+0.3 1-.86 0.18 rowHeight+0.01],'String','Use these buttons to view more doses', 'Style', 'text', 'fontSize',8,'fontWeight','normal', 'BackgroundColor', frameColor);
        end

        set(h, 'userdata', ud);
        set(0, 'CurrentFigure', hOld);

        doseManagementGui('REFRESH_THUMBS')

    case 'PREV_BLOCK'
        ud = get(h, 'userdata');
        if ud.doseBlock == 1
            return;
        end
        ud.doseBlock = ud.doseBlock - 1;
        set(h, 'userdata', ud);
        doseManagementGui('REFRESH_THUMBS')

    case 'NEXT_BLOCK'
        ud = get(h, 'userdata');
        nDoses = length(planC{indexS.dose});
        maxDoseBlocks = ceil(nDoses/9);
        if ud.doseBlock == maxDoseBlocks
            return;
        end
        ud.doseBlock = ud.doseBlock + 1;
        set(h, 'userdata', ud);
        doseManagementGui('REFRESH_THUMBS')

    case 'REFRESH_THUMBS'

        ud = get(h, 'userdata');
        if isfield(ud.handles,'thumbaxis')
            try %To handle double-clicks on next and previous buttons
                delete(ud.handles.thumbaxis);
            catch
                %return
            end
        end

        nDoses      = length(planC{indexS.dose});

        dosesToDisplayV = (ud.doseBlock-1)*9+1:min(nDoses,ud.doseBlock*9);
        nDoses = length(dosesToDisplayV);

        if ud.currentDose > max(dosesToDisplayV)
            ud.currentDose = max(dosesToDisplayV);
        elseif ud.currentDose < min(dosesToDisplayV)
            ud.currentDose = min(dosesToDisplayV);
        end

        %Downsample colormap, redraws much faster.
        cM = CERRColorMap(stateS.optS.doseColormap);
        n  = size(cM, 1);
        newSize = 32;
        interval = (n-1) / 32;
        b = interp1(1:n, cM(:,1), 1:interval:n);
        c = interp1(1:n, cM(:,2), 1:interval:n);
        d = interp1(1:n, cM(:,3), 1:interval:n);
        ud.cM = [b' c' d'];

        %Setup thumbnail pane, with NxN axes.
        dx = 1/x; %pixel width in x,
        dy = 1/y; %pixel width in y, for margins.
        %thumbRegion = [.52 .17 .44 .75];
        thumbRegion = [.52 .23 .44 .70];
        subPlotSize = max(1,ceil(sqrt(nDoses)));
        dh = thumbRegion(4)/subPlotSize;
        dw = thumbRegion(3)/subPlotSize;
        ud.handles.thumbaxis = [];
        for i =1:subPlotSize^2
            row = subPlotSize - ceil(i/subPlotSize) + 1;
            col = mod(i-1,subPlotSize)+1;
            ud.handles.thumbaxis(i) = axes('position', [thumbRegion(1) + dw*(col-1) thumbRegion(2) + dh*(row-1) dw-dx dh-dy], 'box', 'on');
            set(ud.handles.thumbaxis(i), 'ytick',[],'xtick',[], 'color', 'black');
            colormap(ud.handles.thumbaxis(i), ud.cM);
        end

        maxDose = [];
        for i=1:nDoses
            set(ud.handles.thumbaxis(i), 'ButtonDownFcn', ['doseManagementGui(''CHANGEDOSE'', ' num2str(dosesToDisplayV(i)) ');']);
            maxDose{dosesToDisplayV(i)} = num2str(drawThumb(ud.handles.thumbaxis(i), planC, dosesToDisplayV(i), h));
            ud.previewSlice(dosesToDisplayV(i)) = 1; %%%%%%%%%%%%%%%%%
        end

        ud.maxDoses = maxDose;

        try
            set(ud.handles.previewAxis,'nextPlot','add')
            set(ud.handles.thumbaxis,'nextPlot','add')
        end

        set(h,'userdata',ud)

        if ~isempty(maxDose)
            doseManagementGui('REFRESHFIELDS');
        end

    case 'PREVIEWBUTTONDOWN'
        %Button clicked in the preview window.
        ud = get(h, 'userdata');
        ud.previewDown = 1;
        set(h, 'WindowButtonMotionFcn', 'doseManagementGui(''PREVIEWMOTION'')');
        set(h, 'userdata', ud)

    case 'FIGUREBUTTONUP'
        %Mouse up, if in preview window disable motion fcn.
        ud = get(h, 'userdata');
        if ~isfield(ud, 'previewDown') || ud.previewDown == 1
            ud.previewDown = 0;
            set(h, 'WindowButtonMotionFcn', '');
            set(h, 'userdata', ud);
        end

    case 'PREVIEWMOTION'
        %Motion in the preview, with mouse down. Change preview slice.
        ud = get(h, 'userdata');
        cp = get(h, 'currentpoint');
        if isfield(ud, 'previewY')
            if ud.previewY > cp(2)
                ud.previewSlice(ud.currentDose) = ud.previewSlice(ud.currentDose)+1;%min(ud.previewSlice(ud.currentDose)+1, size(getDoseArray(ud.currentDose), 3));
                set(h, 'userdata', ud);
                doseManagementGui('refreshpreviewandfields');
            elseif ud.previewY < cp(2)
                ud.previewSlice(ud.currentDose) = ud.previewSlice(ud.currentDose)-1;%max(ud.previewSlice(ud.currentDose)-1,1);
                set(h, 'userdata', ud);
                doseManagementGui('refreshpreviewandfields');
            end
            ud = get(h, 'userdata');
            ud.previewY = cp(2);
        else
            ud.previewY = cp(2);
        end
        set(h, 'userdata', ud);

    case 'REFRESHPREVIEWANDFIELDS'
        %Refresh both the preview and the fields.  Not modular with RefreshFields
        %to ensure that only one call to getDoseArray is made.
        ud = get(h, 'userdata');
        try
            set(ud.handles.previewAxis,'nextPlot','replace')
            set(ud.handles.thumbaxis,'nextPlot','replace')
        end
        doseNum = ud.currentDose;
        bdf = get(ud.handles.previewAxis, 'buttondownfcn');
        [dA, isCompress, isRemote] = getDoseArray(doseNum, planC);
        minSlice = 1;
        maxSlice = size(dA, 3);
        previewSlice = clip(ud.previewSlice(doseNum), minSlice, maxSlice, 'limits');
        ud.previewSlice(doseNum) = previewSlice;
        set(ud.handles.previewSliceNum, 'string', num2str(previewSlice));
        maxDose = str2num(ud.maxDoses{doseNum});
        cLim =  [0 maxDose];
        if maxDose == 0;
            cLim = [0 1];
        end
        imagesc(dA(:,:,previewSlice), 'parent', ud.handles.previewAxis, 'hittest', 'off');
        set(ud.handles.previewAxis, 'buttondownfcn', bdf, 'box', 'on', 'ytick', [], 'xtick', [], 'CLim', cLim);
        set(h, 'userdata', ud);

        %Refresh fields as well.
        set(ud.handles.name, 'string', planC{indexS.dose}(doseNum).fractionGroupID);
        set(ud.handles.units, 'string', planC{indexS.dose}(doseNum).doseUnits);
        if isCompress
            set(ud.handles.compbutton, 'string', 'Decompress');
            set(ud.handles.storageMethod, 'string', 'Memory, Compressed');
        else
            set(ud.handles.compbutton, 'string', 'Compress');
        end

        if isRemote
            set(ud.handles.remotebutton, 'string', 'Use Memory');
            set(ud.handles.storageMethod, 'string', 'On Disk');
        else
            set(ud.handles.remotebutton, 'string', 'Use Disk');
        end

        if ~isCompress & ~isRemote
            set(ud.handles.storageMethod, 'string', 'Memory');
        end

        %Refresh borders.
        set(ud.handles.thumbaxis, 'xcolor', 'black', 'ycolor', 'black');
        set(ud.handles.thumbaxis(doseNum - (ud.doseBlock-1)*9), 'xcolor', 'yellow', 'ycolor', 'yellow');

        set(ud.handles.maxdose, 'string', ud.maxDoses{doseNum});
        doseSize = getByteSize(planC{indexS.dose}(doseNum));
        set(ud.handles.dosesize, 'string', [num2str(doseSize/(1024*1024), '%6.2f') 'MB']);

        set(ud.handles.previewAxis,'nextPlot','add')
        set(ud.handles.thumbaxis,'nextPlot','add')

    case 'CHANGEDOSE'
        %New dose has been clicked on.
        ud = get(h, 'userdata');
        oldDose = ud.currentDose;
        newDose = varargin{1};
        ud.currentDose = newDose;
        set(h, 'userdata', ud);
        doseManagementGui('refreshpreviewandfields');

    case 'NAMEFIELD'
        %Dose name has changed, update in planC.
        ud = get(h, 'userdata');
        doseNum = ud.currentDose;
        doseUID = planC{indexS.dose}(doseNum).doseUID;
        oldString = planC{indexS.dose}(doseNum).fractionGroupID;
        string = get(gcbo, 'string');

%       Check for duplicate or blank names.
        if (length(string) == 0)
            warndlg('The dose name cannot be blank. Please enter a name.', 'Warning','modal');
            return;
        end
        if ~strcmp(oldString, string)
            for i = 1:length(planC{indexS.dose})
                if ischar(planC{indexS.dose}(1,i).fractionGroupID)
                    idString = planC{indexS.dose}(1,i).fractionGroupID;
                else
                    idString = num2str(planC{indexS.dose}(1,i).fractionGroupID);
                end
                if strcmpi(idString, string) && ~strcmp(idString, oldString)
                    warndlg('That plan name is already in use. Please choose another name.', 'Warning','modal');
                    return;
                end
            end
            for i = 1:length(planC{indexS.DVH})
                if ischar(planC{indexS.DVH}(1,i).planIDOfOrigin)
                    idString = planC{indexS.DVH}(1,i).planIDOfOrigin;
                else
                    idString = num2str(planC{indexS.DVH}(1,i).planIDOfOrigin);
                end
                if strcmpi(idString, string) && ~strcmp(idString, oldString)
                    warndlg('That plan name is already in use. Please choose another name.', 'Warning','modal');
                    return;
                end
            end
        end

        planC{indexS.dose}(doseNum).fractionGroupID = string;
        statusString = ['Renamed dose number ' num2str(doseNum) ' from ''' oldString ''' to ''' string '''.'];
        doseManagementGui('status', statusString);

%       Write new fractionGroupID to beam geometry
        for i = 1:length(planC{indexS.beamGeometry})
            if ischar(planC{indexS.beamGeometry}(1,i).fractionGroupID)
                idString = planC{indexS.beamGeometry}(1,i).fractionGroupID;
            else
                idString = num2str(planC{indexS.beamGeometry}(1,i).fractionGroupID);
            end
            if strcmp(idString, oldString)
                planC{indexS.beamGeometry}(1,i).fractionGroupID = string;
            end
        end

%       Write new fractionGroupID to DVH as fractionIDOfOrigin
        for i = 1:length(planC{indexS.DVH})
            if ischar(planC{indexS.DVH}(1,i).fractionIDOfOrigin)
                idString = planC{indexS.DVH}(1,i).fractionIDOfOrigin;
            else
                idString = num2str(planC{indexS.DVH}(1,i).fractionIDOfOrigin);
            end
            if strcmp(idString, oldString)
                planC{indexS.DVH}(1,i).fractionIDOfOrigin = string;
            end

            try
                if strcmpi(planC{indexS.DVH}(1,i).assocDoseUID, doseUID) && ...
                        ~isempty(planC{indexS.DVH}(1,i).planIDOfOrigin)
                    planC{indexS.DVH}(1,i).planIDOfOrigin = string;
                end
            end
        end

    case 'UNITFIELD'
        %Dose units have changed, update in planC.
        ud = get(h, 'userdata');
        doseNum = ud.currentDose;
        oldString = planC{indexS.dose}(doseNum).doseUnits;
        string = get(gcbo, 'string');
        planC{indexS.dose}(doseNum).doseUnits = string;
        statusString = ['Renamed units of dose number ' num2str(doseNum) ' from ''' oldString ''' to ''' string '''.'];
        doseManagementGui('status', statusString);

    case 'STATUS'
        %Display passed string in status bar.
        statusString = varargin{1};
        ud = get(gcbf, 'userdata');
        h = ud.handles.status;
        set(h, 'string', statusString);
        drawnow

    case 'COMPRESS'
        %Compress/decompress selected dose.
        ud = get(h, 'userdata');
        doseNum = ud.currentDose;
        doseName = planC{indexS.dose}(doseNum).fractionGroupID;

        if ~isCompressed(planC{indexS.dose}(doseNum).doseArray)
            statusString = ['Compressing dose number ' num2str(doseNum) ', ''' doseName ''' please wait...'];
            doseManagementGui('status', statusString);
            planC{indexS.dose}(doseNum).doseArray = compress(getDoseArray(doseNum, planC));
            drawThumb(ud.handles.thumbaxis(doseNum), planC, doseNum, h);
            statusString = ['Compressed dose number ' num2str(doseNum)  ', ''' doseName '''.'];
            doseManagementGui('status', statusString);
            set(ud.handles.compbutton, 'string', 'Decompress');
        else
            statusString = ['Decompressing dose number ' num2str(doseNum) ', ''' doseName ''', please wait...'];
            doseManagementGui('status', statusString);
            %Use getDoseArray and not decompress to use the cached value.
            planC{indexS.dose}(doseNum).doseArray = getDoseArray(doseNum, planC);
            maxDose = drawThumb(ud.handles.thumbaxis(doseNum), planC, doseNum, h);
            %Update the maxDoses value, but be sure to get a fresh ud since
            %a user could have clicked during compression.
            ud = get(h, 'userdata');
            ud.maxDoses{doseNum} = num2str(maxDose);
            set(h, 'userdata', ud);
            statusString = ['Decompressed dose number ' num2str(doseNum)  ', ''' doseName '''.'];
            doseManagementGui('status', statusString);
            set(ud.handles.compbutton, 'string', 'Compress');
        end
        doseManagementGui('refreshpreviewandfields');

    case 'REMOTE'
        %Make/unmake selected dose remote.
        ud = get(h, 'userdata');
        doseNum = ud.currentDose;
        doseName = planC{indexS.dose}(doseNum).fractionGroupID;

        doseUID = planC{indexS.dose}(doseNum).doseUID;

        if isLocal(planC{indexS.dose}(doseNum).doseArray)
            statusString = ['Writing to disk dose number ' num2str(doseNum) ', ''' doseName ''' please wait...'];
            doseManagementGui('status', statusString);
            [fpath,fname] = fileparts(stateS.CERRFile);
            %planC{indexS.dose}(doseNum).doseArray = setRemoteVariable(getDoseArray(doseNum, planC), 'LOCAL',[fpath,'\',fname,'_store'],['doseArray_',doseName,'.mat']);
            planC{indexS.dose}(doseNum).doseArray = setRemoteVariable(getDoseArray(doseNum, planC), 'LOCAL',fullfile(fpath,[fname,'_store']),['doseArray_',doseUID,'.mat']);
            drawThumb(ud.handles.thumbaxis(doseNum), planC, doseNum, h);
            statusString = ['Wrote to disk dose number ' num2str(doseNum)  ', ''' doseName '''.'];
            doseManagementGui('status', statusString);
            set(ud.handles.remotebutton, 'string', 'Use Memory');
            uiwait(msgbox(['doseArray stored in folder ',fullfile(fpath,[fname,'_store']),'Note the Location','modal']))

        else
            statusString = ['Reading from disk dose number ' num2str(doseNum) ', ''' doseName ''', please wait...'];
            doseManagementGui('status', statusString);

            remotePath = planC{indexS.dose}(doseNum).doseArray.remotePath;
            filenam = planC{indexS.dose}(doseNum).doseArray.filename;
            %Use getDoseArray and not decompress to use the cached value.
            planC{indexS.dose}(doseNum).doseArray = getDoseArray(doseNum, planC);

            stateS.reqdRemoteFiles(strcmp(fullfile(remotePath,filenam),stateS.reqdRemoteFiles)) = [];
            delete(fullfile(remotePath,filenam))

            %remove remote storage directory if it is empty
            dirRemoteS = dir(remotePath);
            if ~any(~cellfun('isempty',strfind({dirRemoteS.name},'.mat')))
                rmdir(remotePath)
            end
            maxDose = drawThumb(ud.handles.thumbaxis(doseNum), planC, doseNum, h);
            statusString = ['Read from disk dose number ' num2str(doseNum)  ', ''' doseName '''.'];
            doseManagementGui('status', statusString);
            set(ud.handles.remotebutton, 'string', 'Use Disk');

            %Update the maxDoses value, but be sure to get a fresh ud since
            %a user could have clicked during remote writing to disk.
            ud = get(h, 'userdata');
            ud.maxDoses{doseNum} = num2str(maxDose);
            set(h, 'userdata', ud);
        end
        doseManagementGui('refreshpreviewandfields');

    case 'SAVE'
        %Open dialog to save dose array as .mat file.
        ud = get(h, 'userdata');
        doseNum = ud.currentDose;
        doseName = planC{indexS.dose}(doseNum).fractionGroupID;
        [filename, pathname] = uiputfile('*.mat', ['Save (uncompressed) dose array number ' num2str(doseNum) ' as:']);
        if filename==0
            return;
        end
        dose3D = getDoseArray(doseNum, planC);
        save(fullfile(pathname, filename), 'dose3D');
        statusString = ['Saved dose number ' num2str(doseNum) ', ''' doseName ''' to ' [filename '.mat'] '.'];
        doseManagementGui('status', statusString);

    case 'DELETE'
        %Delete selected dose.  If being displayed, verify deletion with user.
        ud = get(h, 'userdata');
        doseNum = ud.currentDose;
        doseName = planC{indexS.dose}(doseNum).fractionGroupID;

        refreshViewer = 0;
        axesV = checkDisplayedDoses(doseNum);
        if ~isempty(axesV)
            choice = questdlg('One or more CERR axes are currently displaying this dose.  If you delete it, these axes will be set to display no dose.  Proceed?', 'Continue?', 'Continue', 'Abort', 'Continue');
            if strcmpi(choice, 'Abort')
                statusString = ['Delete aborted.'];
                doseManagementGui('status', statusString);
                return;
            else
                %Set the axes dose value to null.
                setAxisDoseToNull(axesV);
            end
        end

        if stateS.doseSet == doseNum
            stateS.doseSet = [];
        end
        del = questdlg(['Are you sure you want to delete dose number' num2str(doseNum) ', ''' doseName '''?'], 'Continue?', 'Continue', 'Abort', 'Continue');
        if strcmpi(del, 'Continue')
            statusString = ['Deleted dose number ' num2str(doseNum) ', ''' doseName '''.'];
            doseManagementGui('status', statusString);
            planC{indexS.dose}(doseNum) = [];
            %If dose below displayed dose deleted, its number has changed.
            if doseNum < stateS.doseSet
                stateS.doseSet = stateS.doseSet - 1;
            end
            refreshViewer = 1;
            updateAxesForDeletedDose(doseNum);
            doseManagementGui;
        else
        end

        % Delete colorbar children if plan contains no dose
        if length(planC{indexS.dose}) == 0
            delete(get(stateS.handle.doseColorbar.trans, 'children'));
        end

        %Refresh CERR axes in the case of changed dose sets.
        if refreshViewer
            sliceCallBack('refresh');
        end

    case 'DOSESCALE'
        ud          = get(h, 'userdata');
        doseNum     = ud.currentDose;
        doseName    = planC{indexS.dose}(doseNum).fractionGroupID;
        newMaxDose  = str2num(get(gcbo, 'string'));

        %Get the dose array and its compression state.
        [dA, isCompress, isRemote] = getDoseArray(doseNum, planC);

        maxDose = max(dA(:));
        if maxDose == 0
            set(gcbo, 'string', '0')
            statusString = ['Unable to rescale all zero dose distribution.'];
            doseManagementGui('status', statusString);
            return;
        end

        doseManagementGui('status', 'Rescaling dose, please wait...');

        if isCompress | isRemote
            statusString = ['Rescaling compressed or remote dose can take a moment, please wait...'];
            doseManagementGui('status', statusString);
        end

        %Perform the rescale.
        dA = dA *  (1/maxDose * newMaxDose);

        oldDoseUID = planC{indexS.dose}(doseNum).doseUID;

        % set new dose UID
        planC{indexS.dose}(doseNum).doseUID = createUID('DOSE');

        % update doseArray in dose struct.
        planC = setDoseArray(doseNum, dA, planC);

        % remove old remote file if it exists
        if ~isLocal(planC{indexS.dose}(doseNum).doseArray)
            if ~ismember(fullfile(planC{indexS.dose}(doseNum).doseArray.remotePath,['doseArray_',oldDoseUID,'.mat']),stateS.reqdRemoteFiles)
                delete(fullfile(planC{indexS.dose}(doseNum).doseArray.remotePath,['doseArray_',oldDoseUID,'.mat']))
            end
        end

        statusString = ['Rescaled dose number ' num2str(doseNum) ', ''' doseName ''' from [0 ' num2str(maxDose) '] Gy to [0 ' num2str(newMaxDose) '] Gy.'];
        ud.maxDoses{doseNum} = num2str(newMaxDose);
        set(h, 'userdata', ud);

        % update dose UID's for CERR axes
        for i = 1:length(stateS.handle.CERRAxis)
            if getAxisInfo(stateS.handle.CERRAxis(i),'doseSets') == doseNum
                setAxisInfo(stateS.handle.CERRAxis(i),'doseSets',doseNum)
            end
        end

        stateS.doseSetChanged = 1;
        stateS.CTDisplayChanged = 1;
        CERRRefresh
        doseManagementGui('status', statusString);

    case 'SCALEBYFACTOR'
        ud          = get(h, 'userdata');
        doseNum     = ud.currentDose;
        doseName    = planC{indexS.dose}(doseNum).fractionGroupID;
        scale       = str2num(get(gcbo, 'string'));
        if isempty(scale) | ~isnumeric(scale)
            return
        end

        %Get the dose array and its compression state.
        [dA, isCompress, isRemote] = getDoseArray(doseNum, planC);

        doseManagementGui('status', 'Scaling dose, please wait...');

        if isCompress | isRemote
            statusString = ['Scaling compressed or remote dose can take a moment, please wait...'];
            doseManagementGui('status', statusString);
        end

        %Perform the rescale.
        dA = dA *  scale;
        set(ud.handles.maxdose,'string',num2str(max(dA(:))))

        statusString = ['Rescaled dose number ' num2str(doseNum) ', ''' doseName ''' by a factor of ' num2str(scale)];

        oldDoseUID = planC{indexS.dose}(doseNum).doseUID;

        % set new dose UID
        planC{indexS.dose}(doseNum).doseUID = createUID('DOSE');

        % update doseArray in dose struct.
        planC = setDoseArray(doseNum, dA, planC);

        % remove old remote file if it exists
        if ~isLocal(planC{indexS.dose}(doseNum).doseArray)
            if ~ismember(fullfile(planC{indexS.dose}(doseNum).doseArray.remotePath,['doseArray_',oldDoseUID,'.mat']),stateS.reqdRemoteFiles)
                delete(fullfile(planC{indexS.dose}(doseNum).doseArray.remotePath,['doseArray_',oldDoseUID,'.mat']))
            end
        end

        ud.maxDoses{doseNum} = num2str(str2num(ud.maxDoses{doseNum})*scale);
        set(h, 'userdata', ud);

        % update dose UID's for CERR axes
        for i = 1:length(stateS.handle.CERRAxis)
            if  getAxisInfo(stateS.handle.CERRAxis(i),'doseSets') == doseNum
                setAxisInfo(stateS.handle.CERRAxis(i),'doseSets',doseNum)
            end
        end

        stateS.doseSetChanged = 1;
        stateS.CTDisplayChanged = 1;
        CERRRefresh
        doseManagementGui('status', statusString);

    case 'QUIT'
        close;
end

return;

function nBytes = getByteSize(data)
%"getByteSize"
%Returns the number of bytes in the passed data
infoStruct = whos('data');
nBytes = infoStruct.bytes;
return;

function maxDose = drawThumb(hAxis, planC, index, hFigure)
%"drawThumb"
%In passed dose array, find slice with highest dose and draw in hAxis.
%Also denote the index in the corner.  If compressed show compressed.
set(hFigure, 'CurrentAxes', hAxis);
toDelete = get(hAxis, 'children');
delete(toDelete);

%Get the dose array and its compression state.
[dA, isCompress, isRemote] = getDoseArray(index, planC);

bdf = get(hAxis, 'buttondownfcn');

maxDose = arrayMax(dA);
maxLoc = find(dA == maxDose);
[r,c,s] = ind2sub(size(dA), maxLoc(1));
thumbImage = dA(:,:,s(1));
imagesc(thumbImage, 'hittest', 'off', 'parent', hAxis);
set(hAxis, 'ytick',[],'xtick',[]);
if isCompress && isRemote
    text(.1, .1, 'Compressed', 'units', 'normalized', 'fontsize', 8, 'color', 'black', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
    text(.1, .2, 'Remote', 'units', 'normalized', 'fontsize', 8, 'color', 'black', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
elseif isRemote
    text(.1, .1, 'Remote', 'units', 'normalized', 'fontsize', 8, 'color', 'black', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
elseif isCompress
    text(.1, .1, 'Compressed', 'units', 'normalized', 'fontsize', 8, 'color', 'black', 'horizontalAlignment', 'left', 'hittest', 'off', 'parent', hAxis);
end

xLim = get(hAxis, 'xlim');
yLim = get(hAxis, 'ylim');
x = (xLim(2) - xLim(1)) * .87 + xLim(1);
y = (yLim(2) - yLim(1)) * .15 + yLim(1);
text(x, y, num2str(index), 'fontsize', 8, 'color', 'white', 'hittest', 'off', 'parent', hAxis);
set(hAxis, 'buttondownfcn', bdf);
axis ij;
drawnow;
return;

function [nAxesV] = checkDisplayedDoses(doseNum)
%"checkDisplayedDoses"
%   Check to see if any CERR axes are displaying the requested dose num.
%   Returns the numbers of any axes that are.

global planC stateS

nAxesV = [];

%Iterate over axes.
for i=1:length(stateS.handle.CERRAxis)
    %Get axis info for this axis.
    %aI = get(stateS.handle.CERRAxis(i), 'userdata');
    aI = stateS.handle.aI(i);
    if ismember(doseNum, aI.doseSets)
        nAxesV = union(nAxesV, i);
    end
end
return;

function setAxisDoseToNull(nAxesV)
%"setAxisDoseToNull"
%   Sets the passed axis number's userdata doseUID fields to null if they
%   were using 'manual'.

global planC stateS

%Iterate over axes.
for i=1:length(nAxesV)
    %Get axis info for this axis.
    %aI = get(stateS.handle.CERRAxis(nAxesV(i)), 'userdata');
    aI = stateS.handle.aI(nAxesV(i));
    if strcmpi(aI.doseSelectMode, 'manual');
        aI.doseSets = [];
        stateS.handle.aI(nAxesV(i)) = aI;
    end
end
return;

function updateAxesForDeletedDose(delIndex)
%"updateAxesForDeletedDose"
%   Shifts all doseNums being displayed in CERR axes to account for deleted
%   dose distributions, given the number of the deleted dose.

global planC stateS
indexS = planC{end};

nAxesV = [];

%Iterate over axes.
for i=1:length(stateS.handle.CERRAxis)
    %Get axis info for this axis.
    %aI = get(stateS.handle.CERRAxis(i), 'userdata');
    aI = stateS.handle.aI(i);

    % DK
    aI.doseSets = aI.doseSets(aI.doseSets > delIndex) - 1;
%     doseSet = getAssociatedDose(aI.doseSets);
%     doseSet = doseSet(doseSet > delIndex)- 1;
%     if doseSet~=0
%         aI.doseSets = planC{indexS.dose}(doseSet).doseSets;
%     else
%         aI.doseSets = [];
%     end
    % DK

    stateS.handle.aI(i) = aI;
end
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function f = cos(f)
%COS   Cosine of a CHEBFUN3T object.
%
%   COS(F) returns the cosine of F.

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

% Check for empty:
if ( isempty(f) )
    return
end

op = @(x,y,z) cos(feval(f, x, y, z));  % Resample.
f = chebfun3t(op, f.domain);           % Call constructor.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [AP,TP,FN,FP,prec,rec,mdAll,levh]=averageprecision(coord_est,coord_true,radiustrue,intens)
% Evaluates average precision for set of true and estimated points. See S407_report for details.
% References: Everingham, M., Gool, L., Williams, C.K.I., Winn, J. & Zisserman, A. The Pascal Visual Object Classes (VOC) Challenge. International Journal of Computer Vision 88, 303-338 (2009).
%
% [AP,TP,FN,FP,prec,rec,mdAll]=averageprecision(coord_est,coord_true,radiustrue,intens)
% input:  coord_est - Nx2 matrix of estimated coordinates
%         coord_true- Nx2 matrix of true coordinates
%         radiustrue - limit on distance between true and estimated coordinates. For smaller distance the estimated point is considered as a true positive (can be set to sigma/2, where sigma is the std of the gaussian approximation of hte PSF.)
%         intens - Nx1 vector of intensity of the estimated source
%
% output: AP - average precision
%         TP - true positives
%         FN - false negatives
%         FP - false poitives
%         prec - precision
%         rec - recall
%         mdAll - localisation precision at the lowest confidence level

% levn = 100; % number of confidence levels for AP (see S407_report.pdf)

levh=sort(unique(sqrt(double(intens))));
levn=length(levh); % number of different intensities

nT=size(coord_true,1);
nE=size(coord_est,1);
minh=min(levh);
maxh=max(levh);
% levh=minh+(maxh-minh)/(levn-1)*[0:levn-1];

TP=zeros(levn,1);
FP=nE*ones(levn,1);
FN=nT*ones(levn,1);
T = coord_true;
mdAll=mindistsep(T,coord_est);
for ii=1:levn
    index = sqrt(intens)>=levh(ii); % above the limit brightness
    E = coord_est(index,:);
    [md,mT,mE]=mindistsep(T,E); % mT and mE are indeces of points connected by a distance md.
    TP(ii) = sum(md<=radiustrue); %true positives
    setDiffE=setdiff(1:size(E,1),mE); % These estimated points have not been assighned to any true point.
    setDiffT=setdiff(1:size(T,1),mT); % These true points have not been assighned to any estimated point.
    FP(ii) = numel(setDiffE)+sum(md>radiustrue); %False positives: points futher from true then md, and estimated pints without any assigned true point
    FN(ii) = numel(setDiffT)+sum(md>radiustrue); %False negatives: points futher from true then md, and true poitions without any assigned esitmated point
end

prec = TP./(TP+FP);
rec = TP./(TP+FN);

precinterpol = zeros(11,1);
for ii=0:10;
    recind=rec>=(.1*ii);
    maxprec = max(prec(recind));
    if ~isempty(maxprec)
        precinterpol(ii+1)=maxprec;
    end
end

AP=1/11*sum(precinterpol);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function F = diff(F, dim, n)
%DIFF   Componentwise derivative of a SPHEREFUNV.
%   DIFF(F) is the tangetial derivative of each component of F in
%   x-direction.
%
%   DIFF(F, DIM) is the first tangential derivative of F along the
%   dimension DIM.
%     DIM = 1 (default) is the derivative in the x-direction.
%     DIM = 2 is the derivative in the y-direction.
%     DIM = 3 is the derivative in the z-direction.
%
%   DIFF(F, DIM, N) is the Nth tangential derivative of each component of F
%   along the dimension specified.
%
%   See also SPHEREFUNV/CURL and SPHEREFUNV/DIV.

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

% Empty check:
if ( isempty(F) )
    return
end

% Defaults:
if ( nargin == 1 || isempty(dim) )
    dim = 1;
end
if ( nargin < 3 )
    n = 1;
end

% Diff each component.
for j = 1:F.nComponents
    F.components{j} = diff(F.components{j}, dim, n);
end

end


text:
% NIFTI Object
%
%   create     - Create a NIFTI-1 file
%   disp       - Disp a NIFTI-1 object
%   display    - Display a NIFTI-1 object
%   fieldnames - Fieldnames of a NIFTI-1 object
%   nifti      - Create a NIFTI-1 object
%   subsasgn   - Subscript assignment
%   subsref    - Subscript referencing
%
%   other operations are unlikely to work.
%
% Example usage.
%
% % Example of creating a simulated .nii file.
% dat            = file_array;
% dat.fname = 'junk.nii';
% dat.dim     = [64 64 32];
% dat.dtype  = 'FLOAT64-BE';
% dat.offset  = ceil(348/8)*8;
%
% % alternatively:
% % dat = file_array( 'junk.nii',dim,dtype,off,scale,inter)
%
% disp(dat)
%
% % Create an empty NIFTI structure
% N = nifti;
%
% fieldnames(N) % Dump fieldnames
%
% % Creating all the NIFTI header stuff
% N.dat = dat;
% N.mat = [2 0 0 -110 ; 0 2 0 -110; 0 0 -2 92; 0 0 0 1];
% N.mat_intent = 'xxx'; % dump possibilities
% N.mat_intent = 'Scanner';
% N.mat0 = N.mat;
% N.mat0_intent = 'Aligned';
%
% N.diminfo.slice = 3;
% N.diminfo.phase = 2;
% N.diminfo.frequency = 2;
% N.diminfo.slice_time.code='xxx'; % dump possibilities
% N.diminfo.slice_time.code = 'sequential_increasing';
% N.diminfo.slice_time.start = 1;
% N.diminfo.slice_time.end = 32;
% N.diminfo.slice_time.duration = 3/32;
%
% N.intent.code='xxx' ; % dump possibilities
% N.intent.code='FTEST'; % or N.intent.code=4;
% N.intent.param = [4 8];
%
% N.timing.toffset = 28800;
% N.timing.tspace=3;
% N.descrip = 'This is a NIFTI-1 file';
% N.aux_file='aux-file-name.txt';
% N.cal = [0 1];
%
% create(N); % Writes hdr info
%
% dat(:,:,:)=0; % Write out the data as all zeros
%
% [i,j,k] = ndgrid(1:64,1:64,1:32);
% dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 30^2))=1; % Write some ones in the file
% dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 15^2))=2;
%
%
% % displaying a slice
% imagesc(dat(:,:,12));colorbar
%
% % get a handle to 'junk.nii';
% M=nifti('junk.nii');
%
% imagesc(M.dat(:,:,12));
%
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: Contents.m 2696 2009-02-05 20:29:48Z guillaume $

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% function spm_induced_demo
% Demo routine for induced responses
%==========================================================================
%
% This demonstration illustrates the generative or forward model used for
% time frequency responses - in other words, a biophysically plausible
% dynamic causal model for induced responses. The basic idea is to
% integrate a neural mass model to obtain expected hidden neuronal
% states produced by an unknown (parameterised) exogenous input. The states
% are used as the expansion point for a linear perturbation analysis of
% the frequency response properties that are local in peristimulus time.
% The ensuing spectra (induced complex cross spectra) are then convolved
% with a wavelet window to generate predictions of a conventional time
% frequency (wavelet) transform. Crucially, these predictions are complex
% and can be used to characterise delays - in terms of cross covariance
% functions. Nonlinearities in the neural mass model mean that the spectral
% responses caused by random neuronal fluctuations are state dependent and
% therefore change with the expected hidden states over peristimulus time.
%
% This routine first creates a simple - two source - generative model using
% a canonical microcircuit architecture and convolution based dynamics. It
% then produces predictions of induced responses to a short and sustained
% input to the first source - as measured by two local field potential
% recordings at each source. Exactly the same model is then integrated in
% time,  using (serially correlated) random fluctuations to drive each
% source (in addition to the exogenous input). This is repeated over 16
% trials and the simulated responses are characterised in terms of a
% wavelet transform - to produce complex cross spectral  data features.
% These are shown graphically with their analytic predictions from the
% generative model.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_induced_demo.m 7679 2019-10-24 15:54:07Z spm $

% Model specification
%==========================================================================
rng('default')
clear spm_fp_cmc_tfm

% number of regions in coupled map lattice
%--------------------------------------------------------------------------
Nc    = 2;                                       % number of channels
Ns    = 2;                                       % number of sources
options.spatial  = 'LFP';
options.model    = 'TFM';
options.analysis = 'TFM';
M.dipfit.model = options.model;
M.dipfit.type  = options.spatial;
M.dipfit.Nc    = Nc;
M.dipfit.Ns    = Ns;

% extrinsic connections (forward an backward)
%--------------------------------------------------------------------------
A{1} = [0 0; 1 0];
A{2} = [0 1; 0 0];
A{3} = [0 0; 0 0];
C    = [1 0;0 1];

% get priors
%--------------------------------------------------------------------------
pE    = spm_dcm_neural_priors(A,{},C,options.model);
pE    = spm_L_priors(M.dipfit,pE);
pE    = spm_ssr_priors(pE);
x     = spm_dcm_x_neural(pE,options.model);

% extrinsic connections
%--------------------------------------------------------------------------
pE.A{1}   = pE.A{1} + 1;
pE.A{2}   = pE.A{2} + 1;
pE.A{3}   = pE.A{3} - 1;
pE.A{4}   = pE.A{4} - 1;

% supress noise
%--------------------------------------------------------------------------
pE.a(1,:) = -2;
pE.b(1,:) = -4;
pE.c(1,:) = -4;

% orders and model
%==========================================================================
nx    = length(spm_vec(x));
nu    = size(pE.C,2);

% create forward model
%--------------------------------------------------------------------------
M.g   = 'spm_gx_erp';
M.f   = 'spm_fx_cmc_tfm';
M.h   = 'spm_fp_cmc_tfm';
M.x   = x;
M.n   = nx;
M.pE  = pE;
M.m   = nu;
M.l   = Nc;
M.Hz  = 2:64;
M.Rft = 4;
M.ds  = 4;

% solve for steady state
%--------------------------------------------------------------------------
M.x   = spm_dcm_neural_x(pE,M);

% Integrate system to see response (time-frequency)
%==========================================================================

% invoke exogenous inputs
%--------------------------------------------------------------------------
N     = 128;
U.dt  = 6/1024;
pst   = ((1:N)' - N/4)*U.dt;
U.u   = sparse(N,M.m);

% exogenous input - a sustained input of about 128 seconds
%--------------------------------------------------------------------------
M.ons    = 64;
U.u(:,1) = spm_erp_u(pst,pE,M);

% integrate generative model to simulate a time frequency response
%--------------------------------------------------------------------------
[csd,erp,CSD,mtf,w,pst,x,dP] = spm_csd_int(pE,M,U);

% plot expected responses
%==========================================================================
spm_figure('GetWin','Simulated time-frequency responses');
pst = 1000*pst;

subplot(4,2,1)
plot(pst,U.u)
xlabel('peristimulus time (ms)')
title('Exogenous input','FontSize',16)
spm_axis tight

% LFP - expectation
%--------------------------------------------------------------------------
subplot(4,2,2)
j   = find(kron(sparse(1,[3 1 5],1,1,8),ones(Ns,1)));
plot(pst,x(j,:))
xlabel('peristimulus time (ms)')
title('Hidden neuronal states','FontSize',16)
spm_axis tight

subplot(4,2,3)
plot(pst,erp{1})
xlabel('peristimulus time (ms)')
title('Evoked response','FontSize',16)
spm_axis tight

% LFP - expectation
%--------------------------------------------------------------------------
subplot(4,2,4)
plot(pst,dP{1})
xlabel('peristimulus time (ms)')
title('intrinsic connectivity','FontSize',16)
spm_axis tight

% expected time frequency response (coherence and cross-covariance)
%--------------------------------------------------------------------------
spm_dcm_tfm_image(CSD{1},pst,w,1)

% expected time frequency response
%--------------------------------------------------------------------------
spm_figure('GetWin','transfer functions');

spm_dcm_tfm_transfer(mtf{1},pst,w)

% simulated responses
%==========================================================================
spm_figure('GetWin','Predicted responses');

% time-frequency
%--------------------------------------------------------------------------
xY.erp = erp;
xY.csd = csd;
spm_dcm_tfm_response(xY,pst,w)

% Integrate system to simulate responses
%==========================================================================
spm_figure('GetWin','Simulated trials');

% get spectral density of neuronal fluctuations
%--------------------------------------------------------------------------
Hz    = 2:128;
Gu    = spm_csd_mtf_gu(pE,Hz);

% simulate Nt trials
%--------------------------------------------------------------------------
M.ds  = 0;                                       % switch to ERP generation
Nt    = 8;
V     = U;
for j = 1:Nt

    fprintf('\nsimulating trial %i',j)
    y    = 32;
    while max(abs(spm_vec(y))) > 1
        V.u  = U.u + spm_rand_power_law(Gu,Hz,U.dt,N)*2;
        y    = spm_csd_int(pE,M,V);
    end
    D(:,:,j) = full(real(y{1}));

    % LFP - random fluctuations
    %----------------------------------------------------------------------
    subplot(2,2,1)
    plot(pst,D(:,:,j)), hold on
    xlabel('time (s)')
    title('simulated response and ERP','FontSize',16)
    spm_axis tight, axis square; drawnow

end
hold off

% LFP - expectations
%--------------------------------------------------------------------------
subplot(2,2,2)
plot(pst*1000,erp{1},'-.'), hold on
plot(pst*1000,mean(D,3)), hold off
xlabel('time (s)')
title('LFP response - expectation','FontSize',16)
spm_axis tight, axis square

% Time frequency response
%--------------------------------------------------------------------------
Nf    = length(M.Hz);
P     = zeros(N,Nf,Nc,Nc);
Q     = zeros(N,Nf,Nc,Nc);
E     = mean(D,3);
for k = 1:Nt

    d     = full(double(D(:,:,k)));
    G     = spm_morlet(d - E,w*U.dt,M.Rft);
    for i = 1:Nc
        for j = 1:Nc
            P(:,:,i,j) = (G(:,:,i).*conj(G(:,:,j)));
        end
    end
    Q = Q + P;
end

% scale induced responsest sample variance
%--------------------------------------------------------------------------
Vm    = mean(mean((var(D))));
Vs    = mean(diag(squeeze(mean(squeeze(mean(Q))))));
Q     = Vm*Q/Vs;

% time-frequency
%--------------------------------------------------------------------------
spm_dcm_tfm_image(Q,pst,w,1)

% simulated responses
%==========================================================================
spm_figure('GetWin','Predicted responses - simulated');

% time-frequency
%--------------------------------------------------------------------------
xY.erp{1} = E;
xY.csd{1} = Q;
spm_dcm_tfm_response(xY,pst,w)

% compare expected and simulated responses
%==========================================================================
spm_figure('GetWin','Expected and simulated responses');

spm_dcm_tfm_image(csd{1},pst,w,0)
spm_dcm_tfm_image(Q,pst,w,1)

author: spm


text:
function Qbar = computeQbar(C, Q, S, ny)
    Qbar_tmp = cell(ny, ny);

    for i=1:ny
       for j=1:ny
          if i==j
              if i == ny
                  Qbar_tmp{ny,ny} = C'*S*C;
              else
                  Qbar_tmp{i,j} = C'*Q*C;
              end
          else
              Qbar_tmp{i,j} = zeros(size(C'*Q*C));
          end
       end
    end

    Qbar = cell2mat(Qbar_tmp);
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function patch = img_patch(img, bb, randomize,p_par)

if nargin == 4 && randomize > 0

    rand('state',randomize);
    randn('state',randomize);

    NOISE = p_par.noise;
    ANGLE = p_par.angle;
    SCALE = p_par.scale;
    SHIFT = p_par.shift;

    cp  = bb_center(bb)-1;
    Sh1 = [1 0 -cp(1); 0 1 -cp(2); 0 0 1];

    sca = 1-SCALE*(rand-0.5);
    Sca = diag([sca sca 1]);

    ang = 2*pi/360*ANGLE*(rand-0.5);
    ca = cos(ang);
    sa = sin(ang);
    Ang = [ca, -sa; sa, ca];
    Ang(end+1,end+1) = 1;

    shR  = SHIFT*bb_height(bb)*(rand-0.5);
    shC  = SHIFT*bb_width(bb)*(rand-0.5);
    Sh2 = [1 0 shC; 0 1 shR; 0 0 1];

    bbW = bb_width(bb)-1;
    bbH = bb_height(bb)-1;
    box = [-bbW/2 bbW/2 -bbH/2 bbH/2];

    H     = Sh2*Ang*Sca*Sh1;
    bbsize = bb_size(bb);
    patch = uint8(warp(img,inv(H),box) + NOISE*randn(bbsize(1),bbsize(2)));

else

    % All coordinates are integers
    if isempty(find((round(bb)-bb) ~= 0, 1)) == 1
        L = max([1 bb(1)]);
        T = max([1 bb(2)]);
        R = min([size(img,2) bb(3)]);
        B = min([size(img,1) bb(4)]);
        patch = img(T:B,L:R);

        % Sub-pixel accuracy
    else

        cp = 0.5 * [bb(1)+bb(3); bb(2)+bb(4)]-1;
        H = [1 0 -cp(1); 0 1 -cp(2); 0 0 1];

        bbW = bb(3,:)-bb(1,:);
        bbH = bb(4,:)-bb(2,:);
        if bbW <= 0 || bbH <= 0
            patch = [];
            return;
        end
        box = [-bbW/2 bbW/2 -bbH/2 bbH/2];

        if size(img,3) == 3
            for i = 1:3
                P = warp(img(:,:,i),inv(H),box);
                patch(:,:,i) = uint8(P);
            end
        else
            patch = warp(img,inv(H),box);
            patch = uint8(patch);
        end

    end
end

author: yuxng


text:
function alpha = neper2db(alpha, y)
%NEPER2DB   Convert nepers to decibels.
%
% DESCRIPTION:
%       neper2db converts an attenuation coefficient in units of
%       Nepers/((rad/s)^y m)to units of dB/(MHz^y cm).
%
% USAGE:
%       alpha = neper2db(alpha)
%       alpha = neper2db(alpha, y)
%
% INPUTS:
%       alpha       - attenuation in Nepers/((rad/s)^y m)
%
% OPTIONAL INPUTS:
%       y           - power law exponent (default = 1)
%
% OUTPUTS:
%       alpha       - attenuation in dB/(MHz^y cm)
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 3rd December 2009
%       last update - 5th December 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also db2neper

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

if nargin == 1
    y = 1;
end

alpha = 20*log10(exp(1))*alpha*(2*pi/1e-6)^y/100;

repo: Aero-matlab


text:
function r = isequal_weak(a,b)
% Compare 2 values for equality under relaxed conditions.
% Equal = isequal_weak(A,B)
%
% In:
%   A : first value
%
%   B : second value
%
% Out:
%   Equal : whether the two values are weakly equal, ignoring:
%           * NaNs
%           * contents of Java objects
%           * variables accessible to anonymous functions that are not referenced in the code
%           * new-style classes that have the same struct() representation
%
% Examples:
%   % test two data sets for equivalence
%   isequal_weak(EEG1,EEG2)
%
% See also:
%   isequal, isequalwithequalnans
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-02-19

% Copyright (C) Christian Kothe, SCCN, 2011, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

if isequal(a,b)
    r = true;
else
    ca = class(a);
    if ~strcmp(ca,class(b))
        r = false;
        % unequal class
        return;
    else
        % same class; check in detail
        switch ca
            case {'double','char','logical'}
                % common non-recursive data structure
                r = isequalwithequalnans(a,b);
            case 'cell'
                % recurse into cells
                if ~isequal(size(a),size(b))
                    r = false;
                else
                    for k=1:numel(a)
                        if ~isequal_weak(a{k},b{k})
                            r = false;
                            return;
                        end
                    end
                    r = true;
                end
            case 'struct'
                % recurse into structs
                if ~isequal(size(a),size(b)) || ~isequal(fieldnames(a),fieldnames(b))
                    r = false;
                else
                    r = isequal_weak(struct2cell(a),struct2cell(b));
                end
            case 'function_handle'
                if isequal(a,b)
                    r = true;
                else
                    % functions are considered different by isequal...
                    sa = char(a);
                    if strcmp(sa,char(b)) && sa(1) == '@'
                        % but are both anonymous functions with identical code:
                        fa = functions(a);
                        fb = functions(b);
                        % equality is determined by referenced variables
                        if ~isempty(fa.workspace) && ~isempty(fb.workspace)
                            r = isequal_weak(fa.workspace{1},fb.workspace{1});
                        elseif isempty(fa.workspace) && isempty(fb.workspace)
                            r = true;
                        elseif (isempty(fa.workspace) && isempty(fieldnames(fb.workspace{1}))) || (isempty(fb.workspace) && isempty(fieldnames(fa.workspace{1})))
                            r = true;
                        else
                            r = false;
                        end
                    else
                        r = false;
                    end
                end
            otherwise
                % misc class
                if isnumeric(a)
                    % misc number format
                    r = isequalwithequalnans(a,b);
                elseif isobject(a)
                    % MATLAB objects with same class: compare as structs
                    r = isequal_weak(struct(a),struct(b));
                elseif isjava(a)
                    % Java objects with same class: ignore contents
                    r = true;
                else
                    try
                        % new-style classes with the same class?
                        r = isequal_weak(struct(a),struct(b));
                    catch
                        % some unknown data type
                        r = isequalwithequalnans(a,b);
                    end
                end
        end
    end
end

repo: REST


text:
function Timefreq = in_bst_timefreq(TimefreqFile, LoadFull, varargin)
% IN_TIMEFREQ_BST: Read a sources file in Brainstorm format.
%
% USAGE:  Timefreq = in_bst_timefreq(TimefreqFile, LoadFull=1, FieldsList=All)
%
% INPUT:
%    - TimefreqFile : Absolute or relative path to the file to read
%    - LoadFull     : If 1, and if file is a "Kernel only" file, load the data, compute Kernel*Data to get the full sources time series
%    - FieldsList   : List of fields to read from the file
% OUTPUT:
%    - Timefreq     : Brainstorm timefreq file structure

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

%% ===== PARSE INPUTS =====
% Get protocol folders
ProtocolInfo = bst_get('ProtocolInfo');
% Filename: Relative / absolute
if ~file_exist(TimefreqFile)
    TimefreqFile = bst_fullfile(ProtocolInfo.STUDIES, TimefreqFile);
    if ~file_exist(TimefreqFile)
        error('Timefreq file not found.');
    end
end
% Load full?
if (nargin < 2) || isempty(LoadFull)
    LoadFull = 1;
end
isKernel = ~isempty(strfind(TimefreqFile, '_KERNEL_'));
% Specific fields
if (nargin < 3)
    % Read all fields
    Timefreq = load(TimefreqFile);
    % Get all the fields
    FieldsToRead = fieldnames(Timefreq);
    FieldsToRead{end+1} = 'ZScore';
else
    % Get fields to read
    FieldsToRead = varargin;
    % Add some conditional fields
    if ismember('TF', FieldsToRead)
        FieldsToRead{end+1} = 'Measure';
        FieldsToRead{end+1} = 'ZScore';
        if LoadFull && isKernel
            FieldsToRead{end+1} = 'DataFile';
        end
    end
    % When reading Leff, make sure nAvg is read as well
    if ismember('Leff', FieldsToRead) && ~ismember('nAvg', FieldsToRead)
        FieldsToRead{end+1} = 'nAvg';
    end
    % Read each field only once
    FieldsToRead = unique(FieldsToRead);
    % Read specified files only
    warning off MATLAB:load:variableNotFound
    Timefreq = load(TimefreqFile, FieldsToRead{:});
    warning on MATLAB:load:variableNotFound
end

%% ===== FILL OTHER MISSING FIELDS =====
for i = 1:length(FieldsToRead)
    if ~isfield(Timefreq, FieldsToRead{i}) || isempty(Timefreq.(FieldsToRead{i}))
        switch(FieldsToRead{i})
            case 'Measure'
                Timefreq.(FieldsToRead{i}) = 'none';
            case 'nAvg'
                Timefreq.(FieldsToRead{i}) = 1;
            case 'Leff'
                if isfield(Timefreq, 'nAvg') && ~isempty(Timefreq.nAvg)
                    Timefreq.Leff = Timefreq.nAvg;
                else
                    Timefreq.Leff = 1;
                end
            case 'Method'
                if ~isempty(strfind(TimefreqFile, '_psd'))
                    Timefreq.(FieldsToRead{i}) = 'psd';
                elseif ~isempty(strfind(TimefreqFile, '_fft'))
                    Timefreq.(FieldsToRead{i}) = 'fft';
                elseif ~isempty(strfind(TimefreqFile, '_hilbert'))
                    Timefreq.(FieldsToRead{i}) = 'hilbert';
                elseif ~isempty(strfind(TimefreqFile, '_mtmconvol'))
                    Timefreq.(FieldsToRead{i}) = 'mtmconvol';
                else
                    Timefreq.(FieldsToRead{i}) = 'morlet';
                end
            otherwise
                Timefreq.(FieldsToRead{i}) = [];
        end
    end
end

%% ===== FILL EMPTY FIELDS =====
if isfield(Timefreq, 'Measure') && isempty(Timefreq.Measure)
    Timefreq.Measure = 'none';
elseif isfield(Timefreq, 'Measure')
    Timefreq.Measure = lower(Timefreq.Measure);
end

%% ===== FIX OLD TIME/FREQUENCY BANDS =====
% Old structure:  {name, fstart, fend}
% New structure:  {name, 'expression', 'stat')
if isfield(Timefreq, 'Freqs') && iscell(Timefreq.Freqs) && ~isempty(Timefreq.Freqs) && ~ischar(Timefreq.Freqs{1,2})
    for iBand = 1:size(Timefreq.Freqs,1)
        Timefreq.Freqs{iBand, 2} = [num2str(Timefreq.Freqs{iBand, 2}), ', ', num2str(Timefreq.Freqs{iBand, 3})];
        Timefreq.Freqs{iBand, 3} = 'mean';
    end
end
if isfield(Timefreq, 'TimeBands') && ~isempty(Timefreq.TimeBands) && ~ischar(Timefreq.TimeBands{1,2})
    for iBand = 1:size(Timefreq.TimeBands,1)
        Timefreq.TimeBands{iBand, 2} = [num2str(Timefreq.TimeBands{iBand, 2}), ', ', num2str(Timefreq.TimeBands{iBand, 3})];
        Timefreq.TimeBands{iBand, 3} = 'mean';
    end
end

%% ===== REBUILD FULL SOURCES =====
if ismember('TF', FieldsToRead) && LoadFull && isKernel && ismember(file_gettype(Timefreq.DataFile), {'link', 'results'})
    % Get imaging kernel
    ResultsMat = in_bst_results(Timefreq.DataFile, 0, 'ImagingKernel');
    % Initialize full matrix
    TFfull = zeros(size(ResultsMat.ImagingKernel,1), size(Timefreq.TF,2), size(Timefreq.TF,3));
    % Multiply the TF values
    for i = 1:size(Timefreq.TF,3)
        TFfull(:,:,i) = ResultsMat.ImagingKernel * Timefreq.TF(:,:,i);
    end
    Timefreq.TF = TFfull;
end

%% ===== APPLY DYNAMIC ZSCORE =====
% DEPRECATED
% Check for structure integrity
if ismember('ZScore', FieldsToRead) && ~isempty(Timefreq.ZScore) && (~isfield(Timefreq.ZScore, 'mean') || ~isfield(Timefreq.ZScore, 'std') || ~isfield(Timefreq.ZScore, 'abs') || ~isfield(Timefreq.ZScore, 'baseline') || isempty(Timefreq.ZScore.abs))
    Timefreq.ZScore = [];
end
% Apply to file values
if ismember('ZScore', FieldsToRead) && ismember('TF', FieldsToRead) && ~isempty(Timefreq.ZScore) && ~isempty(Timefreq.TF)
    % Error in the case of source-based sources
    if isKernel && ~LoadFull
        error('Unsupported operation.');
    end
    Timefreq.TF = process_zscore_dynamic('Compute', Timefreq.TF, Timefreq.ZScore);
    Timefreq = rmfield(Timefreq, 'ZScore');
end

% ===== FIX TRANSPOSED TIME VECTOR =====
if isfield(Timefreq, 'Time') && (size(Timefreq.Time,1) > 1)
    Timefreq.Time = Timefreq.Time';
end

% ===== FIX TRANSPOSED ROWNAMES =====
if isfield(Timefreq, 'RowNames') && (size(Timefreq.RowNames,1) > 1)
    Timefreq.RowNames = Timefreq.RowNames';
end
if isfield(Timefreq, 'RefRowNames') && (size(Timefreq.RefRowNames,1) > 1)
    Timefreq.RefRowNames = Timefreq.RefRowNames';
end

author: brainstorm-tools


text:
function [] = gsp_plot_filter(G, filters, param)
%GSP_PLOT_FILTER  Plot a system of filters
%   Usage:  [filter_data,test_sum]=gsp_plot_filter(G,filters);
%           [filter_data,test_sum]=gsp_plot_filter(G,filters,param);
%
%   Input parameters:
%       G       : Graph object (Or lmax)
%       filters : Cell array of filters (or single filter)
%       param   : Optional variable containing additional parameters
%   Output parameters:
%       none
%
%
%   'gsp_plot_filters(G, filter, param)' plots a system of graph spectral
%   filters.
%
%   Example:::
%
%         Nf = 4;
%         G = gsp_sensor(100);
%         G = gsp_estimate_lmax(G);
%         g = gsp_design_mexican_hat(G, Nf);
%         gsp_plot_filter(G, g);
%
%
%   Additional parameters
%   ---------------------
%
%   * *param.line_width* : Width of the filter plots (default 4).
%   * *param.npoints* : Number of point where the filters are evaluated
%     (default 1000).
%   * *param.x_tic* : Distance between x-tick labels.
%   * *param.y_tic* : Distance between y-tick labels (default 0.25).
%   * *param.minor_tick* : To show minor tick marks (default 1).
%   * *param.plot_eigenvalues* : To plot black X marks at all eigenvalues
%     of the graph (You need to compute the Fourier basis to use this
%     option). By default the eigenvalues are plot if they are contained in
%     the Graph.
%   * *param.lambda_highlights* : To plot red X marks at highlight
%     eigenvalues (default 0).
%   * *param.x_width* : Width of X marks for the eigenvalues (default 3).
%   * *param.x_size* : Size of X marks for the eigenvalues (default 8).
%   * *param.show_sum* : To plot an extra line showing the sum of the
%     squared magnitudes of the filters (default 1 if there is multiple
%     filters).
%   * *param.colors_rgb* : To specify the line colors.
%   * *param.cla* : Clear axis (default 1).
%   * *param.yrange* : To specify a range for the y axis.
%   * *param.verbose* : Verbosity level (1 display the warning - 0 no log)
%     (default 1).
%
%   Demos: gsp_demo
%

% Author : David I Shuman, Nathanael Perraudin
% Testing: test_filter

% Read input parameters
if nargin < 3
   param = struct;
end

if ~isstruct(G)
   G.lmax = G;
   param.plot_eigenvalues = 0;
end

if ~isfield(G,'lmax')
    G = gsp_estimate_lmax(G);
    warning(['GSP_FILTER_ANALYSIS: The variable lmax is not ',...
        'available. The function will compute it for you. ',...
        'However, if you apply many time this function, you ',...
        'should precompute it using the function: ',...
        'gsp_estimate_lmax']);
end

if ~isfield(param,'verbose'), param.verbose = 1; end
if ~isfield(param,'line_width'), param.line_width = 4; end
if ~isfield(param,'x_width'), param.x_width = 3; end
if ~isfield(param,'x_size'), param.x_size = 8; end
if ~isfield(param,'show_sum'), param.show_sum = length(filters)>1; end
if ~isfield(param,'y_tic'), param.y_tic = 0.25; end
if ~isfield(param,'minor_tick'), param.minor_tick = 1; end
if ~isfield(param,'npoints'), param.npoints = 1000; end
if ~isfield(param,'cla'), param.cla = 1; end
if ~isfield(param,'x_tic')
    param.x_tic=max(1,ceil(G.lmax/10));
end
if ~isfield(param,'plot_eigenvalues')
    param.plot_eigenvalues = isfield(G,'e');
end

if ~isfield(G,'lmax')
    if param.verbose
        fprintf('GSP_KERNEL_MEXICAN_HAT has to compute lmax \n')
    end
    G = gsp_estimate_lmax(G);
end

lambdas = linspace(0,G.lmax,param.npoints);

if param.cla
    cla;
end
hold on;
% apply the filter
fd = gsp_filter_evaluate(filters,lambdas);

% plot the filter
plot(lambdas,fd,'LineWidth',param.line_width);

% plot the eigenvalues
if param.plot_eigenvalues
    if isfield(G,'e')
        plot(G.e,zeros(G.N,1),'xk','LineWidth',...
            param.x_width,'MarkerSize',param.x_size);
    else
        if param.verbose
            warning('GSP_PLOT_FILTER: No eigenvalues found in the graph');
        end
    end
end

% plot hightlights eigenvalues
if isfield(param,'lambda_highlights')
    plot(param.lambda_highlights, ...
        zeros(length(param.lambda_highlights),1),...
        'xr','LineWidth',param.x_width,'MarkerSize',param.x_size);
end

% plot the sum
if param.show_sum
    test_sum=sum(fd.^2,2);
    plot(lambdas,test_sum,'k','LineWidth',param.line_width);
end

box on;
% X axis
xlim(full([0 G.lmax]));
% Y axis
if isfield(param,'yrange')
    yrange=param.yrange;
    ylim(yrange);
    set(gca,'YTick',yrange(1):y_tic:yrange(2));
end

% Add marks
if param.minor_tick
    set(gca,'XMinorTick','on','YMinorTick','on');

end

% Change the color
if isfield(param,'colors_rgb');
    set(gca, 'ColorOrder', param.colors_rgb);
end
set(gca,'XTick',0:param.x_tic:G.lmax);

hold off;
end

repo: gspbox


text:
% Genralized Sentiment Analysis
%

clc;
close all;
clear all;

tic;
%
 fid = fopen('data\english.stop');

 stopwords = textscan(fid, '%s');
 stopwords = stopwords{1,1};
 fclose(fid);

% params
% n:minimum appearances of a stem
n=30;

% reading the data file
% xls is easier to read than csv
[num,txt,raw] = xlsread('C:\MatlabNLP\examples\gsa\data\final104.xls');

% reading the description of each shoe
descriptions = raw(2:size(raw,1),2);
style_ratings = num(1:size(num,1),1);
comfort_ratings = num(1:size(num,1),4);
overal_ratings = num(1:size(num,1),5);

% only take m data points

m=50;
m = length(descriptions);
descriptions = descriptions(1:m);
style_ratings = style_ratings(1:m);
comfort_ratings = comfort_ratings(1:m);
overal_ratings = overal_ratings(1:m);
g = containers.Map();
toc;

for i = 1:size(descriptions,1)
    comment = descriptions{i};
    comment = SanitizeComment(comment);
    comment = lower(comment);
    r=regexp(comment,' ','split');
    for j =1:size(r,2)
        word = porterStemmer(cell2mat(r(j)));
        tfflag = ~isStopWord(word, stopwords);
        if isKey(g, word) && tfflag
            g(word) = g(word)+1;
        elseif tfflag
            g(word) = 1;
        end

    end

end
toc;
selectedheaders =containers.Map();
gkeys = keys(g);

for i=1:size(gkeys,2)
    if g(gkeys{i})>=n

        selectedheaders(gkeys{i})=1;

    end

end
headers = keys(selectedheaders);

outputMatrix = zeros(m,length(headers));
for i = 1:size(descriptions,1)
    comment = descriptions{i};
    comment = SanitizeComment(comment);
    comment = lower(comment);

    r=regexp(comment,' ','split');
    comment = [];
    for j =1:size(r,2)
        word = porterStemmer(cell2mat(r(j)));

        comment = [comment,' ',word];
    end
    outputMatrix(i,:) = term_count(comment, headers);

    if mod(i,300)==0
        a = sprintf('%d', i);
        disp(a)
    end

end

toc;
csvwrite('forWeka_featuresonly.csv', outputMatrix);

outputmatrix_all = [style_ratings,comfort_ratings,overal_ratings, outputMatrix];
outputmatrix_all = [ones(1,size(outputmatrix_all,2));outputmatrix_all];
csvwrite('forWeka_everything.csv', outputmatrix_all);

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [result] = disag_fhmm(evaluation_and_training_days, setup, fid)

    % load parameters of algorithm
    dataset = setup.dataset;
    household = setup.household;
    granularity = setup.granularity;
    filteringMethod = setup.filtering;
    filtLength = setup.filtLength;
    evaluation_days = evaluation_and_training_days{1};
    training_days = evaluation_and_training_days{2};
    num_appliances = setup.num_appliances;

    % prepare plug and smart meter data for training
    appliances = {};
    sum_of_plug_data = zeros(1, size(training_days,1) * 24 * 3600);
%     sum_of_plug_data = sum_of_plug_data(1:100);
    for i=1:num_appliances
        appliance = struct;
        eval(['appliance.name = setup.appliance' num2str(i) ';']);
        appliance.id = getApplianceID(appliance.name);
        fprintf('Reading data from appliance %s - id: %d\n', appliance.name, appliance.id);
        appliance.plug_data_training = read_plug_data(dataset, household, appliance.id, training_days, granularity);
%         appliance.plug_data_training = appliance.plug_data_training(1:100);
        sum_of_plug_data = sum_of_plug_data + appliance.plug_data_training;
        eval(['appliance.num_states = setup.num_states_appliance' num2str(i) ';']);

        appliances{i} = appliance;
    end

    % prepare "other" plug
    if setup.states_other > 0
        smart_meter_data_training = read_smartmeter_data(dataset, household, training_days, granularity, 'powerallphases');
%         smart_meter_data_training = smart_meter_data_training(1:100);
        other_data = smart_meter_data_training - sum_of_plug_data;
        other_appliance = struct;
        appliance.name = 'other';
        appliance.plug_data_training = other_data;
        appliance.num_states = setup.states_other;
        appliances{end+1} = appliance;
    end

    smart_meter_data_evaluation = read_smartmeter_data(dataset, household, evaluation_days, granularity, 'powerallphases');
%     smart_meter_data_evaluation = smart_meter_data_evaluation (1:100);

    % perform training
    pyObj = py.nilmtk.disaggregate.chris.fhmm.Fhmm();
    pyObj.train(appliances);

    fprintf('Done with training');

    ret = pyObj.disaggregate(smart_meter_data_evaluation);

    states = cell(ret){1};
    consumption = cell(ret){2};

    result = struct;
    inferred_consumption = zeros(num_appliances, length(smart_meter_data_evaluation));
    appliance_names = {};
    for i=1:num_appliances
        eval(['name = setup.appliance' num2str(i) ';']);
        appliance_names{end+1} = name;
        % eval(['app_result.states = cell2mat(cell(py.list(struct(states).' app_result.name ')));']);
        eval(['cons = cell2mat(cell(py.list(struct(consumption).' name ')));']);
        inferred_consumption(i,:) = cons;
    end
    result.appliance_names = appliance_names;
    result.consumption = inferred_consumption;
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function zipfilename = exporttozip(funcname,abspath,zipfilename)
%% EXPORTTOZIP - Creates a ZIP file containing dependencies of a function
%
% zipfilename = exporttozip(funcname);
% zipfilename = exporttozip(funcname,abspath);
% zipfilename = exporttozip(funcname,abspath,zipfilename);
%
% Finds and Zips all M-files which are required by the specified function(s)
% and which are NOT inside MATLAB's "toolbox" directory.
% Note: this does not include functions called by eval() statements
%
% funcname      is one or more function names, as a string or cell array of strings
% abspath       flag for absolute (rather than relative) paths
%               1 = use absolute paths, write the zipfile to the local directory
%               0 = use relative paths, write the zipfile to the commom rootpath
% zipfilename   the zipfile to export the functions to
%
% abspath is useful when you are using networked drives where a common
% rootpath may not exist or be useful
%
% If no zipfilename is specified, the ZIP file is created in the current
% working directory with name X.zip where X is the first function name in
% the supplied list
%
% The file names in the ZIP file are relative to the
% common root directory directory of all the required files.  If no
% common root directory can be established (e.g. if files are on
% different drives) an error is thrown (sometimes).
%

%% Version
% Copyright 2006-2010 The MathWorks, Inc.
% Original: Malcolm Wood    11 Apr 2006
% Modified: ?               15 Mar 2010     ? on Matlab FileExchange
%           Mark Morehead      Mar 2011     added abspath ability
%                                             relative paths do not work well
%                                             across multiple network drives
%                                           changed strcmp to strcmpi to
%                                             avoid problems with mixing cases
%                                             recommend by Mut Ante
%           Mark Morehead   25 Feb 2013     Added nargin<2 test
%                                           recommended by Ioannis Filippidis
%

%% Start

if ~iscell(funcname)
    funcname = {funcname};
end

if isempty(funcname)
    error('No function names specified');
end

if ~iscellstr(funcname)
    error('Function names must be strings');
end

if nargin>1 && ~isnumeric(abspath)
    error('abspath must be a numeric value');
end

if nargin < 2
    abspath = 0;
end

if nargin<3
    zipfilename = fullfile(pwd,[ funcname{1} '.zip' ]);
end

req = cell(size(funcname));
for i=1:numel(funcname)
    req{i} = mydepfun(funcname{i},1); % recursive
end
req = vertcat(req{:}); % cell arrays of full file names
req = unique(req);

% Use relative rootpath
if ~abspath
    % Find the common root directory
    d = i_root_directory(req);
    % Calculate relative paths for all required files.
    n = numel(d);
    for i=1:numel(req)
        % This is the bit that can't be vectorised
        req{i} = req{i}(n+1:end); % step over last character (which is the separator)
    end
    zip(zipfilename,req,d);
else
    zip(zipfilename,req);
end

fprintf(1,'Created %s \n',zipfilename);
fprintf(1,' with %d entries (exludes functions called by eval())\n',numel(req));

%%%%%%%%%%%%%%%%%%%%%
% Identifies the common root directory of all files in cell array "req"
function d = i_root_directory(req)

d = i_parent(req{1});
for i=1:numel(req)
    t = i_parent(req{i});
    if strncmpi(t,d,numel(d))
        % req{i} is in directory d.  Next file.
        continue;
    end
    % req{i} is not in directory d.  Up one directory.
    count = 1;
    while true
        % Remove trailing separator before calling fileparts.  Add it
        % again afterwards.
        tempd = i_parent(d(1:end-1));
        if strcmp(d,tempd)
            % fileparts didn't find us a higher directory
            error('Failed to find common root directory for %s and %s',req{1},req{i});
        end
        d = tempd;
        if strncmpi(t,d,numel(d))
            % req{i} is in directory d.  Next file.
            break;
        end
        % Safety measure for untested platform.
        count = count+1;
        if count>1000
            error('Bug in i_root_directory.');
        end
    end
end

%%%%%%%%%%%%%%%%%%%
function d = i_parent(d)
% Identifies the parent directory, including a trailing separator

% Include trailing separator in all comparisons so we don't assume that
% file C:\tempX\file.txt is in directory C:\temp
d = fileparts(d);
if d(end)~=filesep
    d = [d filesep];
end

author: Sable

repo: mcbench-benchmarks


text:
function display(b)
% bsarray/display: display of bsarray object
% usage: display(b)
%
% arguments:
%   b - bsarray object
%

% author: Nathan D. Cahill
% email: ndcahill@gmail.com
% date: 18 April 2008

n = b.tensorOrder;

switch n
    case 0
        str = sprintf([ ...
            'DataArray:     [ Empty ] double array\n',...
            'Coefficients:  [ Empty ] double array\n',...
            'Degree:          Empty\n',...
            'Centred:         Empty\n',...
            'Element Spacing: Empty\n',...
            'Lambda:          Empty\n']);
    case 1
        str = sprintf([ ...
            'DataArray:     [ %d x %d ] double array\n',...
            'Coefficients:  [ %d x %d ] double array\n',...
            'Degree:          %d\n',...
            'Centred:         %d\n',...
            'Element Spacing: %g\n',...
            'Lambda:          %g\n'],...
            b.dataSize(1),b.dataSize(2),...
            b.coeffsSize(1),b.coeffsSize(2),...
            b.degree(1),b.centred(1),b.elementSpacing(1),b.lambda(1));
    case 2
        str = sprintf([ ...
            'DataArray:       [ %d x %d ] double array\n',...
            'Coefficients:    [ %d x %d ] double array\n',...
            'Degree:          [ %d %d ]\n',...
            'Centred:         [ %d %d ]\n',...
            'Element Spacing: [ %g %g ]\n',...
            'Lambda:          [ %g %g ]\n'],...
            b.dataSize(1),b.dataSize(2),...
            b.coeffsSize(1),b.coeffsSize(2),...
            b.degree(1),b.degree(2),...
            b.centred(1),b.centred(2),...
            b.elementSpacing(1),b.elementSpacing(2),...
            b.lambda(1),b.lambda(2));
    case 3
        str = sprintf([ ...
            'DataArray:       [ %d x %d x %d ] double array\n',...
            'Coefficients:    [ %d x %d x %d ] double array\n',...
            'Degree:          [ %d %d %d ]\n',...
            'Centred:         [ %d %d %d ]\n',...
            'Element Spacing: [ %g %g %g ]\n',...
            'Lambda:          [ %g %g %g ]\n'],...
            b.dataSize(1),b.dataSize(2),b.dataSize(3),...
            b.coeffsSize(1),b.coeffsSize(2),b.coeffsSize(3),...
            b.degree(1),b.degree(2),b.degree(3),...
            b.centred(1),b.centred(2),b.centred(3),...
            b.elementSpacing(1),b.elementSpacing(2),b.elementSpacing(3),...
            b.lambda(1),b.lambda(2),b.lambda(3));
    otherwise
        str = sprintf([...
            'DataArray:       [ %d - D double array ]\n',...
            'Coefficients:    [ %d - D double array ]\n',...
            'Degree:          [ ',repmat('%d ',1,numel(b.degree)),']\n',...
            'Centred:         [ ',repmat('%d ',1,numel(b.centred)),']\n',...
            'Element Spacing: [ ',repmat('%g ',1,numel(b.elementSpacing)),']\n',...
            'Lambda:          [ ',repmat('%g ',1,numel(b.lambda)),']\n'],...
            n,n,b.degree,b.centred,b.elementSpacing,b.lambda);
end

if isequal(get(0,'FormatSpacing'),'compact')
    disp([inputname(1) ' =']);
    disp(str);
else
    disp(' ');
    disp([inputname(1) ' =']);
    disp(' ');
    disp(str);
end


text:
classdef InputFemFilePrinter < FilePrinter

    properties (Access = private)
        connectivities
        coordinates
        isElemInThisSet
        masterSlave
        corners
        nnode
        npnod
        nelem
        ndim
        resultsDir
        scale
        pdim
        ptype
        type
    end

    methods (Access = public)

        function obj = InputFemFilePrinter(inputData)
            obj.init(inputData);
            obj.createOutPutFolder();
        end

        function print(obj)
            obj.openFile();
            obj.printDataProblem();
            obj.printCoordinates();
            obj.printConnectivities();
            obj.printDirichletData();
            obj.printMasterSlave();
            obj.printMaterialSets();
            obj.closeFile();
        end

    end

    methods (Access = private)

        function init(obj,d)
            obj.connectivities  = d.connec;
            obj.coordinates     = d.coord;
            obj.isElemInThisSet = d.isElemInThisSet;
            obj.masterSlave     = d.masterSlave;
            obj.corners         = d.corners;
            obj.scale           = d.scale;
            obj.pdim            = d.pdim;
            obj.ptype           = d.ptype;
            obj.resultsDir      = d.resultsDir;
            obj.fileName        = d.fileName;
            obj.type            = d.type;
            obj.nnode           = size(d.connec,2);
            obj.nelem           = size(d.connec,1);
            obj.npnod           = size(d.coord,1);
            obj.ndim            = 2;
        end

        function createOutPutFolder(obj)
            dir = obj.resultsDir;
            if ~exist(dir,'dir')
                mkdir(dir)
                addpath(dir)
            end
        end

        function printDataProblem(obj)
            iD = obj.fileID;
            fprintf(iD,'Data_prb = { \n');
            fprintf(iD,['''',obj.type,'''','; \n']);
            fprintf(iD,'''SI''; \n');
            fprintf(iD,['''',obj.pdim,'''','; \n']);
            fprintf(iD,'''Plane_Stress''; \n');
            fprintf(iD,['''',obj.ptype,'''','; \n']);
            fprintf(iD,['''',obj.scale,'''','; \n']);
            fprintf(iD,'}; \n');
        end

        function printCoordinates(obj)
            iD      = obj.fileID;
            coor    = obj.coordinates;
            nodInEl = obj.npnod;
            nDim    = obj.ndim;
            fprintf(iD,'\n');
            fprintf(iD,'coord = [\n');
            printFormat = ['%6.0f ',repmat('%12.5d ',1,nDim),' 0 \n'];
            toPrint = [1:nodInEl;coor'];
            fprintf(iD,printFormat,toPrint);
            fprintf(iD,']; \n');
        end

        function printConnectivities(obj)
            iD     = obj.fileID;
            con    = obj.connectivities;
            nNodes = obj.nnode;
            fprintf(iD,'\n');
            fprintf(iD,'connec = [ \n');
            printFormat = [repmat('%6.0f ',1,nNodes+1),' \n'];
            fprintf(iD,printFormat,[1:obj.nelem;con']);
            fprintf(iD,']; \n');
        end

        function printMasterSlave(obj)
            if ~ isempty(obj.masterSlave)
                iD = obj.fileID;
                fprintf(iD,'\n');
                fprintf(iD,'Master_slave = [\n');
                printFormat = [repmat('%12.0d ',1,2),'\n'];
                fprintf(iD,printFormat,obj.masterSlave');
                fprintf(iD,']; \n');
            end
        end

        function printDirichletData(obj)
            iD = obj.fileID;
            fprintf(iD,'\n');
            fprintf(iD,'dirichlet_data = [\n');
            corn = obj.corners;
            ncorners = size(corn,1);
            repCorn = repmat(corn,1,2)';
            repDir  = repmat((1:obj.ndim)',ncorners,1);
            ToPrint = [repCorn(:)';repDir(:)'];
            printFormat = [repmat('%12.0f ',1,2),'       0 \n'];
            fprintf(iD,printFormat,ToPrint);
            fprintf(iD,']; \n');
        end

        function printMaterialSets(obj)
            iD = obj.fileID;
            fprintf(iD,'\n');
            fprintf(iD,'MaterialSets = [\n');
            printFormat = [repmat('%12.0d ',1,2),'\n'];
            for imatSet = 1:numel(obj.isElemInThisSet)
                elemSet = obj.isElemInThisSet{imatSet};
                n       = length(elemSet);
                toPrint = [elemSet';imatSet*ones(1,n)];
                fprintf(iD,printFormat,toPrint);
            end
            fprintf(iD,']; \n');
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function stop(robot)

    robot.mode('off');

author: petercorke

repo: robotics-toolbox-matlab


text:
% Reed-Solomon Errors and Erasures Decoding
% Decoding is based on the Massey-Berlekamp decoder.  Based on the notes of Adina Matache: http://www.ee.ucla.edu/~matache/rsc/slide.html
% This program was written by Jaco Versfeld. For any questions or remarks email Jaco at jaco.versfeld@gmail.com.
% This code is provided as is and Jaco Versfeld accepts no liability what so ever.
% Feel free to use and modify this code for whatever purpose.

clear
clc

%**************************
%*** RS code Parameters ***
%**************************

% Here we specify the parameters of the (n,k) Reed-Solomon code

m = 4          %Determine the Galois Field, GF(2^m)
n = 2^m - 1    %This is fixed for a Reed-Solomon, the length of the codeword
k = 3          %The number of data symbols, can be anything between 1 to n - 1
h = n-k
t = h/2

%**************************

%*** Generate the Galois Field and Generator polynomial ***

% This step is neccessary for Matlab.  Here we create the Galois Field which is used for
% computations of the Reed-Solomon code

field = gftuple([-1:2^m-2]', m, 2);

%Generator Polynomial:
%Lin + Costello, p.171
%The Generator polynomial is one way of encoding Reed-Solomon codes

%Construct the generator polynomial
c = [1 0];
p(1) = c(1);

for i = 1:h-1
    p(1) = gfmul(p(1),1,field);
    p(2) = 0;
    c = gfconv(c,p,field);
end
g = c;

%**************************

%*** RS Encode ***

%Generate Random Data
DATA_IN = randint(1,k,[-1 n-1]);

%RS encoding
parity = RS_ENC4(DATA_IN,n,k,g,field);
RS_CODE = [parity DATA_IN];

%********************************

%*** Channel ***
RECEIVED = RS_CODE

%I introduce the errors manually here, but any channel can be used here, like an AWGN.
% A maximum of 2*t + e <= (n-k) errors and erasures can be corrected, where t is the number of
% errors and e the number of erasures

%Introduce some errors
RECEIVED(3) = gfadd(RECEIVED(3),randint(1,1,[-1 n-1]),field);
RECEIVED(5) = gfadd(RECEIVED(3),randint(1,1,[-1 n-1]),field);

%Introduce some erasures
erasures = [ 2 3 4 5 6 7 8 9 11 14 15];  %This polynomial contains the positions of the erasures

%RECEIVED(1) = -2
RECEIVED(2) = -2
RECEIVED(3) = -2
RECEIVED(4) = -2
RECEIVED(5) = -2
RECEIVED(6) = -2
RECEIVED(7) = -2
RECEIVED(8) = -2
RECEIVED(9) = -2
RECEIVED(11) = -2
RECEIVED(14) = -2
RECEIVED(15) = -2

%****************

%*** Decoding ***

DECODED = RS_E_E_DEC(RECEIVED, erasures,n,k,t,h,g,field);

%****************

DECODED
RS_CODE

if all(DECODED == RS_CODE)
    disp('Decoding Success')
else
    disp('Decoding Failure')
end


text:
function th=rotateticklabel(h,rot,demo)
%ROTATETICKLABEL rotates tick labels
%   TH=ROTATETICKLABEL(H,ROT) is the calling form where H is a handle to
%   the axis that contains the XTickLabels that are to be rotated. ROT is
%   an optional parameter that specifies the angle of rotation. The default
%   angle is 90. TH is a handle to the text objects created. For long
%   strings such as those produced by datetick, you may have to adjust the
%   position of the axes so the labels don't get cut off.
%
%   Of course, GCA can be substituted for H if desired.
%
%   TH=ROTATETICKLABEL([],[],'demo') shows a demo figure.
%
%   Known deficiencies: if tick labels are raised to a power, the power
%   will be lost after rotation.
%
%   See also datetick.

%   Written Oct 14, 2005 by Andy Bliss
%   Copyright 2005 by Andy Bliss

%DEMO:
if nargin==3
    x=[now-.7 now-.3 now];
    y=[20 35 15];
    figure
    plot(x,y,'.-')
    datetick('x',0,'keepticks')
    h=gca;
    set(h,'position',[0.13 0.35 0.775 0.55])
    rot=90;
end

%set the default rotation if user doesn't specify
if nargin==1
    rot=90;
end
%make sure the rotation is in the range 0:360 (brute force method)
while rot>360
    rot=rot-360;
end
while rot<0
    rot=rot+360;
end
%get current tick labels
a=get(h,'XTickLabel');
%erase current tick labels from figure
set(h,'XTickLabel',[]);
%get tick label positions
b=get(h,'XTick');
c = get(h,'YLim');
ydir = get(h,'YDir');
if strcmpi(ydir,'reverse'),
  c = c([2,1]);
end
%make new tick labels
if rot<180
    th=text(b,repmat(c(1)-.1*(c(2)-c(1)),length(b),1),a,'HorizontalAlignment','right','rotation',rot,'Parent',h);
    for i = 1:numel(th),
      ex = get(th(i),'Extent');
      pos = get(th(i),'Position');
      if strcmpi(ydir,'normal'),
        off = c(1) - (ex(2)+ex(4));
        pos(2) = pos(2) + off;
        set(th(i),'Position',pos);
      else
        off = ex(2)-ex(4) - c(1);
        pos(2) = pos(2) - off;
        set(th(i),'Position',pos);
      end
    end
else
    th=text(b,repmat(c(1)-.1*(c(2)-c(1)),length(b),1),a,'HorizontalAlignment','left','rotation',rot,'Parent',h);
end
set(th,'Interpreter','none');

author: kristinbranson

repo: JAABA