text:
classdef CoreReporting < CoreReportingInterface
    % CoreReporting. Provides error, message and progress reporting.
    %
    %     CoreReporting. Implementation of CoreReportingInterface, which is used by
    %     CoreMat and other libraries for progress and error/message reporting. This
    %     implementation displays warnings and messages on the command window,
    %     and uses Matlab's error() command to process errors. Logging
    %     information is writen to a log file. An object which implements
    %     the CoreProgressInterface can be passed in for progress reporting.
    %
    %     Usage
    %     -----
    %
    %     You should create a single CoreReporting object and pass it into all the
    %     CoreMat routines you use in order to provide error, warning,
    %     message and progress reporting during execution of routines.
    %
    %     If you are not writing a gui application but would like a standard
    %     pop-up progress dialog to appear while waiting for plugins to execute,
    %     consider creating a CoreReportingDefault object instead. Use CoreReporting
    %     if you want to specify your own progress dialog, or specify a gui
    %     viewing panel, or if you want no progress dialog at all.
    %
    %         reporting = CoreReporting(progress_dialog, viewing_panel);
    %
    %             progress_dialog - a CoreProgressInterface object such as
    %                 a CoreProgressDialog for displaying a progress bar. You can omit this argument
    %                 or replace it with [] if you are writing scripts to run
    %                 in the background and do not want progress dialogs popping
    %                 up. Otherwise, you should create a CoreProgressDialog
    %                 or else implement your own progress class if you want
    %                 custom behaviour.
    %
    %     See CoreReportingIntertface.m for details of the methods this class
    %     implements.
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Constant)
        CancelErrorId = 'CoreReporting:UserCancel'
    end

    properties
        ProgressDialog  % Handle to a CoreProgressInterface object
    end

    properties (Access = private)
        LogFileName     % Full path to log file

        % Stack for nested progress reporting
        ProgressStack
        CurrentProgressStackItem
        ParentProgressStackItem
        VerboseMode
    end

    methods
        function obj = CoreReporting(progress_dialog, verbose_mode, log_file_name)
            if nargin > 0
                obj.ProgressDialog = progress_dialog;
            end
            if nargin > 1
                obj.VerboseMode = verbose_mode;
            else
                obj.VerboseMode = false;
            end
            if nargin > 2
                obj.LogFileName = log_file_name;
            else
                obj.LogFileName = fullfile(CoreDiskUtilities.GetUserDirectory, 'corereporting.log');
            end

            obj.ClearProgressStack;
        end

        function Log(obj, message)
            [calling_function, ~] = CoreErrorUtilities.GetCallingFunction(2);

            obj.AppendToLogFile([calling_function ': ' message]);
        end

        function LogVerbose(obj, message)
            if obj.VerboseMode
                [calling_function, ~] = CoreErrorUtilities.GetCallingFunction(2);
                obj.AppendToLogFile([calling_function ': ' message]);
            end
        end

        function ShowMessage(obj, identifier, message)
            [calling_function, ~] = CoreErrorUtilities.GetCallingFunction(2);
            if isempty(calling_function)
                calling_function = 'Command Window';
            end
            disp(message);
            obj.AppendToLogFile([calling_function ': ' identifier ':' message]);
        end

        function ShowMessageFromException(obj, identifier, message, ex)
            [calling_function, ~] = CoreErrorUtilities.GetCallingFunction(2);
            if isempty(calling_function)
                calling_function = 'Command Window';
            end
            message = [calling_function ': ' identifier ':' message '; Exception message:' ex.message ' : ' ex.stack(1).name];
            [summary, stack_text] = obj.GetExceptionSummary(ex);

            obj.AppendToLogFile(message);
            obj.AppendToLogFile(summary);
            obj.AppendToLogFile(stack_text);
            disp(message);
            disp(' ');
            disp('If reporting an error, please include the following information:');
            disp(['  ' summary]);
            disp(['  ' stack_text '\n']);
            disp(' ');
        end

        function ShowWarning(obj, identifier, message, supplementary_info)
            [calling_function, ~] = CoreErrorUtilities.GetCallingFunction(2);

            obj.AppendToLogFile([calling_function ': WARNING: ' identifier ':' message]);
            disp(['WARNING: ' message]);
            if nargin > 3 && ~isempty(supplementary_info)
                disp('Additional information on this warning:');
                disp(supplementary_info);
            end

        end

        function Error(obj, identifier, message)
            [calling_function, stack] = CoreErrorUtilities.GetCallingFunction(2);

            if ischar(calling_function) && length(calling_function) >= 13 && strcmp(calling_function(1:13), 'CoreReporting')
                [calling_function, stack] = CoreErrorUtilities.GetCallingFunction(3);
            end

            msgStruct = [];
            msgStruct.message = ['Error in function ' calling_function ': ' message];
            msgStruct.identifier = identifier;
            msgStruct.stack = stack;
            obj.AppendToLogFile([calling_function ': ERROR: ' identifier ':' message]);
            error(msgStruct);
        end

        function ErrorFromException(obj, identifier, message, ex)
            [calling_function, stack] = CoreErrorUtilities.GetCallingFunction(2);

            msgStruct = [];
            [summary, stack_text] = obj.GetExceptionSummary(ex);
            message_text = ['Error in function ' calling_function ': ' identifier ': ' message ' Exception message:' ex.message ' : ' ex.stack(1).name];
            msgStruct.message = message_text;
            msgStruct.identifier = identifier;
            msgStruct.stack = stack;

            obj.AppendToLogFile(message_text);
            obj.AppendToLogFile(summary);
            obj.AppendToLogFile(stack_text);

            disp(message_text);
            disp(' ');
            disp('If reporting an error, please include the following information');
            disp(['  ' summary]);
            disp(['  ' stack_text]);
            disp(' ');

            error(msgStruct);
        end

        function ShowProgress(obj, text)
            adjusted_text = obj.AdjustProgressText(text);

            if ~isempty(obj.ProgressDialog) && isvalid(obj.ProgressDialog)
                if nargin > 1
                    obj.ProgressDialog.SetProgressText(adjusted_text);
                else
                    obj.ProgressDialog.SetProgressText();
                end
                obj.CurrentProgressStackItem.Visible = true;
            end
        end

        function CompleteProgress(obj)
            if ~isempty(obj.ProgressDialog)
                if isempty(obj.ProgressStack) || ~obj.ParentProgressStackItem.Visible
                    obj.ProgressDialog.Complete;
                    obj.CurrentProgressStackItem.Visible = false;
                end
            end
        end

        function UpdateProgressMessage(obj, text)
            adjusted_text = obj.AdjustProgressText(text);

            if ~isempty(obj.ProgressDialog)
                obj.ProgressDialog.SetProgressText(adjusted_text);
                obj.CurrentProgressStackItem.Visible = true;
            end
        end

        function UpdateProgressValue(obj, progress_value)
            adjusted_value = obj.AdjustProgressValue(progress_value, []);

            if ~isempty(obj.ProgressDialog)
                obj.ProgressDialog.SetProgressValue(adjusted_value);
                obj.CurrentProgressStackItem.Visible = true;
            end
            obj.CheckForCancel;
        end

        function UpdateProgressStage(obj, progress_stage, num_stages)
            progress_value = 100*progress_stage/num_stages;
            value_change = 100/num_stages;
            adjusted_value = obj.AdjustProgressValue(progress_value, value_change);
            if ~isempty(obj.ProgressDialog)
                obj.ProgressDialog.SetProgressValue(adjusted_value);
                obj.CurrentProgressStackItem.Visible = true;
            end
            obj.CheckForCancel;
        end

        function UpdateProgressAndMessage(obj, progress_value, text)
            adjusted_value = obj.AdjustProgressValue(progress_value, []);
            adjusted_text = obj.AdjustProgressText(text);

            if ~isempty(obj.ProgressDialog)
                obj.ProgressDialog.SetProgressAndMessage(adjusted_value, adjusted_text);
                obj.CurrentProgressStackItem.Visible = true;
            end

            obj.CheckForCancel;

        end

        function cancelled = HasBeenCancelled(obj)
            if ~isempty(obj.ProgressDialog)
                cancelled = obj.ProgressDialog.CancelClicked;
            else
                cancelled = false;
            end
        end

        function CheckForCancel(obj)
            if obj.HasBeenCancelled
                obj.Error(CoreReporting.CancelErrorId, 'User cancelled');
            end
        end

        function PushProgress(obj)
            obj.ProgressStack(end + 1) = obj.ParentProgressStackItem;
            obj.ParentProgressStackItem = obj.CurrentProgressStackItem;
            obj.CurrentProgressStackItem = CoreProgressStackItem('', obj.ParentProgressStackItem.MinPosition, obj.ParentProgressStackItem.MaxPosition);
            obj.CurrentProgressStackItem.Visible = obj.ParentProgressStackItem.Visible;

        end

        function PopProgress(obj)
            obj.CurrentProgressStackItem = obj.ParentProgressStackItem;
            if isempty(obj.ProgressStack)
                obj.ParentProgressStackItem = [];
            else
                obj.ParentProgressStackItem = obj.ProgressStack(end);
                obj.ProgressStack(end) = [];
            end
        end

        function ClearProgressStack(obj)
            obj.ProgressStack = CoreProgressStackItem.empty(0);
            obj.CurrentProgressStackItem = CoreProgressStackItem('', 0, 100);
            obj.ParentProgressStackItem = CoreProgressStackItem('', 0, 100);
        end

        function ShowAndClearPendingMessages(obj)
        end

        function OpenPath(obj, file_path, message)
            disp([message, ': ', file_path]);
        end
    end

    methods (Access = private)
        function AppendToLogFile(obj, message)
            file_id = fopen(obj.LogFileName, 'at');
            message = [datestr(now) ': ' message];
            fprintf(file_id, '%s\n', message);
            fclose(file_id);
        end

        function [summary, stack] = GetExceptionSummary(obj, exc)
            stack = 'Stack:';
            for s = exc.stack'
                stack = [stack s.name '(' num2str(s.line) '); '];
            end
            summary = ['Error in ' exc.stack(1).name '(PTK' PTKSoftwareInfo.Version ') : ' exc.message];
        end

        function adjusted_text = AdjustProgressText(obj, text)
            adjusted_text = text;
            obj.CurrentProgressStackItem.ProgressText = text;
        end

        function adjusted_value = AdjustProgressValue(obj, value, value_change)
            if isempty(value_change)
                value_change = value - obj.CurrentProgressStackItem.LastProgressValue;
            end
            obj.CurrentProgressStackItem.LastProgressValue = value;

            scale = (obj.ParentProgressStackItem.MaxPosition - obj.ParentProgressStackItem.MinPosition)/100;
            adjusted_value = obj.ParentProgressStackItem.MinPosition + scale*value;
            obj.CurrentProgressStackItem.MinPosition = adjusted_value;
            if value_change > 0
                obj.CurrentProgressStackItem.MaxPosition = adjusted_value + scale*value_change;
            end
        end

        function SetValueChange(obj, value_change)
            value = obj.CurrentProgressStackItem.LastProgressValue;
            scale = (obj.ParentProgressStackItem.MaxPosition - obj.ParentProgressStackItem.MinPosition)/100;
            adjusted_value = obj.ParentProgressStackItem.MinPosition + scale*value;
            obj.CurrentProgressStackItem.MinPosition = adjusted_value;
            if value_change > 0
                obj.CurrentProgressStackItem.MaxPosition = adjusted_value + scale*value_change;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%kgcreate2dtex 'Create a 2d texture from an image.'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros gcreate2dtex.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% Double: umin 'Umin', default: 0: 'Minimum u-coordinate value associated with origin of image data.'
