text:
function varargout = metCompareGUI(varargin)
% Creates a comparison GUI for deciding which metabolite are
% in a new or existing reaction.
% Called by `metCompare`, calls `findMetMatch`, `colText`, `addMetInfo`.
%
% USAGE:
%
%    [RxnInfo, stopFlag] = metCompareGUI(RxnInfo)
%
% INPUTS
%    RxnInfo:     Structure containing relevent info. See `metCompare` function
%                 `fillRxnInfo`.
%
% OPTIONAL OUTPUTS:
%    RxnInfo:     Update structure.
%    stopFlag:    Indicates if more metabolites need to be reviewed
%    CMODEL:      global input
%    TMODEL:      global input
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

gui_Singleton = 1; % Begin initialization code - DO NOT EDIT
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @metCompareGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @metCompareGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
end
% End initialization code - DO NOT EDIT

% --- Executes just before metCompareGUI is made visible.
function metCompareGUI_OpeningFcn(hObject, eventdata, handles, varargin)

% Declare globals
global CMODEL TMODEL

% workaround to enable java handles
axes(handles.axes1)
set(handles.axes1, 'Visible', 'off')

% Choose default command line output for metCompareGUI
handles.output = hObject;

% Load rxn, metabolite, and match data.
handles.RxnInfo = varargin{1};

% Populate reaction info.
cModelName = CMODEL.description ;
set(handles.staticCmodel, 'String', cModelName) ;
Tmodels = fieldnames(TMODEL.Models) ;
if length(Tmodels) == 1
    set(handles.staticTmodel, 'String', Tmodels{1})
end
rxnNum = num2str(handles.RxnInfo.rxnIndex) ;

% Set the reaction number if the metabolites all come from a reaction.
if handles.RxnInfo.rxnIndex
    set(handles.CRxnNumField, 'String', rxnNum) ;
else
    % Otherwise change how GUI appears.
    set(handles.staticCRxnNum, 'String', 'Met Number:') ;
    set(handles.staticRxnName, 'String', 'In Crxn(s):') ;
    set(handles.staticRxnEquation, 'String', 'In Trxn(s):') ;
    set(handles.staticTRxnNum, 'String', 'Current match in T.') ;
    set(handles.staticNewRxnEquation, 'Visible', 'Off') ;
    set(handles.newRxnEquationField, 'Visible', 'Off') ;
    % And disable the ability to skip checking the mets.
    set(handles.buttonSkip, 'Enable', 'off') ;
end
rxnName = handles.RxnInfo.rxnName ;
set(handles.rxnNameField, 'String', rxnName) ;
rxnEq = handles.RxnInfo.matchRxnEquation ;
set(handles.rxnEquationField, 'String', rxnEq) ;

rxnMatch = handles.RxnInfo.rxnMatch ;

% Indicate if the reaction is new, has already been matched, or N/A.
if ~rxnMatch && handles.RxnInfo.rxnIndex
    set(handles.TRxnNumField, 'String', 'New Rxn') ;
elseif ~rxnMatch && ~handles.RxnInfo.rxnIndex
    set(handles.TRxnNumField, 'String', 'N/A') ;
else
    set(handles.TRxnNumField, 'String', num2str(rxnMatch)) ;
    % You cannot choose new metabolite for reactions with matches.
    set(handles.chooseNew, 'Enable', 'off') ;
end

% List unseen metabolites.
if handles.RxnInfo.rxnIndex
  handles.RxnInfo.unseen = ~handles.RxnInfo.goodMatch ;
else
    % If just comparing metabolites, do not mark any as seen.
    handles.RxnInfo.unseen = ones(handles.RxnInfo.nMets, 1) ;
end

% Put metabolite names in the popup menu.
set(handles.popup_met, 'String', handles.RxnInfo.metData(1, :)) ;

% Populate first metabolite info and matches
% nowMet is the index of the metabolite in RxnInfo. It is selected for by
% the first metabolite that has not been seen.
handles.RxnInfo.nowMet = find(handles.RxnInfo.unseen == 1, 1) ;
if isempty(handles.RxnInfo.nowMet)
    handles.RxnInfo.nowMet = 1 ;
end
set(handles.popup_met, 'Value', handles.RxnInfo.nowMet) ;
popup_met_Callback(handles.popup_met, eventdata, handles)

% Create reaction equation based of best match suggestions.
handles = metMatchRxnEquation(hObject, handles) ;

% Update stats
fillStats(handles)

% Sets the cell selection callback feature for the uitables.
set(handles.uitable_matches, 'CellSelectionCallback', {@matchMetCallback, ...
                                                     handles})
set(handles.uitable_met, 'CellSelectionCallback', {@metTableCallback, ...
                                                 handles})

% Update handles structure ;
guidata(hObject, handles) ;

% UIWAIT makes metCompareGUI wait for user response (see UIRESUME)
uiwait(handles.figure1);
end

% --- Outputs from this function are returned to the command line.
function varargout = metCompareGUI_OutputFcn(hObject, eventdata, handles)
% Get default command line output from handles structure
handles.output = handles.RxnInfo ;
varargout{1} = handles.output ;
varargout{2} = handles.stopFlag ;
close
end

%% Button and Input Change Functions
% --- Executes on button press in buttonChoose.
function handles = buttonChoose_Callback(hObject, eventdata, handles)
global CMODEL
RxnInfo = handles.RxnInfo ;
% nowMetIndex is the index of the metabolite in CMODEL
nowMetIndex = RxnInfo.metIndex(RxnInfo.nowMet) ;

choice = get(get(handles.chooseMatch, 'SelectedObject'), 'Tag') ;

matchScores = RxnInfo.matchScores{RxnInfo.nowMet} ;
matchIndex  = RxnInfo.matchIndex{RxnInfo.nowMet} ;
% [matchScores, matchIndex] = findMetMatch(nowMetIndex,RxnInfo.rxnMatch) ;

switch choice
    case 'choose1'
        RxnInfo.matches(RxnInfo.nowMet) = matchIndex(1) ;
    case 'choose2'
        RxnInfo.matches(RxnInfo.nowMet) = matchIndex(2) ;
    case 'choose3'
        RxnInfo.matches(RxnInfo.nowMet) = matchIndex(3) ;
    case 'choose4'
        RxnInfo.matches(RxnInfo.nowMet) = matchIndex(4) ;
    case 'choose5'
        RxnInfo.matches(RxnInfo.nowMet) = matchIndex(5) ;
    case 'chooseNew'
        % if metabolite has been declared new, ensure that all reactions it is
        % involved in are also declared new.
        nowRxnIndex = find(CMODEL.S(nowMetIndex,:) ~= 0) ;
        if sum(RxnInfo.rxnList(nowRxnIndex) > 0) > 0

            set(handles.errorField, 'String',...
                [CMODEL.mets{nowMetIndex} 'Cannot be new, it is in matched Rxns.' ...
                num2str(nowRxnIndex(RxnInfo.rxnList(nowRxnIndex) > 0)) ])
            return
        end
        RxnInfo.matches(RxnInfo.nowMet) = 0 ;
    case 'chooseOther'
         newRxnIndex = str2double(get(handles.chooseOtherNo, 'String')) ;
         if newRxnIndex > 0
             RxnInfo.matches(RxnInfo.nowMet)
         else
             set(handles.errorField, 'String',...
                'Please enter a valid Rxn number.')
             return
         end
    otherwise
end

% Add data back to handle.
RxnInfo.matches = RxnInfo.matches ;

% Indicate the user has reviewed this match.
RxnInfo.goodMatch(RxnInfo.nowMet) = 1 ;

% Update metabolite as seen.
RxnInfo.unseen(RxnInfo.nowMet) = 0 ;

% Go to next metabolite in the list that has not been seen.
val = find(RxnInfo.unseen == 1, 1) ;
% If all matches are good then just go to the first metabolites.
if isempty(val)
    val = 1 ;
end
set(handles.popup_met, 'Value', val) ;

% Update handles.
handles.RxnInfo = RxnInfo ;
guidata(hObject, handles) ;

% Update new reaction equation. This also checks if the decision was good.
handles = metMatchRxnEquation(hObject, handles) ;

% Go to next met.
popup_met_Callback(handles.popup_met, eventdata, handles)
end

% --- Executes on selection change in popup_met.
function popup_met_Callback(hObject, eventdata, handles)
% Get index of choice and associate with met, save to handles.
val = get(hObject, 'Value') ;
handles.RxnInfo.nowMet = val ;
guidata(hObject, handles) ;
% Populate table with met and match information.
populateMetTables(handles) ;
end

% --- Executes on selection of item in uitable_matches
function matchMetCallback(hObject, eventdata, handles)
global TMODEL
if ~isempty(eventdata.Indices)
    % Declare variables
    RxnInfo = handles.RxnInfo ; % For convenience.
    RxnInfo.nowMet = get(handles.popup_met, 'Value') ;
    nowMetIndex = RxnInfo.metIndex(RxnInfo.nowMet) ;
    matchScores = RxnInfo.matchScores{RxnInfo.nowMet} ;
    matchIndex  = RxnInfo.matchIndex{RxnInfo.nowMet} ;
%     [matchScores,matchIndex] = findMetMatch(nowMetIndex,RxnInfo.rxnMatch) ;

    % Set match reaction number based on column clicked.
    set(handles.chooseOtherNo,'String', 'Met #')
    matchCol = eventdata.Indices(2) ;
    if matchCol == 1
        set(handles.chooseMatch, 'SelectedObject', handles.choose1) ;
    elseif matchCol == 2
        set(handles.chooseMatch, 'SelectedObject', handles.choose2) ;
    elseif matchCol == 3
        set(handles.chooseMatch, 'SelectedObject', handles.choose3) ;
    elseif matchCol == 4
        set(handles.chooseMatch, 'SelectedObject', handles.choose4) ;
    elseif matchCol == 5
        set(handles.chooseMatch, 'SelectedObject', handles.choose5) ;
    else
        set(handles.chooseMatch, 'SelectedObject', handles.chooseOther) ;
        set(handles.chooseOtherNo, 'String', num2str(matchIndex(matchCol)))
    end

    % Opens up KEGG ID site if KEGGID is selected
    matchRow = eventdata.Indices(1) ;
    % If we are in the right row
    if matchRow == 7
       % Grab the ID.
       cID = TMODEL.metKEGGID{matchIndex(matchCol)} ;
       if ~isempty(cID)
           % Enforce the cID is the right length.
           cID = cID(1:6) ;
           % Launch website.
           KEGGIDurl = ['http://www.genome.jp/dbget-bin/www_bget?' cID] ;
           web(KEGGIDurl, '-new', '-noaddressbox', '-notoolbar')
       end
    end
end
end

% --- Executes on selection of item in uitable_met
function metTableCallback(hObject, eventdata, handles)
if ~isempty(eventdata.Indices)
    RxnInfo.nowMet = get(handles.popup_met, 'Value') ;
    % Opens up KEGG ID site if KEGGID is selected
    matchRow = eventdata.Indices(1) ;
    % If we are in the right row.
    if matchRow == 7
       % Grab the ID.
       cID = handles.RxnInfo.metData{5, RxnInfo.nowMet} ;
       if ~isempty(cID)
           % Enforce the cID is the right length.
           cID = cID(1:6) ;
           % Launch browser.
           KEGGIDurl = ['http://www.genome.jp/dbget-bin/www_bget?' cID] ;
           web(KEGGIDurl, '-new', '-noaddressbox', '-notoolbar')
       end
    end
end
end

% --- Executes on button press in buttonAddMets.
function buttonAddMets_Callback(hObject, eventdata, handles)
% Do not stop the next iteration of calling metCompare.
handles.stopFlag = 0 ;
guidata(hObject, handles) ;
% Allows the wait command to be suspended, closing the GUI.
uiresume(handles.figure1)
end

% --- Executes on button press in buttonSkip.
function buttonSkip_Callback(hObject, eventdata, handles)
% Stop the next iteration of calling metCompare.
handles.stopFlag = 1 ;
guidata(hObject, handles) ;
% Allows the wait command to be suspended, closing the GUI.
uiresume(handles.figure1)
end

% --- Executes when selected object is changed in chooseMatch. NOT USED.
function chooseMatch_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in chooseMatch
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
end

% --- Executes on change in buttonChoose. NOT USED.
function chooseOtherNo_Callback(hObject, eventdata, handles)
% hObject    handle to chooseOtherNo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
end

