text:
function outputWaveforms = load_antelope(request, specificDatabase)
   % load a waveform from antelope
   %  w = load_antelope(request, specificDatabase)
   %   CRITERIA is the search criteria, created with buildAntelopeCriteria
   %   request.sTime is the startTimes
   %   endTimes is the endTimes
   %   combineWaves is a logical value:
   %     Should segmented waveforms be combined,(within requested timerange)?
   %   database is the antelope database

   % AUTHOR: Celso Reyes
   % MODIFICATIONS: Glenn Thompson, Carl Tape

   %TODO: maybe order of operations can be changed to avoid unpacking
   %request when unnecessary.

   [~, chanInfo, startTimes, endTimes, combineWaves] = unpackDataRequest(request);
   assert(numel(startTimes) == numel(endTimes),...
      'Waveform:load_antelope:startEndMismatch', 'Unequal number of start and end times');

   %call this routine for each database, then return the waveforms.
   if exist('specificDatabase','var')
      database = specificDatabase;
      outputWaveforms = emptyWaveform();

      [criteria, nCriteria] = buildAntelopeCriteria(chanInfo);
      for i = 1:nCriteria
         %if multiple traces will result, then there may be multiple records for tr
         [tr, database, fdb] = get_antelope_traces(...
            startTimes, endTimes, criteria(i).group, database);
         w = cycleThroughTraces(tr, combineWaves);
         outputWaveforms = [outputWaveforms; w]; %#ok<AGROW>
      end
      dbclose(fdb);
   else
      outputWaveforms = RecursivelyLoadFromEachDatabase(request);
   end
end

function outputWaveforms = RecursivelyLoadFromEachDatabase(request)
   TRY_MULTIDAY = false;
   database = getDatabase(request, TRY_MULTIDAY);
   database = unique(database, 'stable'); %  avoid changing requested order
   outputWaveforms = emptyWaveform();
   for n = 1 : numel(database)
      w = load_antelope(request, database(n));
      outputWaveforms = [outputWaveforms; w(:)]; %#ok<AGROW>
   end
end

function w = emptyWaveform()
   w = waveform;
   w = w([]);
end

function database = getDatabase(request, TRY_MULTIDAY)
   [ds, chanInfo, startTimes, endTimes] = unpackDataRequest(request);
   if TRY_MULTIDAY %good luck splicing
      dbDatesToCheck = subdivide_files_by_date(ds,startTimes, endTimes);
      database =  getfilename(ds,chanInfo, dbDatesToCheck);
   else
      database =  getfilename(ds,chanInfo, startTimes);
   end
end

function w = cycleThroughTraces(tr, COMBINE_WAVEFORMS)
   %cycleThroughTraces  converts each trace into one or more waveforms
   %   returns a Nx1 waveform
   w(numel(tr)).waves = waveform;
   for traceidx = 1:numel(tr)
      w(traceidx) = wavesfromtraces(tr, traceidx, COMBINE_WAVEFORMS);
   end
   w = [w.waves]; %horizontal cat
   w = w(:);
end

function w_scnl = wavesfromtraces(tr, traceidx, COMBINE_WAVEFORMS)
   if ~isstruct(tr{traceidx}) % marker for no data
      w_scnl.waves = emptyWaveform();
   else
      w_scnl.waves = traceToWaveform(tr{traceidx}); %create waveform list
      trdestroy(tr{traceidx});
   end
   if COMBINE_WAVEFORMS %combine all of this trace's records
      w_scnl.waves = combine(w_scnl.waves);
   end;
   w_scnl.waves = reshape(w_scnl.waves,  numel(w_scnl.waves), 1);
end

%% helper functions

function [critList, nCrit] = buildAntelopeCriteria(chanTag)
   % builds a criteria list from SCNL (Station-Channel-Network-Location)
   % returns a structure array of cells with criteria for use with get_antelope_trace
   % structure returned as critList(n).group, where group is a cell containing
   % the criteria, such as {'sta == "OKCF"','chan=="EHZ"'}
   %
   % ex. myscnls = scnlobject({'OKCF','OKFG','OKSO'},{'EHZ','BHZ','BHZ'});
   % ex.  cl =
   % builtAntelopeCriteria(myscnls)
   % ... CL will then contain
   %
   % CL(1).group(1).field = 'sta'
   % CL(1).group(1).relation = '=='
   % CL(1).group(1).data = 'OKCF'
   % CL(1).group(2).field = 'chan'
   % CL(1).group(2).relation = '=='
   % CL(1).group(2).data = 'EHZ'
   % ...
   % CL(N).group(4).field = 'loc'
   % CL(N).group(4).relation = '=='
   % CL(N).group(4).data = '--'
   %
   % so...    getCritListExpression(CL(1).group(1)) --> 'sta == "OKCF"'
   % It was done this way so that the dabase could be checked to ensure
   % that the view actually supports each search parameter.
   %
   % either multiple stations OR multiple channels may work... not both.
   % also, if any of the SCNL parameters are left empty, they will not be used
   % in subsetting the antelope database.
   %

   for N = numel(chanTag) : -1: 1
      if ~isempty(chanTag(N).station)
         critList(N).group(1) = grabCritList('sta', chanTag(N).station);
      end
      if ~isempty(chanTag(N).network)
         critList(N).group(2) = grabCritList('net', chanTag(N).network);
      end
      if ~isempty(chanTag(N).channel)
         critList(N).group(3) = grabCritList('chan', chanTag(N).channel);
      end
      if ~isempty(chanTag(N).location);
         critList(N).group(4) = grabCritList('loc', chanTag(N).location);
      end
   end

   critList = critList(:);
   nCrit = numel(critList);
end

function A = grabCritList(key, value)
   if ~isempty(value)
      value = regexprep(value,'(?<!\.)\*','\.\*'); %replace * with .*, but leave existing .* alone
      value = makecell(value);
      for i = numel(value) : -1 : 1
         A(i).field = key;
         A(i).relationship = '=~';
         A(i).value = value{i};
      end
   elseif ismember(key,{'sta','cha'})  %ew. get rid of special treatment
      erid = sprintf('Waveform:load_antelope:no%sTerm',key);
      error(erid,'No %s  requested. To retrieve all %ss use ''*''', key, key);
   else
      A.field='';
      A.relationship = '';
      A.value = '';
   end
end

function [tr, rawDb, filteredDb] =  get_antelope_traces(startdates, enddates, criteriaList, database)
   % GET_ANTELOPE_TRACE gets interesting data from database, and returns the tracebuf object
   % [tr(1:end)] =  get_antelope_trace(startdates, enddates, criteriaList, database)
   %    STARTDATE is a matlab datenum
   %    ENDDATE is also in matlab datenum
   %    CRITERIALIST is a cell of statements used to filter the antelope
   %        database.  example: {'sta=="OKCF"','chan=="BHZ"'}
   %    DATABASE can either be a database name, or an open antelope database
   %    pointer.
   %
   %    TR is a tracebuf object containing the data requested.
   %
   %  In this usage, with a single output argument, TR, the database is opened
   %  and then closed before the function returns.
   %
   % You can opt to keep the database open, which can make subsequent calls to
   % get_antelope_trace much faster.  To do this, ask for the database pointer
   % as one of the output arguments.  The database will be closed by any
   % function call that uses the database pointer as the 'database' parameter,
   % and does not ask for a database pointer back.
   %
   % [tr, mydb] = get_antelope_trace(...)
   %    if mydb is asked for, then the database will not be closed, allowing for
   %    more efficient reuse.  mydb is the unfiltered (and open) wfdisc table
   %
   % [tr, mydb, filteredDB] = get_antelope_trace(...)
   %   returns wfdisc table in mydb, and a filtered wfdisc table in filteredDB
   %
   % if no records are found, then tr will be set to []
   %
   % will generate error message if no records found, so consider using in a
   % TRY-CATCH loop.
   %
   %
   % Example:
   %  % set up the data
   %  sDate = datenum('7/12/2008 05:00:00.00');
   %  eDate = datenum('7/12/2008 05:10:00.00');
   %  stationsToGet = {'OKFG','OKSO','OMG'};
   %  channelsToGet = {'BHZ', 'BHE'};
   %  mydb = '/mydir/my_antelope_database';
   %
   %  % get our criteria together
   %  for st = 1:numel(stationsToGet)
   %    staCriteria(st) = {sprintf('sta=="%s"',stationsToGet(st)};
   %  end
   %  for ch = 1:numel(channelsToGet)
   %    chanCriteria(ch) = {sprintf('chan=="%s"',channelsToGet(ch)};
   %  end
   %
   %  % open the database once, and read in each tracebuf. Then translate into
   %  % waveforms.
   %  n = 1;
   %  for sta = 1:numel(staCriteria)
   %    for ch = 1:numel(chanCriteria)
   %      critera = [chanCriteria(ch),staCriteria(sta)];
   %      [tr, mydb] = get_antelope_trace(sDate,eDate, criteria, mydb)
   %      w(n) = traceToWaveform(tr);
   %      n = n+1;
   %    end
   %  end
   %  % finally, close up shop.
   %  dbclose(mydb)

   % Modifications
   %%% Glenn Thompson 2012/02/06: Occasionally the C program trload_css cannot even load the trace data. Added try...catch..end to handle this.

   useExistingDatabasePtr =  isAntelopeDatabasePtr(database);

   if useExistingDatabasePtr
      try
         dbnrecs(database);
      catch
         warning('Waveform:load_antelope:databaseNotOpen', ...
            'a Database Pointer was passed to trace, but the database was not open');
         %tr = trnew; %return a new object, forcing the ability to destroy it later
         tr = { -1 };
         return;
      end
   end

   % do not close the database if a pointer is asked for in the return arguments
   onFinishCloseDB = nargout < 2;

   antelope_starts = mep2dep(startdates);
   antelope_ends = mep2dep(enddates);

   %if the database isn't already open, then open it for reading
   if ~useExistingDatabasePtr
      mydb = dbopen(database,'r');
      mydb = dblookup_table(mydb,'wfdisc');
   else
      mydb = database;
   end
   %check to ensure wfdisk table exists and is populated

   if safe_dbnrecs(mydb) == 0,
      cleanUpFail('Waveform:load_antelope:databaseNotFound', ...
         'Database not found: %s', dbquery(mydb,'dbTABLE_FILENAME'));
      return;
   end;

   rawDb = mydb;   % keep a copy of the pre-subset (raw) database

   % subset the data based upon the desired criteria
   criteriaList = keepCriteriaThatMatchDatabaseFields(criteriaList, mydb);
   allExp = getAsExpressions(criteriaList);

   %subset the database based on this particular criterion
   mydb = dbsubset(mydb,allExp);
   if safe_dbnrecs(mydb) == 0
      cleanUpFail('Waveform:load_antelope:dataNotFound', 'No records found for criteria [%s].', allExp);
      return;
   end;

   filteredDb = mydb;

   [st, ed] = dbgetv(mydb,'time','endtime');
   %% Get the tracebuf object for this starttime, endtime
   % Loop through all times.  Result is tr(1:numel(starttimes) of all tracebuffers.
   for mytimeIDX = 1:numel(antelope_starts)
      someDataExists = any(antelope_starts(mytimeIDX)<= (ed) & antelope_ends(mytimeIDX) >= (st));
      if someDataExists
         %%% Glenn Thompson 2012/02/06: Occasionally the C program trload_css cannot even load the trace data.
         % This error needs to be handled. So adding a try..catch..end around the original instruction.
         try
            tr{mytimeIDX} = trload_css(mydb, antelope_starts(mytimeIDX), antelope_ends(mytimeIDX)); %#ok<AGROW>
         catch
            cleanUpFail('Waveform:load_antelope:trload_css failed', ...
               'Database not found: %s', dbquery(mydb,'dbTABLE_FILENAME'));
            disp(database)
            disp(allExp)
            starttimes = antelope_starts(mytimeIDX);
            endtimes = antelope_ends(mytimeIDX);
            fprintf('%.0f %.0f\n ',starttimes,endtimes);
            disp(['starttimes: ', datestr(epoch2datenum(starttimes))]);
            disp(['endtimes: ', datestr(epoch2datenum(endtimes))]);
            fprintf('trload_css(mydb, starttimes, endtimes))\n');
            return
         end
         trsplice(tr{mytimeIDX},20);
      else
         tr{mytimeIDX} = -1; %#ok<AGROW>
      end
   end %mytimeIDX
   closeIfAppropriate(mydb, onFinishCloseDB);

   function cleanUpFail(varargin)
      %cleanUpFail   tidy up database and records and display warning
      closeIfAppropriate(mydb, onFinishCloseDB);
      tr = { -1 };
      filteredDb = dbinvalid;
      warning(varargin{:});
   end
end

function  critList = keepCriteriaThatMatchDatabaseFields(critList, mydb)
   dbFields = dbquery(mydb, 'dbTABLE_FIELDS');
   critFields = {critList.field};  % all criteria we've parsed
   validCritIdx = true(size(critFields)); % preallocation
   for n = 1:numel(critFields)
      validCritIdx(n) = ...                %make sure:
         ~isempty(critFields{n}) && ...    % isn't empty and
         ismember(critFields{n}, dbFields);%  is a database field
   end
   critList = critList(validCritIdx);   %keep only good ones
end
function allExp = getAsExpressions(criteria)
   for n=1:numel(criteria)
      eachExp(n) = crit2expression(criteria(n)); %#ok<AGROW>
   end
   allExp = eachExp{1};
   for i= 2 : (numel(eachExp))
      allExp = [allExp,' && ', eachExp{i}]; %#ok<AGROW>
   end
end

function cle = crit2expression(cl)
   if ischar(cl.value)
      cle = {sprintf('%s%s/%s/',cl.field, cl.relationship, cl.value)};
   else
      cle =  {sprintf('%s %s %s',cl.field,cl.relationship,num2str(cl.value))};
   end
end

function n = safe_dbnrecs(mydb)
   %safe_dbnrecs   return either the number of records or 0
   %   n = safe_dbnrecs(mydb)  when accessing dbnrecs fails, this returns 0
   try
      n = dbnrecs(mydb);
   catch
      n = 0;
   end
end

function closeIfAppropriate(mydb, closeDatabase)
   if closeDatabase
      dbclose(mydb);
   end
end

function result = isAntelopeDatabasePtr(database)
   result = isa(database,'struct');
   result = result && any(strcmpi(fieldnames(database),'database'));
end

function X = makecell(X)
   if ~iscell(X)
      X = {X};
   end
end

function allw = traceToWaveform(tr)
   %traceToWaveform converts traceobjects to waveforms
   %    w = traceToWaveform(blankwaveform, all_tracobjects)
   %
   % Note: this may return multiple waveform objects, depending upon
   % how many segments and/or scnl's.

   persistent wBlank           % cannot test if wBlank is empty because
   persistent blankWaveExists  %   isempty(wBlank) is always true.

   if isempty(blankWaveExists)
      wBlank = waveform;
      wBlank =  addfield(wBlank,'calib', 0);
      wBlank = addfield(wBlank, 'calibration_applied', 'NO');
      blankWaveExists = true;
   end

   try
      traceCount = dbnrecs(tr);
      assert(traceCount > 0);
   catch
      allw = wBlank([]);
      return
   end

   % preallocations
   spikes(traceCount).mask = false; %used to track data spikes & inf values
   allw = repmat(wBlank,traceCount,1);

   maxAllowableSignal = (realmax('single') * 1e-2);

   % LOOP twice through segments represented by this trace object
   % 1st loop: find signal spikes, and get header info.
   for seg = 1:traceCount
      tr.record = seg - 1;
      s = db2struct(tr); %do once, get one for each segment
      [sunit, ~] = segtype2units(s.segtype);
      allw(seg) = set(wBlank, ...
         'station', s.sta, ...
         'channel', s.chan, ...
         'network', s.net, ...
         ... 'location', s.loc, ... % unknown if 'loc' really is a field
         'start', dep2mep(s.time), ...
         'freq', s.samprate, ...
         'units', sunit, ...
         'calib', s.calib);

      % data spikes must be known PRIOR to applying calibration
      data = trextract_data(tr);
      spikes(seg).mask =(abs(data) >= maxAllowableSignal) | isinf(data);
   end

   % now, apply calibrations to all traces at once
   trapply_calib(tr);
   hasCalibs = get(allw,'calib') ~= 0;
   allw(hasCalibs) = set(allw(hasCalibs),'calibration_applied','YES');

   % 2nd loop: assign data to the waveforms.
   replaceWithNan = @(W,BAD) setsamples(W, BAD.mask, nan);
   for seg = 1:traceCount
      tr.record = seg - 1;
      allw(seg) = set(allw(seg), 'data', trextract_data(tr));
      allw(seg) = replaceWithNan(allw(seg), spikes(seg));
   end
end

function [units, data_description] = segtype2units(unitCode)
   %segtype2units   retrieves unit and description based on the segtype code
   %'segtype' in antelope datasets indicate the natural units of the detector
   persistent codeKey
   if isempty(codeKey)
      codeKey = createCodeKey();
   end
   if codeKey.isKey(unitCode)
      details = codeKey(unitCode);
      units = details{1};
      data_description = details{2};
   else
      [units, data_description] = deal('null');
   end

   function SU = createCodeKey()
      %createSegUnits   creates a map  where SU(char) = {units, data_description}
      SU = containers.Map;
      SU('A') = {'nm / sec / sec','acceleration'};
      SU('B') = {'25 mw / m / m','UV (sunburn) index(NOAA)'};
      SU('D') = {'nm', 'displacement'};
      SU('H') = {'Pa','hydroacoustic'};
      SU('I') = {'Pa','infrasound'};
      SU('J') = {'watts','power (Joulses/sec) (UCSD)'};
      SU('K') = {'kPa','generic pressure (UCSB)'};
      SU('M') = {'mm','Wood-Anderson drum recorder'};
      SU('P') = {'mb','barometric pressure'};
      SU('R') = {'mm','rain fall (UCSD)'};
      SU('S') = {'nm / m','strain'};
      SU('T') = {'sec','time'};
      SU('V') = {'nm / sec','velocity'};
      SU('W') = {'watts / m / m', 'insolation'};
      SU('a') = {'deg', 'azimuth'};
      SU('b') = {'bits/ sec', 'bit rate'};
      SU('c') = {'counts', 'dimensionless integer'};
      SU('d') = {'m', 'depth or height (e.g., water)'};
      SU('f') = {'micromoles / sec / m /m', 'photoactive radiation flux'};
      SU('h') = {'pH','hydrogen ion concentration'};
      SU('i') = {'amp','electric curent'};
      SU('m') = {'bitmap','dimensionless bitmap'};
      SU('n') = {'nanoradians','angle (tilt)'};
      SU('o') = {'mg/l','diliution of oxygen (Mark VanScoy)'};
      SU('p') = {'percent','percentage'};
      SU('r') = {'in','rainfall (UCSD)'};
      SU('s') = {'m / sec', 'speed (e.g., wind)'};
      SU('t') = {'C','temperature'};
      SU('u') = {'microsiemens/cm','conductivity'};
      SU('v') = {'volts','electric potential'};
      SU('w') = {'rad / sec', 'rotation rate'};
      SU('-') = {'null','null'};
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function sfun_keyboard_input_v1_2(block)
% modified from sfun_keyboard_input_v1_01.m of Marc Compere by Emanuele Ruffaldi
% created : 17 June 2003
% modified: 20 June 2003
% created: 19 May 2009 => Level 2 and terminate
%

% Level-2 M file S-Function with keyboard
%   Copyright 1990-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.1 $

  setup(block);

%endfunction

function setup(block)

  block.NumDialogPrms  = 0;
  block.NumInputPorts  = 0;
  block.NumOutputPorts = 2;

  block.OutputPort(1).Dimensions       = 1;
  block.OutputPort(1).SamplingMode = 'Sample';
  block.OutputPort(1).DatatypeID  = 0;
  block.OutputPort(2).Dimensions       = 1;
  block.OutputPort(2).SamplingMode = 'Sample';
  block.OutputPort(2).DatatypeID  = 0;
  block.SampleTimes = [-1, 0];
%end

  %% Register methods
%  block.RegBlockMethod('CheckParameters', @CheckPrms);
  block.RegBlockMethod('PostPropagationSetup',    @DoPostPropSetup);
  block.RegBlockMethod('InitializeConditions',    @InitConditions);
  block.RegBlockMethod('Outputs',                 @Output);
  block.RegBlockMethod('Update',                  @Update);
  block.RegBlockMethod('Terminate',                 @Terminate);
%endfunction

function DoPostPropSetup(block)

  %% Setup Dwork
  block.NumDworks = 3;
  block.Dwork(1).Name = 'key';
  block.Dwork(1).Dimensions      = 1;
  block.Dwork(1).DatatypeID      = 0;
  block.Dwork(1).Complexity      = 'Real';
  block.Dwork(1).UsedAsDiscState = true;

  block.Dwork(2).Name = 'fig';
  block.Dwork(2).Dimensions      = 1;
  block.Dwork(2).DatatypeID      = 0;
  block.Dwork(2).Complexity      = 'Real';
  block.Dwork(2).UsedAsDiscState = false;

  block.Dwork(3).Name = 'new';
  block.Dwork(3).Dimensions      = 1;
  block.Dwork(3).DatatypeID      = 0;
  block.Dwork(3).Complexity      = 'Real';
  block.Dwork(3).UsedAsDiscState = true;
%endfunction

function InitConditions(block)

  %% Initialize Dwork
  block.Dwork(1).Data = 0;
   handle.figure=findobj('type','figure','Tag','keyboard input figure');

   if isempty(handle.figure)
      % 'position' args -> [left, bottom, width, height]
      handle.figure=figure('position',[100 100 400 200],...
                           'WindowStyle','Modal',...
                           'Name','Keyboard Input',...
                           'Color',get(0,'DefaultUicontrolBackgroundColor')); %,...
                           %'HandleVisibility','callback');
      %handle.figure=figure('position',[800 620 400 300]);
      %handle.figure=figure('position',[800 620 400 300],'WindowButtonDownFcn',@myCallback)
      %handle.figure=figure('position',[800 620 400 300],'WindowButtonMoveFcn',@myCallback_move,'WindowButtonDownFcn',@myCallback_clickdown)
      set(handle.figure,'Tag','keyboard input figure');

      % Make the OFF button (position args->[left bottom width height])
      handle.offbutton = uicontrol(handle.figure,...
          'Style','pushbutton',...
          'Units','characters',...
          'Position',[5 5 46 2],...
          'String','Disable exclusive figure-keyboard input',...
          'Callback',{@turn_modal_off,handle});

      % Make the ON button (position args->[left bottom width height])
      handle.onbutton = uicontrol(handle.figure,...
          'Style','pushbutton',...
          'Units','characters',...
          'Position',[5 1 46 2],...
          'String','Re-enable exclusive figure-keyboard input',...
          'Callback',{@turn_modal_on,handle});

   else, % reset the figure to 'modal' to continue accepting keyboard input
      set(handle.figure,'WindowStyle','Modal')
   end
  block.Dwork(2).Data = handle.figure;

%endfunction

function Output(block)

  block.OutputPort(1).Data = block.Dwork(1).Data;
  block.OutputPort(2).Data = block.Dwork(3).Data;

%endfunction

function Update(block)

   handle.figure = block.Dwork(2).Data;
   %handle = get(handle.figure,'userdata');

   current_char=get(handle.figure,'CurrentCharacter'); % a single character, like 'b'

   % update the grahics object
   %set(handle.point,'Xdata',[x(1) x(1)],'Ydata',[x(2) x(2)],'Zdata',[-1 +1]);

   % conditionally update the (numeric) state
   if ~isempty(current_char)
      block.Dwork(1).Data =double(current_char);
      block.Dwork(3).Data = 1;
	   % reset 'CurrentCharacter' so if user lifts up from key, this is noticed
	   set(handle.figure,'CurrentCharacter',char(0)) % the plus key is the only key that may be
   else,
      block.Dwork(3).Data = 0;
   end

                                   % pressed, but when the user stops, is not noticed

   % notes:
   %    - use sprintf() to convert string -> number contained in a string (or character) variable
   %    - use char() to convert floating point number -> ascii character
   %    - use str2num() to convert the string-number into a (double) floating point number
   %    - use str2num() with char() to convert a number-string into a char-string
   %         char(97) == char(str2num('97')) --> 'a'
   %    - use num2str() to convert a (double) number into the same number but contained
   %      in a string variable
   % For example:
   %    tmp='a';                        % assign a string
   %    tmp_num=sprintf('%i',tmp)       % convert that string into a number contained in a string variable, tmp_num
   %        (tmp_num is the string containing the characters '97')
   %    tmp_char=char(str2num(tmp_num)) % convert that string variable back into a number, then into the original string

%endfunction

% Callback for turning Modal OFF
function turn_modal_off(obj,eventdata,handle)
%disp('turn_modal_off:')
%handle
set(handle.figure,'WindowStyle','Normal')
%end

% Callback for turning Modal ON
function turn_modal_on(obj,eventdata,handle)
%disp('turn_modal_on:')
%handle
set(handle.figure,'WindowStyle','Modal')
%end

% % Callback for 'WindowButtonMoveFcn' in figure
% function myCallback_move(obj,eventdata)
% str=sprintf('\tWindowButtonMoveFcn callback executing');disp(str)
% end
%
% % Callback for 'WindowButtonDownFcn' in figure
% function myCallback_clickdown(obj,eventdata)
% str=sprintf('\t\tWindowButtonDownFcn callback executing');disp(str)
% end

function Terminate(block)
   handle.figure = block.Dwork(2).Data;
   if handle.figure ~= 0
       close (handle.figure);
   end

author: Sable

repo: mcbench-benchmarks


text:
function [CTTraceS, RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
   getPBRayData(edgeS, sourceS, numSamplePts, xySampleRate,planC, xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);
%JOD
% global planC
%
% indexS=planC{end};

water = 1000; %Assumes water equals 1000.

%-----------Get CT scan---------------------%

[CTUniform3D, CTUniformInfoS] = getUniformizedCTScan(1,planC);

xOffset = CTUniformInfoS.xOffset;
yOffset = CTUniformInfoS.yOffset;

%-----------Fix source characteristics---------------------%

orgV   = [sourceS.isocenter.x, sourceS.isocenter.y, sourceS.isocenter.z];

%-----------Get ray parameters---------------------%

% [RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV] = ...
%       getPBRays(edgeS, sourceS, xySampleRate);
[RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
      getPBRays(xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);

rayLength = 500;  %in cm.  This is the length of the rays which are passed
                  %through the CT matrix to determine cumulative CT densities.

numSlices = size(CTUniform3D,3);

zFirst = CTUniformInfoS.firstZValue;

sliceThickness = CTUniformInfoS.sliceThickness;

%minBox an maxBox containt coordinates of the corners of the CT box as used
%by the ray intersection routine.

maxBoxS.z = (numSlices - 1) * sliceThickness + zFirst + sliceThickness/2;

minBoxS.z = zFirst - sliceThickness/2;

delta_xy = CTUniformInfoS.grid1Units;

imageWidth = CTUniformInfoS.size(1);

minBoxS.x = - imageWidth/2 * delta_xy + xOffset;
maxBoxS.x =   imageWidth/2 * delta_xy + xOffset;

minBoxS.y = - imageWidth/2 * delta_xy + yOffset;
maxBoxS.y =   imageWidth/2 * delta_xy + yOffset;

imageWidth = CTUniformInfoS.sizeOfDimension1;

CTTraceS = struct('CTNumsRay',[],'CTCumNumsRay',[],'distSamplePts',[]);

for i = 1 : size(RTOGPBVectorsM,1)

  rayDeltaS.x = RTOGPBVectorsM(i,1) * rayLength;
  rayDeltaS.y = RTOGPBVectorsM(i,2) * rayLength;
  rayDeltaS.z = RTOGPBVectorsM(i,3) * rayLength;

  %are the components of the ray's direction and maximum length.

  deltaV = [rayDeltaS.x, rayDeltaS.y, rayDeltaS.z];

  t_entrance = rayBoxIntersection(sourceS,rayDeltaS,minBoxS,maxBoxS);

  %The entrance point is therefore
  entranceV = orgV + t_entrance * deltaV;

  %Reflect to find exit point (assume length of ray is long enough that ray does exit):

  if t_entrance ~= -1

    %find exit point
    %get end of ray
    rayOrgS2.xRel = sourceS.xRel + rayDeltaS.x;  %reflected source positions
    rayOrgS2.yRel = sourceS.yRel + rayDeltaS.y;
    rayOrgS2.zRel = sourceS.zRel + rayDeltaS.z;

    rayDeltaS2.x = - rayDeltaS.x;
    rayDeltaS2.y = - rayDeltaS.y;
    rayDeltaS2.z = - rayDeltaS.z;

    rayOrgS2.isocenter = sourceS.isocenter;

    t = rayBoxIntersection(rayOrgS2,rayDeltaS2,minBoxS,maxBoxS);

    t_exit = 1 - t;

    exitV = orgV + t_exit * deltaV;

else
    error('PB Ray does not intersect CT scan.');
end

  %Now produce a set of sampling points between the entrance and the exit

  nV = 1 : numSamplePts;

  delta_t  = (t_exit - t_entrance)/(numSamplePts - 1);

  tV = t_entrance + (nV - 1) * delta_t;

  CTTraceS(i).distSamplePts = tV * sum(deltaV.^2).^0.5;

  sampleV.x =  sourceS.xRel + tV * deltaV(1);

  sampleV.y =  sourceS.yRel + tV * deltaV(2);

  sampleV.z =  sourceS.zRel + tV * deltaV(3);

  sampleRTOGV.x = sampleV.x + sourceS.isocenter.x;
  sampleRTOGV.y = sampleV.y + sourceS.isocenter.y;
  sampleRTOGV.z = sampleV.z + sourceS.isocenter.z;

  %---------Sample CT densities----------%

  %To go from sample points in RTOG system to CT densities, we convert as follows:

  %What is the slice number?
  sliceV = 1 + (sampleRTOGV.z - zFirst)/sliceThickness;

  %Now do 3-D interpolation:
  zFieldV = [minBoxS.z + 0.5 * sliceThickness:sliceThickness:maxBoxS.z - 0.5 * sliceThickness];
  xFieldV = [minBoxS.x + 0.5 * delta_xy, delta_xy, maxBoxS.x - 0.5 * delta_xy];
  yFieldV = [minBoxS.y + 0.5 * delta_xy, delta_xy, maxBoxS.y - 0.5 * delta_xy];

  [CTNumsV] = finterp3(sampleRTOGV.x, sampleRTOGV.y, sampleRTOGV.z, CTUniform3D, xFieldV, yFieldV, zFieldV, 0);

%   [CTNumsV] = finterp3(sampleRTOGV.x, sampleRTOGV.y, sampleRTOGV.z, CTUniform3D, xFieldV, yFieldV, zFieldV, 0);

  CTTraceS(i).densityRay = delta_t * norm(deltaV) * CTNumsV/water;
  CTTraceS(i).cumDensityRay = cumsum(CTTraceS(i).densityRay);  %Account for sampling rate to convert to g/cm^2.

  if any(CTTraceS(i).cumDensityRay > 50)
      warning('Cumulative density ray appears to exceed maximum length.');
  end

end

%-----------fini---------------------%

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function new_struct = struct_strrep(old_struct, old_string, new_string, depth)
% Traverses depth-first through an entire structure, replacing
% old_string with new_string everywhere it goes
%
% :Usage:
% ::
%
%     new_struct = struct_strrep(old_struct, old_string, new_string)

    MAX_STRUCT_DEPTH = 200;
    new_struct = old_struct;

    if(~exist('depth', 'var') || isempty(depth))
        depth = 1;
    elseif(depth > MAX_STRUCT_DEPTH)
        error('Exceeding max structure depth of %d. If this is legitimate, please edit this function.\n', depth);
    end

    fnames = fieldnames(new_struct);
    for i = 1:length(fnames)
        % Below code does not work as Matlab R2007B has some massive bug in the class() function.
        % It works when executed manualy, but not as part of a func... wtf?
        %         field_class = class(new_struct.(fnames{i}));
        %         switch(field_class)
        %             case 'struct'
        %                 new_struct.(fnames{i}) = struct_strrep(new_struct.(fnames{i}), old_string, new_string);
        %             case 'cell'
        %                 if(iscellstr(new_struct.(fnames{i})))
        %                     new_struct.(fnames{i}) = strrep(new_struct.(fnames{i}), old_string, new_string);
        %                 end
        %             case 'char'
        %                 new_struct.(fnames{i}) = char(strrep(cellstr(new_struct.(fnames{i})), old_string, new_string));
        %         end

        if(isstruct(new_struct.(fnames{i})))
            for j=1:length(new_struct.(fnames{i}))
                new_struct.(fnames{i})(j) = struct_strrep(new_struct.(fnames{i})(j), old_string, new_string, depth + 1);
            end
        elseif(iscellstr(new_struct.(fnames{i})))
            new_struct.(fnames{i}) = strrep(new_struct.(fnames{i}), old_string, new_string);
        elseif(ischar(new_struct.(fnames{i})))
            new_struct.(fnames{i}) = char(strrep(cellstr(new_struct.(fnames{i})), old_string, new_string));
        elseif(iscell(new_struct.(fnames{i})))
            for j=1:length(new_struct.(fnames{i}))
                new_struct.(fnames{i}){j} = struct_strrep(new_struct.(fnames{i}){j}, old_string, new_string, depth + 1);
            end
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function result_payload = my_finish_mean(payload, ECG_header)

result_payload.mean = payload.the_sum ./ payload.the_size;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function points=Cart2Ellipse(cartPoints,algorithm,a,f)
%%CART2ELLIPSE Convert Cartesian coordinates to ellipsoidal (latitude,
%              longitude, and altitude) coordinates.
%
%INPUTS: cartPoints A matrix of the points in ECEF Cartesian coordinates
%                   that are to be transformed into ellipsoidal
%                   coordinates. Each column of cartPoints is of the
%                   format [x;y;z].
%         algorithm This specified the algorithm to use for the conversion.
%                   Note that none work at the origin. Possible values are:
%                   0 (The default if this parameter is omitted or an empty
%                     matrix is passed and f<0.01) Use the algorithm of
%                     Olson in [1].
%                   1 Use the Algorithm of Sofair in [2], which is a
%                     modification of [3]. This will not work close to the
%                     center of the Earth.
%                   2 (The default if this parameter is omitted or an empty
%                     matrix is passed and f>=0.01) Use the algorithm of
%                     Fukushima in [4]. This should work close to the
%                     center of the Earth.
%                 a The semi-major axis of the reference ellipsoid. If this
%                   argument is omitted, the value in
%                   Constants.WGS84SemiMajorAxis is used.
%                 f The flattening factor of the reference ellipsoid. If
%                   this argument is omitted, the value in
%                   Constants.WGS84Flattening is used.
%
%OUTPUTS: points A matrix of the converted points. Each column of the
%                matrix has the format [latitude;longitude;altitude], with
%                latitude and longitude given in radians.
%
%The algorithm of Olson in [1] appears to be the most precise non-iterative
%method available for targets far above the Earth. The method of Sofair in
%[2] and [3] is also a non-iterative algorithm, but tends to have
%significantly worse accuracy for such targets. Fukushima's algorithm in
%[4] is iterative and typically converges in six or fewer iterations. It is
%set to run for a maximum number of 500 iterations. More than 6 iterations
%can be necessary when a large flattening is used and the point in question
%is near the center of the Earth. Its accuracy appears to be marginally
%better than [1] for things not near the surface of the Earth/ reference
%ellipsoid, but it is slower.
%
%REFERENCES:
%[1] D. K. Olson, "Converting Earth-centered, Earth-fixed coordinates to
%    geodetic coordinates," IEEE Transactions on Aerospace and Electronic
%    Systems, vol. 32, no. 1, pp. 473-476, Jan. 1996.
%[2] I. Sofair "Improved method for calculating exact geodetic latitude and
%    altitude revisited," Journal of Guidance, Control, and Dynamics, vol.
%    23, no. 2, p. 369, Mar. 2000.
%[3] I. Sofair, "Improved method for calculating exact geodetic latitude
%    and altitude," Journal of Guidance, Control, and Dynamics, vol. 20,
%    no. 4, pp. 824-826, Jul.-Aug. 1997.
%[4] Fukushima, T., "Transformation from Cartesian to geodetic coordinates
%    accelerated by Halley's method", Journal of Geodesy, vol. 79, no. 12,
%    pp. 689-693, Mar. 2006.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(f))
    f=Constants.WGS84Flattening;
end

if(nargin<3||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

if(nargin<2||isempty(algorithm))
    if(f<0.01)
        algorithm=0;
    else
        algorithm=2;
    end
end

switch(algorithm)
    case 0%Olson's algorithm
        [lambda,phi,h]=OlsonAlg(cartPoints,a,f);

        if(any(imag(lambda)~=0)||any(imag(phi)~=0)||any(imag(phi)~=0))
            error('The point given is too close to the center of the Earth for the algorithm of Olson.')
        end

    case 1%Sofair's algorithm
        [lambda,phi,h]=SofairAlg(cartPoints,a,f);
    case 2%Fukushima's algorithm
        [lambda,phi,h]=FukishimaAlg(cartPoints,a,f);
    otherwise
        error('Unknown algorithm specified.')
end

points=[phi;lambda;h];

end

function [lambda,phi,h]=SofairAlg(cartPoints,a,f)
%%SOFAIRALG This implements the algorithm of [1], which is a modified
%           version of the algorithm of [2]. Both techniques will fail if
%           the point in question is too close to the origin (deep within
%           the Earth). If the algorithm fails, then this function wil have
%           an error.
%
%REFERENCES:
%[1] I. Sofair "Improved method for calculating exact geodetic latitude and
%    altitude revisited," Journal of Guidance, Control, and Dynamics, vol.
%    23, no. 2, p. 369, Mar. 2000.
%[2] I. Sofair, "Improved method for calculating exact geodetic latitude
%    and altitude," Journal of Guidance, Control, and Dynamics, vol. 20,
%    no. 4, pp. 824-826, Jul.-Aug. 1997.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numPoints=size(cartPoints,2);

%The semi-minor axis of the reference ellipsoid.
b=a*(1-f);

%The square of the first numerical eccentricity.
e2=2*f-f^2;

%The square of the second numerical eccentricity.
eps2=a^2/b^2-1;

%Allocate space for the results.
phi=zeros(1,numPoints);
lambda=zeros(1,numPoints);
h=zeros(1,numPoints);
for curPoint=1:numPoints
    %Extract the coordinates
    x0=cartPoints(1,curPoint);
    y0=cartPoints(2,curPoint);
    z0=cartPoints(3,curPoint);

    r0=sqrt(x0.^2+y0.^2);
    p=abs(z0)/eps2;
    s=r0.^2/(e2*eps2);
    q=p.^2-b.^2+s;

    lambda(curPoint)=atan2(y0,x0);

    if(q<0)
        error('The point given is too close to the center of the Earth for the algorithm of Sofair.')
    end

    u=p./sqrt(q);
    v=b^2*u.^2./q;
    P=27*v.*s./q;
    Q=(sqrt(P+1)+sqrt(P)).^(2/3);
    t=(1+Q+1./Q)/6;
    %The max command prevents finite precision problems due to
    %subtraction within the square root.
    c=max(0,u.^2-1+2*t);
    c=sqrt(c);
    w=(c-u)/2;

    %The z coordinate of the closest point projected on the ellipsoid.
    %The max command deals with precision problems when the argument
    %is nearly zero. The problems arise due to the subtraction within
    %the square root.
    z=max(0,sqrt(t.^2+v)-u.*w-t/2-1/4);
    z=sign(z0).*sqrt(q).*(w+sqrt(z));
    Ne=a*sqrt(1+eps2.*z.^2/b^2);

    %The min and max terms deals with finite precision problems.
    val=min(z*(eps2+1)./Ne,1);
    val=max(val,-1);
    phi(curPoint)=asin(val);
    h(curPoint)=r0.*cos(phi(curPoint))+z0.*sin(phi(curPoint))-a^2./Ne;
end

end

function [lambda,phi,h]=FukishimaAlg(cartPoints,a,f)
%%FUKUSHIMAALG This function implements the algorithm of [1] with minor
%              modifications.
%
%If one lets the algorithm run for an arbitrary number of iterations, there
%will generally be underflows, since the ratio of S and C matter, but both
%terms can drift by a constant factor during the iterations. Thus, after
%each iterative step, the values are normalized so that C=1 (it takes one
%division). Convergence of Fukushima's method is assumed to occur after
%six iterations.
%
%REFERENCES:
%[1] Fukushima, T., "Transformation from Cartesian to geodetic coordinates
%    accelerated by Halley's method", J.Geodesy (2006) 79: 689-693.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numPoints=size(cartPoints,2);

%The semi-minor axis of the reference ellipsoid.
b=a*(1-f);

%The square of the first numerical eccentricity.
e2=2*f-f^2;

ec=sqrt(1-e2);

%Allocate space for the results.
phi=zeros(1,numPoints);
lambda=zeros(1,numPoints);
h=zeros(1,numPoints);
for curPoint=1:numPoints
    %Extract the coordinates
    x0=cartPoints(1,curPoint);
    y0=cartPoints(2,curPoint);
    z0=cartPoints(3,curPoint);

    lambda(curPoint)=atan2(y0,x0);

    p=sqrt(x0^2+y0^2);
    P=p/a;
    Z=(ec/a)*abs(z0);

    S=Z;
    C=ec*P;

    %Loop until convergence. Normally, only 6 iterations or fewer is
    %required. When some points are near the surface of the Earth and f is
    %close to 1, the required number of iterations can be higher.
    for curIter=1:500
        A=sqrt(S^2+C^2);
        B=1.5*e2*S*C^2*((P*S-Z*C)*A-e2*S*C);
        F=P*A^3-e2*C^3;
        D=Z*A^3+e2*S^3;

        SNew=D*F-B*S;
        CNew=F^2-B*C;

        SOld=S;
        COld=C;

        SNew=SNew/CNew;
        if(~isfinite(SNew))
            S=SNew;
            break;
        else
            S=SNew;
            C=1;
        end

        if(S==SOld&&C==COld)
           break;
        end
    end
    Cc=ec*C;

    %If the point is along the z-axis, then SNew and CNew will
    %both be zero, leading to a non-finite result.
    if(~isfinite(S))
        phi(curPoint)=sign(z0)*(pi/2);
        h(curPoint)=abs(z0)-b;
    else
        phi(curPoint)=sign(z0)*atan(S/Cc);
        h(curPoint)=(p*Cc+abs(z0)*S-b*A)/sqrt(Cc^2+S^2);
    end
end

end

function [lambda,phi,h]=OlsonAlg(cartPoints,a,f)
%%OLSONALG This function implements the algorithm of [1], removing a test
%          for a radius being too small as the function will still work at
%          smaller radii.
%
%REFERENCES:
%[1] D. K. Olson, "Converting Earth-centered, Earth-fixed coordinates to
%   geodetic coordinates," IEEE Transactions on Aerospace and Electronic
%   Systems, vol. 32, no. 1, pp. 473-476, Jan. 1996.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numPoints=size(cartPoints,2);

%The square of the eccentricity.
e2=2*f-f^2;

a1=a*e2;
a2=a1^2;
a3=a1*e2/2;
a4=(5/2)*a2;
a5=a1+a3;
a6=1-e2;

%Allocate space for the results.
phi=zeros(1,numPoints);
lambda=zeros(1,numPoints);
h=zeros(1,numPoints);
for curPoint=1:numPoints
    %Extract the coordinates
    x=cartPoints(1,curPoint);
    y=cartPoints(2,curPoint);
    z=cartPoints(3,curPoint);

    zp=abs(z);
    w2=x^2+y^2;
    w=sqrt(w2);
    z2=z^2;
    r2=w2+z2;
    %The algorithm will work with points close to the origin. Thus, there
    %is no need to have a test for r being too small as is the case in [1].
    r=sqrt(r2);

    lambda(curPoint)=atan2(y,x);
    s2=z2/r2;
    c2=w2/r2;
    u=a2/r;
    v=a3-a4/r;
    if(c2>0.3)
        s=(zp/r)*(1+c2*(a1+u+s2*v)/r);
        phi(curPoint)=asin(s);
        ss=s^2;
        c=sqrt(1-ss);
    else
        c=(w/r)*(1-s2*(a5-u-c2*v)/r);
        phi(curPoint)=acos(c);
        ss=1-c^2;
        s=sqrt(ss);
    end

    g=1-e2*ss;
    rg=a/sqrt(g);
    rf=a6*rg;
    u=w-rg*c;
    v=zp-rf*s;
    f=c*u+s*v;
    m=c*v-s*u;
    p=m/(rf/g+f);
    phi(curPoint)=phi(curPoint)+p;
    h(curPoint)=f+m*p/2;
    if(z<0)
        phi(curPoint)=-phi(curPoint);
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function cadj=ref_spreadadj_5(coef)
%REF_SPREADADJ_5  Symbol of adjoint spreading function.
%   Usage: cadj=ref_spreadadj_5(c,number);
%
%   Development version by FJ for comparison of different implementations
%   cadj=SPREADADJ(c) will compute the symbol cadj of the spreading
%   operator that is the adjoint of the spreading operator with symbol c.
%
%   Implementation for sparse matrix using loop

L=size(coef,1);

[row,col,val]=find(coef);

% Optimization note : As said in note of case 3, we are computing
% the Lth root of unity which have special properties and symetries
% taht could be exploited to highly reduce this computation.
% Furthermore here we precompute every possible exponential
% term even if some are unneeded. But there is no simple way to
% know which one should be computed and the computation needed to
% know it could be worst than this computation. Nevertheless, it
% could be used in a different implementation
temp=exp((-i*2*pi/L)*(0:L-1));

cadj=spalloc(L,L,length(val));
for k=1:length(val)
  ii=mod(L-row(k)+1, L);
  jj=mod(L-col(k)+1, L);
  cadj(ii+1,jj+1)=conj(val(k))*temp(mod(ii*jj,L)+1);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [c, ceq] = ConstraintFCN_models(u,uold,x,N,LBo,UBo,LBdu,UBdu,p,select_model)
%% Constraint function of nonlinear MPC for F8 system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon length
%   uold:   latest applied control input
%   LBo:    Lower bound of output x
%   UBo:    Upper bound of output x
%   LBdu:   Lower bound for input difference uk - uk-1
%   UBdu:   Upper bound for input difference uk - uk-1
%   p:      Parameters for model
%   select_model: Selects model future-state prediction
%
% Output:
%   c:      inequality constraints applied across prediction horizon
%   ceq:    equality constraints  (empty)
%

Nvar = length(x);
%% Nonlinear MPC design parameters
% Ensure that all cell populations are positive
zMin = LBo;

%% Integrate system
if strcmp(select_model,'DelayDMDc')
    Uinput = getHankelMatrix_MV([p.udelay(2:end)';u],p.Ndelay)';
    [xk,~] = lsim(p.sys,Uinput,[0:N-1].*p.dt,p.xdelay);
    xk = xk(:,end-Nvar+1:end);
    xk = xk + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'DMDc')
    [xk,~] = lsim(p.sys,[u',0],[0:N].*p.dt,x-p.xmean);
    xk = xk(2:end,:) + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'eDMDc')
    Y0 = poolData((x-p.xmean)',Nvar,p.polyorder,p.usesine)';
    [xk,~] = lsim(p.sys,[u' 0],[0:N].*p.dt,Y0(2:end));
    xk = xk(2:end,1:Nvar) + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'SINDYc')
    Ns = size(x,1);
    xk = zeros(Ns,N+1); xk(:,1) = x;
    for ct=1:N
        % Obtain plant state at next prediction step.
        xk(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xk(:,ct),u(ct),p.dt,1,[],p);
    end
    xk = xk(:,2:N+1);
elseif strcmp(select_model,'partialSINDYc')
    Ns = size(x,1);
    xk = zeros(Ns,N+1); xk(:,1) = x;
    for ct=1:N
        % Obtain plant state at next prediction step.
        xk_tmp = rk4u(@sparseGalerkinControl_Discrete,xk(p.SelectVars,ct),u(ct),p.dt,1,[],p);
        xk(p.SelectVars,ct+1) = xk_tmp;
    end
    xk = xk(:,2:N+1);
elseif strcmp(select_model,'NARX')
    Hu = [u',0];
    Hx = zeros(Nvar,length(Hu)); Hx(:,1) = x;
    if p.TRANSFORM_LOG == 1
        Hx = log(Hx);
    end
    [Us,Ui,Si] = preparets(p.net,con2seq(Hu),{},con2seq(Hx));
    xk = p.net(Us,Ui,Si);
    xk = cell2mat(xk);
    if p.TRANSFORM_LOG == 1
        xk = exp(xk);
    end
end

%% Inequality constraints calculation
c = zeros(N,1);
% Apply N population size constraints across prediction horizon, from time
% k+1 to k+N

for ct=1:N
    % -z + zMin < 0 % lower bound
    c(ct) = -xk(1,ct)+zMin;

end

%% No equality constraints
ceq = [];

author: eurika-kaiser

repo: SINDY-MPC


text:
classdef plotlyfig < handle

    %----CLASS PROPERTIES----%
    properties
        data; % data of the plot
        layout;  % layout of the plot
        url; % url response of making post request
        error; % error response of making post request
        warning; % warning response of making post request
        message; % message response of making post request

    end

    properties (SetObservable)
        UserData;% credentials/configuration/verbose
        PlotOptions; % filename,fileopt,world_readable
    end

    properties (Hidden = true)
        PlotlyDefaults; % plotly specific conversion defualts
        State; % state of plot (FIGURE/AXIS/PLOTS)
    end

    properties (Access = private)
        PlotlyReference; % load the plotly reference
        InitialState; % inital userdata
    end

    %----CLASS METHODS----%
    methods

        %----CONSTRUCTOR---%
        function obj = plotlyfig(varargin)

            %-Core-%
            obj.data = {};
            obj.layout = struct();
            obj.url = '';

            %-UserData-%
            try
                [obj.UserData.Username,...
                    obj.UserData.ApiKey,...
                    obj.UserData.PlotlyDomain] = signin;
            catch
                errkey = 'plotlyfigConstructor:notSignedIn';
                error(errkey, plotlymsg(errkey));
            end

            obj.UserData.Verbose = true;

            %-PlotOptions-%
            obj.PlotOptions.CleanFeedTitle = true;
            obj.PlotOptions.FileName = '';
            obj.PlotOptions.FileOpt = 'new';
            obj.PlotOptions.WorldReadable = true;
            obj.PlotOptions.ShowURL = true;
            obj.PlotOptions.OpenURL = true;
            obj.PlotOptions.Strip = true;
            obj.PlotOptions.Visible = 'on';
            obj.PlotOptions.TriangulatePatch = false;

            % offline options
            obj.PlotOptions.Offline = false;
            obj.PlotOptions.ShowLinkText = true;
            obj.PlotOptions.LinkText = obj.get_link_text;
            obj.PlotOptions.IncludePlotlyjs = true;

            %-PlotlyDefaults-%
            obj.PlotlyDefaults.MinTitleMargin = 80;
            obj.PlotlyDefaults.TitleHeight = 0.01;
            obj.PlotlyDefaults.TitleFontSizeIncrease = 40;
            obj.PlotlyDefaults.FigureIncreaseFactor = 1.5;
            obj.PlotlyDefaults.AxisLineIncreaseFactor = 1.5;
            obj.PlotlyDefaults.MarginPad = 0;
            obj.PlotlyDefaults.MaxTickLength = 20;
            obj.PlotlyDefaults.ExponentFormat = 'none';
            obj.PlotlyDefaults.ErrorbarWidth = 6;
            obj.PlotlyDefaults.ShowBaselineLegend = false;
            obj.PlotlyDefaults.Bargap = 0;
            obj.PlotlyDefaults.CaptionMarginIncreaseFactor = 1.2;
            obj.PlotlyDefaults.MinCaptionMargin = 80;

            %-State-%
            obj.State.Figure = [];
            obj.State.Axis = [];
            obj.State.Plot = [];
            obj.State.Text = [];
            obj.State.Legend = [];
            obj.State.Colorbar = [];

            % figure object management
            obj.State.Figure.NumAxes = 0;
            obj.State.Figure.NumPlots = 0;
            obj.State.Figure.NumLegends = 0;
            obj.State.Figure.NumColorbars = 0;
            obj.State.Figure.NumTexts = 0;

            %-PlotlyReference-%
            obj.PlotlyReference = [];

            %-InitialState-%
            obj.InitialState.Username = obj.UserData.Username;
            obj.InitialState.ApiKey = obj.UserData.ApiKey;
            obj.InitialState.PlotlyDomain = obj.UserData.PlotlyDomain;

            % initialize figure handle
            fig_han = [];

            % initialize autoupdate key
            updatekey = false;

            % parse inputs
            switch nargin

                case 0
                case 1
                    % check for figure handle
                    if ishandle(varargin{1})
                        if strcmp(get(varargin{1},'type'),'figure')
                            fig_han = varargin{1};
                            updatekey = true;
                        end
                    else
                        errkey = 'plotlyfigConstructor:invalidInputs';
                        error(errkey , plotlymsg(errkey));
                    end

                otherwise

                    % check for figure handle
                    if ishandle(varargin{1})
                        if strcmp(get(varargin{1},'type'),'figure')
                            fig_han = varargin{1};
                            updatekey = true;
                            parseinit = 2;
                        end
                    else
                        parseinit = 1;
                    end

                    % check for proper property/value structure
                    if mod(length(parseinit:nargin),2) ~= 0
                        errkey = 'plotlyfigConstructor:invalidInputs';
                        error(errkey , plotlymsg(errkey));
                    end

                    % parse property/values
                    for a = parseinit:2:nargin
                        if(strcmpi(varargin{a},'filename'))
                            obj.PlotOptions.FileName = varargin{a+1};
                            % overwrite if filename provided
                            obj.PlotOptions.FileOpt = 'overwrite';
                        end
                        if(strcmpi(varargin{a},'fileopt'))
                            obj.PlotOptions.FileOpt = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'world_readable'))
                            obj.PlotOptions.WorldReadable = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'link'))
                            obj.PlotOptions.ShowURL = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'open'))
                            obj.PlotOptions.OpenURL = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'strip'))
                            obj.PlotOptions.Strip = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'visible'))
                            obj.PlotOptions.Visible = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'offline'))
                            obj.PlotOptions.Offline = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'showlink'))
                            obj.PlotOptions.ShowLinkText = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'linktext'))
                            obj.PlotOptions.LinkText = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'include_plotlyjs'))
                            obj.PlotOptions.IncludePlotlyjs = varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'layout'))
                            obj.layout= varargin{a+1};
                        end
                        if(strcmpi(varargin{a},'data'))
                            obj.data = varargin{a+1};
                        end
                    end
            end

            % create figure/axes if empty
            if isempty(fig_han)
                fig_han = figure;
                axes;
            end

            % plotly figure default style
            set(fig_han,'Name',obj.PlotOptions.FileName,'Color',[1 1 1],'NumberTitle','off', 'Visible', obj.PlotOptions.Visible);

            % figure state
            obj.State.Figure.Handle = fig_han;

            % update
            if updatekey
                obj.update;
            end

            % add figure listeners
            addlistener(obj.State.Figure.Handle,'Visible','PostSet',@(src,event)updateFigureVisible(obj,src,event));
            addlistener(obj.State.Figure.Handle,'Name','PostSet',@(src,event)updateFigureName(obj,src,event));

            % add plot options listeners
            addlistener(obj,'PlotOptions','PostSet',@(src,event)updatePlotOptions(obj,src,event));

            % add user data listeners
            addlistener(obj,'UserData','PostSet',@(src,event)updateUserData(obj,src,event));
        end

        %-------------------------USER METHODS----------------------------%

        %----GET OBJ.STATE.FIGURE.HANDLE ----%
        function plotlyFigureHandle = gpf(obj)
            plotlyFigureHandle = obj.State.Figure.Handle;
            set(0,'CurrentFigure', plotlyFigureHandle);
        end

        %----LOAD PLOTLY REFERENCE-----%
        function obj = loadplotlyref(obj)
            if isempty(obj.PlotlyReference)

                % plotly reference
                plotlyref = load('plotly_reference.mat');

                % update the PlotlyRef property
                obj.PlotlyReference = plotlyref.pr;
            end
        end

        %----KEEP THE MATLAB STYLE DEFAULTS----%
        function obj = revert(obj)

            % set the PlotOptions.Strip property
            obj.PlotOptions.Strip = false;

            % update the object
            obj.update;

        end

        %----STRIP THE STYLE DEFAULTS----%
        function obj = strip(obj)

            % set the PlotOptions.Strip property
            obj.PlotOptions.Strip = true;

            % load plotly reference
            obj.loadplotlyref;

            % strip the style keys from data
            for d = 1:length(obj.data)
                obj.data{d} = obj.stripkeys(obj.data{d}, obj.data{d}.type, 'style');
            end

            % strip the style keys from layout
            obj.layout = obj.stripkeys(obj.layout, 'layout', 'style');

        end

        %----GET THE FIELDS OF TYPE DATA----%
        function data = getdata(obj)

            % load plotly reference
            obj.loadplotlyref;

            % initialize output
            data = cell(1,length(obj.data));

            % remove style / plot_info types in data
            for d = 1:length(obj.data)
                data{d} = obj.stripkeys(obj.data{d}, obj.data{d}.type, {'style','plot_info'});
            end

        end

        %----VALIDATE THE DATA/LAYOUT KEYS----%
        function validate(obj)

            % load plotly reference
            obj.loadplotlyref;

            % validate data fields
            for d = 1:length(obj.data)
                obj.stripkeys(obj.data{d}, obj.data{d}.type, {'style','plot_info'});
            end

            % validate layout fields
            obj.stripkeys(obj.layout, 'layout', 'style');

        end

        %----GET PLOTLY FIGURE-----%
        function obj = download(obj, file_owner, file_id)
            plotlyfig = plotlygetfile(file_owner, file_id);
            obj.data = plotlyfig.data;
            obj.layout = plotlyfig.layout;
        end

        %---------------------STATIC IMAGE METHODS------------------------%

        %----SAVE STATIC PLOTLY IMAGE-----%
        function obj = saveas(obj, filename, varargin)

            % strip keys
            if obj.PlotOptions.Strip
                obj.strip;
            end

            % handle title
            if isempty(filename)
                handleFileName(obj);
                filename = obj.PlotOptions.FileName;
            end

            % create image figure
            imgfig.data = obj.data;
            imgfig.layout = obj.layout;

            % save image
            plotlygenimage(imgfig, filename, varargin{:});
        end

        %----SAVE STATIC JPEG IMAGE-----%
        function obj = jpeg(obj, filename)
            if nargin > 1
                obj.saveas(filename,'jpeg');
            else
                obj.saveas(obj.PlotOptions.FileName,'jpeg');
            end
        end

        %----SAVE STATIC PDF IMAGE-----%
        function obj = pdf(obj, filename)
            if nargin > 1
                obj.saveas(filename,'pdf');
            else
                obj.saveas(obj.PlotOptions.FileName,'pdf');
            end
        end

        %----SAVE STATIC PNG IMAGE-----%
        function obj = png(obj, filename)
            if nargin > 1
                obj.saveas(filename,'png');
            else
                obj.saveas(obj.PlotOptions.FileName,'png');
            end
        end

        %----SAVE STATIC SVG IMAGE-----%
        function obj = svg(obj, filename)
            if nargin > 1
                obj.saveas(filename,'svg');
            else
                obj.saveas(obj.PlotOptions.FileName,'svg');
            end
        end

        %----ADD A CUSTOM CAPTION-----%
        function obj = add_caption(obj, caption_string, varargin)

            caption.text = caption_string;

            % defaults
            caption.xref = 'paper';
            caption.yref = 'paper';
            caption.xanchor = 'left';
            caption.yanchor = 'top';
            caption.x = 0.1;
            caption.y = -0.05;
            caption.showarrow = false;

            % inject any custom annotation specs
            for n = 1:2:length(varargin)
                caption = setfield(caption, varargin{n}, varargin{n+1});
            end

            % adjust the bottom margin
            obj.layout.margin.b = max(obj.layout.margin.b, ...
                    obj.PlotlyDefaults.MinCaptionMargin);

            % add the new caption to the figure
            obj.State.Figure.NumTexts = obj.State.Figure.NumTexts + 1;
            obj.layout.annotations{obj.State.Figure.NumTexts} = caption;

            % update the figure state
            obj.State.Text(obj.State.Figure.NumTexts).Handle = NaN;
            obj.State.Text(obj.State.Figure.NumTexts).AssociatedAxis = gca;
            obj.State.Text(obj.State.Figure.NumTexts).Title = false;

        end

        %------------------------REST API CALL----------------------------%

        %----SEND PLOT REQUEST (NO UPDATE)----%
        function obj = plotly(obj)

            % strip keys
            if obj.PlotOptions.Strip
                obj.strip;
            end

            % validate keys
            validate(obj);

            % handle filename
            handleFileName(obj);

            % handle title (for feed)
            if obj.PlotOptions.CleanFeedTitle
                cleanFeedTitle(obj);
            end

            %args
            args.filename = obj.PlotOptions.FileName;
            args.fileopt = obj.PlotOptions.FileOpt;
            args.world_readable = obj.PlotOptions.WorldReadable;

            %layout
            args.layout = obj.layout;

            %send to plotly
            if ~obj.PlotOptions.Offline
                response = plotly(obj.data, args);

                %update response
                obj.url = response.url;
                obj.error = response.error;
                obj.warning = response.warning;
                obj.message = response.message;

                %open url in browser
                if obj.PlotOptions.OpenURL
                    web(response.url, '-browser');
                end
            else
                obj.url = plotlyoffline(obj);
                if obj.PlotOptions.OpenURL
                    web(obj.url, '-browser');
                end
            end

        end

        %-----------------------FIGURE CONVERSION-------------------------%

        %automatic figure conversion
        function obj = update(obj)

            % reset figure object count
            obj.State.Figure.NumAxes = 0;
            obj.State.Figure.NumPlots = 0;
            obj.State.Figure.NumLegends = 0;
            obj.State.Figure.NumColorbars = 0;
            obj.State.Figure.NumTexts = 0;

            % find axes of figure
            ax = findobj(obj.State.Figure.Handle,'Type','axes','-and',{'Tag','','-or','Tag','PlotMatrixBigAx','-or','Tag','PlotMatrixScatterAx'});
            obj.State.Figure.NumAxes = length(ax);

            % update number of annotations (one title per axis)
            obj.State.Figure.NumTexts = length(ax);

            % find children of figure axes
            for a = 1:length(ax)

                % reverse axes
                axrev = length(ax) - a + 1;

                % set axis handle field
                obj.State.Axis(a).Handle = ax(axrev);

                % add title
                obj.State.Text(a).Handle = get(ax(axrev),'Title');
                obj.State.Text(a).AssociatedAxis = handle(ax(axrev));
                obj.State.Text(a).Title = true;

                % find plots of figure
                plots = findobj(ax(axrev),'-not','Type','Text','-not','Type','axes','-depth',1);

                % add baseline objects
                baselines = findobj(ax(axrev),'-property','BaseLine');

                for np = 1:length(plots)

                    % reverse plots
                    nprev = length(plots) - np + 1;

                    % update the plot fields
                    obj.State.Figure.NumPlots = obj.State.Figure.NumPlots + 1;
                    obj.State.Plot(obj.State.Figure.NumPlots).Handle = handle(plots(nprev));
                    obj.State.Plot(obj.State.Figure.NumPlots).AssociatedAxis = handle(ax(axrev));
                    obj.State.Plot(obj.State.Figure.NumPlots).Class = getGraphClass(plots(nprev));
                end

                % find text of figure
                texts = findobj(ax(axrev),'Type','text','-depth',1);

                for t = 1:length(texts)
                    obj.State.Text(obj.State.Figure.NumTexts + t).Handle = handle(texts(t));
                    obj.State.Text(obj.State.Figure.NumTexts + t).Title = false;
                    obj.State.Text(obj.State.Figure.NumTexts + t).AssociatedAxis = handle(ax(axrev));
                end

                % update number of annotations
                obj.State.Figure.NumTexts = obj.State.Figure.NumTexts + length(texts);

            end

            % find legends of figure
            if isHG2
                legs = findobj(obj.State.Figure.Handle,'Type','Legend');
            else
                legs = findobj(obj.State.Figure.Handle,'Type','axes','-and','Tag','legend');
            end

            obj.State.Figure.NumLegends = length(legs);

            for g = 1:length(legs)
                obj.State.Legend(g).Handle = handle(legs(g));

                % find associated axis
                legendAxis = findLegendAxis(obj, handle(legs(g)));

                % update colorbar associated axis
                obj.State.Legend(g).AssociatedAxis = legendAxis;
            end

            % find colorbar of figure
            if isHG2
                cols = findobj(obj.State.Figure.Handle,'Type','Colorbar');
            else
                cols = findobj(obj.State.Figure.Handle,'Type','axes','-and','Tag','Colorbar');
            end

            obj.State.Figure.NumColorbars = length(cols);

            for c = 1:length(cols)

                % update colorbar handle
                obj.State.Colorbar(c).Handle = handle(cols(c));

                % find associated axis
                colorbarAxis = findColorbarAxis(obj, handle(cols(c)));

                % update colorbar associated axis
                obj.State.Colorbar(c).AssociatedAxis = colorbarAxis;

            end

            %--------------------UPDATE PLOTLY FIGURE---------------------%

            % reset dataget(obj.State.Figure.Handle,'Children')
            obj.data = {};

            % reset layout
            obj.layout = struct();

            % update figure
            updateFigure(obj);

            % update axes
            for n = 1:obj.State.Figure.NumAxes
                updateAxis(obj,n);
            end

            % update plots
            for n = 1:obj.State.Figure.NumPlots
                updateData(obj,n);
            end

            % update annotations
            for n = 1:obj.State.Figure.NumTexts
                updateAnnotation(obj,n);
            end

            % update legends
            for n = 1:obj.State.Figure.NumLegends
                updateLegend(obj,n);
            end

            % update colorbars
            for n = 1:obj.State.Figure.NumColorbars
                updateColorbar(obj,n);
            end

        end

        %----------------------EXTRACT PLOTLY INDICES---------------------%

        %----GET CURRENT AXIS INDEX ----%
        function currentAxisIndex = getAxisIndex(obj,axishan)
            currentAxisIndex = find(arrayfun(@(x)(eq(x.Handle,axishan)),obj.State.Axis));
        end

        %----GET CURRENT DATA INDEX ----%
        function currentDataIndex = getDataIndex(obj,plothan)
            currentDataIndex = find(arrayfun(@(x)(eq(x.Handle,plothan)),obj.State.Plot));
        end

        %----GET CURRENT ANNOTATION INDEX ----%
        function currentAnnotationIndex = getAnnotationIndex(obj,annothan)
            currentAnnotationIndex = find(arrayfun(@(x)(eq(x.Handle,annothan)),obj.State.Text));
        end

        %-------------------CALLBACK FUNCTIONS--------------------------%

        %----UPDATE FIGURE OPTIONS----%
        function obj = updateFigureVisible(obj,src,event)
            % update PlotOptions.Visible
            obj.PlotOptions.Visible = get(obj.State.Figure.Handle,'Visible');
        end

        function obj = updateFigureName(obj,src,event)
            % update PlotOptions.Name
            obj.PlotOptions.FileName = get(obj.State.Figure.Handle,'Name');
        end

        %----UPDATE PLOT OPTIONS----%
        function obj = updatePlotOptions(obj,src,event)
            if iscell(obj.PlotOptions.FileName)
                for i = 1 : numel(obj.PlotOptions.FileName)
                    obj.PlotOptions.FileName{i} = [obj.PlotOptions.FileName{i} ' - '];
                end
                obj.PlotOptions.FileName = cell2mat(obj.PlotOptions.FileName');
            end
            set(obj.State.Figure.Handle, 'Name', obj.PlotOptions.FileName, 'Visible', obj.PlotOptions.Visible);
        end

        %----UPDATE USER DATA----%
        function obj = updateUserData(obj,src,event)
            signin(obj.UserData.Username,...
                obj.UserData.ApiKey,...
                obj.UserData.PlotlyDomain);
        end

        %-------------------OVERLOADED FUNCTIONS--------------------------%

        %----PLOT FUNCTIONS----%
        function han = image(obj, varargin)
            %call image
            han = image(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = imagesc(obj, varargin)
            %call imagesc
            han = imagesc(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = line(obj, varargin)
            %call line
            han = line(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = patch(obj, varargin)
            %call patch
            han = patch(varargin{:});
            %update object
            obj.update;
        end

        function han = rectangle(obj, varargin)
            %call rectangle
            han = rectangle(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = area(obj,varargin)
            %call area
            han = area(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = bar(obj,varargin)
            %call bar
            han = bar(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = contour(obj,varargin)
            %call contour
            han = contour(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = plot(obj,varargin)
            %call plot
            han = plot(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = errorbar(obj,varargin)
            %call errorbar
            han = errorbar(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = quiver(obj,varargin)
            %call quiver
            han = quiver(varargin{:});
            %update object
            obj.update;
        end

        function han = scatter(obj, varargin)
            %call bar
            han = scatter(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = stairs(obj,varargin)
            %call stairs
            han = stairs(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = stem(obj,varargin)
            %call stem
            han = stem(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = boxplot(obj,varargin)
            %call boxplot
            han = boxplot(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function han = mesh(obj,varargin)
            %call mesh
            han = mesh(varargin{:});
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        function [y,t,x] = initial(obj,varargin)
            % plot initial
            initial(varargin{:});
            % call initial
            [y,t,x] = initial(varargin{:});
            % fake output by calling plot
            plot(t,y);
            %update object
            obj.update;
            %send to plotly
            obj.plotly;
        end

        %----OTHER----%
        function delete(obj)
            % reset persistent USERNAME, KEY, and DOMAIN
            % of signin to original state
            if isfield(obj.InitialState,'Username')
                signin(obj.InitialState.Username, ...
                    obj.InitialState.ApiKey,...
                    obj.InitialState.PlotlyDomain);
            end
        end
    end

    methods (Access=private)
        %----STRIP THE FIELDS OF A SPECIFIED KEY-----%
        function stripped = stripkeys(obj, fields, fieldname, key)

            %plorlt reference
            pr = obj.PlotlyReference;

            % initialize output
            stripped = fields;

            % get fieldnames
            fn = fieldnames(stripped);
            fnmod = fn;

            try
                for d = 1:length(fn);

                    % clean up axis keys
                    if any(strfind(fn{d},'xaxis')) || any(strfind(fn{d},'yaxis'))
                        fnmod{d} = fn{d}(1:length('_axis'));
                    end

                    % keys:(object, style, plot_info, data)
                    keytype = getfield(pr,fieldname,fnmod{d},'key_type');

                    % check for objects
                    if strcmp(keytype,'object')

                        % clean up font keys
                        if any(strfind(fn{d},'font'))
                            fnmod{d} = 'font';
                        end

                        % handle annotations
                        if strcmp(fn{d},'annotations')
                            annot = stripped.(fn{d});
                            fnmod{d} = 'annotation';
                            for a = 1:length(annot)
                                %recursive call to stripkeys
                                stripped.annotations{a} = obj.stripkeys(annot{a}, fnmod{d}, key);
                            end
                        else
                            %recursive call to stripkeys
                            stripped.(fn{d}) = obj.stripkeys(stripped.(fn{d}), fnmod{d}, key);
                        end

                        % look for desired key and strip if not an exception
                    elseif any(strcmp(keytype, key))
                        if ~isExceptionStrip(stripped,fn{d})
                            stripped = rmfield(stripped, fn{d});
                        end
                    end

                end

                %----CLEAN UP----%
                remfn = fieldnames(stripped);

                for n = 1:length(remfn)
                    if isstruct(stripped.(remfn{n}))
                        if isempty(fieldnames(stripped.(remfn{n})))
                            stripped = rmfield(stripped,remfn{n});
                        end
                    elseif isempty(stripped.(remfn{n}))
                        stripped = rmfield(stripped,remfn{n});
                    end
                end

                %---HANDLE ERRORS---%
            catch exception
                if obj.UserData.Verbose
                    % catch 3D output until integrated into graphref
                    if ~(strcmpi(fieldname,'surface') || strcmpi(fieldname,'scatter3d'))
                        fprintf(['\nWhoops! ' exception.message(1:end-1) ' in ' fieldname '\n\n']);
                    end
                end
            end
        end

        function link_text = get_link_text(obj)
           plotly_domain = obj.UserData.PlotlyDomain;
           link_domain = strrep(plotly_domain, 'https://', '');
           link_domain = strrep(link_domain, 'http://', '');
           link_text = ['Export to ' link_domain];
        end
    end
end

author: goGPS-Project