% Double: vmin 'Vmin', default: 0: 'Mininum v-coordinate assigned to the origin of the texture.'
% Double: umax 'Umax', default: 1: 'Maximum u-coordinate to use for the non-origin corner of the texture image.'
% Double: vmax 'Vmax', default: 1: 'Maximum v-coordinate to use for the non-origin corner of the texture image.'
%
% Example: o = kgcreate2dtex(i, {'i','';'o','';'umin',0;'vmin',0;'umax',1;'vmax',1})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% gcreate2dtex - Create a 2d texture from an image.
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% gbindtexture
%
%  RESTRICTIONS
% Textures are created so as to contain 3 channels of floating point RGB data,
% so space efficiency is not optimal.  Only RGB textures are created.  Future
% work includes generating luminance textures, color-mapped textures and other
% lower-color-resolution texture formats.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996,1997 , The Regents of the University of California.  All rights reserved.
%

function varargout = kgcreate2dtex(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kgcreate2dtex(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'umin', 0;'vmin', 0;'umax', 1;'vmax', 1};
maxval={0,0,1,1,1,1};
minval={0,0,0,0,0,0};
istoggle=[0,0,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Double','Double','Double','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'gcreate2dtex"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function varargout = gui_newapproach(varargin)
% Open a dialog to select a new BCI approach
% Approach = gui_newapproch()
%
% Out:
%   Approach : the newly created approach, or [] if the user cancelled the process
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-10-25

% Last Modified by GUIDE v2.5 10-Dec-2011 18:46:11

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @gui_newapproach_OpeningFcn, ...
                   'gui_OutputFcn',  @gui_newapproach_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 initialization code - DO NOT EDIT

% --- Executes during object creation, after setting all properties.
function figure1_CreateFcn(hObject, eventdata, handles)
% retrieve and attach the known approaches...
[handles.approaches,flat,handles.indexable_approaches] = gui_listapproaches;
guidata(hObject,handles);

% --- Executes just before gui_newapproach is made visible.
function gui_newapproach_OpeningFcn(hObject, eventdata, handles, varargin)
% invoke the popupmenu at least once
popupmenu1_Callback(handles.popupmenu1,{},guidata(hObject));
% UIWAIT makes gui_newapproach wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = gui_newapproach_OutputFcn(hObject, eventdata, handles)
try
    varargout = handles.output;
    full_edit = get(handles.checkbox1,'Value') == get(handles.checkbox1,'Max');
    % ... and close the figure
    close(handles.figure1);
catch
    varargout = {[],'Cancel'};
end

if strcmp(varargout{2},'OK')
    if full_edit
        % full review
        [varargout{1:2}] = gui_reviewapproach(varargout{1},true);
    else
        % quick edit
        [varargout{1:2}] = gui_configapproach(varargout{1},true);
    end
else
    varargout = {[],'Cancel'};
end

% --- Executes on selection change in popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
idx = gui_select_bracketed(hObject);
% store the current approach in handles
handles.approach = handles.indexable_approaches{idx};
% update the description
set(handles.edit1,'String',handles.approach.description);
% write back
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)
% populate with pretty-printed list
set(hObject,'String',gui_print_grouplist(handles.approaches,@(x)x.name));

function pushbutton1_Callback(hObject, eventdata, handles)
handles.output = {handles.approach,'OK'};
guidata(hObject,handles);
uiresume(handles.figure1);

function pushbutton3_Callback(hObject, eventdata, handles)
handles.output = {[],'Cancel'};
guidata(hObject,handles);
uiresume(handles.figure1);

function pushbutton4_Callback(hObject, eventdata, handles)
env_doc code/paradigms

function edit1_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)

function checkbox1_Callback(hObject, eventdata, handles)

function pushbutton1_KeyPressFcn(hObject, eventdata, handles)

function figure1_KeyPressFcn(hObject, eventdata, handles)
if strcmp(eventdata.Key,'return')
    pushbutton1_Callback(handles.pushbutton1, eventdata, handles); end
if strcmp(eventdata.Key,'escape')
    pushbutton3_Callback(handles.pushbutton3, eventdata, handles); end

function popupmenu1_KeyPressFcn(hObject, eventdata, handles)
if strcmp(eventdata.Key,'return')
    pushbutton1_Callback(handles.pushbutton1, eventdata, handles); end
if strcmp(eventdata.Key,'escape')
    pushbutton3_Callback(handles.pushbutton3, eventdata, handles); end

repo: REST


text:
function [y, dh_struct] = findgrad(X_struct, f, e, varargin)

[X_vec, X_template] = struct2vector(X_struct);

[y] = feval(f, X_struct, varargin{:});

dh = nan(length(X_vec),1) ;

for j = 1:length(X_vec)

  dx = zeros(length(X_vec),1);
  dx(j) = dx(j) + e;

  [y2] = feval(f, vector2struct(X_vec + dx, X_template), varargin{:});
  [y1] = feval(f, vector2struct(X_vec - dx, X_template), varargin{:});

  dh(j) = (y2 - y1)/(2*e);
end

dh_struct = vector2struct(dh, X_template);

repo: CategoryShapes


text:
function [betas] = esvm_perform_platt_calibration(grid, imageset, models, ...
                                                  params)
% Perform calibration by learning the sigmoid parameters (linear
% transformation of svm scores) for each model independently. This
% type of SVM classifier calibration is due to John Platt who used
% this trick to convert SVM output scors into probabilities for
% comparison.  If we perform an operation such as NMS, we will now
% have "comparable" scores.  This is performed on the 'trainval' set
% for PASCAL VOC.
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if length(grid) == 0 || length(models) == 0
  betas = [];
  ALL_bboxes = [];
  return;
end

if length(params.dataset_params.localdir) > 0
  CACHE_FILES = 1;
else
  CACHE_FILES = 0;
end

%if enabled, do NMS, if disabled return raw detections
DO_NMS = 0;

% if enabled, display images
display = params.dataset_params.display;

% if display is enabled and dump_images is enabled, then dump images
% into DUMPDIR
dump_images = 0;

models_name = '';
if length(models)>=1 && isfield(models{1},'models_name') && ...
      isstr(models{1}.models_name)
  models_name = models{1}.models_name;
end

DUMPDIR = sprintf('%s/www/calib/',params.dataset_params.localdir, ...
                  models_name);

if dump_images==1 && ~exist(DUMPDIR,'dir')
  mkdir(DUMPDIR);
end
% show A NxN grid of top detections (if display is turned on)
SHOW_TOP_N_SVS = 10;

if nargin < 1
  fprintf(1,'Not enough arguments, need at least the grid\n');
  return;
end

final_dir = ...
    sprintf('%s/models',params.dataset_params.localdir);

if CACHE_FILES == 1 && ~exist(final_dir','dir')
  mkdir(final_dir);
end

final_file = ...
    sprintf('%s/%s-betas.mat',...
            final_dir, ...
            models_name);

if CACHE_FILES == 1
  lockfile = [final_file '.lock'];
  if fileexists(final_file) || (mymkdir_dist(lockfile)==0)

    %wait until lockfiles are gone
    wait_until_all_present({lockfile},5,1);
    fprintf(1,'Loading final file %s\n',final_file);
    res = load_keep_trying(final_file);
    betas = res.betas;
    return;
  end
end

for i = 1:length(models)
  if ~isfield(models{i},'curid')
    models{i}.curid = '-1';
  end
end

model_ids = cellfun2(@(x)x.curid,models);
targets = 1:length(models);

cls = models{1}.cls;

%if isfield(params.dataset_params,'classes')
%  targetc = find(ismember(params.dataset_params.classes, ...
%                          models{1}.cls));
%else
targetc = models{1}.cls;
%end

%fprintf(1,'Preparing boxes for calibration\n');
for i = 1:length(grid)
  if mod(i,100)==0
    fprintf(1,'.');
  end
  cur = grid{i};

  %do not process grids with no bboxes
  if size(cur.bboxes,1) == 0
    continue;
  end

  if size(cur.bboxes,1) >= 1
    cur.bboxes(:,5) = 1:size(cur.bboxes,1);
    cur.coarse_boxes(:,5) = 1:size(cur.bboxes,1);
    if DO_NMS == 1
      cur.bboxes = esvm_nms_within_exemplars(cur.bboxes,.5);
      cur.coarse_boxes = cur.coarse_boxes(cur.bboxes(:,5),:);
    end

    if length(cur.extras)>0
      cur.extras.os = cur.extras.maxos(cur.bboxes(:,5));
      try
      cur.extras.os = cur.extras.os.* ...
          reshape(double(ismember(cur.extras.maxclass,targetc)),...
                  size(cur.extras.os));
      catch
        keyboard
      end

      % cur.extras.os = cur.extras.os.* ...
      %     reshape((cur.extras.maxclass(cur.bboxes(:,5)) == ...
      %              targetc),size(cur.extras.os));

    end
  end

  cur.bboxes(:,5) = grid{i}.index;
  cur.coarse_boxes(:,5) = grid{i}.index;

  coarse_boxes{i} = cur.coarse_boxes;
  bboxes{i} = cur.bboxes;

  %if we have overlaps, collect them
  if length(cur.extras) > 0

    %use all objects as ground truth
    %goods = 1:length(cur.extras.cats);

    %% find the ground truth examples of the right category
    %goods = find(ismember(cur.extras.cats,cls));

    exids = cur.bboxes(:,6);

    %if length(goods) == 0
    %  os{i} = zeros(size(bboxes{i},1),1);
    %else
    %  curos = cur.extras.os(:,goods);
    os{i} = cur.extras.maxos; %max(curos,[],2);
    %end
  else
    os{i} = zeros(size(bboxes{i},1),1);
  end

  scores{i} = cur.bboxes(:,7)';
end

ALL_bboxes = cat(1,bboxes{:});
ALL_coarse_boxes = cat(1,coarse_boxes{:});
ALL_os = cat(1,os{:});

curids = cellfun2(@(x)x.curid,grid);
% Pre-processing models for calibration

for exid = 1:length(models)
  fprintf(1,'.');
  sourcegrid = find(ismember(curids,models{exid}.curid));
  if length(sourcegrid) == 0
    sourcegrid = -1;
  end

  hits = find((ALL_bboxes(:,6)==exid));
  all_scores = ALL_bboxes(hits,end);
  all_os = ALL_os(hits,:);

  good_scores = all_scores(all_os>=.5);
  good_os = all_os(all_os>=.5);

  bad_scores = all_scores(all_os<.5);
  bad_os = all_os(all_os<.5);

  %add virtual sample at os=1.0, score=1.0
  %good_os = cat(1,good_os,1.0);
  %good_scores = cat(1,good_scores,1.0);

  if length(good_os) <= 1 || (length(bad_os) ==0)
    beta = [.1 100];
  else

    [aa,bb] = sort(bad_scores,'descend');
    curlen = min(length(bb),10000*length(good_scores));
    bb = bb(round(linspace(1,length(bb),curlen)));

    bad_scores = bad_scores(bb);
    bad_os = bad_os(bb);
    all_scores = [good_scores; bad_scores];
    all_os = [good_os; bad_os];

    [aaa,bbb] = sort(all_scores, 'descend');
    %bbb = bbb(1:min(1000,length(bbb)));
    %all_scores = all_scores(bbb);
    %all_os = all_os(bbb);
    beta = esvm_learn_sigmoid(all_scores, all_os);
  end

  if beta(1)<.001
    fprintf(1,['warning[esvm_perform_platt_calibration.m]: beta(1)' ...
               ' is low']);
    % beta(1) = .001;
  end

  betas(exid,:) = beta;

  if (sum(ismember(exid,targets))==0)
    continue
  end

  if display == 1

    figure(1)
    clf
    subplot(1,2,1)

    xs = linspace(min(all_scores),max(all_scores),1000);
    fx = @(x)(1./(1+exp(-beta(1)*(x-beta(2)))));

    plot(xs,fx(xs),'b','LineWidth',2)
    hold on
    plot(all_scores,all_os,'r.','MarkerSize',14)

    axis([min(xs) max(xs) 0 1])
    xlabel('SVM score')
    ylabel(sprintf('Max Overlap Score with %s',models{exid}.cls))

    title(sprintf('Learned Sigmoid \\beta=[%.3f %.3f]',beta(1), ...
                  beta(2)))

    subplot(1,2,2)
    Iex = convert_to_I(models{exid}.I);
    imagesc(Iex)
    plot_bbox(models{exid}.gt_box)
    axis image
    axis off

    title(sprintf('Exemplar %s.%d.%s',...
                  models{exid}.curid,...
                  models{exid}.objectid, ...
                  models{exid}.cls))
    drawnow
    snapnow

    bbs=ALL_coarse_boxes(hits,:);
    bbs_os = ALL_os(hits,:);
    [aa,bb] = sort(bbs(:,end),'descend');
    bbs_show = bbs(bb,:);

    %models{exid}.model.svids = {};
    %m = try_reshape(models{exid},bbs_show,100);

    %[models{exid}.model.svids,models{exid}.model.nsv] = ...
    %    esvm_reconstruct_features(bbs_show,100);%,'trainval','');

    models{exid}.model.svbbs = bbs_show;
    m2 = models(exid);
    if length(imageset) > 0
      m2{1}.train_set = imageset;
      m2{1}.model.svbbs(:,6) = 1;
      m2{1}.model.svxs = [];
      figure(445)
      clf
      imagesc(esvm_show_det_stack(m2{1},8))
      axis image
      axis off
      drawnow
      title(sprintf('Topdets Ex %s.%d.%s',...
                    models{exid}.curid,...
                    models{exid}.objectid, ...
                    models{exid}.cls))
      drawnow
      snapnow
    end
  end

  if (display == 0)
    continue
  end

  if dump_images == 1
    figure(2)
    filer = sprintf('%s/result.%d.%s.%s.png', DUMPDIR, ...
                    exid,models{exid}.cls,models{exid}.models_name);
    set(gcf,'PaperPosition',[0 0 20 20]);
    print(gcf,filer,'-dpng');

  else
    pause(.01)
  end
end

if CACHE_FILES == 1
  fprintf(1,['\nLoaded calibration parameters "betas", saving to' ...
             ' %s\n'],final_file);
  save(final_file,'betas');
  rmdir(lockfile);
end

author: quantombone

repo: exemplarsvm


text:
function y = atanh(x)
%ATANH        Implements  atanh(x)  for intervals
%
%   y = atanh(x)
%
%interval standard function implementation
%

% written  12/30/98     S.M. Rump
% modified 08/31/99     S.M. Rump  complex allowed, sparse input,
%                                  major revision, improved accuracy
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 12/04/05     S.M. Rump  'realstdfctsexcptnignore' added and some
%                                     improvements, tocmplx replaced by cintval
% modified 09/06/07     S.M. Rump  approximate std fcts removed, exceptional arguments
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/18/08     S.M. Rump  StdFctsException ignore/NaN
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  if issparse(x)
    [ix,jx,sx] = find(x);
    [m,n] = size(x);
    y = sparse(ix,jx,atanh(full(sx)),m,n);
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if x.complex
    y = log( 2./(1-x) - 1 )/2;
    if rndold
      setround(rndold)
    end
    return
  end

  % input x real and full
  % real range of definition:  [-1,1]
  INTLAB_STDFCTS_EXCPTN = getappdata(0,'INTLAB_STDFCTS_EXCPTN');
  indexneg = ( x.inf<-1 );               % (partially) exceptional indices
  indexpos = ( x.sup>1 );                % (partially) exceptional indices
  if ~isempty(find(indexneg)) | ~isempty(find(indexpos)) % handle input out-of-range
    if INTLAB_STDFCTS_EXCPTN<=1   % out-of-range input handled as complex
      if INTLAB_STDFCTS_EXCPTN==1
        warning('ATANH: Real interval input out of range changed to be complex')
      end
      y = x;
      index = indexneg | indexpos;
      %VVVV  y(index) = atanh(cintval(x(index)));
      s.type = '()'; s.subs = {index}; y = subsasgn(y,s,atanh(cintval(subsref(x,s))));
      %AAAA  Matlab bug fix
      index = ~index;
      if any(index(:))
        %VVVV  y(index) = atanh(x(index));
        s.type = '()'; s.subs = {index}; y = subsasgn(y,s,atanh(subsref(x,s)));
        %AAAA  Matlab bug fix
      end
      if rndold
        setround(rndold)
      end
      return
    end
    setappdata(0,'INTLAB_STDFCTS_EXCPTN_',1);
    if INTLAB_STDFCTS_EXCPTN==3    % ignore input out of range (ignore-mode)
      x.inf(indexneg) = -1;               % completely exceptional indices treated below
      x.sup(indexpos) = 1;
      index = ( x.sup<-1 ) | ( x.inf>1 ); % completely exceptional indices
    end
  else
    index = [];                           % make sure Index is not undefined
  end

  % input x real and full
  y = x;
  wng = warning;
  warning off

  % treat non-exceptional arguments
  xinf = x.inf(:);
  xsup = x.sup(:);

  IndexInfPos = ( xinf>=0 );
  len1 = sum(IndexInfPos);
  IndexSupNeg = ( xsup<=0 );
  len2 = sum(IndexSupNeg);

  Y = atanh_pos( [ xinf(IndexInfPos) ; -xsup(IndexSupNeg) ] , -1 );
  y.inf(IndexInfPos) = Y(1:len1);
  y.sup(IndexSupNeg) = -Y( len1+1 : end );

  IndexInfNeg = ( xinf<0 );
  len1 = sum(IndexInfNeg);
  IndexSupPos = ( xsup>0 );
  len2 = sum(IndexSupPos);

  Y = atanh_pos( [ -xinf(IndexInfNeg) ; xsup(IndexSupPos) ] , 1 );
  y.inf(IndexInfNeg) = -Y(1:len1);
  y.sup(IndexSupPos) = Y( len1+1 : len1+len2 );

  y.inf(xinf==-1) = -inf;
  y.sup(xsup==-1) = -inf;
  y.inf(xinf==1) = inf;
  y.sup(xsup==1) = inf;

  if INTLAB_STDFCTS_EXCPTN==3      % ignore input out of range (ignore-mode)
    if ~isempty(find(index))              % completely exceptional arguments to NaN
      y.inf(index) = NaN;
      y.sup(index) = NaN;
    end
  else                                    % any input out of range to NaN (NaN-mode)
    index = indexneg | indexpos;
    if ~isempty(find(index))              % exceptional arguments to NaN
      y.inf(index) = NaN;
      y.sup(index) = NaN;
    end
  end

  setround(rndold)
  warning(wng)

repo: OpenMAS


text:
% This make.m is for MATLAB and OCTAVE under Windows, Mac, and Unix
function make()
try
	% This part is for OCTAVE
	if (exist ('OCTAVE_VERSION', 'builtin'))
		mex libsvmread.c
		mex libsvmwrite.c
		mex -I.. svmtrain.c ../svm.cpp svm_model_matlab.c
		mex -I.. svmpredict.c ../svm.cpp svm_model_matlab.c
	% This part is for MATLAB
	% Add -largeArrayDims on 64-bit machines of MATLAB
	else
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmread.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmwrite.c

        % Activate this if we want OMP implementation
% 		mex CFLAGS="\$CFLAGS -std=c99" CXXFLAGS="\$CXXFLAGS -fopenmp" -largeArrayDims -I.. -lgomp svmtrain.c ../svm.cpp svm_model_matlab.c
% 		mex CFLAGS="\$CFLAGS -std=c99" CXXFLAGS="\$CXXFLAGS -fopenmp" -largeArrayDims -I.. -lgomp svmpredict.c ../svm.cpp svm_model_matlab.c

        % Normal CFlage without OMP
        mex CFLAGS="\$CFLAGS -std=c99"  -largeArrayDims -I.. svmtrain.c ../svm.cpp svm_model_matlab.c
		mex CFLAGS="\$CFLAGS -std=c99"  -largeArrayDims -I.. svmpredict.c ../svm.cpp svm_model_matlab.c

	end
catch err
	fprintf('Error: %s failed (line %d)\n', err.stack(1).file, err.stack(1).line);
	disp(err.message);
	fprintf('=> Please check README for detailed instructions.\n');
end

author: epfl-lasa


text:
%
% References:
%
% C. Lu. A Library of ADMM for Sparse and Low-rank Optimization. National University of Singapore, June 2016.
% https://github.com/canyilu/LibADMM.
% C. Lu, J. Feng, S. Yan, Z. Lin. A Unified Alternating Direction Method of Multipliers by Majorization
% Minimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 40, pp. 527-541, 2018
%

addpath(genpath(cd))
clear

%% Examples for testing the low-rank tensor models
% For detailed description of the sparse models, please refer to the Manual.

opts.mu = 1e-6;
opts.rho = 1.1;
opts.max_iter = 500;
opts.DEBUG = 1;

%% Tensor RRPCA based on sum of nuclear norm minimization (rpca_snn)
n1 = 50;
n2 = n1;
n3 = n1;
r = 5
L = rand(r,r,r);
U1 = rand(n1,r);
U2 = rand(n2,r);
U3 = rand(n3,r);
L = nmodeproduct(L,U1,1);
L = nmodeproduct(L,U2,2);
L = nmodeproduct(L,U3,3); % low rank part

p = 0.05;
m = p*n1*n2*n3;
temp = rand(n1*n2*n3,1);
[~,I] = sort(temp);
I = I(1:m);
Omega = zeros(n1,n2,n3);
Omega(I) = 1;
E = sign(rand(n1,n2,n3)-0.5);
S = Omega.*E; % sparse part, S = P_Omega(E)

Xn = L+S;

lambda = sqrt([max(n1,n2*n3), max(n2,n1*n3), max(n3,n1*n2)]);
lambda = [1 1 1]
[Lhat,Shat,err,iter] = trpca_snn(Xn,lambda,opts);

err
iter

%% low-rank tensor completion based on sum of nuclear norm minimization (lrtc_snn)
n1 = 50;
n2 = n1;
n3 = n1;
r = 5;
X = rand(r,r,r);
U1 = rand(n1,r);
U2 = rand(n2,r);
U3 = rand(n3,r);
X = nmodeproduct(X,U1,1);
X = nmodeproduct(X,U2,2);
X = nmodeproduct(X,U3,3);
p = 0.5;
omega = find(rand(n1*n2*n3,1)<p);
M = zeros(n1,n2,n3);
M(omega) = X(omega);

lambda = [1 1 1];
[Xhat,err,iter] = lrtc_snn(M,omega,lambda,opts);
err
iter
RSE = norm(X(:)-Xhat(:))/norm(X(:))

%% regularized low-rank tensor completion based on sum of nuclear norm minimization (lrtcR_snn)
n1 = 50;
n2 = n1;
n3 = n1;
r = 5;
X = rand(r,r,r);
U1 = rand(n1,r);
U2 = rand(n2,r);
U3 = rand(n3,r);
X = nmodeproduct(X,U1,1);
X = nmodeproduct(X,U2,2);
X = nmodeproduct(X,U3,3);
p = 0.5;
omega = find(rand(n1*n2*n3,1)<p);
M = zeros(n1,n2,n3);
M(omega) = X(omega);
lambda = [1 1 1];
[Xhat,err,iter] = lrtcR_snn(M,omega,lambda,opts);
err
iter

%% Tensor RRPCA based on tensor nuclear norm minimization (rpca_tnn)
n1 = 50;
n2 = n1;
n3 = n1;
r = 0.1*n1 % tubal rank
L1 = randn(n1,r,n3)/n1;
L2 = randn(r,n2,n3)/n2;
L = tprod(L1,L2); % low rank part

p = 0.1;
m = p*n1*n2*n3;
temp = rand(n1*n2*n3,1);
[~,I] = sort(temp);
I = I(1:m);
Omega = zeros(n1,n2,n3);
Omega(I) = 1;
E = sign(rand(n1,n2,n3)-0.5);
S = Omega.*E; % sparse part, S = P_Omega(E)

Xn = L+S;
lambda = 1/sqrt(n3*max(n1,n2));

tic
[Lhat,Shat] = trpca_tnn(Xn,lambda,opts);

RES_L = norm(L(:)-Lhat(:))/norm(L(:))
RES_S = norm(S(:)-Shat(:))/norm(S(:))
trank = tubalrank(Lhat)

%% low-rank tensor completion based on tensor nuclear norm minimization (lrtc_tnn)
n1 = 50;
n2 = n1;
n3 = n1;
r = 0.1*n1 % tubal rank
L1 = randn(n1,r,n3)/n1;
L2 = randn(r,n2,n3)/n2;
X = tprod(L1,L2); % low rank part
p = 0.5;
omega = find(rand(n1*n2*n3,1)<p);
M = zeros(n1,n2,n3);
M(omega) = X(omega);

[Xhat,obj,err,iter] = lrtc_tnn(M,omega,opts);

err
iter
RSE = norm(X(:)-Xhat(:))/norm(X(:))
trank = tubalrank(Xhat)

%% regularized low-rank tensor completion based on tensor nuclear norm minimization (lrtcR_tnn)
n1 = 50;
n2 = n1;
n3 = n1;
r = 0.1*n1 % tubal rank
L1 = randn(n1,r,n3)/n1;
L2 = randn(r,n2,n3)/n2;
X = tprod(L1,L2); % low rank part
p = 0.5;
omega = find(rand(n1*n2*n3,1)<p);
M = zeros(n1,n2,n3);
M(omega) = X(omega);

lambda = 0.5;
[Xhat,Ehat,obj,err,iter] = lrtcR_tnn(M,omega,lambda,opts);
err
iter

%% low-rank tensor recovery from Gaussian measurements based on tensor nuclear norm minimization (lrtr_Gaussian_tnn)
n1 = 30;
n2 = n1;
n3 = 5;
r = 0.2*n1; % tubal rank
X = tprod(randn(n1,r,n3),randn(r,n2,n3)); % size: n1*n2*n3

m = 3*r*(n1+n2-r)*n3+1; % number of measurements
n = n1*n2*n3;
A = randn(m,n)/sqrt(m);

b = A*X(:);
Xsize.n1 = n1;
Xsize.n2 = n2;
Xsize.n3 = n3;

opts.DEBUG = 1;
[Xhat,obj,err,iter]  = lrtr_Gaussian_tnn(A,b,Xsize,opts);

RSE = norm(Xhat(:)-X(:))/norm(X(:))
trank = tubalrank(Xhat)

author: canyilu

repo: LibADMM-toolbox


text:
function varargout = pdeSolve(pdeFun, tIn, u0, bc, varargin)
%PDESOLVE   Solve PDEs using Chebfun.
%
%   PDESOLVE() solves solves an initial-boundary value problem via a method of
%   lines approach. See @CHEBFUN/PDE15S() for documentation.
%
% See also PDE15S, PDE23T, PDESET.

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

global DIFFORDER SYSSIZE
DIFFORDER = 0;
SYSSIZE = 0;

% Default options:
tol = 1e-6;             % 'eps' in Chebfun terminology
% The default behaviour with no outputs is to plot. If the method is called with
% outputs, by default, we don't plot.
doPlot = ( nargout == 0 );
doHold = false;         % Hold plot?
plotOpts = {'-'};       % Plotting style
adjustBCs = true;       % Adjust inconsistent BCs
throwBCwarning = true;  % Throw a warning for inconsistent BCs
timeChunks = 51;        % Default number of time slices if not specified

% Parse the variable inputs:
if ( numel(varargin) == 2 )
    opt = varargin{1};
    opt.N = varargin{2};
elseif ( numel(varargin) == 1 )
    if ( isstruct(varargin{1}) )
        opt = varargin{1};
    else
        opt = pdeset;
        opt.N = varargin{1};
    end
else
    opt = pdeset;
end
optN = opt.N;
if ( isempty(optN) )
    optN = NaN;
end

% Get the domain:
DOMAIN = domain(u0, 'ends');

if ( isfield(opt, 'AdjustBCs') && ~isempty(opt.AdjustBCs) && ~opt.AdjustBCs )
    throwBCwarning = false;
    adjustBCs = false;
end

% PDE solver options:
if ( ~isempty(opt.Eps) )
    tol = opt.Eps;
end
if ( ~isempty(opt.Plot) )
    doPlot = strcmpi(opt.Plot, 'on');
end
if ( ~isempty(opt.HoldPlot) )
    doHold = strcmpi(opt.HoldPlot, 'on') || strcmp(opt.HoldPlot, '1');
end
if ( ~isempty(opt.PlotStyle) )
    plotOpts = opt.PlotStyle;
end
if ( ~isempty(opt.ODESolver) )
    ODESOLVER = opt.ODESolver;
else
    ODESOLVER = @ode15s;
end

% Set time chunks:
if ( length(tIn) == 2 )
    tt = linspace(tIn(1), tIn(2), timeChunks);
else
    tt = tIn;
end

userMassSet = false;
ISPERIODIC = false;
DONE = false;
COUNTER = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%  PARSE INPUTS TO PDEFUN  %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Determine the size of the system, i.e., number of dependent variables.
SYSSIZE = min(size(u0));
[pdeFun, varNamesParsed] = parseFun(pdeFun);
if ( isfield(opt, 'difforder') )
    DIFFORDER = opt.difforder;
else
    getDIFFORDER(pdeFun, DOMAIN);
end

if ( (max(DIFFORDER) < 2) && isequal(ODESOLVER, @ode15s) )
    warning('CHEBFUN:PdeSolver', ...
        'PDE23T() is recommended for non-diffusive problems.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  EVENT SETUP  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( ~isnan(optN) )
    opt.OutputFcn = @nonAdaptiveEvent;
else
    opt.OutputFcn = @adaptiveEvent;
end

    function status = nonAdaptiveEvent(t, U, flag)
        % This event is called at the end of each chunk in non-adaptive mode.
        status = false;
        if ( ~isempty(flag) )
            return
        end

        % Sometimes we get given more than one time slice.
        for kk = 1:numel(t)

            if ( ~any(abs(tSpan - t(kk)) < 1e-6) )
                % This is not a designated time slice!
                continue
            end

            % Reshape solution:
            Uk = reshape(U(:,kk), n, SYSSIZE);
            uCurrent = chebfun(Uk, DOMAIN, 'tech', techHandle);
            tCurrent = t(kk);
            % Store for output:
            COUNTER = COUNTER + 1;
            uOut{COUNTER} = uCurrent;

            % Plot current solution:
            if ( doPlot )
                plotFun(uCurrent, t(kk));
            end
        end

        if ( guiFlag )
            status = guiEvent(status);
        end

    end

    function status = adaptiveEvent(t, U, flag)
        % This event is called at the end of each chunk in adaptive mode.

        status = false;
        if ( ~isempty(flag) )
            return
        end

        % Sometimes we get given more than one time slice.
        for kk = 1:numel(t)
            % Reshape solution:
            Uk = reshape(U(:,kk), currentLength, SYSSIZE);

            % Happiness check:
            c = (1+sin(1:SYSSIZE)).'; % Arbitrarily linear combination.
            Uk2 = (Uk*c/sum(c));
            uk2 = tech.make(Uk2, pref);
            [ishappy, cutoff] = happinessCheck(uk2, Uk2, [], [], pref); %#ok<ASGLU>

            if ( ishappy )

                if ( ~any(abs(tSpan - t(kk)) < 1e-6) )
                    % This is not a designated time slice!
                    continue
                end

                % Store these values:
                tCurrent = t(kk);
                uCurrent = chebfun(Uk, DOMAIN, 'tech', techHandle);
                uCurrent = simplify(uCurrent);

                COUNTER = COUNTER + 1;
                uOut{COUNTER} = uCurrent;
                tOut(COUNTER) = tCurrent;

                % Plot current solution (if plotting is on)
                if ( doPlot )
                    plotFun(uCurrent, tCurrent);
                end
                % TODO: Re-insert this?
%                 % If we have 2.5 times as many coefficients as we need then
%                 % shorten the representation and cause the integrator to stop.
%                 if ( cutoff < 0.4*n && n > 17)
%                     currentLength = round(1.25*cutoff)';
%                     %currentLength = floor( currentLength / 1.5  );
%                     currentLength = currentLength + 1 - rem(currentLength,2);
%                     currentLength = max(currentLength, 17);
%                     status = true;
%                     return
%                 end

            else

                % Increase length and bail out:
                currentLength = 2*currentLength-1;
                status = true;
                break

            end

            if ( guiFlag )
                status = guiEvent(status);
            end

        end

    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%  PLOTTING SETUP  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Determine which figure to plot to (for CHEBGUI) and set default display values
% for variables.
YLim = opt.YLim;
gridOn = 0;
guiFlag = false;
if ( isfield(opt, 'handles') )
    if ( opt.handles.gui )
        guiFlag = true;
        axesSol = opt.handles.fig_sol;
        axesNorm = opt.handles.fig_norm;
        gridOn = opt.handles.guifile.options.grid;
        solveButton = opt.handles.button_solve;
        clearButton = opt.handles.button_clear;
        panelSol = opt.handles.panel_figSol;
    end
    varNames = opt.handles.varnames;
    xLabel = opt.handles.indVarName{1};
    tlabel = opt.handles.indVarName{2};
    fontsize = opt.handles.fontsizePanels;
    set(axesSol, 'fontsize', fontsize);
else
    % Obtain the correct variable names from the parsing of the operator
    if ( length(varNamesParsed) == numColumns(u0) )
        % Space and time variables did not get passed explicitly to operator
        tlabel = 't';
        xLabel = 'x';
    else
        tlabel = varNamesParsed{1};
        xLabel = varNamesParsed{2};
    end
    varNames = varNamesParsed(end-numColumns(u0)+1:end);
    if ( doPlot )    % Set up ploting in non-gui mode
        axesSol = gca;
        cla(axesSol)
        xlabel(axesSol, xLabel);
        % Bring attention to the figure
        shg
    end
end

% Set plot options
if ( doPlot )

    % Use linear scale. Always.
    if ( ~ishold )
        set(gca, 'XScale', 'Linear', 'YScale', 'Linear');
    end

    set(axesSol, 'NextPlot', 'replacechildren');

    % Fix x limits
    set(axesSol, 'xLim', u0.domain);

    % Fix y limits if they were specified
    if ( ~isempty(YLim) )
        set(axesSol, 'ylim', YLim);
    end

    % Initialize lines for plotting
    hLines = plot(axesSol, DOMAIN, NaN(length(DOMAIN), numColumns(u0)), ...
        plotOpts{:});

    % Grid on?
    if ( gridOn )
        grid(axesSol, 'on')
    else
        grid(axesSol, 'off')
    end

    % Hold on?
    if ( doHold )
        hold(axesSol, 'on')
    else
        hold(axesSol, 'off')
    end

    % For plotting, it's useful to know whether we're running in old or new
    % Matlab graphics mode
    if ( ~verLessThan('matlab', '8.4') )
        newMatlabVersion = true;
    else
        newMatlabVersion = false;
    end
end

    function status = guiEvent(status)
        %GUIEVENT   Deal with GUI events ('stop', 'pause', etc).
        % OUTPUTS:
        %   status = true exits the current time chunk.
        %   DONE = true exits PDE solver. (Note DONE is a GLOBAL variable).

        % Interrupt computation if stop or pause button is pressed in the GUI.
        if ( strcmp(get(solveButton, 'String'), 'Solve') )
            % Stop.
            tt = tt( tt <= tCurrent );
            uOut(COUNTER+1:end) = [];
            status = true;
            DONE = true;
        elseif ( strcmp(get(clearButton, 'String'), 'Continue') )
            % Wait, pause.

            % Plot a waterfall plot to the bottom figure window:
            axes(axesNorm)
            uuTmp = prepare4output(uOut(1:COUNTER));
            if ( SYSSIZE == 1 )
                waterfall(uuTmp, tt(tt<=tCurrent), 'linewidth', 2);
            else
                cols = get(0, 'DefaultAxesColorOrder');
                waterfall(uuTmp, tt(tt<=tCurrent), 'linewidth', 2, ...
                    'EdgeColors', cols);
            end
            xlabel(xLabel), ylabel(tlabel), zlabel(varNames)
            if ( gridOn )
                grid on
            end
            view([322.5 30])
            box off
            axes(axesSol)
            % Hang around until 'CONTINUE' or 'STOP' is presed.
            waitfor(clearButton, 'String');
            % Call again to see if 'STOP' was pressed.
            status = guiEvent(status);

            set(axesNorm, 'fontsize', fontsize);
        end
    end

    function plotFun(U, t)
        %PLOTFUN    Plot current solution U at a time t.

        % Obtain plotting data for the solution
        Udata = plotData(U);
        UdataX = Udata.xLine;
        UdataY = Udata.yLine;

        % Update the YData of the lines plotted
        if ( ~doHold )
            for hCounter = 1:length(hLines)
                set(hLines(hCounter), 'XData', UdataX);
                set(hLines(hCounter), 'YData', UdataY(:, hCounter)');
            end
        else
            % Reset color cycle prior to point plot if running R2014b.
            if ( newMatlabVersion )
                set(axesSol, 'ColorOrderIndex', 1);
            end
            hLines = plot(axesSol, UdataX, UdataY, plotOpts{:});
        end

        % Update the title, either of the GUI panel or the figure
        if ( guiFlag )
            set(panelSol, 'Title', sprintf('%s = %.3f', tlabel, t))
        elseif ( nargin > 1 )
            title(sprintf('%s = %.3f', tlabel, t))
        end

        % Update the plot
        drawnow

    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%% PARSE BOUNDARY CONDITIONS %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( ischar(bc) && (strcmpi(bc, 'neumann') || strcmpi(bc, 'dirichlet')) )
    bc = struct( 'left', bc, 'right', bc);
elseif ( iscell(bc) && numel(bc) == 2 )
    bc = struct( 'left', bc{1}, 'right', bc{2});
end
if ( isstruct(bc) && ~isfield(bc, 'middle') )
    bc.middle.op = [];
end

% Initialise some global variables.
numLinBCs = 0; % linear BCS (total)
numLBCs = 0;   % nonlinear LBCs
numMBCs = 0;   % nonlinear other constraints
numRBCs = 0;   % nonlinear RBCs
BCRHS = {};    % RHS values for BCs

    % Differentiation matrices at left and right endpoints:
    function D = neumannLeft(n)
        x = -cos(pi*(1:n-2)/(n-1));
        D = [0, 2./(1+x), .5];
        D(1:2:end) = -D(1:2:end);
        D(1) = -sum(D);
    end
    function D = neumannRight(n)
        x = -cos(pi*(1:n-2)/(n-1));
        D = [-.5, 2./(x-1), 0];
        D(end:-2:1) = -D(end:-2:1);
        D(end) = -sum(D);
    end

if ( ischar(bc) && any(strcmpi(bc, {'periodic', 'trig'})) )
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%% PERIODIC BCS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    ISPERIODIC = true;

    % One can still use a Chebyshev basis and enforce periodic conditions by
    % enforcing suitable constraints on the derivative. However, using a
    % periodic basis (for now, trigtech) is much more efficient.
%     r = cell(sum(DIFFORDER), 1);
%     count = 1;
%     for j = 1:SYSSIZE
%         for k = 0:DIFFORDER(j)-1
%             c = (diff(DOMAIN)/2)^k;
%             A = @(n) [1 zeros(1, n-2) -1]*chebcolloc2.diffmat(n, k)*c;
%             r{count} = @(n) [zeros(1, (j-1)*n) A(n) zeros(1,(SYSSIZE-j)*n)];
%             count = count + 1;
%         end
%     end
%     bc = struct( 'left', [], 'right', []);
%     bc.left.op = r(1:2:end);
%     bc.right.op = r(2:2:end);
%     BCRHS = num2cell(zeros(1, numel(r)));

else
    %%%%%%%%%%%%%%%%%%%%%%%%%%% NONPERIODIC BCS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ( isfield(bc, 'left') && ~isfield(bc, 'right') )
        bc.right = [];
    elseif ( isfield(bc, 'right') && ~isfield(bc, 'left') )
        bc.left = [];
    elseif ( ~isfield(bc, 'left') && ~isfield(bc, 'right') )
        bc = struct('left', [], 'right', [], 'middle', bc);
    end

    % Deal with struct and numeric input:
    bc.left = dealWithStructInput(bc.left);
    bc.right = dealWithStructInput(bc.right);

    % Temporary chebdouble for evaluating BC function handles:
    uTmp = chebdouble(ones(1, SYSSIZE), DOMAIN);

    %% LEFT BCS
    if ( isempty(bc.left) )
        bc.left = struct('op', []);

    elseif ( ischar(bc.left) || (iscell(bc.left) && ischar(bc.left{1})) )
        %%%%%%%%%%%%%%%%%%% DIRICHLET AND NEUMANN BCS (LEFT) %%%%%%%%%%%%%%%%%
        if ( iscell(bc.left) )
            v = bc.left{2};
            bc.left = bc.left{1};
        else
            v = 0;
        end
        if ( ~isnumeric(v) )
            error('CHEBFUN:CHEBFUN:pde15s:nonNumericVal1', ...
                'For BCs of the form {char, val} val must be numeric.')
        end
        BCRHS = num2cell(repmat(v, SYSSIZE, 1));
        numLinBCs = numel(v)*SYSSIZE;

        if ( strcmpi(bc.left, 'dirichlet') )
            A = @(n) [1, zeros(1, n-1)];
        elseif ( strcmpi(bc.left, 'neumann') )
            A = @(n) neumannLeft(n)*(diff(DOMAIN)/2);
        else
            error('CHEBFUN:CHEBFUN:pde15s:bcSyntax1', 'Unknown BC syntax');
        end
        bc.left = struct('op', []);
        bc.left.op = cell(SYSSIZE, 1);
        for k = 1:SYSSIZE
            bc.left.op{k} = ...
                @(n) [zeros(1,(k-1)*n),  A(n),  zeros(1, (SYSSIZE-k)*n)];
        end

    elseif ( numel(bc.left) == 1 && isa(bc.left, 'function_handle') )
        %%%%%%%%%%%%%%%%%%%%%%%  GENERAL BCS (LEFT)  %%%%%%%%%%%%%%%%%%%%%%%%%%
        opLBC = parseFun(bc.left, 'lbc');
        numLBCs = max(size(opLBC(0, DOMAIN(1), uTmp)));
        bc.left = struct('op', []);

    else
        error('CHEBFUN:CHEBFUN:pde15s:bcSyntax2', 'Unknown BC syntax');

    end

    %% MIDDLE BCS / OTHER CONSTRAINTS
    if ( isfield(bc, 'middle') && isa(bc.middle, 'function_handle') )
        %%%%%%%%%%%%%%%%%%%%%  GENERAL BCS (MIDDLE)  %%%%%%%%%%%%%%%%%%%%%%%%%%
        opMBC = parseFun(bc.middle, 'bc');
        numMBCs = max(size(opMBC(0, mean(DOMAIN), uTmp)));
        bc.middle = [];

    end

    %% RIGHT BCS
    if ( isempty(bc.right) )
        bc.right = struct('op', []);

    elseif ( ischar(bc.right) || (iscell(bc.right) && ischar(bc.right{1})) )
        %%%%%%%%%%%%%%%%%%% DIRICHLET AND NEUMANN BCS (RIGHT) %%%%%%%%%%%%%%%%%
        if ( iscell(bc.right) )
            v = bc.right{2};
            bc.right = bc.right{1};
        else
            v = 0;
        end
        if ( ~isnumeric(v) )
            error('CHEBFUN:CHEBFUN:pde15s:nonNumericVal1', ...
                'For BCs of the form {char, val} val must be numeric.')
        end
        BCRHS = [BCRHS, num2cell(repmat(v, SYSSIZE, 1))];
        numLinBCs = numLinBCs + numel(v)*SYSSIZE;

        if ( strcmpi(bc.right, 'dirichlet') )
            A = @(n) [zeros(1, n-1), 1];
        elseif ( strcmpi(bc.right, 'neumann') )
            A = @(n) neumannRight(n)*(diff(DOMAIN)/2);
        else
            error('CHEBFUN:CHEBFUN:pde15s:bcSyntax3', 'Unknown BC syntax');
        end
        bc.right = struct('op', []);
        bc.right.op = cell(SYSSIZE, 1);
        for k = 1:SYSSIZE
            bc.right.op{k} = ...
                @(n) [zeros(1,(k-1)*n),  A(n),  zeros(1,(SYSSIZE-k)*n)];
        end

    elseif ( numel(bc.right) == 1 && isa(bc.right, 'function_handle') )
        %%%%%%%%%%%%%%%%%%%%%%  GENERAL BCS (RIGHT)  %%%%%%%%%%%%%%%%%%%%%%%%%%
        opRBC = parseFun(bc.right, 'rbc');
        numRBCs = max(size(opRBC(0, DOMAIN(end), uTmp)));
        bc.right = struct('op', []);

    else
        error('CHEBFUN:CHEBFUN:pde15s:bcSyntax4', 'Unknown BC syntax');

    end

end

% Total number of BCs:
numBCs = numLinBCs + numLBCs + numMBCs + numRBCs;

if ( isstruct(bc) && ~isfield(bc, 'middle') )
    bc.middle.op = [];
end

%% %%%%%%%%%%%%%%%%%%%%%%% SUPPORT FOR COUPLED BVP-PDES! %%%%%%%%%%%%%%%%%%%%%%%
% Experimental feature for coupled ode/pde systems: (An entry equal to 1 denotes
% that the corresponding variable appears with a time derivative. 0 otherwise.)
if ( isfield(opt, 'PDEflag') && ~isempty(opt.PDEflag) )
    pdeFlag = opt.PDEflag;
else
    pdeFlag = true;
end
if ( numel(pdeFlag) == 1 )
    pdeFlag = repmat(pdeFlag, 1, SYSSIZE);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PERIODIC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( ~ISPERIODIC )
    techHandle = @chebtech2;
    points = @chebpts;
    mydouble = @chebdouble;
else
    techHandle = @trigtech;
    points = @trigpts;
    mydouble = @trigdouble;
    u0 = chebfun(u0, 'trig', 'eps', tol);
    u0 = simplify(u0, tol);
end
tech = techHandle();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%% SETUP TOLERANCES AND INITIAL CONDITION %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ODE tolerances:
aTol = odeget(opt, 'AbsTol', tol);
rTol = odeget(opt, 'RelTol', tol);
if ( isnan(optN) )
    aTol = min(aTol, tol);
    rTol = min(rTol, tol);
end
opt.AbsTol = aTol;
opt.RelTol = rTol;

% Check for (and try to remove) piecewise initial conditions:
u0Trans = u0(1).isTransposed;
if ( u0Trans )
    u0 = transpose(u0);
end
for k = 1:numel(u0)
    if ( numel(u0(k).funs) > 1 )
        u0(k) = merge(u0(k), 'all', 1025, tol);
        if ( u0(k).nfuns > 1 )
            error('CHEBFUN:CHEBFUN:pde15s:piecewise', ...
                'Piecewise initial conditions are not supported.');
        end
    end
end

% Simplify initial condition to tolerance or fixed size in optN:
if ( isnan(optN) )
    u0 = simplify(u0);
else
    for k = 1:numel(u0)
        u0(k).funs{1}.onefun = prolong(u0(k).funs{1}.onefun, optN);
    end
end

% Initial condition:
uCurrent = u0;
% Storage:
uOut = cell(1, numel(tt));
uOut{1} = uCurrent;
tOut(1) = tt(1);

% Initialise variables for ONESTEP():
M = []; P = []; n = []; B = 0;

% Set the preferences:
pref = tech.techPref();
pref.chebfuneps = tol;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% TIME CHUNKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plot initial condition:
currentLength = length(u0);
if ( doPlot )
    plotFun(u0, tt(1));
    leg = legend(axesSol, varNames);
    set(leg, 'HandleVisibility', 'off')
    drawnow
end
DONE = false;

if ( ~isnan(optN) )
    % Non-adaptive in space:
    currentLength = optN;
    tSpan = tt;
    solvePDE(tSpan); % Do all chunks at once!

else
    % Adaptive in space

    % Min-length is 9
    currentLength = max(currentLength, 9);
    tCurrent = tt(1);
    tSpan = tt;
    while ( tCurrent < tt(end) && ~DONE )
        tSpan(tSpan < tCurrent) = [];
        solvePDE(tSpan);
    end

end

if ( doPlot && ~doHold )
    hold off
end

% Ensure tOut is a column vector
tOut = tOut(:);

uOut = prepare4output(uOut);

    function uOut = prepare4output(uIn)
        % If we only had one dependent variable, return an array valued CHEBFUN
        % instead of a QUASIMATRIX.
        if ( SYSSIZE == 1 )
            uOut = horzcat(uIn{:});
        else
            blocks = cell(SYSSIZE, numel(uIn));
            for kk = 1:SYSSIZE
                blocks(kk,:) = cellfun(@(u) extractColumns(u, kk), uIn, ...
                    'UniformOutput', false);
            end
            uOut = chebmatrix(blocks); % CHEBMATRIX
        end
    end

switch nargout
    case 0
    case 1
        if ( length(tIn) == 2 )     % Only T_start and T_end were given.
            uOut = uOut(:,[1,end]); % Strip intermediate steps.
        end
        varargout{1} = uOut;
    case 2
        varargout{1} = tOut;
        varargout{2} = uOut;
    case 1 + size(uOut, 1)
        varargout{1} = tOut;
        for varCounter = 1:size(uOut, 1)
            varargout{varCounter + 1} = chebfun(uOut(varCounter,:)); %#ok<AGROW>
        end
    otherwise
        error('CHEBFUN:CHEBFUN:pdeSolve:nargout', ...
            'Incorrect number of output arguments.');
end

clear global

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  ONESTEP  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function solvePDE(tSpan)
        % Constructs the result of one time chunk at fixed discretization.

        % Evaluate the chebfun at discrete points:
        x = points(currentLength, DOMAIN);
        U0 = feval(uCurrent, x);

        if ( ISPERIODIC )

            % The discretisation length
            n = currentLength;

        elseif ( isempty(n) || (n ~= currentLength) )
            % This depends only on the size of n. If this is the same, reuse!

            % The new discretisation length
            n = currentLength;

            % Linear constraints:
            if ( numLinBCs > 0 )
                bcop = [bc.left.op ; bc.middle.op ; bc.right.op];
                B = cell2mat(cellfun(@(f) feval(f, n), bcop, 'UniformOutput', false));
            end

            % Zero matrix for algebraic (boundary) terms:
            Z = zeros(numBCs, SYSSIZE*n);
            % Initialise evaluated BC rhs vector:
            q = zeros(numBCs, 1);

            % Project / mass matrix.
            P = cell(1, SYSSIZE);
            M = cell(1, SYSSIZE);

            for kk = 1:SYSSIZE
%                 xk = chebpts(n-DIFFORDER(kk), DOMAIN, 1);
                xk = chebpts(n-DIFFORDER(kk)+2, DOMAIN, 2); xk([1,end]) = [];
                P{kk} = barymat(xk, x);
                M{kk} = pdeFlag(kk)*P{kk};
            end
            P = [ Z ; blkdiag(P{:}) ];
            M = [ Z ; blkdiag(M{:}) ];

            % Multiply by user-defined mass matrix
            if ( userMassSet )
                M = feval(userMass, n)*M;
            end

            % ODE options: (mass matrix)
            opt = odeset(opt, 'Mass', M, 'MassSingular', 'yes', ...
                'MStateDependence', 'none');

        end

        if ( ~ISPERIODIC )
            % We have to ensure the starting condition satisfies the boundary
            % conditions, or else we will get a singularity in time. We do this
            % by tweaking the BC values to match the reality. Changing the IC
            % itself is much trickier. Find out what the BC deviance from
            % nominal really is:
            BCVALOFFSET = 0;             % recover nominal value in next call
            F = odeFun(tSpan(1), U0(:)); % also assigns values to "q"

            % If this is for the initial chunk, check whether the initial
            % condition nearly satisfies the BCs. We're quite lax about this,
            % because discretization at low N can cause derivatives to look
            % fairly bad.
            BCrows = 1:numBCs;
            if ( throwBCwarning && (length(uOut) > 1) && ...
                    (norm(F(BCrows)) > 0.05*norm(F)) )
                warning('CHEBFUN:CHEBFUN:pde15s:BadIC',...
                    'Initial state may not satisfy the boundary conditions.')
                throwBCwarning = false;
            end
            if ( adjustBCs )
                BCVALOFFSET = F(BCrows) - q;
            end
        end

        % Solve ODE over time chunk with the selected solver:
        try
            [~, ~] = ODESOLVER(@odeFun, tSpan, U0, opt);
        catch ME
            if ( strcmp(ME.identifier, 'MATLAB:odearguments:SizeIC') )
                error('CHEBFUN:CHEBFUN:pde15s:dims', ...
                    'Dimension mismatch. Check boundary conditions.');
            else
                rethrow(ME)
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%  ODEFUN  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function F = odeFun(t, U)
            % This is what the ODESOLVER() calls.

            % Reshape to n by SYSSIZE:
            U = reshape(U, n, SYSSIZE);

            % Evaluate the PDEFUN:
            myU = mydouble(U, DOMAIN);
            F = double(pdeFun(t, x, myU));
            F = F(:);

            if ( ISPERIODIC )
                return
            end

            % Enforce boundary constraints:

            % Project:
            F = P*F;

            % Replacements for the BC algebraic conditions:
            if ( numLinBCs > 0 )
                % Get the right-hand sides: (may be time-dependent)
                for l = 1:numLinBCs
                    if ( isa(BCRHS{l}, 'function_handle') )
                        q(l,1) = feval(BCRHS{l}, t);
                    else
                        q(l,1) = BCRHS{l};
                    end
                end
                indx = 1:numLinBCs;
                F(indx) = B*U(:) - q(indx);
            else
                indx = 0;
            end

            % Replacements for the nonlinear BC conditions:
            if ( numLBCs > 0 )
                indx = indx(end) + (1:numLBCs);
                fTmp = double(opLBC(t, x, myU));
                F(indx) = fTmp(1, :);
            end
            if ( numMBCs > 0 )
                % TODO: This won't work if there are also left nonlin BCs
                indx = indx(end) + (1:numMBCs);
                F(indx) = double(opMBC(t, x, myU));
            end
            if ( numRBCs > 0 )
                indx = indx(end) + (1:numRBCs);
                fTmp = double(opRBC(t, x, myU));
                F(indx) = fTmp(end,:);
            end

            % Adjust BC rows by the needed offset from original time:
            indx = 1:numBCs;
            F(indx) = F(indx) - BCVALOFFSET;

        end
    end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%  MISC  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [outFun, varNamesParsed] = parseFun(inFun, bcFlag)
% Rewrites the input function handle to call the right DIFF, SUM, methods, etc,
% and convert the input @(t, x, u, v, w, ...) to @(t, x, U).
%
% For PARSEFUN(INFUN), if the INFUN has only one independent variable as an
% input, we assume this is x. PARSEFUN(INFUN, 'lbc') or PARSEFUN(INFUN, 'rbc')
% will assume the variable is t. PARSEFUN(INFUN, 'bc') will assume space, but
% throw a warning.

global SYSSIZE

% Ensure backwards compatibility by processing the input function.
[inFun, varNamesParsed] = backCompat(inFun);

% V4: We don't accept only time or space as input args (i.e., both or nothing).
% V5: Actually, we now accept @(u, x) diff(u, 2) + sin(x).*u (for CHEBOPs).
Nind = nargin(inFun) - SYSSIZE;
if ( Nind == 0 )
    outFun = @(t, x, u) conv2cell(inFun, u);
elseif ( Nind == 1 )
    if ( nargin > 1 )
        if ( strcmp(bcFlag, 'bc') )
            warning('CHEBFUN:CHEBFUN:pde15s:bcInput', ...
                ['Please input both time and space independent variable to ' ...
                 '.BC function handle inputs.\nAssuming given variable is ' ...
                 'time.']);
             outFun = @(t, x, u) conv2cell(inFun, u, x);
        else
            outFun = @(t, x, u) conv2cell(inFun, u, t);
        end
    else
        outFun = @(t, x, u) conv2cell(inFun, u, x);
    end
elseif ( Nind == 2 )
    outFun = @(t, x, u) conv2cell(inFun, u, t, x);
else
    error('CHEBFUN:CHEBFUN:pde15s:inputs_ind', ...
        ['Incorrect number of independent variables in input function. ' ...
         '(Must be 0, 1, or 2).']);
end

end

function newFun = conv2cell(oldFun, u, varargin)
% This function allows the use of different variables in the anonymous function,
% rather than using the clunky quasi-matrix notation.
%
% Inputs:
%   oldFun in the original function handle
%   u will be the solution we evaluate at
%   varargin possibly contains the time and space variables
% Output
%   newFun the modified function handle

global SYSSIZE
% Note. This is faster than mat2cell!
tmpCell = cell(1, SYSSIZE);
for qk = 1:SYSSIZE
    tmpCell{qk} = extractColumns(u, qk);
end
newFun = oldFun(varargin{:}, tmpCell{:});
end

function getDIFFORDER(pdeFun, DOMAIN)
% Extract the DIFFORDER by evaluating the operator at some NaN values.
global DIFFORDER SYSSIZE
% Find the order of each of the variables:
tmp = chebdouble(zeros(SYSSIZE), DOMAIN);
v = pdeFun(0, 0, tmp);
DIFFORDER = get(v, 'diffOrder');
DIFFORDER = max(DIFFORDER, 0);
end

function out = dealWithStructInput(in)
% Parse the inputs for BCs.
% Throw a warning or an error for a struct. Convert numeric value to 'dirichet'.
if ( isstruct(in) )
    if ( numel(in) == 1)
        warning('CHEBFUN:CHEBFUN:pde15s:bcstruct', ...
            'PDE15S no longer supports struct inputs for bc.left and bc.right.')
        if ( isfield(in, 'val') )
            out = {in.op, in.val};
        else
            out = in.op;
        end
    else
        error('CHEBFUN:CHEBFUN:pde15s:bcstruct', ...
            'PDE15S no longer supports struct inputs for bc.left and bc.right.')
    end
elseif ( isempty(in) )
    out = [];
elseif ( isnumeric(in) )
    out = {'dirichlet', in};
else
    out = in;
end
end

function [outFun, varNamesOut] = backCompat(inFun)
% In V4 PDE15S required the user to pass in dummy function handles for the
% differential operators. For example, u'' + u' + sum(u) would have needed
%  pdefun = @(u, t, x, diff, sum) diff(u, 2) + diff(u) + sum(u).
% V5 deals with this in a different way (by using the CHEBDOUBLE class), but we
% would still like to support the old syntax (at least for now).
%
% To do this, we parse the input function for the strings 'diff', 'Diff', and
% 'D'. If we find any of these, we assume the old syntax is being used, and
% redefine the function handle appropriately. Of course, this is not fool-proof,
% but it should work most of the time.

global SYSSIZE

% Determine the variable names:
str = func2str(inFun);
parLoc = strfind(str, ')');
varList = str(3:parLoc(1)-1);
commaLoc = strfind(varList, ',');
varNames = cell(1, numel(commaLoc)+1);
for k = 1:numel(commaLoc)
    varNames{k} = varList(1:commaLoc(k)-1);
    varList(1:commaLoc(k)) = [];
    commaLoc = commaLoc - commaLoc(k);
end
varNames{numel(commaLoc)+1} = varList;

% Look for 'diff', 'Diff', or 'D':
diffLoc = cellfun(@(v) any(strcmp(v, {'diff', 'Diff', 'D'})), varNames);
idx = find(diffLoc);
if ( isempty(idx) )
    % None present - new syntax!
    outFun = inFun;
    varNamesOut = varNames;
    return
end

warning('CHEBFUN:CHEBFUN:pde15s:oldSyntax', ...
    ['The syntax\n ' str '\nis deprecated and may not be supported in ', ...
     'future releases.\n Please see the PDE15S documentation for details.'] )

% Switch the order for V5 syntax:
varNamesNewOrder = varNames([(1+SYSSIZE):(idx-1), 1:SYSSIZE]);

% Get the new list of variables:
varList1 = varNamesNewOrder{1};
varList2 = varNames{1};
for k = 2:idx-1
    varList1 = [varList1, ',', varNamesNewOrder{k}]; %#ok<AGROW>
    varList2 = [varList2, ',', varNames{k}]; %#ok<AGROW>
end

% Compile the new function string:
funList = {'@diff', '@sum', '@cumsum', '@fred'};
newStr = ['@(', varList1, ')' 'inFun(' varList2];
for k = 0:(nargin(inFun) - idx)
    newStr = [newStr, ',', funList{k+1}]; %#ok<AGROW>
end
newStr = [newStr, ')'];

% Make the new function handle:
outFun = eval(newStr);

% Return the obtained variable names
varNamesOut = varNamesNewOrder;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%
%
% DIFFERENTIAL SEARCH ALGORITHM (DSA) (in MATLAB)
% STANDARD VERSION of DSA (16.July.2013)
%
%
% usage : > ds(method,fnc,mydata,popsize,dim,low,up,maxcycle)
%
% method
%--------------------------------------------------------------
% 1: Bijective DSA (B-DSA)
% 2: Surjective DSA (S-DSA)
% 3: Elitist DSA (strategy 1) (E1-DSA)
% 4: Elitist DSA (strategy 2) (E2-DSA)
% if method=[. . . ...],   Hybrid-DSA (H-DSA)
%--------------------------------------------------------------
% example :
% ds(1,'circlefit',mydata,10,3,-10,10,2000) ; % B-DSA
% ds(2,'circlefit',mydata,10,3,-10,10,2000) ; % S-DSA
% ds(3,'circlefit',mydata,10,3,-10,10,2000) ; % E1-DSA
% ds(4,'circlefit',mydata,10,3,-10,10,2000) ; % E2-DSA
% ds([1 2],'circlefit',mydata,10,3,-10,10,2000) ; % Hybrid-DSA, in this case B-DSA and S-DSA are hybridized.
%--------------------------------------------------------------
% Please cite this article as;
% P.Civicioglu, "Transforming geocentric cartesian coordinates to geodetic coordinates by using differential search algorithm",  Computers & Geosciences, 46 (2012), 229-247.
% P.Civicioglu, "Understanding the nature of evolutionary search algorithms", Additional technical report for the project of 110Y309-Tubitak,2013, Ankara, Turkey.
%
%
%
%
%--------------------------------------------------------------
%{
19.March.2013
Copyright Notice
Copyright (c) 2012, Pinar Civicioglu
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 copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution

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

function ds(method,fnc,mydata,size_of_superorganism,size_of_one_clan,low_habitat_limit,up_habitat_limit,epoch)

% size_of_superorganism ; size of population.
% size_of_one_clan ; size of problem dimension (1,2,3,...,d), where each clan (i.e. sub-superorganism) includes d-individuals.
% mydata ; additional parameters for objective function, use mydata=[], if it is not needed. See Best-Fit Circle example (circlefit.m) for usage of mydata.

%Initialization

% control of habitat limits
if numel(low_habitat_limit)==1,
    low_habitat_limit=low_habitat_limit*ones(1,size_of_one_clan);
    up_habitat_limit=up_habitat_limit*ones(1,size_of_one_clan);
end

% generate initial individuals, clans and superorganism.
superorganism=genpop(size_of_superorganism,size_of_one_clan,low_habitat_limit,up_habitat_limit);
% success of clans/superorganism
fit_superorganism=feval(fnc,superorganism,mydata);

for epk=1:epoch

    % SETTING OF ALGORITHMIC CONTROL PARAMETERS
    % Trial-pattern generation strategy for morphogenesis; 'one-by-one morphogenesis'.
    % p1=0.0*rand;  % i.e.,  0.0 <= p1 <= 0.0
    % p2=0.0*rand;  % i.e.,  0.0 <= p2 <= 0.0

    % Trial-pattern generation strategy for morphogenesis; 'one-or-more morphogenesis'. (DEFAULT)
    p1=0.3*rand;  % i.e.,  0.0 <= p1 <= 0.3
    p2=0.3*rand;  % i.e.,  0.0 <= p2 <= 0.3

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

   [direction,msg]=generate_direction(method(randi(numel(method))),superorganism,size_of_superorganism,fit_superorganism);

   map=generate_map_of_active_individuals(size_of_superorganism,size_of_one_clan,p1,p2);

  %-------------------------------------------------------------------
    % Recommended Methods for generation of Scale-Factor; R
    % R=4*randn;  % brownian walk
    % R=4*randg;  % brownian walk
    % R=lognrnd(rand,5*rand);  % brownian walk
     R=1./gamrnd(1,0.5);   % pseudo-stable walk
    % R=1/normrnd(0,5);    % pseudo-stable walk

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

    % bio-interaction (morphogenesis)
    stopover=superorganism+(R.*map).*(direction-superorganism);

   % Boundary Control
    stopover=update(stopover,low_habitat_limit,up_habitat_limit);

    % Selection-II
    fit_stopover=feval(fnc,stopover,mydata);
    ind=fit_stopover<fit_superorganism;
    fit_superorganism(ind)=fit_stopover(ind);
    superorganism(ind,:)=stopover(ind,:);

   % update results
    [globalminimum,indexbest]=min(fit_superorganism);
    globalminimizer=superorganism(indexbest,:);

    % export results
    assignin('base','globalminimum',globalminimum);
    assignin('base','globalminimizer',globalminimizer);
    fprintf('%s  | %5.0f   --->   %10.16f\n',msg,epk,globalminimum)

end

function pop=genpop(a,b,low,up)
pop=ones(a,b);
for i=1:a
    for j=1:b
        pop(i,j)=rand*(up(j)-low(j))+low(j);
    end
end

function p=update(p,low,up)
[popsize,dim]=size(p);
for i=1:popsize
    for j=1:dim
        % first (standard)-method
        if p(i,j)<low(j), if rand<rand, p(i,j)=rand*(up(j)-low(j))+low(j); else p(i,j)=low(j); end, end
        if p(i,j)>up(j),  if rand<rand, p(i,j)=rand*(up(j)-low(j))+low(j); else p(i,j)=up(j); end, end

        %{
       %  second-method
        if rand<rand,
            if p(i,j)<low(j) || p(i,j)>up(j), p(i,j)=rand*(up(j)-low(j))+low(j); end
        else
            if p(i,j)<low(j), p(i,j)=low(j); end
            if p(i,j)>up(j),  p(i,j)=up(j); end
        end
        %}
    end
end

function [direction,msg]=generate_direction(method,superorganism,size_of_superorganism,fit_superorganism);
 switch method
        case 1,
            % BIJECTIVE DSA  (B-DSA) (i.e., go-to-rnd DSA);
            % philosophy: evolve the superorganism (i.e.,population) towards to "permuted-superorganism (i.e., random directions)"
            direction=superorganism(randperm(size_of_superorganism),:); msg=' B-DSA';
        case 2,
            % SURJECTIVE DSA (S-DSA) (i.e., go-to-good DSA)
            % philosophy: evolve the superorganism (i.e.,population) towards to "some of the random top-best" solutions
            ind=ones(size_of_superorganism,1);
            [null_,B]=sort(fit_superorganism);
            for i=1:size_of_superorganism, ind(i)=B(randi(ceil(rand*size_of_superorganism),1)); end;
            direction=superorganism(ind,:);  msg=' S-DSA';
        case 3,
            % ELITIST DSA #1 (E1-DSA) (i.e., go-to-best DSA)
            % philosophy: evolve the superorganism (i.e.,population) towards to "one of the random top-best" solution
            [null,jind]=sort(fit_superorganism); ibest=jind(ceil(rand*size_of_superorganism)); msg='E1-DSA';
            direction=repmat(superorganism(ibest,:),[size_of_superorganism 1]);
        case 4,
            % ELITIST DSA #2 (E2-DSA) (i.e., go-to-best DSA)
            % philosophy: evolve the superorganism (i.e.,population) towards to "the best" solution
            [null_,ibest]=min(fit_superorganism); msg='E2-DSA';
            direction=repmat(superorganism(ibest,:),[size_of_superorganism 1]);
 end
return

function map=generate_map_of_active_individuals(size_of_superorganism,size_of_one_clan,p1,p2);
    % strategy-selection of active/passive individuals
    map=zeros(size_of_superorganism,size_of_one_clan);
        if rand<rand,
            if rand<p1,
                % Random-mutation #1 strategy
                for i=1:size_of_superorganism
                    map(i,:)=rand(1,size_of_one_clan) < rand;
                end
            else
                % Differential-mutation strategy
                for i=1:size_of_superorganism
                    map(i,randi(size_of_one_clan))=1;
                end
            end
        else
             % Random-mutation #2 strategy
            for i=1:size_of_superorganism
                map(i,randi(size_of_one_clan,1,ceil(p2*size_of_one_clan)))=1;
            end
        end
return

author: Sable

repo: mcbench-benchmarks


text:
function [U2] = mle_lam1(X,phi2,v1,PP,PnP,n,T);
%
%
% this function calculates the maximum likelihood function U2
% (lambda) with 1 Fourier harmonic
%
b=phi2+X(1)*sin(n/T+X(2));
mb=kron(b,v1);   % produces a matrix with identical lines
A=PnP.*log(mb);
B=PP.*mb;
%
U2=sum(sum(A))-sum(sum(B));
U2=-U2;


text:
function VBA_setup()
% VBA_setup()
% Cleanly add the VBA_toolbox to the Matlab path.
% Please prefer this function to addpath() or addpath(genpath()) to ensure
% all subfolders are included with the notable exception of the .git tracker.

fprintf('\n');
fprintf(' ======================================================================================\n')
fprintf(' VBA-toolbox installation \n')
fprintf(' ======================================================================================\n')
fprintf('\n');

%%

isWin = ~isempty(strfind(computer, 'PCWIN')) || strcmp(computer, 'i686-pc-mingw32');
isOSX = ~isempty(strfind(computer, 'MAC')) || ~isempty(strfind(computer, 'apple-darwin'));
isLinux = strcmp(computer,'GLNX86') || strcmp(computer,'GLNXA64') || ~isempty(strfind(computer, 'linux-gnu'));

%%
% Locate ourselves:
targetdirectory=fileparts(mfilename('fullpath'));
if ~strcmpi(targetdirectory, pwd)
    fprintf(' *** You need to change your working directory to the VBA-toolbox folder before running this routine!\n');
    return
else
    fprintf(' Will setup %s as the working copy of the VBA-toolbox folder.\n',targetdirectory);

end
fprintf('\n');

%%
% Does SAVEPATH work?
err=savepath;

if err
    % warn user the new path wont be saved:
    p=fullfile(matlabroot,'toolbox','local','pathdef.m');

    fprintf(' *** savepath() failed. Probably the pathdef.m file lacks write permission. \n\n');

    fprintf('     You will have to re-run %s again after each Matlab startup untill you fix this permission issue.\n',mfilename);

    fprintf('     Please ask a user with administrator privileges to enable write by everyone for the file:\n')
    fprintf('     ''%s''\n\n',p);

    fprintf(' // Installation will continue but your changes will not be saved.\n')

end

%%
% Remove old "VBA-toolbox" from path:
if isOSX || isLinux
    p = strsplit(path,':');
elseif isWin
    p = strsplit(path,';');
end

r = regexp(p, '^.*VBA-toolbox.*$','match');
r = r(~cell2mat(cellfun(@isempty,r,'UniformOutput',false)));

if ~isempty(r)
    try
        % try to get version
        ver = VBA_version();
    catch
        % find root path
        pl = cell2mat(cellfun(@(x) numel(x{1}),r,'UniformOutput',false));
        ridx = find(pl == min(pl));

        ver.path = r{ridx};
        ver.version = 'unknown';
    end

    fprintf(' *** The VBA-toolbox already appears in your path:\n');
    fprintf('     - installation path: %s\n',ver.path);
    fprintf('     - release: %s\n',ver.version);

    fprintf('     Installation will remove this toolbox from your path (without deleting the files).\n');

    cont = NaN;
    while isnan(cont)
        answer = input('\n     Are you sure you want to continue (yes/no)? ','s');
        switch answer
            case 'yes'
                fprintf('\n //  Alright, let''s go...\n')
                cont = true;
            case 'no'
                fprintf('\n //  You''re the boss! \n');
                cont = false;
            otherwise
                fprintf('\n //  I didn''t understand, please try again.\n');
        end
    end

    if ~cont
        fprintf('\n Run %s again if you change your mind.\n\n',mfilename);
        return
    end

    % remove all old path

    for iP = 1:numel(r)
        rmpath(str2mat(r{iP}))
    end
    try
        savepath
    end
end

%% find folder to install
p = strsplit(genpath(pwd),':');
ref_path = ['refs' filesep 'heads'];
r = regexp(p, ['^.*VBA-toolbox' filesep '.git' filesep '(?!' ref_path ').*$'],'match');
p = p(cell2mat(cellfun(@isempty,r,'UniformOutput',false)));

for iP = 1:numel(p)
    addpath(str2mat(p{iP}))
end
try
    savepath
end

fprintf('\n')
try
    ver = VBA_version();
    fprintf(' The VBA-toolbox has been successfully installed in %s\n\n',ver.path);
    fprintf(' Type demo_Qlearning to give it a try\n\n');
    fprintf(' ======================================================================================\n');

catch err
    fprintf(' *** Something wrong happened:\n')
    fprintf('     -  %s \n\n',err.message);
end

end

author: MBB-team

repo: VBA-toolbox