% --- Executes on change in editNMatches. NOT USED.
function editNMatches_Callback(hObject, eventdata, handles)
% hObject    handle to editNMatches (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
end

%% Subfunctions
% Populates both met and match tables.
function populateMetTables(handles)
% Pull out variables for easier referencing. Declare globals.
RxnInfo = handles.RxnInfo ;
global CMODEL TMODEL

% resize Row headers uitable_met
jscroll = findjobj(handles.uitable_met) ;
rowHeaderViewport = jscroll.getComponent(4) ;
rowHeader = rowHeaderViewport.getComponent(0);
newWidth = 125 ;
rowHeaderViewport.setPreferredSize(java.awt.Dimension(newWidth, 0));
height = rowHeader.getHeight;
rowHeader.setPreferredSize(java.awt.Dimension(newWidth, height));
rowHeader.setSize(newWidth,height);
% format left
rend = rowHeader.getCellRenderer(1, 0);
rend.setHorizontalAlignment(javax.swing.SwingConstants.LEFT);
jscroll.repaint %apply changes

nowMetIndex = RxnInfo.metIndex(RxnInfo.nowMet) ;
metDataTable = RxnInfo.metData(:, RxnInfo.nowMet) ;
nowMetDataTable = [ {''} ; ...
                     metDataTable{1} ; ...
                     CMODEL.mets{nowMetIndex}(end-1) ; ...
                     metDataTable(2:end)] ;

% Index of the met in CMODEL.
nMatches = str2double(get(handles.editNMatches, 'String')) ;
metMatchTable = cell(7, nMatches) ;
if nMatches <= RxnInfo.nMets
    matchScores = RxnInfo.matchScores{RxnInfo.nowMet} ;
    matchIndex  = RxnInfo.matchIndex{RxnInfo.nowMet} ;
else
    [matchScores, matchIndex] = findMetMatch(nowMetIndex, RxnInfo.rxnMatch) ;
end
for iMatch = 1:nMatches
    metMatchTable{1, iMatch} = num2str(matchScores(iMatch) / ...
                                      sum(matchScores(1:nMatches)) ) ;
    metMatchTable{2, iMatch} = [TMODEL.mets{matchIndex(iMatch)}, ', ', ...
                               num2str(matchIndex(iMatch))] ;
    metMatchTable{3, iMatch} = TMODEL.mets{matchIndex(iMatch)}(end-1) ;
    metMatchTable{4, iMatch} = TMODEL.metNames{matchIndex(iMatch)} ;
    metMatchTable{5, iMatch} = TMODEL.metFormulas{matchIndex(iMatch)} ;
    metMatchTable{6, iMatch} = num2str(TMODEL.metCharge(matchIndex(iMatch))) ;
    metMatchTable{7, iMatch} = TMODEL.metKEGGID{matchIndex(iMatch)} ;
    metMatchTable{8, iMatch} = TMODEL.metSEEDID{matchIndex(iMatch)} ;
    metMatchTable{9, iMatch} = TMODEL.metID{matchIndex(iMatch)} ;
end

% weed out accidently remaining cells in cells
for iac = 1:numel(nowMetDataTable)
    if iscell(nowMetDataTable{iac})
        nowMetDataTable{iac} = [nowMetDataTable{iac}{:}] ;
    end
end
for iac = 1:numel(metMatchTable)
    if iscell(metMatchTable{iac})
        metMatchTable{iac} = [metMatchTable{iac}{:}] ;
    end
end

% If compartment doesn't match or if score is below 0.3
if ~strcmpi(metMatchTable{3, 1}, nowMetDataTable{3}) ...
        || str2double(metMatchTable{1, 1}) < 0.3
    set(handles.chooseNew, 'Value', 1)
else
    set(handles.choose1, 'Value', 1)
end

% Set color of good and bad correspondence in table.
colOptions = {'blue', 'red'} ;
for ic = 1:size(metMatchTable, 2)
    metMatchTable{1, ic} = ['<text>  ' metMatchTable{1,ic} ...
        '</text>             '] ;
    metMatchTable{2, ic} = colText(['  ' metMatchTable{2, ic}], colOptions{ ...
        2 - strcmpi(metMatchTable{2, ic}(1:strfind(metMatchTable{2, ic}, ...
                                                ', ')-1), ...
        nowMetDataTable{2, 1}(1:strfind(nowMetDataTable{2, 1}, ',') - 1))}) ;
    for ir = 3:size(metMatchTable, 1)
        if ~isempty(metMatchTable{ir, ic}) && ...
                ~isempty(nowMetDataTable{ir, 1})
            metMatchTable{ir, ic} = colText(['  ' metMatchTable{ir, ic}], ...
                colOptions{1 + isempty(strfind(metMatchTable{ir, ic}, nowMetDataTable{ir, 1}))}) ;
        else
            metMatchTable{ir, ic} = ['  ' metMatchTable{ir, ic}] ;
        end
    end
end

% If we are just comparing metabolites, put in information about what
% reactions that metabolite is involved in.
if ~RxnInfo.rxnIndex
    % Metabolite number in C.
    set(handles.CRxnNumField, 'String', nowMetIndex)
    % Current match in T.
    set(handles.TRxnNumField, 'String', RxnInfo.metList(nowMetIndex)) ;
    % Reactions it is involved in C.
    involvedRxns = find(CMODEL.S(nowMetIndex, :)) ;
    set(handles.rxnNameField, 'String', involvedRxns)
    % Reactions that it is matched to in T.
    set(handles.rxnEquationField, 'String', ...
        num2str(RxnInfo.rxnList(involvedRxns))) ;

    % Also give a hint as how best to proceed.
    if isempty(involvedRxns)
        hintString = 'Metabolite not in C, trivial problem, declare new.' ;
        set(handles.chooseNew, 'Value', 1)
    elseif ~RxnInfo.rxnList(involvedRxns)
    	hintString = ...
            'Metabolite not in a reaction matched to T, declare new' ;
        set(handles.chooseNew, 'Value', 1)
    else
        hintString = ...
            'Metabolite matched to a reaction in T. Best to keep match.' ;
        set(handles.choose1, 'Value', 1)
    end
    set(handles.errorField, 'String', hintString)
end

set(handles.uitable_matches, 'Data', metMatchTable, 'Fontsize', 10) ;

for ic = 1:size(nowMetDataTable, 2)
    for ir = 1:size(nowMetDataTable, 1)
        % add leading whitespace to all cells
        nowMetDataTable{ir, ic} = ['  ' nowMetDataTable{ir, ic}] ;
    end
end
set(handles.uitable_met, 'Data', nowMetDataTable, 'Fontsize', 10) ;

% Update unseen mets field
unseenMets = '' ;
for iMet = 1:RxnInfo.nMets
    if handles.RxnInfo.unseen(iMet) == 1
        unseenMets = [unseenMets RxnInfo.metData{1, iMet} '; '] ;
    end
end
set(handles.unseenMetsField, 'String', unseenMets) ;

end

% Makes reaction equation based on current best matches and checks the
% validity of the choices.
function handles = metMatchRxnEquation(hObject, handles)
% Declare some variables just to not have to type handles so much.
global CMODEL TMODEL
RxnInfo = handles.RxnInfo ;

% If metabolites come from a reaction, make the new equation.
% Disable choosing mets until checks have passed.
set(handles.buttonAddMets, 'Enable', 'off') ;
if RxnInfo.rxnIndex
    % Build the equation.
    reactants = cell(1, 1) ;
    products = cell(1, 1) ;
    for iMet = 1:RxnInfo.nMets
        % Reactants
        if RxnInfo.metStoichs(iMet) < 0
            % First reactant.
            if isempty(reactants{1})
               % Add Stoich
               reactants{1} = num2str(abs(RxnInfo.metStoichs(iMet))) ;
               % If a match has been found in TMODEL, choose that name.
               if RxnInfo.matches(iMet) ~= 0 && ...
                       (RxnInfo.matches(iMet) < length(TMODEL.mets))
                   reactants{2} = [' ' ...
                       TMODEL.mets{RxnInfo.matches(iMet)} ' '];
               else % Add met name from CMODEL.
                   reactants{2} = [' ' ...
                       CMODEL.mets{RxnInfo.metIndex(iMet)} ' '] ;
               end
            % nth reactants.
            else
               reactants{end + 1} = ' + ' ;
               reactants{end + 1} = num2str(abs(RxnInfo.metStoichs(iMet))) ;
               if RxnInfo.matches(iMet) ~= 0 && ...
                       (RxnInfo.matches(iMet) < length(TMODEL.mets))
                   reactants{end + 1} = [' ' ...
                                   TMODEL.mets{RxnInfo.matches(iMet)} ' '];
               else
                   reactants{end + 1} = [' ' ...
                             CMODEL.mets{RxnInfo.metIndex(iMet)} ' '] ;
               end
            end
        % Products
        else
            % First product
            if isempty(products{1})
                if RxnInfo.rev == 1
                    products{1} = '<==> ' ;
                else
                    products{1} = '--> ' ;
                end
                % Add stoich.
                products{2} = num2str(abs(RxnInfo.metStoichs(iMet))) ;
                % Add name.
                if RxnInfo.matches(iMet) ~= 0 && ...
                        (RxnInfo.matches(iMet) < length(TMODEL.mets))
                   products{3} = [' ' TMODEL.mets{RxnInfo.matches(iMet)}] ;
                else
                   products{3} = [' ' ...
                             CMODEL.mets{RxnInfo.metIndex(iMet)} ' '] ;
                end
            else % nth product.
                products{end + 1} = ' + ' ;
                products{end + 1} = num2str(abs(RxnInfo.metStoichs(iMet))) ;
                if RxnInfo.matches(iMet) ~= 0 && ...
                       (RxnInfo.matches(iMet) < length(TMODEL.mets))
                   products{end + 1} = [' ' ...
                                  TMODEL.mets{RxnInfo.matches(iMet)} ' '] ;
                else
                   products{end + 1} = [' ' ...
                             CMODEL.mets{RxnInfo.metIndex(iMet)} ' '] ;
                end
            end
        end
    end
    % If there were no products, must be exchange reaction.
    if isempty(products{1})
        products{1} = '-->' ;
    end
    % Combine both halves of the equation and print to field.
    matchRxnEquation = [reactants{1, :} products{1, :} ] ;
    set(handles.newRxnEquationField, 'String', matchRxnEquation) ;
    handles.RxnInfo.matchRxnEquation = matchRxnEquation ;

    % NaN in metList create problem, so remove them (really, they shouldn't
    % be there in the first place)
    RxnInfo.metList(isnan(RxnInfo.metList)) = 0 ;

    % Set addReaction button off, then check if all the mets in the
    % reaction have good matches, then allowing the user to add the
    % reaction. Also, check for errors.
    errorString = 'Review required.' ;
    bgcolor = [0.7 0.7 0.7] ;
    % All metabolites must be good matches and new metabolites must be
    % reviewed.
    if ~ismember(1, RxnInfo.unseen) && ~ismember(0, RxnInfo.goodMatch)
        errorString = 'All clear.' ;
        set(handles.buttonAddMets, 'Enable', 'on') ;
    else
        for iMet = 1:RxnInfo.nMets
            if RxnInfo.metList(RxnInfo.metIndex(iMet)) ~= ...
                    RxnInfo.matches(iMet) && ...
                    RxnInfo.metList(RxnInfo.metIndex(iMet))
                errorString = ['ERROR: cMet ' ...
                    num2str(RxnInfo.metIndex(iMet)) ...
                    ' already assigned to tMet ' ...
                    num2str(RxnInfo.metList(RxnInfo.metIndex(iMet))) ...
                    '.'] ;
                bgcolor = [1 0 0] ;
                set(handles.buttonAddMets, 'Enable', 'off') ;
            elseif ~RxnInfo.metList(RxnInfo.metIndex(iMet)) ...
                    && ~isempty(find(RxnInfo.metList ...
                    == RxnInfo.matches(iMet), 1)) ...
                    && RxnInfo.matches(iMet) ~= 0
                errorString = ['ERROR: Metabolite ' ...
                    num2str(RxnInfo.metIndex(iMet)) ...
                    's match (tMet ' ...
                    num2str(RxnInfo.matches(iMet)) ...
                    ') already matched to cMet ' ...
                    num2str(find(RxnInfo.metList ...
                    == RxnInfo.matches(iMet), 1)) ...
                    '.'] ;
                bgcolor = [1 0 0] ;
                set(handles.buttonAddMets, 'Enable', 'off') ;
            end
        end
    end
% Set the error bar color to indicate an error or not.
set(handles.errorField, 'String', errorString, 'BackgroundColor', bgcolor) ;
else
    if ~RxnInfo.unseen
        set(handles.buttonAddMets, 'Enable', 'on') ;
    end
end

guidata(hObject, handles) ;
end

% Update stats.
function fillStats(handles)
RxnInfo = handles.RxnInfo ;
rxnString = [num2str(length(find(RxnInfo.rxnList >= 0))) ' / ' ...
             num2str(length(RxnInfo.rxnList))] ;
metString = [num2str(length(find(RxnInfo.metList))) ' / ' ...
             num2str(length(RxnInfo.metList))] ;
set(handles.text_metsDeclared, 'String', metString)
set(handles.text_rxnsAdded, 'String', rxnString)
end

%% Object Creation Functions
function popup_met_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popup_met (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

if ispc && isequal(get(hObject, 'BackgroundColor'), get(0, 'defaultUicontrolBackgroundColor'))
    set(hObject, 'BackgroundColor', 'white');
end
end

function chooseOtherNo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to chooseOtherNo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

if ispc && isequal(get(hObject, 'BackgroundColor'), get(0, 'defaultUicontrolBackgroundColor'))
    set(hObject, 'BackgroundColor', 'white');
end
end

function editNMatches_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editNMatches (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

if ispc && isequal(get(hObject, 'BackgroundColor'), get(0, 'defaultUicontrolBackgroundColor'))
    set(hObject, 'BackgroundColor', 'white');
end
end

function pushbutton_addinfo_Callback(hObject, eventdata, handles)
nowmet = get(handles.popup_met, 'String') ;
nowmet = nowmet{get(handles.popup_met, 'Value')} ;
nowmet = regexp(nowmet, ',', 'split') ;
nowmet = str2double(nowmet{end}) ;
addMetInfo(nowmet)

popup_met_Callback(handles.popup_met, [], handles)

end

author: opencobra

repo: cobratoolbox


text:
function [all_data, units,apt_info] = compute_apt(trx,n,fn)

flies = trx.exp2flies{n};
nflies = numel(flies);
all_data = cell(1,nflies);
all_units = cell(1,nflies);
parts = strsplit(fn,'_');
view_str = regexp(parts{1},'view(\d*)','tokens');
view = str2double(view_str{1}{1});
fn_type = parts{2};
comp_type = parts{3};

apt_info = trx(flies(1)).trkInfo;
if isfield(apt_info,'crop_loc') && ~isempty(apt_info.crop_loc) && numel(apt_info.crop_loc)==4
  crop_loc = apt_info.crop_loc;
  x_center = mean(crop_loc(1:2));
  y_center = mean(crop_loc(3:4));
  global_center = [y_center x_center];
else
  global_center = double(cell2mat(apt_info.params.imsz))/2;
end

if startsWith(fn_type,'social')
  [all_data,units] = compute_apt_social(trx,n,fn,apt_info);
  return;
end

for fndx = 1:nflies
  if view > 1
    x = trx(flies(fndx)).(sprintf('x_view%d',view));
    y = trx(flies(fndx)).(sprintf('y_view%d',view));
    theta = trx(flies(fndx)).(sprintf('theta_view%d',view));
  else
    x = trx(flies(fndx)).x;
    y = trx(flies(fndx)).y;
    theta = trx(flies(fndx)).theta;
  end
  dt = trx(flies(fndx)).dt;
  nframes = trx(flies(fndx)).nframes;
  pxpermm = trx.pxpermm;
  n_parts = apt_info.params.n_classes;
  mod_apt_data = trx(flies(fndx)).kpts;
  mod_apt_data = reshape(mod_apt_data,n_parts,[],nframes);

  switch fn_type
    case 'global'
      part_num = str2double(parts{end});
      [data,cur_units] = compute_global(mod_apt_data,dt,theta,comp_type,part_num,pxpermm,global_center);
    case 'body'
      part_num = str2double(parts{end});
      [data,cur_units] = compute_body(mod_apt_data,x,y,dt,theta,comp_type,part_num,pxpermm);
    case 'pair'
      part2 = str2double(parts{end});
      part1 = str2double(parts{end-1});
      x1 = permute(mod_apt_data(part1,1,:),[1,3,2]);
      y1 = permute(mod_apt_data(part1,2,:),[1,3,2]);
      x2 = permute(mod_apt_data(part2,1,:),[1,3,2]);
      y2 = permute(mod_apt_data(part2,2,:),[1,3,2]);
      [x1,y1] = convert_to_body(x1,y1,x,y,theta);
      [x2,y2] = convert_to_body(x2,y2,x,y,theta);
      [data,cur_units] = pair_fn(x1,y1,x2,y2,dt,comp_type,pxpermm);
    case 'triad'
      part3 = str2double(parts{end});
      part2 = str2double(parts{end-1});
      part1 = str2double(parts{end-2});
      x1 = permute(mod_apt_data(part1,1,:),[1,3,2]);
      y1 = permute(mod_apt_data(part1,2,:),[1,3,2]);
      x2 = permute(mod_apt_data(part2,1,:),[1,3,2]);
      y2 = permute(mod_apt_data(part2,2,:),[1,3,2]);
      x3 = permute(mod_apt_data(part3,1,:),[1,3,2]);
      y3 = permute(mod_apt_data(part3,2,:),[1,3,2]);
      [data,cur_units] = triad_fn(x1,x2,x3,y1,y2,y3,dt,comp_type,pxpermm);

    otherwise
      error('Undefined computation type type %s',comp_type);
  end
  all_data{fndx} = data;
  all_units{fndx} = cur_units;
end
units = all_units{1};

function [data,units] = compute_global(apt_data, dt, theta, comp_type, part,pxpermm,global_center)
x = permute(apt_data(part,1,:),[1,3,2]);
y = permute(apt_data(part,2,:),[1,3,2]);
relative_fns = {'x','y','dx','dy','sin','cos','dtheta'};
if any(strcmp(relative_fns,comp_type))
  i_x = global_center(2);
  i_y = global_center(1);
  [data,units] = relative(x,y,i_x,i_y,dt,fn,pxpermm);
else
  switch comp_type
    case 'velmag'
     [data,units] = velmag(x, y, dt,pxpermm);
    case 'dphi'
     [data,units] = dphi(x, y, dt, theta);
    case 'dvelmag'
     [data,units] = dvelmag(x, y, dt,pxpermm);
    otherwise
      error('Undefined computation type type %s',comp_type);
  end
end

function [b_x,b_y] = convert_to_body(x,y,trx_x,trx_y,theta)
% convert global coordinates into body coordinates
b_x = nan(size(x)); b_y = nan(size(y));
for t = 1:numel(x)
  T = [1,0,0
    0,1,0
    -trx_x(t),-trx_y(t),1];
  R = [cos(theta(t)+pi/2),-sin(theta(t)+pi/2),0
    sin(theta(t)+pi/2),cos(theta(t)+pi/2),0
    0,0,1];
  A = T*R;
  M = [x(t) y(t) 1]*A;
  b_x(t) = M(1); b_y(t) = M(2);
end

function [data,units] = compute_body(apt_data, trx_x, trx_y, dt, theta, comp_type, part,pxpermm)
relative_fns = {'x','y','dx','dy','sin','cos','dtheta'};
x = permute(apt_data(part,1,:),[1,3,2]);
y = permute(apt_data(part,2,:),[1,3,2]);
[x,y] = convert_to_body(x,y,trx_x,trx_y,theta);

switch comp_type
  case 'velmag'
   [data,units] = velmag(x, y, dt, pxpermm);
  case 'dphi'
   [data,units] = dphi(x, y, dt, theta);
  case 'dvelmag'
   [data,units] = dvelmag(x, y, dt, pxpermm);
  case relative_fns
    [data,units] = relative(x,y,zeros(size(x)),zeros(size(y)),dt, comp_type, pxpermm);
  case 'distcenter'
    [data,units] = dist_center(x,y,zeros(size(x)),zeros(size(y)), pxpermm);
  case 'ddistcenter'
    [data,units] = ddist_center(x,y,zeros(size(x)),zeros(size(y)),dt, pxpermm);
  otherwise
    error('Undefined computation type type %s',comp_type);
end

function [data, units] = velmag(x,y,dt,pxpermm)
if numel(x) == 1
  data = 0;
else
  dx = diff(x,1,2);
  dy = diff(y,1,2);
  data = sqrt(dx.^2 + dy.^2)./dt;
end
data = data/pxpermm;
units = parseunits('mm/s');

function [data, units] = dvelmag(x,y,dt,pxpermm)
if numel(x) == 1
  data = [];
elseif numel(x) == 2
  data = 0;
else
  dx = diff(x,1,2);
  dy = diff(y,1,2);
  tmp = sqrt(diff(dx./dt,1,2).^2 + diff(dy./dt,1,2).^2)./dt(2:end);
  data = [0,tmp];
end
data = data/pxpermm;
units = parseunits('mm/s/s');

function [data, units] = dphi(x,y,dt,theta)
if numel(x) > 1
  dy1 = [y(2)-y(1),(y(3:end)-y(1:end-2))/2,y(end)-y(end-1)];
  dx1 = [x(2)-x(1),(x(3:end)-x(1:end-2))/2,x(end)-x(end-1)];
  badidx = dy1 == 0 & dx1 == 0;
  phi = atan2(dy1,dx1);
  phi(badidx) = theta(badidx);
  data = modrange(diff(phi),-pi,pi)./dt;
else
  data = [];
end
units = parseunits('rad/s');

function [data, units] = dist_center(x,y,x_body,y_body,pxpermm)
data = sqrt((x-x_body).^2 + (y-y_body).^2);
data = data/pxpermm;
units = parseunits('mm');

function [data, units] = ddist_center(x,y, x_body, y_body,dt,pxpermm)
if numel(x)>1
  dist = sqrt((x-x_body).^2 + (y-y_body).^2);
  data = diff(dist)./dt;
else
  data = 0;
end
data = data/pxpermm;
units = parseunits('mm/s');

function [data, units] = relative(x1, y1, x2, y2, dt, fn,pxpermm)

switch fn
  case 'x'
    data = x1-x2;
    data = data/pxpermm;
    units = parseunits('mm');
  case 'y'
    data = y1-y2;
    data = data/pxpermm;
    units = parseunits('mm');
  case 'dx'
    data = diff(x1-x2)./dt;
    data = data/pxpermm;
    units = parseunits('mm/s');
  case 'dy'
    data = diff(y1-y2)./dt;
    data = data/pxpermm;
    units = parseunits('mm');
  case 'sin'
    len = sqrt((x1-x2).^2 + (y1-y2).^2);
    data = (x1-x2)./len;
    data(len==0) = 0;
    units = parseunits('');
  case 'cos'
    len = sqrt((x1-x2).^2 + (y1-y2).^2);
    data = (y1-y2)./len;
    data(len==0) = 0;
    units = parseunits('');
  case 'dtheta'
    if numel(x1) > 2
      theta = atan2(y1-y2, x1-x2);
      data = modrange(diff(theta,1,2),-pi,pi)./dt;
    else
      data = [];
    end
    units = parseunits('rad/s');
end

function [data,units] = pair_fn(x1,y1,x2,y2,dt, fn,pxpermm)

relative_fns = {'x','y','dx','dy','sin','cos','dtheta'};
if any(strcmp(relative_fns,fn))
  [data,units] = relative(x1,y1,x2,y2,dt,fn,pxpermm);
elseif strcmp(fn,'velmag')
  [data,units] = velmag(x1-x2,y1-y2,dt,pxpermm);
elseif strcmp(fn,'dvelmag')
  [data,units] = dvelmag(x1-x2,y1-y2,dt,pxpermm);
elseif strcmp(fn,'dist')
  [data,units] = dist_center(x1,y1,x2,y2,pxpermm);
elseif strcmp(fn,'ddist')
  [data,units] = ddist_center(x1,y1,x2,y2,dt,pxpermm);
elseif strcmp(fn,'dphi')
  theta = atan2(y1-y2,x1-x2);
  [data,units] = dphi(x1-x2, y1-y2, dt, theta);
elseif strcmp(fn,'u1')
  theta = atan2(y1-y2,x1-x2);
  dx = diff(x1,1,2)./dt;
  dy = diff(y1,1,2)./dt;
  phi = [atan2(dy,dx)];
  angle_vel = pi - theta(1:end-1) + phi;
  % This is the angle between velocity vector and the vector from x2 to x1.
  data = sqrt(dx.^2 + dy.^2).*cos(angle_vel);
  data = data/pxpermm;
  units = parseunits('mm/s');
elseif strcmp(fn,'v1')
  theta = atan2(y1-y2,x1-x2);
  dx = diff(x1,1,2)./dt;
  dy = diff(y1,1,2)./dt;
  phi = [atan2(dy,dx)];
  angle_vel = pi - theta(1:end-1) + phi;
  % This is the angle between velocity vector and the vector from x2 to x1.
  data = sqrt(dx.^2 + dy.^2).*sin(angle_vel);
  data = data/pxpermm;
  units = parseunits('mm/s');
elseif strcmp(fn,'u2')
  theta = atan2(y2-y1,x2-x1);
  dx = diff(x2,1,2)./dt;
  dy = diff(y2,1,2)./dt;
  phi = [atan2(dy,dx)];
  angle_vel = pi - theta(1:end-1) + phi;
  % This is the angle between velocity vector and the vector from x2 to x1.
  data = sqrt(dx.^2 + dy.^2).*cos(angle_vel);
  data = data/pxpermm;
  units = parseunits('mm/s');
elseif strcmp(fn,'v2')
  theta = atan2(y2-y1,x2-x1);
  dx = diff(x2,1,2)./dt;
  dy = diff(y2,1,2)./dt;
  phi = atan2(dy,dx);
  angle_vel = pi - theta(1:end-1) + phi;
  % This is the angle between velocity vector and the vector from x2 to x1.
  data = sqrt(dx.^2 + dy.^2).*sin(angle_vel);
  data = data/pxpermm;
  units = parseunits('mm/s');
elseif strcmp(fn,'areaswept')
  area = zeros(1,numel(x1));
  for ndx = 1:numel(x1)-1
    curx = [x1(ndx),x1(ndx+1),x2(ndx+1),x2(ndx)];
    cury = [y1(ndx), y1(ndx+1), y2(ndx+1), y2(ndx)];
    area(ndx) = polyarea(curx,cury);
  end
  data = area/pxpermm/pxpermm;
  units = parseunits('mm*mm');
else
  error('Unknown function for APT pair %s',fn);
end

function [data,units] = triad_fn(x1,x2,x3, y1,y2, y3, dt,fn,pxpermm)
theta1 = atan2(y1-y2,x1-x2);
theta2 = atan2(y3-y2,x3-x2);
theta = modrange(theta1-theta2,-pi,pi);
switch fn
  case 'cos'
    data = cos(theta);
    units = parseunits('');
  case 'sin'
    data = sin(theta);
    units = parseunits('');
  case 'dangle'
    data = diff(theta)./dt;
    units = parseunits('rad/s');
  case 'area'
    area = zeros(1,numel(x1));
    for ndx = 1:numel(x1)
      area(ndx) = polyarea([x1(ndx),x2(ndx),x3(ndx)],[y1(ndx),y2(ndx),y3(ndx)]);
    end
    data = area;
    data = data/pxpermm/pxpermm;
    units = parseunits('mm*mm');
  case 'darea'
    area = zeros(1,numel(x1));
    for ndx = 1:numel(x1)
      area(ndx) = polyarea([x1(ndx),x2(ndx),x3(ndx)],[y1(ndx),y2(ndx),y3(ndx)]);
    end
    data = diff(area);
    data = data/pxpermm/pxpermm;
    units = parseunits('mm*mm/s');
  case 'dlen'
    len1 = sqrt((x1-x2).^2 + (y1-y2).^2);
    len2 = sqrt((x3-x2).^2 + (y3-y2).^2);
    data = len1-len2;
    data = data/pxpermm;
    units = parseunits('mm');
end

author: kristinbranson

repo: JAABA


text:
function gmt_Colormap_ConvertPovRayToGMT(sPovRayFile, sGMTFile, fMin, fMax, bFlip)

mMatlabColormap = gui_Colormap_ReadPovRay(sPovRayFile, 256);
if exist('bFlip', 'var')
  if bFlip
    mMatlabColormap = flipud(mMatlabColormap);
  end
end
nRow = length(mMatlabColormap(:,1));
fDiff = fMax - fMin;
vMinRow = fMin:(fDiff/nRow):(fMax-(fDiff/nRow));
vMaxRow = (fMin+(fDiff/nRow)):(fDiff/nRow):fMax;
mGMTColormap = [vMinRow' mMatlabColormap.*255 vMaxRow' mMatlabColormap.*255];
save(sGMTFile, 'mGMTColormap', '-ascii');

author: CelsoReyes


text:
function [d] = readraw_Gannet(filename, type, ds, endian,offset)

if( nargin == 3 )
	endian = 'b';
end

fp=fopen(filename, 'rb', endian);
%fseek(fp, 0, -1);
%fseek(fp, 262144, -1);
fseek(fp, offset, -1);
d = fread(fp, prod(ds), type);
fclose(fp);
d = reshape(d, ds);

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [st] = getTimeStr()
    c = clock();
    MM = num2str(c(2) + 100); MM = MM(2:end);
    dd = num2str(c(3) + 100); dd = dd(2:end);
    hh = num2str(c(4) + 100); hh = hh(2:end);
    mm = num2str(c(5) + 100); mm = mm(2:end);
    ss = num2str(round(c(6)) + 100); ss = ss(2:end);

    st = strcat(MM,dd, '_', hh, mm, ss);
    fprintf('Time now: %s/%s/%4d - %s:%s:%s\n', MM, dd, c(1), hh, mm, ss);
    % if nargin == 2

end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function fx = p44_fun ( n, x )

%*****************************************************************************80
%
%% P44_FUN evaluates the integrand for problem 57.
%
%  Discussion:
%
%    The problem has a parameter ALPHA that can be set by calling
%    P44_PARAM_SET.
%
%    The suggested parameter range is 0.0 <= ALPHA <= 9.0.
%
%    As ALPHA increases, the integrand becomes more oscillatory.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    exp ( 20 * ( x - 1 ) ) * sin ( 2^alpha * x )
%
%  Exact Integral:
%
%    ( 20 sin ( 2^alpha ) - 2^alpha cos ( 2^alpha )
%    + 2^alpha exp ( -20 ) ) / ( 400 + 4^alpha )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 84.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  alpha = p44_param_get ( );

  fx = exp ( 20.0 * ( x(1:n) - 1.0 ) ) .* sin ( 2.0^alpha * x(1:n) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function []=draw(f,a,b)

%
% draws f in the intervall [a,b]
%

x=linspace(a,b,100);
y=feval(f,x);
plot(x,y,'k-');

print -deps p

repo: Aero-matlab


text:
function linpack_z_test18 ( )

%*****************************************************************************80
%
%% TEST18 tests ZPBDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  m = 1;
  lda = m+1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST18\n' );
  fprintf ( 1, '  For a double precision complex (C)\n' );
  fprintf ( 1, '  positive definite hermitian band matrix (PB),\n' );
  fprintf ( 1, '  ZPBDI computes the determinant as\n' );
  fprintf ( 1, '    det = MANTISSA * 10**EXPONENT\n' );
  fprintf ( 1, '\n' );
%
%  Set the value of the superdiagonal and diagonal.
%
  a(1,1) = complex ( 0.0000,  0.0000 );
  a(1,2) = complex ( 2.1341, -0.2147 );
  a(1,3) = complex ( 2.0905,  1.1505 );

  a(2,1) = complex ( 4.5281,  0.0000 );
  a(2,2) = complex ( 5.0371,  0.0000 );
  a(2,3) = complex ( 4.7638,  0.0000 );

  [ a, info ] = zpbfa ( a, lda, n, m );

  if ( info ~= 0 )
    fprintf ( 1, '  Error!  ZPBFA returns INFO = %d\n', info );
    return
  end

  det = zpbdi ( a, lda, n, m );

  fprintf ( 1, '  Determinant = %f * 10^(%f)\n', det(1), det(2) );

  return
end


text:
function [hp] = Btuph2hp(Btuph)
% Convert power from British thermal units per hour to mechanical horsepower.
% Chad A. Greene 2012
hp = Btuph*0.000393015;

author: Sable

repo: mcbench-benchmarks


text:
function [ x, y ] = hand_acquire ( )

%*****************************************************************************80
%
%% HAND_ACQUIRE prompts the user to enter data about the hand shape.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 April 2011
%
%  Author:
%
%    Cleve Moler
%    Modifications by John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625.
%
%  Parameters:
%
%    Output, real X(*), Y(*), the screen coordinates of points selected by the
%    user using mouse clicks.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HAND_ACQUIRE:\n' );
  fprintf ( 1, '  This function will help you record the shape of your hand\n' );
  fprintf ( 1, '  as a pair of X, Y data vectors.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In a moment, a blank screen will appear.\n' );
  fprintf ( 1, '  Place your hand on the screen.\n' );
  fprintf ( 1, '  Using the mouse, trace the outline of your hand, clicking\n' );
  fprintf ( 1, '  with the mouse at each point you want to record.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Once you have outlined your hand, press RETURN to save the data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Press RETURN to begin recording your hand!\n' );

  screen = get ( 0, 'ScreenSize' );
%
%  Specify the location and size of a figure window
%  that takes up the whole screen.
%
  figure ( 'Position', screen );
%
%  Create a graphical coordinate system on the figure.
%
  axes ( 'Position', [ 0, 0, 1, 1 ] );
%
%  Get input from the user.
%
  [ x, y ] = ginput ( );

  n = length ( x );

  printf ( 1, '\n' );
  printf ( 1, '  %d points were recorded.\n' );

  return
end


text:
function [ObjHandle, varargout] = guipanel(varargin)

%% Parse Child -----------------------------------------------------------%
child = inputParser;
child.StructExpand = true;
child.KeepUnmatched = true;

%Interuptible: Defaults to off
child.addParamValue('Interruptible','off')

%BusyAction: Defaults to cancel
child.addParamValue('BusyAction','queue')

%Child's Parent: Defaults to new fig if Param/Value pair is missing
if (nargin >= 1) && ~ischar(varargin{1})
    assert(ishandle(varargin{1}),'guipanel:badHandle',...
        ['When creating a guipanel inside an existing figure/panel, ',...
        'the first argument must be a currently existing/valid graphics handle'])
    child.addOptional('Parent',varargin{1})
else
    %-% NOTICE %-%
    % The figure handle is generated by evaluating a figure command at the
    % time the addParamValue method is called.
    child.addParamValue('Parent',figure);
end

%% Finish Parsing and Instantiate ----------------------------------------%
%   We have only slightly differed from the calling syntax of uipanel.
%   Specifically, guipanel may be called with a single argument, the handle
%   of a pre-existing figure or pre-exisiting panel (this allows easy
%   nesting of panels).

child.parse(varargin{:});

if ~any(strcmp('BusyAction',child.UsingDefaults));
    warning('guipanel:defaultOverWrite',...
        'Setting ''BusyAction'' to ''cancel'' may cause unexpected gui behavior')
end

if ~any(strcmp('Interruptible',child.UsingDefaults));
    warning('guipanel:defaultOverWrite',...
        'Setting ''Interruptible'' to ''on'' may cause unexpected gui behavior')
end

assert(ishandle(child.Results.Parent),'guipanel:badHandle',...
    '%s is not a valid graphics handle',num2str(child.Results.Parent))
    %Create a panel in the specified container (Figure or Panel)
    if isempty(fieldnames(child.Unmatched))
        ObjHandle = uipanel(child.Results,...
            'ButtonDownFcn',@PanelButtonDownFcn);
    else
        ObjHandle = uipanel(child.Results,...
            child.Unmatched,...
            'ButtonDownFcn',@PanelButtonDownFcn);
    end
   %inputParser class is poorly behaved, manually free mem
    clear child

    %Create togglebutton in the panel.  Togglebutton will lock panel in
    %place when activated, and release the panel when deactivated.
    uicontrol(ObjHandle,...
        'Style','togglebutton',...
        'Units','pixels',...
        'BackgroundColor','red',...
        'Position',[4 2 9 9],...
        'ToolTipString','Lock',...
        'Tag','Lock',...
        'Callback',@LockPanelCallback);

    %% Begin callbacks and guipanel subfcns

    function PanelButtonDownFcn(src,event)

        uistack(src,'top') %just in case
        if ~get(findobj('Parent',src,...
                'Type','uicontrol',...
                'Style','togglebutton',...
                'Tag','Lock'),'Value')

            %Get parent container of src
            parent = get(src,'Parent');

            %Get parent figure of src
            fig = ancestor(src,'figure');

            %Sync units
            old_units = SyncUnits([0,src,parent,fig],'pixels');

            %Position of src relative to parent figure
            pos_r2fig = getpixelposition(src,true);

            %Position of src relative to parent container
            pos_r2par = getpixelposition(src);

            %Position of pointer relative to parent figure
            point_r2fig = Root2Fig(get(0,'PointerLocation'),fig);

            [epsX,epsY] = MakeEpsilon(point_r2fig,pos_r2fig);
            if any(abs(epsX) <= 10) || any(abs(epsY) <= 10)
                %RESIZE
                set(fig,'Pointer','fleur')
                if fig == parent
                    position = rbbox(pos_r2par);
                    setpixelposition(src,position)
                else
                    position = rbbox(pos_r2fig);
                    setpixelposition(src,position,true)
                end
            else
                %MOVE
                set(fig,'Pointer','fleur')
                if fig == parent
                    position = dragrect(pos_r2fig);
                    setpixelposition(src,position)
                else
                    position = dragrect(pos_r2fig);
                    setpixelposition(src,position,true)
                end
            end
            set(fig,'Pointer','arrow')
            UnSyncUnits(old_units)
            refresh(fig)
        end
    end

    function LockPanelCallback(src,event)
        if get(src,'Value')
            set(src,'BackgroundColor','white')
        else
            set(src,'BackgroundColor','red')
        end
    end

    varargout{1} = struct(...
        'Root2Fig',@Root2Fig,...
        'Parent2Child',@Parent2Child,...
        'MakeEpsilon',@MakeEpsilon,...
        'MakeRect',@MakeRect,...
        'SyncUnits',@SyncUnits,...
        'UnSyncUnits',@UnSyncUnits,...
        'OverlapingAt',@OverlappingAt);
end

%% Utility Functions -----------------------------------------------------%
%

function coords = Root2Fig(point, fig)
% Transforms point from root-centric to child-centric coordinates.
% Point is measured in the parent's reference frame

if fig == 0
    coords = point;
    return
end

oldUnits = SyncUnits([0,fig],'Pixels');

%Get position of user specified point relative to the Root
PxRoot = point(1); PyRoot = point(2);

%Get position of Figure relative to the Root
figloc = get(fig,'Position');
FigOriginX = figloc(1); FigOriginY = figloc(2);

%Send cursor coordinates from Root to Figure reference frame:
PxRoot2Fig = (PxRoot - FigOriginX); PyRoot2Fig = (PyRoot - FigOriginY);

%Clean up and return result
coords = [PxRoot2Fig,PyRoot2Fig];
UnSyncUnits(oldUnits)

end
function coords = Parent2Child(point, child)
% Transforms point from parent-centric to child-centric coordinates.
% Point is measured in the parent's reference frame
parent = ancestor(child,{'figure','uipanel'});
fig = ancestor(child,'figure');
oldUnits = SyncUnits([0,fig,parent,child],'pixels');

%Root->Figure
r2f = Root2Fig(point, fig);

%Parent's position in Figure frame
if (parent ~= 0) && (parent ~= fig)
    parent_position = getpixelposition(parent,true);
    x = r2f(1) - parent_position(1);
    y = r2f(2) - parent_position(2);
    coords = [x,y];
else
    coords = r2f;
end

UnSyncUnits(oldUnits)

end
function varargout = MakeEpsilon(point, position_vector)
%Determine a position vector's left/right & top/bottom distances from a point
%
%   Calling syntax:
%       [epsX,epsY] = MakeEpsilon(point, vector)
%       [X_right,X_left,Y_bottom,Y_top] = MakeEpsilon(point,vector)
%       Where point is 1x2
%       or
%
%       [X_right,X_left,Y_bottom,Y_top] = MakeEpsilon(vector1,vector2)
%       Where vector1 is 1x4xN, and vector2 is 1x4.  In the case that
%       N > 1, each of X_right, X_left, Y_bottom, and Y_top will be
%       1x1xN.
%
%   If epsX(1) <= 0 &&  epsX(2) >= 0 then px is inside the
%   horizontal boundaries.
%
%   If epsY(1) <= 0 &&  epsY(2) >= 0 then py is inside the vertical
%   boundaries.
%
    assert(nargin == 2,'Too few/Too many inputs, MakeEpsilon requires two input args')
    if nargout
        %assert((nargout == 2) || (nargout == 4), 'MakeEpsilon takes either two or four output args')
    end

    sz = size(point);

    switch mat2str(sz)
        case '[1 2]'
            [epsX,epsY] = epsilonPointSubFcn(point, position_vector);
            varargout{1} = epsX; varargout{2} = epsY;
            return
        case '[1 4]'
            [X_right,X_left,Y_bottom,Y_top] = epsilonVectorSubFcn(point,position_vector);
        otherwise
            try
                [X_right,X_left,Y_bottom,Y_top] = epsilonVectorSubFcn(point,position_vector);
            catch ME
                ME.message
                error('guiSubFcn:MakeEpsilon','Cannot reconcile input')
            end
    end

    varargout{1} = X_right; varargout{2} = X_left;
    varargout{3} = Y_bottom; varargout{4} = Y_top;

    function varargout = epsilonPointSubFcn(pnt, pos)
        %Do horz
        leftXlim = pos(1)-pnt(:,1);
        rightXlim = pos(1)+pos(3)-pnt(:,1);
        epsX = [leftXlim, rightXlim];
        %Do vert
        bottomYlim = pos(2)-pnt(:,2);
        topYlim = pos(2)+pos(4)-pnt(:,2);
        epsY = [bottomYlim, topYlim];
        %Parse outputs
        if nargout ==2
            varargout{1} = epsX; varargout{2} = epsY;
        else
            varargout{1} = epsX(:,1); varargout{2} = epsX(:,2);
            varargout{3} = epsY(:,1); varargout{4} = epsY(:,2);
        end
    end
    function varargout = epsilonVectorSubFcn(pos1,pos2)
       [rectX,rectY] = MakeRect(pos1); %todo, overload MakeRect to accept pos vecs
       %rectX, rectY will each be 4x1
       pnt_array = [rectX,rectY];
       [X_right,X_left,Y_bottom,Y_top] = epsilonPointSubFcn(pnt_array, pos2);
       varargout{1} = X_right; varargout{2} = X_left;
       varargout{3} = Y_bottom; varargout{4} = Y_top;
    end

end
function [rectX,rectY] = MakeRect(position_vector)
%MakeRect takes a position vector and creates X,Y vertices for the
%corresponding rectangle. The position vector can be a (1,4,N)
%multidimensional array.
bottomLeft = [position_vector(:,1,:),position_vector(:,2,:)];
bottomRight = [position_vector(:,1,:)+position_vector(:,3,:),position_vector(:,2,:)];
topRight = [position_vector(:,1,:)+position_vector(:,3,:),position_vector(:,2,:)+position_vector(:,4,:)];
topLeft = [position_vector(:,1,:),position_vector(:,2,:)+position_vector(:,4,:)];
coords = [bottomLeft; bottomRight; topRight; topLeft];
rectX = coords(:,1,:); rectY = coords(:,2,:);
end
function oldHandleUnits = SyncUnits(handles, trgt_units)
% oldUnits is a n x 2 cell array where oldUnits{:,1} == handles
% and oldUnits{:,2} == "Units"
if nargout < 1
    warning('utilityFcns:reconcileUnits',...
        'old units not recorded, new units still set')
elseif nargout > 1
    error('Too many output arguments')
elseif nargin < 2 || ~ischar(trgt_units)
    error('Target units not specified')
end
[m,n] = size(handles);
handles = reshape(handles,m*n,1); % Turn handles into column vector
oldUnits = get(handles,'Units');

handles_cell = mat2cell(handles, ones(1,m*n), 1); % Turn handles into cell
oldHandleUnits = [handles_cell, oldUnits];

%Set new units last as a precautionary -- if exception arises at
%least we might recover the old units.
try
    set(handles, 'Units', trgt_units)
catch ME
    ME.message
    error('Unable to set units')
end

end
function UnSyncUnits(old_handle_units)
    arrayfun(@(k)set(old_handle_units{k,1},'Units',old_handle_units{k,2}),...
        1:1:size(old_handle_units,1))
end
function [adjacencyMatrix,varargout] = OverlappingAt(position_vector)
%Returns the adjacency matrix indicating the overlap of any of N rectangles
%as specified by the Nx4 array of position vectors, position_vector.  If
%called with two outputs, the second output is matrix where the (i,j)th
%element is the area of overlap between position_vector(i,:) and position_vector(j,:).
%

area = rectint(position_vector,position_vector);
density = nnz(area) / numel(area);
if density <= 0.4
    area = sparse(area);
end
area = area - diag(diag(area)); %removes self-reference
if nargout == 2
    varargout{1} = area;
end
adjacencyMatrix = (area ~= 0);
end


text:
function fd1d_display ( prefix )

%*****************************************************************************80
%
%% FD1D_DISPLAY displays data for a 1D finite difference function.
%
%  Discussion:
%
%    The function is represented as the sum of piecewise linear functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Usage:
%
%    fd1d_display ( 'prefix' )
%
%    where
%
%    'prefix_nodes.txt'    contains the point coordinates;
%    'prefix_values.txt'   contains the point values.
%
%  Parameters:
%
%    Input, string PREFIX, the prefix of the input files.
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read files defining a finite difference computatation,\n' );
  fprintf ( 1, '  involving a set of nodes, and values of one or more\n' );
  fprintf ( 1, '  functions at those nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot the functions in a MATLAB graphics window.\n' );
%
%  First argument is the file prefix.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FD1D_DISPLAY:\n' );
    prefix = input ( '  Enter the common prefix of the files (in ''quotes''!).\n' );
  end
%
%  Get the node coordinates.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  node_x = load ( node_filename );
%
%  Get the function data.
%
  value_filename = strcat ( prefix, '_values.txt' );
  node_value = load ( value_filename );
%
%  Make the plot.
%
  clf

  plot ( node_x, node_value )
  hold on
  plot ( node_x, node_value, '*' )

  grid

  title_string = s_escape_tex ( prefix );
  title ( title_string );
  xlabel ( 'X axis' );
  ylabel ( 'Function values' );

  hold off
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_DISPLAY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' | ...
         s1(s1_pos) == '_' | ...
         s1(s1_pos) == '^' | ...
         s1(s1_pos) == '{' | ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function new_x = monolist(x,dmax,dmin)
% MONOLIST Generate monomials
%
% y = MONOLIST(x,dmax,dmin)
%
% Returns the monomials [1 x(1) x(1)^2 ... x(1)^dmax(1) x(2) x(1)x(2)  etc...]
%
% >>sdpvar x y z
% >>sdisplay(monolist([x y z],4))
%
%  Input
%     x      : Vector with SDPVAR variables
%     dmax   : Integers > 0
%
%  See also POLYNOMIAL, DEGREE

% Flatten
x = reshape(x,1,length(x));
x_orig = x;

if nargin == 3
    if length(dmin)>1 || any(dmin > dmax) || ~isequal(dmin,fix(dmin))
        error('dmin has to be an integer scalar larger than dmax');
    end
elseif nargin == 2
    dmin = 0;
end

if (length(dmax)==1 | all(dmax(1)==dmax)) & islinear(x) & ~isa(x,'ncvar')
    dmax = dmax(1);
    % Powers in monomials
    powers = monpowers(length(x),dmax);

    powers = powers(find(sum(powers,2)>=dmin),:);

    % Use fast method for x^alpha
    if isequal(getbase(x),[zeros(length(x),1) eye(length(x))])
        new_x = recovermonoms(powers,x);
        return
    end

    % Monolist on dense vectors is currently extremely slow, but also
    % needed in some applications (stability analysis using SOS) For
    % performance issue, the code below is hard-coded for special cases
    % FIX : Urgent, find underlying indexing...

    % Vectorize quadratic and quadrtic case
    if dmax==2 & length(x)>1
        V=x.'*[1 x];
        ind=funkyindicies(length(x));
        new_x = [1 V(ind(:)).'].';
        return
    elseif (length(x)==4 & dmax==6)

         ind =[    1           2           3           4           5           6           7           8,

           9          10          11          12          13          14          15          16,

          17          18          19          20          21          22          23          24,

          25          26          27          28          29          30          31          32,

          33          34          49          50          51          52          86          53,

          54          55          89          56          57          91          58          92,

         126          59          60          61          95          62          63          97,

          64          98         132          65          66         100          67         101,

         135          68         102         136         170         185         186         187,

         188         222         256         189         190         191         225         259,

         192         193         227         261         194         228         262         296,

         330         364         195         196         197         231         265         198,

         199         233         267         200         234         268         302         336,

         370         201         202         236         270         203         237         271,

         305         339         373         204         238         272         306         340,

         374         408         442         476         510         525         526         527,

         528         562         596         630         529         530         531         565,

         599         633         532         533         567         601         635         534,

         568         602         636         670         704         738         772         806,

         840         535         536         537         571         605         639         538,

         539         573         607         641         540         574         608         642,

         676         710         744         778         812         846         541         542,

         576         610         644         543         577         611         645         679,

         713         747         781         815         849         544         578         612,

         646         680         714         748         782         816         850         884,

         918         952         986        1020        1054        1088        1122        1156,

        1190          0            0          0          0             0           0            0];
        ind = ind';
        ind = ind(find(ind));
        v=monolist(x,3);
        V=v(2:end)*v.';
        new_x = [1;V(ind(:))];
        return
    elseif dmax==4 & (1<=length(x)) & length(x)<=4 %& length(x)>1
        v=monolist(x,2);
        V=v(2:end)*v.';

        % Cone to generate indicies
        %p = sdpvar(n,1);
        %v = monolist(p,2);
        %V = v(2:end)*v';V=V(:);
        %m = monolist(p,4)
        %ind = [];
        %for i = 2:length(m)
        % ind = [ind min(find(~any(V-m(i))))];
        %end

        switch length(x)
            case 1
                new_x = [1; V([1 2 4 6]')];
                return
            case 2
                new_x = [1;V([1 2 3 4 5 8 9 10 15 18 19 20 25 30]')];
                return;
            case 3
                new_x=[1;V([1 2 3 4 5 6 7 8 9 13 14 15 24 16 17 26 18 27 36 40 41 42 51 60 43 44 53 62 45 54 63 72 81 90]')];
                return
            case 4
                new_x=[1;V([    1     2     3     4     5     6     7     8     9    10    11    12    13    14    19    20    21    35    22    23    37    24    38    52    25    26    40    27    41    55    28    42 56    70    75    76    77    91   105    78    79    93   107    80    94   108   122   136  150    81    82    96   110    83    97   111   125   139   153    84    98   112   126   140  154   168   182   196   210]')];
                return
            otherwise
        end
    end

    % Na, we have things like (c'x)^alpha
    % precalc x^p
    for i = 1:length(x)
        temp = x(i);
        precalc{i,1} = temp;
        for j = 2:1:dmax
            temp = temp*x(i);
            precalc{i,j} = temp;
        end
    end

    new_x = [];

    for i = 1:size(powers,1) % All monomials
        temp = 1;

        for j = 1:size(powers,2) % All variables
            if powers(i,j)>0
                temp = temp*precalc{j,powers(i,j)};
            end
        end
        new_x = [new_x temp];

    end

else

    dmax = dmax(:)*ones(1,length(x_orig));

    x = [1 x];

    % Lame loop to generate all combinations
    new_x = 1;
    for j = 1:1:max(dmax)
        temp = [];
        for i = 1:length(x)
            temp = [temp x(i)*new_x];
        end
        new_x = temp;
        new_x = fliplr(unique(new_x));
        new_degrees = degree(new_x,x(2:end));
        remv = [];
        for i = 1:length(dmax);
            if new_degrees(i)>=dmax(i)
                x = recover(setdiff(getvariables(x),getvariables(x_orig(i))));
                x = [1;x(:)];
                remv = [remv i];
            end
        end
        dmax = dmax(setdiff(1:length(dmax),remv));
    end
end
new_x = reshape(new_x(:),length(new_x),1);
if dmin > 0
    for i = 1:length(new_x)
        if sum(powers(i,:)) < dmin
            keep(i) = 0;
        else
            keep(i) = 1;
        end
    end
    if any(keep==0)
        new_x = new_x(find(keep));
    end
end

function ind = funkyindicies(n)

M=reshape(1:n*(n+1),n,n+1);
ind = M(:,1)';
for i = 1:n
    ind = [ind M(i,2:i+1)];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function circle()
    %   This subroutine "circle"  selects the Ni closest earthquakes
    %   around a interactively selected point.  Resets ZG.newcat and ZG.newt2
    %   Operates on "primeCatalog".
    %
    % axis: h1
    % plots to: plos1 as xk
    % inCatalog: a
    % outCatalog: newt2, newcat
    % mouse controlled
    % closest events
    %
    % turned into function by Celso G Reyes 2017
    %
%
%
    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();

    ShapeGeneral.clearplot(); % was tag plos1

    % interactively get the circle of interest
    shape=ShapeCircle();

    [ZG.newt2, max_km] = selectCircle(newa, shape.toStruct());

    fprintf('Radius of selected Circle: %s km\n', num2str(max_km) );

    % plot Ni clostest events on map as 'x':

    set(gca,'NextPlot','add')
    plot(ZG.newt2.Longitude,ZG.newt2.Latitude,'xk','Tag','plos1');

    %
    ZG.newcat = ZG.newt2;                   % resets ZG.newcat and ZG.newt2

    timeplot(2)

end


text:
function s6_resave_annot_after_remaing()
%image chnaged
% stb = load('annot_test.mat');
% base_annot = stb.videos;
% clear stb;clear all
close all
stb = load('annotV5.mat');
%stb = load('annot.mat');
base_annot = stb.videos;
clear stb;

stp = load('annot_full_Philippe.mat'); % train+test
phill_annot = stp.videos;
clear stp;

stp = load('annot_apt.mat'); % train+test
apt_annot = stp.videos;
clear stp;

baseList = {};
for i=1:length(base_annot)
    baseList{i,1} = base_annot(i).name; % store all the cor vid in a cellarray
end

philList = {};
for i=1:length(phill_annot)
    philList{i,1} = phill_annot(i).name; % store all the cor vid in a cellarray
end

aptList = {};
for i=1:length(apt_annot)
    aptList{i,1} = apt_annot(i).name; % store all the cor vid in a cellarray
end

images_dir = '/mnt/sun-alpha/datasets/UCF101/images';
save_path = './saved/';
%for v = 403 % : length(base_annot)
% rvid = randperm(length(base_annot));

% FigHandle = figure(1);
% FigHandle = figure('Position', [100, 100, 1049, 895]);
% hold on
count = 0;
load('lastdata.mat')
% merged_annot = struct();
load('merged_annot_temporal');
load('checkpoint_improved_temporal.mat')
v = 1;vc =0;
while v <= length(baseList)

    videoName = baseList{v,1};

    base_tubes = base_annot(v).tubes;

    index = find(strcmp(philList,videoName));
    phil_tubes = phill_annot(index).tubes;

    index = find(strcmp(aptList,videoName));
    apt_tubes = apt_annot(index).tubes;

    class = base_annot(v).tubes(1).class;
    num_imgs = base_annot(v).num_images;

    merged_annot(v).num_imgs = num_imgs;
    merged_annot(v).name = base_annot(v).name;

    if ~isgood(v) % && strcmp(inputs{v},'t')
        vc = vc +1;
        newcell = tdata{v};
        tline = newcell{1};
        stl = strsplit(tline,' ');

        if length(stl{1})>1
            inputs{v} = 'ra';
            isgood(v) = 0;
            merged_annot(v).tubes = base_annot(v).tubes;
        else
            merged_annot(v).tubes = extract_bounds(newcell, base_tubes, phil_tubes, apt_tubes, num_imgs,class);
            inputs{v} = 'y';
            isgood(v) = 1;
        end

    end
    v=v+1;

end
% vc
fprintf('%d\n',sum(isgood));
save('checkpoint_improved_after_remaining.mat','isgood','inputs');
save('merged_annot_remianing.mat','merged_annot')

function tubes = extract_bounds(newcell, base_tubes, phil_tubes, apt_tubes, num_imgs,class)
tubes = struct();
alltubes = {base_tubes,phil_tubes,apt_tubes};
for t = 1 : length(newcell)
    tline = newcell{t};
    stl = strsplit(tline,' ');

    if length(stl)<6
        [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,0);
    else
        [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,1);
    end
    tubes(t).sf = sf;
    tubes(t).ef = ef;
    tubes(t).class = class;
    tubes(t).boxes = boxes;
end

function [sf,ef,boxes] = gettube(alltubes,stl,num_imgs,isone)
tube  = [];
annotType = {'b','p','a'};
ind = find(strcmp(annotType,stl{1}));

annot = alltubes{ind};

tid = str2num(stl{2});

sf = max(str2num(stl{3}),annot(tid).sf);

ef = min(str2num(stl{4}),annot(tid).ef);
st = sf-annot(tid).sf+1;
et = ef-sf+st; ct = 1;
boxes = annot(tid).boxes(st:et,:);
if isone
    ind = find(strcmp(annotType,stl{5}));
    annot = alltubes{ind};
    tid = str2num(stl{6});
    sfa = str2num(stl{7});
    if (sfa-ef)>1
        for kk = 1:(sfa-ef-1)
            boxes = [boxes;boxes(end,:)];
        end
    end
    ef = str2num(stl{8});
    st = sfa-annot(tid).sf+1;
    et = ef-sfa+st; ct = 1;
    boxes = [boxes;annot(tid).boxes(st:et,:)];
end

author: gurkirt


text:
function tests = test_spm_dcm_loo
% Unit Tests for test_spm_dcm_peb
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_dcm_loo.m 7102 2017-06-08 12:56:06Z peter $

tests = functiontests(localfunctions);

% -------------------------------------------------------------------------
function test_loo_group(testCase)
% Tests LOO cross-validation in the presence of a binary group effect

data_path = get_data_path();

% Reduce number of subjects to increase performance
s = [1:6 24:30];

% Load first level DCMs
GCM = load(fullfile(data_path,'models','GCM_simulated.mat'));
GCM = GCM.GCM(s,:);

% Prepare group level design matrix
X  = load(fullfile(data_path,'design_matrix.mat'));
X  = X.X(s,:);
ns = size(X,1);
X  = [ones(ns,1) X];

% Run
M.X = X;
[qE,qC,Q] = spm_dcm_loo(GCM(:,2),M,{'B'});

% Test that classical p-value is significant
[T,df] = spm_ancova(M.X(:,1:2),[],qE(:),[0;1]);
p = 1 - spm_Tcdf(T,df(2));
testCase.assertTrue(p < 0.05);

% -------------------------------------------------------------------------
function test_loo_group_null(testCase)
% Tests LOO cross-validation in the absence of a binary group effect

data_path = get_data_path();

% Reduce number of subjects to increase performance
s = [1:6 24:30];

% Load first level DCMs
GCM = load(fullfile(data_path,'models','GCM_simulated.mat'));
GCM = GCM.GCM(s,:);

% Prepare group level design matrix
X  = load(fullfile(data_path,'design_matrix.mat'));
X  = X.X(s,:);
ns = size(X,1);
X  = [ones(ns,1) X];

% Run
M.X = X;
[qE,qC,Q] = spm_dcm_loo(GCM(:,2),M,{'A'});

% Test that classical p-value is significant
[T,df] = spm_ancova(M.X(:,1:2),[],qE(:),[0;1]);
p = 1 - spm_Tcdf(T,df(2));
testCase.assertTrue(p > 0.05);

% -------------------------------------------------------------------------
function test_loo_continuous_null(testCase)
% Tests LOO cross-validation in the absence of a continuous group effect

data_path = get_data_path();

% Reduce number of subjects to increase performance
s = [1:6 24:30];

% Load first level DCMs
GCM = load(fullfile(data_path,'models','GCM_simulated.mat'));
GCM = GCM.GCM(s,:);

% Prepare group level design matrix
X  = load(fullfile(data_path,'design_matrix.mat'));
X  = X.X(s,:);
ns = size(X,1);
X  = [ones(ns,1) X];

% Re-order to test random age covariate
X = X(:,[1 3 2]);

% Run
M.X = X;
[qE,qC,Q] = spm_dcm_loo(GCM(:,2),M,{'B'});

% Test that classical p-value is significant
[T,df] = spm_ancova(M.X(:,1:2),[],qE(:),[0;1]);
p = 1 - spm_Tcdf(T,df(2));
testCase.assertTrue(p >= 0.05);

% -------------------------------------------------------------------------
function data_path = get_data_path()

data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'fMRI', 'simulated_2region');

author: spm


text:
% Rotate position samples
%
% This function rotates provided position samples to specified angle.
%
%  USAGE
%   rpos = rotatePositions(pos, degAngle)
%   pos         Matrix with position samples.
%   degAngle    Rotation angle in degrees.
%   rpos        Rotated position samples. Same dimension as pos.
%
function rpos = rotatePositions(pos, degAngle)
    rpos(:, 2) = pos(:, 2)*cosd(degAngle) - pos(:, 3)*sind(degAngle);
    rpos(:, 3) = pos(:, 2)*sind(degAngle) + pos(:, 3)*cosd(degAngle);
    if size(pos, 2) > 3
        rpos(:, 4) = pos(:, 4)*cosd(degAngle) - pos(:, 5)*sind(degAngle);
        rpos(:, 5) = pos(:, 4)*sind(degAngle) + pos(:, 5)*cosd(degAngle);
    end
end

author: kavli-ntnu


text:
% ==============================================================================
%
%    Software License Agreement (BSD License)
%    Copyright (c) 2019
%    (www.aimlab.wpi.edu)
%
%    All rights reserved.
%
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions
%    are met:
%
%    * Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
%
%    * Redistributions in binary form must reproduce the above
%    copyright notice, this list of conditions and the following
%    disclaimer in the documentation and/or other materials provided
%    with the distribution.
%
%    * Neither the name of authors nor the names of its contributors may
%    be used to endorse or promote products derived from this software
%    without specific prior written permission.
%
%    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
%    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
%    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
%    FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
%    COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
%    INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
%    BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
%    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
%    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
%    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
%    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%    POSSIBILITY OF SUCH DAMAGE.
%
%    \author:    <http://www.aimlab.wpi.edu>
%    \author:    <amunawar@wpi.edu>
%    \author:    Adnan Munawar
%    \version:   0.1$
% ==============================================================================

function [nIterations,sizePath,run_time] =  RRTconnect3D(dim,segmentLength,random_world,show_output)
% dim = 2;
% segmentLength = 5;
% random_world = 0;
% standard length of path segments
if dim ==2
    start_cord = [5,5];
    goal_cord = [95,95];

else

    start_cord = [5,5,5];
    goal_cord = [95,95,95];
end

% create random world
Size = 100;
NumObstacles = 100;

if random_world ==1
    world = createWorld(NumObstacles,ones(1,dim)*Size,zeros(1,dim),dim);
else
    [world NumObstacles] = createKnownWorld(ones(1,dim)*Size,[0;0;0],dim);
end
% randomly select start and end nodes
%start_node = generateRandomNode(world,dim)
%end_node   = generateRandomNode(world,dim)
start_node = [start_cord,0,0,0];
end_node = [goal_cord,0,0,0];
% establish tree starting with the start node
tree = start_node;

a = clock;

% check to see if start_node connects directly to end_node
if ( (norm(start_node(1:dim)-end_node(1:dim))<segmentLength )...
        &&(collision(start_node,end_node,world,dim)==0) )
    path = [start_node; end_node];
else
    nIterations = 0;
    numPaths = 0;
    flag = 0;
    while numPaths<1,
        [tree,flag] = extendTree(tree,end_node,segmentLength,world,dim);
        numPaths = numPaths + flag;
        nIterations = nIterations+1;
    end
end

path = findMinimumPath(tree,end_node,dim);
sizePath = size(path,1);

b = clock;

run_time = 3600*(b(4)-a(4)) + 60 * (b(5)-a(5)) + (b(6) - a(6));

if show_output == 1
    % find path with minimum cost to end_node
    figure;
    plotExpandedTree(world,tree,dim);
    plotWorld(world,path,dim);
    % figure(2);
    % plotWorld(world,path,dim);
    %plotExpandedTree(world,tree,dim);
end
end

function world = createWorld(NumObstacles, endcorner, origincorner,dim)

if dim == 2

    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2))
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = min(endcorner(1)- origincorner(1), endcorner(2)-origincorner(2));
        maxRadius = 5*maxRadius/NumObstacles/2;
        for i=1:NumObstacles,
            % randomly pick radius
            world.radius(i) = maxRadius*rand;
            % randomly pick center of obstacles
            cx = origincorner(1) + world.radius(i)...
                + (endcorner(1)-origincorner(1)-2*world.radius(i))*rand;
            cy = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            world.cx(i) = cx;
            world.cy(i) = cy;
        end
    end

elseif dim ==3;
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) || (endcorner(2) <= origincorner(2)) || (endcorner(3) <= origincorner(3))
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        bounds = [endcorner(1)- origincorner(1), endcorner(2)-origincorner(2), endcorner(3)-origincorner(3)];
        maxRadius = min(bounds);
        maxRadius = 5*maxRadius/NumObstacles;
        for i=1:NumObstacles,
            % randomly pick radius
            world.radius(i) = maxRadius*rand;
            % randomly pick center of obstacles
            cx = origincorner(1) + world.radius(i)...
                + (endcorner(1)-origincorner(1)-2*world.radius(i))*rand;
            cy = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            cz = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            world.cx(i) = cx;
            world.cy(i) = cy;
            world.cz(i) = cz;
        end
    end
end
end

function [world NumObstacles] = createKnownWorld(endcorner, origincorner,dim)
NumObstacles = 5;
if dim == 2
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2)),
        disp('Not valid corner specifications!')
        world=[];
        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = 10;

        world.radius(1) = maxRadius;
        cx = 50;
        cy = 50;
        world.cx(1) = cx;
        world.cy(1) = cy;

        world.radius(2) = maxRadius;
        cx = 75;
        cy = 25;
        world.cx(2) = cx;
        world.cy(2) = cy;

        world.radius(3) = maxRadius;
        cx = 25;
        cy = 75;
        world.cx(3) = cx;
        world.cy(3) = cy;

        world.radius(4) = maxRadius;
        cx = 25;
        cy = 25;
        world.cx(4) = cx;
        world.cy(4) = cy;

        world.radius(5) = maxRadius;
        cx = 75;
        cy = 75;
        world.cx(5) = cx;
        world.cy(5) = cy;
    end

elseif dim == 3

    NumObstacles = 9;
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2)) | (endcorner(3) <= origincorner(3)),
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = 10;

        world.radius(1) = maxRadius;
        cx = 50;
        cy = 50;
        cz = 50;
        world.cx(1) = cx;
        world.cy(1) = cy;
        world.cz(1) = cz;

        world.radius(2) = maxRadius;
        cx = 25;
        cy = 25;
        cz = 25;
        world.cx(2) = cx;
        world.cy(2) = cy;
        world.cz(2) = cz;

        world.radius(3) = maxRadius;
        cx = 75;
        cy = 75;
        cz = 75;
        world.cx(3) = cx;
        world.cy(3) = cy;
        world.cz(3) = cz;

        world.radius(4) = maxRadius;
        cx = 25;
        cy = 25;
        cz = 75;
        world.cx(4) = cx;
        world.cy(4) = cy;
        world.cz(4) = cz;

        world.radius(5) = maxRadius;
        cx = 75;
        cy = 75;
        cz = 25;
        world.cx(5) = cx;
        world.cy(5) = cy;
        world.cz(5) = cz;

        world.radius(6) = maxRadius;
        cx = 25;
        cy = 75;
        cz = 25;
        world.cx(6) = cx;
        world.cy(6) = cy;
        world.cz(6) = cz;

        world.radius(7) = maxRadius;
        cx = 75;
        cy = 25;
        cz = 25;
        world.cx(7) = cx;
        world.cy(7) = cy;
        world.cz(7) = cz;

        world.radius(8) = maxRadius;
        cx = 75;
        cy = 25;
        cz = 75;
        world.cx(8) = cx;
        world.cy(8) = cy;
        world.cz(8) = cz;

        world.radius(9) = maxRadius;
        cx = 25;
        cy = 75;
        cz = 75;
        world.cx(9) = cx;
        world.cy(9) = cy;
        world.cz(9) = cz;
    end
end
end

function node=generateRandomNode(world,dim)

if dim ==2;
    % randomly pick configuration
    px       = (world.endcorner(1)-world.origincorner(1))*rand;
    py       = (world.endcorner(2)-world.origincorner(2))*rand;

    chi      = 0;
    cost     = 0;
    node     = [px, py, chi, cost, 0];

    % check collision with obstacle
    while collision(node, node, world,dim),
        px       = (world.endcorner(1)-world.origincorner(1))*rand;
        py       = (world.endcorner(2)-world.origincorner(2))*rand;

        chi      = 0;
        cost     = 0;
        node     = [px, py, chi, cost, 0];
    end

elseif dim ==3;
    % randomly pick configuration
    px       = (world.endcorner(1)-world.origincorner(1))*rand;
    py       = (world.endcorner(2)-world.origincorner(2))*rand;
    pz       = (world.endcorner(3)-world.origincorner(3))*rand;

    chi      = 0;
    cost     = 0;
    node     = [px, py, pz, chi, cost, 0];

    % check collision with obstacle
    while collision(node, node, world,dim),
        px       = (world.endcorner(1)-world.origincorner(1))*rand;
        py       = (world.endcorner(2)-world.origincorner(2))*rand;
        pz       = (world.endcorner(3)-world.origincorner(3))*rand;

        chi      = 0;
        cost     = 0;
        node     = [px, py, pz, chi, cost, 0];
    end

end

end

function collision_flag = collision(node, parent, world,dim)

collision_flag = 0;

for i=1:dim
    if (node(i)>world.endcorner(i))|(node(i)<world.origincorner(i))
        collision_flag = 1;
    end
end

if collision_flag == 0 && dim ==2
    for sigma = 0:.2:1,
        p = sigma*node(1:dim) + (1-sigma)*parent(1:dim);
        % check each obstacle
        for i=1:world.NumObstacles,
            if (norm([p(1);p(2)]-[world.cx(i); world.cy(i)])<=1*world.radius(i)),
                collision_flag = 1;
                break;
            end
        end
    end

elseif collision_flag == 0 && dim ==3
    for sigma = 0:.2:1,
        p = sigma*node(1:dim) + (1-sigma)*parent(1:dim);
        % check each obstacle
        for i=1:world.NumObstacles,
            if (norm([p(1);p(2);p(3)]-[world.cx(i); world.cy(i); world.cz(i)])<=1*world.radius(i)),
                collision_flag = 1;
                break;
            end
        end
    end
end
end

function collision_flag = is_point_valid(point, world,dim)

collision_flag = 0;

for i=1:dim
    if (point(i)>world.endcorner(i))||(point(i)<world.origincorner(i))
        collision_flag = 1;
    end
end

if collision_flag == 0 && dim ==2

    p = point(1:dim);
    % check each obstacle
    for i=1:world.NumObstacles,
        if (norm([p(1);p(2)]-[world.cx(i); world.cy(i)])<=1*world.radius(i)),
            collision_flag = 1;
            break;
        end
    end

elseif collision_flag == 0 && dim ==3
    p = point(1:dim);
    % check each obstacle
    for i=1:world.NumObstacles,
        if (norm([p(1);p(2);p(3)]-[world.cx(i); world.cy(i); world.cz(i)])<=1*world.radius(i)),
            collision_flag = 1;
            break;
        end
    end
end
end

function flag = canEndConnectToTree(tree,end_node,minDist,world,dim)
flag = 0;
% check only last node added to tree since others have been checked
if ( (norm(tree(end,1:dim)-end_node(1:dim))<minDist)...
        & (collision(tree(end,1:dim), end_node(1:dim), world,dim)==0) ),
    flag = 1;
end

end

function [new_tree,flag] = extendTree(tree,end_node,segmentLength,world,dim)

flag = 0;
% select a random point
randomPoint = zeros(1,dim);
for i=1:dim
    randomPoint(1,i) = (world.endcorner(i)-world.origincorner(i))*rand;
end

% find leaf on node that is closest to randomPoint
tmp = tree(:,1:dim)-ones(size(tree,1),1)*randomPoint;
sqrd_dist = sqr_eucl_dist(tmp,dim);
[min_dist,idx] = min(sqrd_dist);

min_parent_idx = idx;

new_point = tree(idx,1:dim);
new_node = tree(idx,:);

pflag = 0;

while norm(new_point-randomPoint)>0 && pflag==0

    if norm(new_point-randomPoint)<segmentLength

        pflag = collision(randomPoint,tree(min_parent_idx,:),world,dim);
        if pflag == 0

            new_point = randomPoint;
            min_cost = cost_np(tree(min_parent_idx,:),new_point,dim);
            new_node = [new_point,0,min_cost,min_parent_idx];
            tree = [tree;new_node];
            pflag = 1;
            goal_flag = is_goal(new_node,end_node,segmentLength,world,dim);

            if goal_flag == 1;
                tree(end,dim+1)=1;
                flag = 1;
            end
        end

    else

        new_point = (randomPoint-tree(min_parent_idx,1:dim));
        new_point = tree(min_parent_idx,1:dim)+(new_point/norm(new_point))*segmentLength;

        min_cost  = cost_np(tree(min_parent_idx,:),new_point,dim);
        new_node  = [new_point, 0, min_cost, min_parent_idx];

        pflag = collision(new_node,tree(min_parent_idx,:),world,dim);

        if pflag == 0
            tree = [tree ; new_node];
            min_parent_idx = size(tree,1);

            goal_flag = is_goal(new_node,end_node,segmentLength,world,dim);

            if goal_flag == 1;
                tree(end,dim+1)=1;  % mark node as connecting to end.
                pflag = 1;
                flag = 1;
            end
        end
    end
end

new_tree = tree;

end

function goal_flag = is_goal(node,end_node,segmentLength,world,dim)

goal_flag = 0;
if (norm(node(1:dim)-end_node(1:dim))<segmentLength )...
        && (collision(node,end_node,world,dim)==0)
    goal_flag = 1;
end

end

function e_dist = sqr_eucl_dist(array,dim)

sqr_e_dist = zeros(size(array,1),dim);
for i=1:dim

    sqr_e_dist(:,i) = array(:,i).*array(:,i);

end
e_dist = zeros(size(array,1),1);
for i=1:dim

    e_dist = e_dist+sqr_e_dist(:,i);

end

end

%calculate the cost from a node to a point
function [cost] = cost_np(from_node,to_point,dim)

diff = from_node(:,1:dim) - to_point;
eucl_dist = norm(diff);
cost = from_node(:,dim+2) + eucl_dist;

end

%calculate the cost from a node to a node
function [cost] = cost_nn(from_node,to_node,dim)

diff = from_node(:,1:dim) - to_node(:,1:dim);
eucl_dist = norm(diff);
cost = from_node(:,dim+2) + eucl_dist;

end

function [cost] = line_cost(from_node,to_point,dim)
diff = from_node(:,1:dim) - to_point;
cost = norm(diff);
end

function path = findMinimumPath(tree,end_node,dim)

% find nodes that connect to end_node
connectingNodes = [];
for i=1:size(tree,1),
    if tree(i,dim+1)==1,
        connectingNodes = [connectingNodes ; tree(i,:)];
    end
end

% find minimum cost last node
[tmp,idx] = min(connectingNodes(:,dim+2));

% construct lowest cost path
path = [connectingNodes(idx,:); end_node];
parent_node = connectingNodes(idx,dim+3);
while parent_node>1,
    parent_node = tree(parent_node,dim+3);
    path = [tree(parent_node,:); path];
end

end

function plotExpandedTree(world,tree,dim)
ind = size(tree,1);
while ind>0
    size(tree);
    branch = [];
    node = tree(ind,:);
    branch = [ branch ; node ];
    parent_node = node(dim+3);
    while parent_node > 1
        cur_parent = parent_node;
        branch = [branch; tree(parent_node,:)];
        parent_node = tree(parent_node,dim+3);
    end
    ind = ind - 1;

    if dim == 2
        X = branch(:,1);
        Y = branch(:,2);

        p = plot(X,Y);
        set(p,'Color','r','LineWidth',0.5,'Marker','.','MarkerEdgeColor','g');
        hold on;

    elseif dim == 3
        X = branch(:,1);
        Y = branch(:,2);
        Z = branch(:,3);

        p = plot3(X,Y,Z);
        set(p,'Color','r','LineWidth',0.5,'Marker','.','MarkerEdgeColor','g');
        hold on;
    end
end
end

function plotWorld(world,path,dim)
% the first element is the north coordinate
% the second element is the south coordinate
if dim ==2

    N = 10;
    th = 0:2*pi/N:2*pi;
    axis([world.origincorner(1),world.endcorner(1),...
        world.origincorner(2), world.endcorner(2)]);
    hold on

    for i=1:world.NumObstacles,
        X = world.radius(i)*sin(th) + world.cx(i);
        Y = world.radius(i)*cos(th) + world.cy(i);
        fill(X,Y,'blue');
    end

    X = path(:,1);
    Y = path(:,2);
    p = plot(X,Y);

elseif dim ==3
    axis([world.origincorner(1),world.endcorner(1),...
        world.origincorner(2), world.endcorner(2),...
        world.origincorner(3), world.endcorner(3)]);
    hold on

    for i=1:world.NumObstacles,
        [X Y Z] = sphere(10);
        X = (X*world.radius(i));
        Y = (Y*world.radius(i));
        Z = (Z*world.radius(i));
        surf(X+world.cx(i),Y+world.cy(i),Z+world.cz(i));
        colormap([0.5 0.2 0.3]);
    end

    X = path(:,1);
    Y = path(:,2);
    Z = path(:,3);
    p = plot3(X,Y,Z);
end
set(p,'Color','black','LineWidth',3)
xlabel('X axis');
ylabel('Y axis');
zlabel('Z axis');
title('RRT Connect Algorithm');
end

author: adnanmunawar

repo: matlab-rrt-variants


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LMS Algorithm
% Author : Arsal Butt
% Contact information :
% School of Electrical & Electronic Engineering,
% University of Bradford,
% Bradford, BD7 1DP,
% United Kingdom .
% email : a.butt7@bradford.ac.uk
% Date    : 03-8-2010
% Version : 1.0.2
% Reference : Statistical and Adaptive Signal Processing By Dimitris G. Manolakis, Vinay K. Ingle & Stephen M. Kogon.
% Note : The author doesn't take any responsibility for any harm caused by
% the use of this code;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use N iterations of an adapting filter system to identify an unknown
% Mth-order FIR filter.
function lmsalgorithm = lms()
sysorder = input('Enter the System Order? ')
N = input('Enter the number of Itterations? ')
if ((sysorder >= 2) && (N >= 2))
    x  = randn(N,1);     % Input to the filter
b  = fir1(sysorder-1,0.5);  % FIR system to be identified
n  = 0.1*randn(N,1); % Uncorrelated noise signal
d  = filter(b,1,x)+n;  % Desired signal = output of H + Uncorrelated noise signal
w = zeros (sysorder, 1) ; % Initially filter weights are zero
for n = sysorder : N
	u = x(n:-1:n-sysorder+1) ;
    y(n)= w' * u; % output of the adaptive filter
    e(n) = d(n) - y(n) ; % error signal = desired signal - adaptive filter output
    mu=0.008;
    w = w + mu * u * e(n) ; % filter weights update
end
hold on
plot(d,'g')
plot(y,'r');
semilogy((abs(e)),'m') ;
title('System output') ;
xlabel('Samples')
ylabel('True and estimated output')
legend('Desired','Output','Error')
axis([0 N -2 2.5])
figure
plot(b', 'k+')
hold on
plot(w, 'r*')
legend('Actual weights','Estimated weights')
title('Comparison of the actual weights and the estimated weights') ;
else ('System Order and numbers of Itterations should be greater than 1')
end

end


text:
% [PYR, INDICES, STEERMTX, HARMONICS] = buildSFpyr(IM, HEIGHT, ORDER, TWIDTH)
%
% Construct a steerable pyramid on matrix IM, in the Fourier domain.
% This is similar to buildSpyr, except that:
%
%    + Reconstruction is exact (within floating point errors)
%    + It can produce any number of orientation bands.
%    - Typically slower, especially for non-power-of-two sizes.
%    - Boundary-handling is circular.
%
% HEIGHT (optional) specifies the number of pyramid levels to build. Default
% is maxPyrHt(size(IM),size(FILT));
%
% The squared radial functions tile the Fourier plane, with a raised-cosine
% falloff.  Angular functions are cos(theta-k\pi/(K+1))^K, where K is
% the ORDER (one less than the number of orientation bands, default= 3).
%
% TWIDTH is the width of the transition region of the radial lowpass
% function, in octaves (default = 1, which gives a raised cosine for
% the bandpass filters).
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.  This is compatible with the MatLab Wavelet toolbox.
% See the function STEER for a description of STEERMTX and HARMONICS.

% Eero Simoncelli, 5/97.
% See http://www.cis.upenn.edu/~eero/steerpyr.html for more
% information about the Steerable Pyramid image decomposition.

function [pyr,pind,steermtx,harmonics] = buildSFpyr(im, ht, order, twidth)

%-----------------------------------------------------------------
%% DEFAULTS:

max_ht = floor(log2(min(size(im)))+2);

if (exist('ht') ~= 1)
  ht = max_ht;
else
  if (ht > max_ht)
    error(sprintf('Cannot build pyramid higher than %d levels.',max_ht));
  end
end

if (exist('order') ~= 1)
  order = 3;
elseif ((order > 15)  | (order < 0))
  fprintf(1,'Warning: ORDER must be an integer in the range [0,15]. Truncating.\n');
  order = min(max(order,0),15);
else
  order = round(order);
end
nbands = order+1;

if (exist('twidth') ~= 1)
  twidth = 1;
elseif (twidth <= 0)
  fprintf(1,'Warning: TWIDTH must be positive.  Setting to 1.\n');
  twidth = 1;
end

%-----------------------------------------------------------------
%% Steering stuff:

if (mod((nbands),2) == 0)
  harmonics = [0:(nbands/2)-1]'*2 + 1;
else
  harmonics = [0:(nbands-1)/2]'*2;
end

steermtx = steer2HarmMtx(harmonics, pi*[0:nbands-1]/nbands, 'even');

%-----------------------------------------------------------------

dims = size(im);
ctr = ceil((dims+0.5)/2);

[xramp,yramp] = meshgrid( ([1:dims(2)]-ctr(2))./(dims(2)/2), ...
    ([1:dims(1)]-ctr(1))./(dims(1)/2) );
angle = atan2(yramp,xramp);
log_rad = sqrt(xramp.^2 + yramp.^2);
log_rad(ctr(1),ctr(2)) =  log_rad(ctr(1),ctr(2)-1);
log_rad  = log2(log_rad);

%% Radial transition function (a raised cosine in log-frequency):
[Xrcos,Yrcos] = rcosFn(twidth,(-twidth/2),[0 1]);
Yrcos = sqrt(Yrcos);

YIrcos = sqrt(1.0 - Yrcos.^2);
lo0mask = pointOp(log_rad, YIrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);
imdft = fftshift(fft2(im));
lo0dft =  imdft .* lo0mask;

[pyr,pind] = buildSFpyrLevs(lo0dft, log_rad, Xrcos, Yrcos, angle, ht, nbands);

hi0mask = pointOp(log_rad, Yrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);
hi0dft =  imdft .* hi0mask;
hi0 = ifft2(ifftshift(hi0dft));

pyr = [real(hi0(:)) ; pyr];
pind = [size(hi0); pind];

author: ricedsp


text:
function rec=iEWT2D_Ridgelet(ewt,mfb)

%==========================================================================
% function rec=iEWT2D_Ridgelet(ewt,mfb)
%
% This function performs the inverse 2D Empirical Ridgelet of ewt
% accordingly to the filter bank mfb .
%
% TO RUN THIS FUNCTION YOU NEED TO HAVE THE MATLAB POLARLAB TOOLBOX OF
% MICHAEL ELAD: http://www.cs.technion.ac.il/~elad/Various/PolarLab.zip
%
% Inputs:
%   -ewt: cell containing the 2D EWT components
%   -mfb: filter bank used during the EWT
%
% Output:
%   -rec: reconstructed image
%
% Author: Jerome Gilles - Giang Tran
% Institution: UCLA - Department of Mathematics
% Year: 2013
% Version: 1.0
%==========================================================================

PseudoFFT=zeros(size(ewt{1}));

% We perform the reconstruction of the Pseudo-Polar FFT domain by performing
% the adjoint operator of the EWT1D
for k=1:length(ewt)
    PseudoFFT=PseudoFFT+fftshift(fft(ewt{k}).*mfb{k},1);
end

rec=real(IPPFFT(PseudoFFT,1e-6));


text:
function [inputImages, inputSignals, valid, imageSizes, imgFeatures] = filterImages(inputImages, inputSignals, varargin)
	% Filters large and small objects in an set of images, returns filtered matricies along with vector with decisions and sizes.
	% Biafra Ahanonu
	% 2013.10.31
	% based on SpikeE code
	% inputs
		%
	% outputs
		%

	[inputImages, inputSignals, valid, imageSizes, imgFeatures] = ciapkg.image.filterImages(inputImages, inputSignals,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah