text:
function [d,si]=abfload(fn,varargin)
% ** function [d,si]=abfload(fn,varargin)
% loads and returns data in the Axon abf format PRIOR TO VERSION 2.0.
% Data may have been acquired in the following modes:
% (1) event-driven variable-length
% (2) event-driven fixed-length
% (3) gap-free
% Information about scaling, the time base and the number of channels and
% episodes is extracted from the header of the abf file (see also abfinfo.m).
% All optional input parameters listed below (= all except the file name)
% must be specified as parameter/value pairs, e.g. as in
%          d=abfload('d:\data01.abf','start',100,'stop','e');
%
%                    >>> INPUT VARIABLES >>>
%
% NAME        TYPE, DEFAULT      DESCRIPTION
% fn          char array         abf data file name
% start       scalar, 0          only gap-free-data: start of cutout to be read (unit: sec)
% stop        scalar or char,    only gap-free-data: end of cutout to be read (unit: sec).
%             'e'                 May be set to 'e' (end of file).
% sweeps      1d-array or char,  only episodic data: sweep numbers to be read. By default,
%             'a'                 all sweeps will be read ('a')
% channels    cell array         names of channels to be read, like {'IN 0','IN 8'};
%              or char, 'a'       ** make sure spelling is 100% correct (including blanks) **
%                                 if set to 'a', all channels will be read
% chunk       scalar, 0.05       only gap-free-data: the elementary chunk size (Megabytes)
%                                 to be used for the 'discontinuous' mode of reading data
%                                 (when less channels shall be read than exist in file)
% machineF    char array,        the 'machineformat' input parameter of the
%              'ieee-le'          matlab fopen function. 'ieee-le' is the correct
%                                 option for windows; depending on the
%                                 platform the data were recorded/shall be read
%                                 by abfload 'ieee-be' is the alternative.
%
%                    <<< OUTPUT VARIABLES <<<
%
% NAME        TYPE            DESCRIPTION
% d                           the data read, the format depending on the recording mode
%             1. GAP-FREE:
%             2d array        2d array of size
%                              '<data pts>' by '<number of chans>'
%                              Examples of access:
%                              d(:,2)       data from channel 2 at full length
%                              d(1:100,:)   first 100 data points from all channels
%             2. EPISODIC FIXED-LENGTH:
%             3d array        3d array of size
%                              '<data pts per sweep>' by '<number of chans>' by '<number of sweeps>'
%                              Examples of access:
%                              d(:,2,:)            is a matrix which contains all events (at full length)
%                                                  of channel 2 in its columns
%                              d(1:200,:,[1 11])   contains first 200 data points of events #1
%                                                  and #11 of all channels
%             3. EPISODIC VARIABLE-LENGTH:
%             cell array      cell array whose elements correspond to single sweeps. Each element is
%                              a (regular) array of size
%                              '<data pts per sweep>' by '<number of chans>'
%                              Examples of access:
%                              d{1}                 a 2d-array which contains sweep #1 (all of it, all channels)
%                              d{2}(1:100,2)        a 1d-array containing first 100 data points of of channel 2 in sweep #1
%
% si          scalar          the sampling interval in usec
%
% Date of this version: June 17, 2011

% future improvements:
% - handle expansion of header in transition to file version 1.65 better

% --- defaults
% gap-free
start=0.0;
stop='e';
% episodic
sweeps='a';
% general
channels='a';
% the size of data chunks (see above) in Mb. 0.05 Mb is an empirical value
% which works well for abf with 6-16 channels and recording durations of
% 5-30 min
chunk=0.05;
machineF='ieee-le';
verbose=1;

% assign values of optional input parameters, if any were given
pvpmod(varargin);

if verbose, disp(['** ' mfilename]); end
d=[];
si=[];
if ischar(stop)
  if ~strcmpi(stop,'e')
    error('input parameter ''stop'' must be specified as ''e'' (=end of recording) or as a scalar');
  end
end

% --- obtain vital header parameters and initialize them with -1
% temporary initializing var
tmp=repmat(-1,1,16);
% for the sake of typing economy set up a cell array (and convert it to a struct below)
% column order is
%        name, position in header in bytes, type, value)
headPar={
  'fFileVersionNumber',4,'float',-1;
  'nOperationMode',8,'int16',-1;
  'lActualAcqLength',10,'int32',-1;
  'nNumPointsIgnored',14,'int16',-1;
  'lActualEpisodes',16,'int32',-1;
  'lFileStartTime',24,'int32',-1;
  'lDataSectionPtr',40,'int32',-1;
  'lSynchArrayPtr',92,'int32',-1;
  'lSynchArraySize',96,'int32',-1;
  'nDataFormat',100,'int16',-1;
  'nADCNumChannels', 120, 'int16', -1;
  'fADCSampleInterval',122,'float', -1;
  'fSynchTimeUnit',130,'float',-1;
  'lNumSamplesPerEpisode',138,'int32',-1;
  'lPreTriggerSamples',142,'int32',-1;
  'lEpisodesPerRun',146,'int32',-1;
  'fADCRange', 244, 'float', -1;
  'lADCResolution', 252, 'int32', -1;
  'nFileStartMillisecs', 366, 'int16', -1;
  'nADCPtoLChannelMap', 378, 'int16', tmp;
  'nADCSamplingSeq', 410, 'int16',  tmp;
  'sADCChannelName',442, 'uchar', repmat(tmp,1,10);
  'fADCProgrammableGain', 730, 'float', tmp;
  'fInstrumentScaleFactor', 922, 'float', tmp;
  'fInstrumentOffset', 986, 'float', tmp;
  'fSignalGain', 1050, 'float', tmp;
  'fSignalOffset', 1114, 'float', tmp;
  'nTelegraphEnable',4512,'int16',tmp;
  'fTelegraphAdditGain',4576,'float',tmp
};

fields={'name','offs','numType','value'};
s=cell2struct(headPar,fields,2);
numOfParams=size(s,1);
clear tmp headPar;

if ~exist(fn,'file'), error(['could not find file ' fn]); end
if verbose, disp(['opening ' fn '..']); end
[fid,messg]=fopen(fn,'r',machineF);

if fid == -1,error(messg);end;
% determine absolute file size
fseek(fid,0,'eof');
fileSz=ftell(fid);
fseek(fid,0,'bof');

% read all vital information in header
% convert names in structure to variables and read value from header
for g=1:numOfParams,
  if fseek(fid, s(g).offs,'bof')~=0,
    fclose(fid);
    error(['something went wrong locating ' s(g).name]);
  end;
  sz=length(s(g).value);
  eval(['[' s(g).name ',n]=fread(fid,sz,''' s(g).numType ''');']);
  if n~=sz,
    fclose(fid);
    error(['something went wrong reading value(s) for ' s(g).name]);
  end;
end;

if lActualAcqLength<nADCNumChannels,
  fclose(fid);
  error('less data points than sampled channels in file');
end
% the numerical value of all recorded channels (numbers 0..15)
recChIdx=nADCSamplingSeq(1:nADCNumChannels);
% the corresponding indices into loaded data d
recChInd=1:length(recChIdx);
% the channel names, e.g. 'IN 8'
recChNames=(reshape(char(sADCChannelName),10,16))';
recChNames=recChNames(recChIdx+1,:);

chInd=[];
eflag=0;
if ischar(channels)
  if strcmp(channels,'a')
    chInd=recChInd;
  else
    fclose(fid);
    error('input parameter ''channels'' must either be a cell array holding channel names or the single character ''a'' (=all channels)');
  end
else
  for i=1:length(channels)
    tmpChInd=strmatch(channels{i},recChNames,'exact');
    if ~isempty(tmpChInd)
      chInd=[chInd tmpChInd];
    else
      % set error flag to 1
      eflag=1;
    end
  end;
end
if eflag
  fclose(fid);
  disp('**** available channels:');
  disp(recChNames);
  disp(' ');
  disp('**** requested channels:');
  disp(strvcat(channels));
  error('at least one of the requested channels does not exist in data file (see above)');
end

% fFileVersionNumber needs a fix - for whatever reason its value is sometimes
% a little less than what it should be (e.g. 1.6499999xxxx instead of 1.65)
fFileVersionNumber=.001*round(fFileVersionNumber*1000);

% gain of telegraphed instruments, if any
if fFileVersionNumber>=1.65
  addGain=nTelegraphEnable.*fTelegraphAdditGain;
  addGain(addGain==0)=1;
else
  addGain=ones(size(fTelegraphAdditGain));
end

% tell me where the data start
blockSz=512;
switch nDataFormat
  case 0
    dataSz=2;  % bytes/point
    precision='int16';
  case 1
    dataSz=4;  % bytes/point
    precision='float32';
  otherwise
    fclose(fid);
    error('invalid number format');
end;
headOffset=lDataSectionPtr*blockSz+nNumPointsIgnored*dataSz;
% fADCSampleInterval is the TOTAL sampling interval
si=fADCSampleInterval*nADCNumChannels;

if ischar(sweeps) && sweeps=='a'
  nSweeps=lActualEpisodes;
  sweeps=1:lActualEpisodes;
else
  nSweeps=length(sweeps);
end;

switch nOperationMode
  case 1
    if verbose,
      disp('data were acquired in event-driven variable-length mode');
    end
    warndlg('function abfload has not yet been thorougly tested for data in event-driven variable-length mode - please double-check that the data loaded is correct','Just a second, please');
    if (lSynchArrayPtr<=0 || lSynchArraySize<=0),
      fclose(fid);
      error('internal variables ''lSynchArraynnn'' are zero or negative');
    end
    switch fSynchTimeUnit
      case 0  % time information in synch array section is in terms of ticks
        synchArrTimeBase=1;
      otherwise % time information in synch array section is in terms of usec
        synchArrTimeBase=fSynchTimeUnit;
    end;
    % the byte offset at which the SynchArraySection starts
    lSynchArrayPtrByte=blockSz*lSynchArrayPtr;
    % before reading Synch Arr parameters check if file is big enough to hold them
    % 4 bytes/long, 2 values per episode (start and length)
    if lSynchArrayPtrByte+2*4*lSynchArraySize<fileSz,
      fclose(fid);
      error('file seems not to contain complete Synch Array Section');
    end
    if fseek(fid,lSynchArrayPtrByte,'bof')~=0,
      fclose(fid);
      error('something went wrong positioning file pointer to Synch Array Section');
    end
    [synchArr,n]=fread(fid,lSynchArraySize*2,'int32');
    if n~=lSynchArraySize*2,
      fclose(fid);
      error('something went wrong reading synch array section');
    end
    % make synchArr a lSynchArraySize x 2 matrix
    synchArr=permute(reshape(synchArr',2,lSynchArraySize),[2 1]);
    % the length of episodes in sample points
    segLengthInPts=synchArr(:,2)/synchArrTimeBase;
    % the starting ticks of episodes in sample points WITHIN THE DATA FILE
    segStartInPts=cumsum([0 (segLengthInPts(1:end-1))']*dataSz)+headOffset;
    % start time (synchArr(:,1)) has to be divided by nADCNumChannels to get true value
    % go to data portion
    if fseek(fid,headOffset,'bof')~=0,
      fclose(fid);
      error('something went wrong positioning file pointer (too few data points ?)');
    end
    for i=1:nSweeps,
      % if selected sweeps are to be read, seek correct position
      if ~isequal(nSweeps,lActualEpisodes),
        fseek(fid,segStartInPts(sweeps(i)),'bof');
      end;
      [tmpd,n]=fread(fid,segLengthInPts(sweeps(i)),precision);
      if n~=segLengthInPts(sweeps(i)),
        warning(['something went wrong reading episode ' int2str(sweeps(i)) ': ' segLengthInPts(sweeps(i)) ' points should have been read, ' int2str(n) ' points actually read']);
      end;
      dataPtsPerChan=n/nADCNumChannels;
      if rem(n,nADCNumChannels)>0,
        fclose(fid);
        error('number of data points in episode not OK');
      end
      % separate channels..
      tmpd=reshape(tmpd,nADCNumChannels,dataPtsPerChan);
      % retain only requested channels
      tmpd=tmpd(chInd,:);
      tmpd=tmpd';
      % if data format is integer, scale appropriately; if it's float, tmpd is fine
      if ~nDataFormat
        for j=1:length(chInd),
          ch=recChIdx(chInd(j))+1;
          tmpd(:,j)=tmpd(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))...
            *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch);
        end;
      end
      % now place in cell array, an element consisting of one sweep with channels in columns
      d{i}=tmpd;
    end;
  case {2,5}
    if nOperationMode==2
      if verbose
        disp('data were acquired in event-driven fixed-length mode');
      end
    else
      if verbose
        disp('data were acquired in waveform fixed-length mode (clampex only)');
      end
    end
    % determine first point and number of points to be read
    startPt=0;
    dataPts=lActualAcqLength;
    dataPtsPerChan=dataPts/nADCNumChannels;
    if rem(dataPts,nADCNumChannels)>0,
      fclose(fid);
      error('number of data points not OK');
    end
    dataPtsPerChanPerSweep=dataPtsPerChan/lActualEpisodes;
    if rem(dataPtsPerChan,lActualEpisodes)>0
      fclose(fid);
      error('number of data points not OK');
    end
    dataPtsPerSweep=dataPtsPerChanPerSweep*nADCNumChannels;
    if fseek(fid,startPt*dataSz+headOffset,'bof')~=0
      fclose(fid);
      error('something went wrong positioning file pointer (too few data points ?)');
    end
    d=zeros(dataPtsPerChanPerSweep,length(chInd),nSweeps);
    % the starting ticks of episodes in sample points WITHIN THE DATA FILE
    selectedSegStartInPts=((sweeps-1)*dataPtsPerSweep)*dataSz+headOffset;
    for i=1:nSweeps,
      fseek(fid,selectedSegStartInPts(i),'bof');
      [tmpd,n]=fread(fid,dataPtsPerSweep,precision);
      if n~=dataPtsPerSweep,
        fclose(fid);
        error(['something went wrong reading episode ' int2str(sweeps(i)) ': ' dataPtsPerSweep ' points should have been read, ' int2str(n) ' points actually read']);
      end;
      dataPtsPerChan=n/nADCNumChannels;
      if rem(n,nADCNumChannels)>0
        fclose(fid);
        error('number of data points in episode not OK');
      end;
      % separate channels..
      tmpd=reshape(tmpd,nADCNumChannels,dataPtsPerChan);
      % retain only requested channels
      tmpd=tmpd(chInd,:);
      tmpd=tmpd';
      % if data format is integer, scale appropriately; if it's float, d is fine
      if ~nDataFormat
        for j=1:length(chInd),
          ch=recChIdx(chInd(j))+1;
          tmpd(:,j)=tmpd(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))...
            *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch);
        end;
      end
      % now fill 3d array
      d(:,:,i)=tmpd;
    end;

  case 3
    if verbose, disp('data were acquired in gap-free mode'); end
    % from start, stop, headOffset and fADCSampleInterval calculate first point to be read
    %  and - unless stop is given as 'e' - number of points
    startPt=floor(1e6*start*(1/fADCSampleInterval));
    % this corrects undesired shifts in the reading frame due to rounding errors in the previous calculation
    startPt=floor(startPt/nADCNumChannels)*nADCNumChannels;
    % if stop is a char array, it can only be 'e' at this point (other values would have
    % been caught above)
    if ischar(stop),
      dataPtsPerChan=lActualAcqLength/nADCNumChannels-floor(1e6*start/si);
      dataPts=dataPtsPerChan*nADCNumChannels;
    else
      dataPtsPerChan=floor(1e6*(stop-start)*(1/si));
      dataPts=dataPtsPerChan*nADCNumChannels;
      if dataPts<=0
        fclose(fid);
        error('start is larger than or equal to stop');
      end
    end
    if rem(dataPts,nADCNumChannels)>0
      fclose(fid);
      error('number of data points not OK');
    end
    if fseek(fid,startPt*dataSz+headOffset,'bof')~=0,
      fclose(fid);
      error('something went wrong positioning file pointer (too few data points ?)');
    end
    % *** decide on the most efficient way to read data:
    % (i) all (of one or several) channels requested: read, done
    % (ii) one (of several) channels requested: use the 'skip' feature of
    % fread
    % (iii) more than one but not all (of several) channels requested:
    % 'discontinuous' mode of reading data. Read a reasonable chunk of data
    % (all channels), separate channels, discard non-requested ones (if
    % any), place data in preallocated array, repeat until done. This is
    % faster than reading the data in one big lump, separating channels and
    % discarding the ones not requested
    if length(chInd)==1 && nADCNumChannels>1
      % --- situation (ii)
      % jump to proper reading frame position in file
      if fseek(fid,(chInd-1)*dataSz,'cof')~=0
        fclose(fid);
        error('something went wrong positioning file pointer (too few data points ?)');
      end
      % read, skipping nADCNumChannels-1 data points after each read
      dataPtsPerChan=dataPts/nADCNumChannels;
      [d,n]=fread(fid,dataPtsPerChan,precision,dataSz*(nADCNumChannels-1));
      if n~=dataPtsPerChan,
        fclose(fid);
        error(['something went wrong reading file (' int2str(dataPtsPerChan) ' points should have been read, ' int2str(n) ' points actually read']);
      end
    elseif length(chInd)/nADCNumChannels<1
      % --- situation (iii)
      % prepare chunkwise upload:
      % preallocate d
      d=repmat(nan,dataPtsPerChan,length(chInd));
      % the number of data points corresponding to the maximal chunk size,
      % rounded off such that from each channel the same number of points is
      % read (do not forget that each data point will by default be made a
      % double of 8 bytes, no matter what the original data format is)
      chunkPtsPerChan=floor(chunk*2^20/8/nADCNumChannels);
      chunkPts=chunkPtsPerChan*nADCNumChannels;
      % the number of those chunks..
      nChunk=floor(dataPts/chunkPts);
      % ..and the remainder
      restPts=dataPts-nChunk*chunkPts;
      restPtsPerChan=restPts/nADCNumChannels;
      % chunkwise row indices into d
      dix=(1:chunkPtsPerChan:dataPtsPerChan)';
      dix(:,2)=dix(:,1)+chunkPtsPerChan-1;
      dix(end,2)=dataPtsPerChan;
      if verbose && nChunk
        disp(['reading file in ' int2str(nChunk) ' chunks of ~' num2str(chunk) ' Mb']);
      end
      % do it: if no remainder exists loop through all rows of dix,
      % otherwise spare last row for the lines below (starting with
      % 'if restPts')
      for ci=1:size(dix,1)-(restPts>0)
        [tmpd,n]=fread(fid,chunkPts,precision);
        if n~=chunkPts
          fclose(fid);
          error(['something went wrong reading chunk #' int2str(ci) ' (' ...
            int2str(chunkPts) ' points should have been read, ' int2str(n) ' points actually read']);
        end
        % separate channels..
        tmpd=reshape(tmpd,nADCNumChannels,chunkPtsPerChan);
        d(dix(ci,1):dix(ci,2),:)=tmpd(chInd,:)';
      end
      % collect the rest
      [tmpd,n]=fread(fid,restPts,precision);
      if n~=restPts
        fclose(fid);
        error(['something went wrong reading last chunk (' ...
          int2str(restPts) ' points should have been read, ' int2str(n) ' points actually read']);
      end
      % separate channels..
      tmpd=reshape(tmpd,nADCNumChannels,restPtsPerChan);
      d(dix(end,1):dix(end,2),:)=tmpd(chInd,:)';
    else
      % --- situation (i)
      [d,n]=fread(fid,dataPts,precision);
      if n~=dataPts,
        fclose(fid);
        error(['something went wrong reading file (' int2str(dataPts) ' points should have been read, ' int2str(n) ' points actually read']);
      end;
      % separate channels..
      d=reshape(d,nADCNumChannels,dataPtsPerChan);
      d=d';
    end
    % if data format is integer, scale appropriately; if it's float, d is fine
    if ~nDataFormat
      for j=1:length(chInd),
        ch=recChIdx(chInd(j))+1;
        d(:,j)=d(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))...
          *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch);
      end;
    end
  otherwise
    disp('recording mode of data must be event-driven variable-length (1), event-driven fixed-length (2) or gap-free (3) -- returning empty matrix');
    d=[];
    si=[];
end;

fclose(fid);


text:
%Demo for the Fractional-Order Chaotic Systems (FOChS) - functions:

%1. Chen's system:
[t, y]=FOChen([35 3 28 -7], [0.9 0.9 0.9], 100, [-9 -5 14]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%2. CNN - 3 cells net:
[t, y]=FO3CNN([1.24 1.1 1.0 4.4 3.21], [0.99 0.99 0.99], 100, [0.1 0.1 0.1]);
figure
plot3(y(:,1), y(:,2), y(:,3), 'k');
xlabel('x_1(t)'); ylabel('x_2(t)'); zlabel('x_3(t)'); grid;

%3. Arneodo's system:
[t,y]=FOArneodo([-5.5 3.5 0.8 -1.0], [0.97 0.97 0.96], 200, [-0.2 0.5 0.2]);
figure
plot3(y(:,1), y(:,2), y(:,3), 'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%4. Genesio-Tesi's system:
[t, y]=FOGenTesi([1.1 1.1 0.45 1.0], [1 1 0.95], 200, [-0.1 0.5 0.2]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%5. Lorenz's system:
[t, y]=FOLorenz([10 28 8/3],[0.993 0.993 0.993],100,[0.1 0.1 0.1]);
figure
plot3(y(:,1), y(:,2), y(:,3), 'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%6. Newton-Leipnik's system:
[t, y]=FONewLeipnik([0.4 0.175], [0.95 0.95 0.95], 200, [0.19 0 -0.18]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%7. Rossler's system:
[t, y]=FORossler([0.5 0.2 10], [0.9 0.85 0.95], 120, [0.5 1.5 0.1]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%8. Lotka-Volterra system:
[t, y]=FOLotkaVolterra([1 1 1 1 2 3 2.7], [0.95 0.95 0.95], 200, [1 1.4 1]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%9. Duffing's system:
[t, y]=FODuffing([0.15 0.3 1], [0.9 1], 200, [0.21 0.31]);
figure
plot(y(:,1), y(:,2), 'k');
xlabel('x(t)'); ylabel('y(t)'); grid;

%10. Van der Pol's oscillator:
[t, y]=FOvanDerPol(1, [1.2 0.8], 60, [0.2 -0.2]);
figure
plot(y(:,1), y(:,2), 'k');
xlabel('y_1(t)'); ylabel('y_2(t)'); grid;

%11. Volta's system:
[t, y]=FOVolta([19 11 0.73],[0.99 0.99 0.99], 20, [8 2 1]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%12. Lu's system:
[t, y]=FOLu([36 3 20], [0.985 0.99 0.98], 60, [0.2 0.5 0.3]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%13. Liu's system:
[t, y]=FOLiu([1 2.5 5 1 4 4], [0.95 0.95 0.95], 100, [0.2 0 0.5]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

%14. Chua's systems:
[t, y]=FOChuaNR([10.725 10.593 0.268 -0.7872 -1.1726], [0.93 0.99 0.92], 60, [0.6 0.1 -0.6]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;
[t, y]=FOChuaM([10 13 0.1 1.5 0.3 0.8], [0.97 0.97 0.97 0.97], 200, [0.8 0.05 0.007 0.6]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;
figure
plot3(y(:,4), y(:,1), y(:,2),'k');
xlabel('w(t)'); ylabel('x(t)'); zlabel('y(t)'); grid;

%15. Financial system:
[t, y]=FOFinanc([1 0.1 1],[1 0.95 0.99],200, [2 -1 1]);
figure
plot3(y(:,1), y(:,2), y(:,3),'k');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)'); grid;

% Author: Dr. Ivo Petras (ivo.petras@tuke.sk), 2010.

author: Sable

repo: mcbench-benchmarks


text:
function [c,ceq] = snd_nl_terminalconstraints(t, x, varargin)
    if nargin >=3
        battery = cell2mat(varargin(1));
    else
    end

    c(1) = x(2)-battery.range(2);
    c(2) = battery.range(1)-x(2);
    c(3) = -x(3);
    c(4) = x(3)-100;
    ceq = [];
end

author: juchengquan


text:
%%  Performing Gaussian blur
%   I: input image
%   s: standard deviation
%   GI: blurred image

function GI = gaussianBlur(I,s)

M = gaussianMask(1,s);
if max(size(M))==0
    GI=I;
    return;
end
M = M/sum(sum(M));   % normalize the gaussian mask
GI=I;
for i=1:(2*s+1)
    GI=BoundMirrorExpand(GI);
end
GI = conv2(GI,M,'same');
for i=1:(2*s+1)
    GI=BoundMirrorShrink(GI);
end


text:
function [STATS]=scatter_stats(D)

%Fitting Gaussian mixture (1) model
GM=gmdistribution.fit(D,1);
MU=GM.mu; %Means
MU_total=rms(MU); %Total mean
VAR=GM.Sigma; %Variance
STD=sqrt(VAR); %Standard deviations
[VAR_vec,VAR_eig] = eig(VAR); %Get principal components
STD_total=sqrt(trace(VAR_eig)./3); %Standard deviations

%Creating output structure
STATS.GM=GM;
STATS.MU=MU;
STATS.MU_total=MU_total;
STATS.VAR=VAR;
STATS.STD=STD;
STATS.STD_total=STD_total;

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% MPC applied to HIV system using a SINDYc model.
% Function for TESTING

clear all, close all, clc
figpath = '../FIGURES/HIV/'; mkdir(figpath)
datapath = '../DATA/HIV/'; mkdir(datapath)
addpath('../utils');

SystemModel = 'HIV';
ModelName = 'SINDYc';

%% Load Model
Nvar = 5;
InputSignalTypeModel = 'prbs'; % prbs; chirp; noise; sine2; sphs; mixed
load(fullfile(datapath,['EX_',SystemModel,'_SI_SINDYc','_',InputSignalTypeModel,'.mat']))

%% TRUE SYSTEM PARAMETERS
run_HIV_params
yout = poolDataLIST({'x1','x2','x3','x4','x5','u'},Models.SINDYc.Xi,6,3,0);

% ZURAKOWSKI
xi_truth.term{1} = {'1', 'x1', 'x1x2', 'x1x2u'};
xi_truth.coeff{1} = [lambda1, -d, -alpha1, eta*alpha1];
xi_truth.term{2} = {'x2', 'x1x2', 'x2x4', 'x2x5', 'x1x2u'};
xi_truth.coeff{2} = [-a, alpha1, -p1, -p2, -eta*alpha1];
xi_truth.term{3} = {'x3', 'x2x3', 'x1x2x3'};
xi_truth.coeff{3} = [-b2, -c2*q, c2];
xi_truth.term{4} = {'x4', 'x2x4'};
xi_truth.coeff{4} = [-b1, c1];
xi_truth.term{5} = {'x5', 'x2x3'};
xi_truth.coeff{5} = [-h, c2*q];

% Construct true Xi
Xi0 = zeros(size(Models.SINDYc.Xi));
for i = 1:Nvar
    for j = 1:length(xi_truth.term{i})
        idx = find(strcmp(yout,xi_truth.term{i}(j)));
        Xi0(idx,i) = xi_truth.coeff{i}(j);
    end
end

%% Apply Model predictive control to system using SINDYc model
select_model = 'SINDYc';
pest.ahat = Model.Xi(:,1:Nvar); % Use Xi0(:,1:Nvar) to execute true model parameters
pest.polyorder = Model.polyorder;
pest.usesine = Model.usesine;
pest.dt = 1/12; %Model.dt; % Don't need to use the same time step model was trained on

% Parameters MPC
options = optimoptions('fmincon','Algorithm','sqp','Display','none', ...
    'MaxIterations',100);
Nweeks = 50;
Duration = Nweeks*7;             % Run for 'Duration' time units
Ton = 0;                         % Time units when control is turned on
x0n = [10, 0.1, 0.1, 0.1, 0.1]'; % Initial condition
Ts  = pest.dt;                   % Sampling time [1/hr]
Tcontrol = 0;
getMPCparams

% Reference state, which shall be achieved
xref = zeros(1,5);
xref(2) = ((c2*(lambda1-d*q)-b2*alpha1) - sqrt((c2*(lambda1-d*q)-b2*alpha1)^2 - 4*alpha1*c2*q*d*b2))/(2*alpha1*c2*q);
xref(1) = lambda1/(d+alpha1*xref(2));
xref(4) = 0;
xref(5) = (xref(2)*c2*(alpha1*q-a) + b2*alpha1)/(c2*p2*xref(2));
xref(3) = h*xref(5)/(c2*q*xref(2));

% Initialize variables
Nt = (Duration/Ts)+1;
uopt0    = 0;
xhat     = x0n;
uopt     = uopt0.*ones(Nu,1);
xHistory = zeros(Nvar,Nt); xHistory(:,1) = xhat;
uHistory = zeros(1,Nt); uHistory(1)   = uopt(1);
tHistory = zeros(1,Nt); tHistory(1)   = Tcontrol;
rHistory = zeros(Nvar,Nt);

%%
% Start simulation
fprintf('Simulation started.  It might take a while...\n')
tic
for ct = 1:(Duration/Ts)   % For each iteration: take measurements & optimize control input & apply control input

    if mod(ct*Ts,7) == 0 % Update once a week
        % NMPC with full-state feedback
        COSTFUN = @(u) ObjectiveFCN_models(u,xhat,N,Nu,xref,uHistory(:,ct),pest,diag(Q),R,Ru,select_model);
        CONSFUN = @(u) ConstraintFCN_models(u,uHistory(:,ct),xhat,N,LBo,UBo,LBdu,UBdu,pest,select_model);
        uopt = fmincon(COSTFUN,uopt,[],[],[],[],LB,UB,CONSFUN,options);
    end

    % Integrate system
    xhat = rk4u(@HIVsys_ZURAKOWSKI,xhat,uopt(1),Ts/2,2,[],0);
    xHistory(:,ct+1) = xhat;
    uHistory(:,ct+1) = uopt(1);
    tHistory(:,ct+1) = ct*Ts+Tcontrol;
    rHistory(:,ct+1) = xref;

    if mod(ct,1000) == 0
        disp(['PROGRESS: ',num2str(100*ct/(Duration/Ts)),'%'])
    end

end
fprintf('Simulation finished!\n')
toc

%%
clear ph
figure;box on,hold on
ccolors = get(gca,'colororder');
plot(tHistory(1:ct)/7,xref(1)*ones(length(tHistory(1:ct)),1)./max(xHistory(1,1:ct)),'--','Color',ccolors(1,:),'LineWidth',1), hold on
plot(tHistory(1:ct)/7,xref(2)*ones(length(tHistory(1:ct)),1)./max(xHistory(2,1:ct)),'--','Color',ccolors(2,:),'LineWidth',1)
plot(tHistory(1:ct)/7,xref(3)*ones(length(tHistory(1:ct)),1)./max(xHistory(3,1:ct)),'--','Color',ccolors(3,:),'LineWidth',1)
ph(1) = plot(tHistory(1:ct)/7,xHistory(1,1:ct)./max(xHistory(1,1:ct)),'-','Color',ccolors(1,:),'LineWidth',1.5);
ph(2) = plot(tHistory(1:ct)/7,xHistory(2,1:ct)./max(xHistory(2,1:ct)),'-','Color',ccolors(2,:),'LineWidth',1.5);
ph(3) = plot(tHistory(1:ct)/7,xHistory(3,1:ct)./max(xHistory(3,1:ct)),'-','Color',ccolors(3,:),'LineWidth',1.5);
ph(4) = plot(tHistory(1:ct)/7,uHistory(1:ct),'-k','LineWidth',1.5);
ylim([0 1])
legend(ph,'x1','x2','x3','u')
%% Show results
clear ph

figure;hold on, box on,
ccolors = get(gca,'colororder');
plot(tHistory,xref1(1)*ones(length(tHistory),1),'--','Color',ccolors(1,:),'LineWidth',1)
plot(tHistory,xref1(2)*ones(length(tHistory),1),'--','Color',ccolors(2,:),'LineWidth',1)
plot(tHistory,xref1(3)*ones(length(tHistory),1),'--','Color',ccolors(3,:),'LineWidth',1)
ph(1) = plot(tHistory,xHistory(1,:),'-','Color',ccolors(1,:),'LineWidth',1.5);
ph(2) = plot(tHistory,xHistory(2,:),'-','Color',ccolors(2,:),'LineWidth',1.5);
ph(3) = plot(tHistory,xHistory(3,:),'-','Color',ccolors(3,:),'LineWidth',1.5);
ph(4) = plot(tHistory,uHistory,'-k','LineWidth',1.5);
xlabel('Time')
ylabel('xi')
legend(ph,'x1','x2','x3','Control')
axis tight
set(gca,'xtick',[50,100,150,200])
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')

%% Save Results
Results.t = tHistory;
Results.x = xHistory;
Results.u = uHistory;
Results.J = evalObjectiveFCN(uHistory,xHistory,rHistory,diag(Q),R,Ru);

author: eurika-kaiser

repo: SINDY-MPC


text:
%CONFMAT Construct confusion matrix
%
%  [C,NE,LABLIST] = CONFMAT(LAB1,LAB2,METHOD,FID)
%
% INPUT
%  LAB1        Set of labels
%  LAB2        Set of labels
%  METHOD      'count' (default) to count number of co-occurences in
%	             LAB1 and LAB2, 'disagreement' to count relative
%		           non-co-occurrence.
%  FID         Write text result to file
%
% OUTPUT
%  C           Confusion matrix
%  NE          Total number of errors (empty labels are neglected)
%  LABLIST     Unique labels in LAB1 and LAB2
%
% DESCRIPTION
% Constructs a confusion matrix C between two sets of labels LAB1
% (corresponding to the rows in C) and LAB2 (the columns in C). The order of
% the rows and columns is returned in LABLIST. NE is the total number of
% errors (sum of non-diagonal elements in C).
%
% When METHOD = 'count' (default), co-occurences in LAB1 and LAB2 are counted
% and returned in C. When METHOD = 'disagreement', the relative disagreement
% is returned in NE, and is split over all combinations of labels in C
% (such that the rows sum to 1). (The total disagreement for a class equals
% one minus the sensitivity for that class as computed by TESTC).
%
%   [C,NE,LABLIST] = CONFMAT(D,METHOD)
%
% If D is a classification result D = A*W, the labels LAB1 and LAB2 are
% internally retrieved by CONFMAT before computing the confusion matrix.
%
% When no output argument is specified, or when FID is given, the
% confusion matrix is displayed or written a to a text file. It is assumed
% that LAB1 contains true labels and LAB2 stores estimated labels.
%
% EXAMPLE
% Typical use of CONFMAT is the comparison of true and and estimated labels
% of a testset A by application to a trained classifier W:
% LAB1 = GETLABELS(A); LAB2 = A*W*LABELD.
% More examples can be found in PREX_CONFMAT, PREX_MATCHLAB.
%
% SEE ALSO
% MAPPINGS, DATASETS, GETLABELS, LABELD

% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: confmat.m,v 1.7 2008/10/14 21:34:32 duin Exp $

function [CC,ne,lablist,lablist_true] = confmat_new (arg1,arg2,arg3,fid,lablist_true)

	prtrace(mfilename);

	% Check arguments.
  if nargin < 5, lablist_true = []; end
  if nargin < 4, fid = 1; end
	if nargin < 3 | isempty(arg3)
		if isdataset(arg1)
			lab1 = getlabels(arg1); lab2 = arg1*labeld;
			if nargin < 2| isempty(arg2)
				method = 'count';
				prwarning(4,'no method supplied, assuming count');
			else
				method = arg2;
			end
		else
			method = 'count';
			prwarning(4,'no method supplied, assuming count');
			lab1 = arg1;
			if (nargin < 2 | isempty(arg2))
				error('prtools_addin','Second label list not supplied')
			end
			lab2 = arg2;
		end
	else
		lab1 = arg1;
		lab2 = arg2;

        method = arg3;
	end
	if nargin < 2
		if ~isdataset(arg1)
			error('prtools_addin','two labellists or one dataset should be supplied')
		end
	end

	% Renumber LAB1 and LAB2 and find number of unique labels.

	m = size(lab1,1);
	if (m~=size(lab2,1))
		error('prtools_addin','LAB1 and LAB2 have to have the same lengths.');
    end

    if( isempty(lablist_true) )

        [nlab1,nlab2,lablist] = renumlab(lab1,lab2);
    % 	n = max(nlab1);
    % 	n = max(nlab2);
    % 	n = max(n,n);
        n = size(lablist,1);

    else
        nlab1 = renumlab(lab1,lablist_true);
        nlab2 = renumlab(lab2,lablist_true);
        lablist = lablist_true;
        n = size(lablist,1);
    end

	% Construct matrix of co-occurences (confusion matrix).

	C = zeros(n+1,n+1);
	for i = 0:n
		K = find(nlab1==i);
		if (isempty(K))
			C(i+1,:) = zeros(1,n+1);
		else
			for j = 0:n
				C(i+1,j+1) = length(find(nlab2(K)==j));
			end
		end
	end

	% position rejects and unlabeled object at the end of the matrix

	D = C;
	D(1:end-1,1:end-1) = C(2:end,2:end);
	D(end,:) = [C(1,2:end) C(1,1)];
	D(1:end-1,end) = C(2:end,1);
	C = D;
	D = D(1:end-1,1:end-1);
	DD = D(1:min(n,n),1:min(n,n));
	% Calculate number of errors ('count') or disagreement ('disagreement').
	% Neglect rejects

  if nargout > 1
    J = find(nlab1~=0 & nlab2~=0);
    ne = nlabcmp(lab1(J,:),lab2(J,:));
  end

	switch (method)
		case 'count'
			ne = sum(sum(D)) - sum(diag(DD));       % Diagonal entries are correctly
                                                % classified, so all off-diagonal
                                                % entries denote wrong ones.

		case 'disagreement'
			ne = (sum(sum(D)) - sum(diag(DD)))/m;   % Relative sum of off-diagonal
                                             % entries.
			E = repmat(sum(D,2),1,n);              % Disagreement = 1 -
			D = ones(n,n)-D./E;                    % relative co-occurence.
			D = D / (n-1);
		otherwise
			error('prtools_addin','unknown method');
	end

	%Distinguish 'rejects / no_labels' from 'non_rejects / fully labeled'
	if (any(C(:,end) ~= 0) | any(C(end,:)~=0)) & strcmp(method,'count')
		n = n+1;
        n_real = n;
		labch = char(strlab(lablist),'reject');
		labcv = char(strlab(lablist),'No');
    else
        if( isempty(lablist_true) )
            labcv = strlab(lablist);
            n_real = n;
        else
            labcv = lablist_true;
            n_real = size(lablist_true,1);
        end
        labch = strlab(lablist);
		%labcv = labch;
		C = D;
	end

    % If no output argument is specified, pretty-print C.

	if (nargout == 0) | nargin == 4

		if nargin < 4, fid = 1; end

		% Make sure labels are stored in LABC as matrix of characters,
    % max. 6 per label.

		if (size(labch,2) > 6)
			labch = labch(:,1:6);
			%labcv = labcv(:,1:6);
		end
		if (size(labch,2) < 5)
			labch = [labch repmat(' ',n,ceil((5-size(labch,2))/2))];
% 			labcv = [labcv repmat(' ',n,ceil((5-size(labcv,2))/2))];
		end

		%C = round(1000*C./repmat(sum(C,2),1,size(C,2)));

		nspace = max(size(labcv,2)-7,0);
		cspace = repmat(' ',1,nspace);
		%fprintf(fid,['\n' cspace '        | Estimated Labels']);
		fprintf(fid,['\n  True   ' cspace '| Estimated Labels']);
		fprintf(fid,['\n  Labels ' cspace '|']);
		for j = 1:n, fprintf(fid,'%7s',labch(j,:)); end
		fprintf(fid,'|');
		fprintf(fid,' Totals');
		fprintf(fid,'\n ');
		fprintf(fid,repmat('-',1,8+nspace));
		fprintf(fid,'|%s',repmat('-',1,7*n));
		fprintf(fid,'|-------');
		fprintf(fid,'\n ');

		for j = 1:min(n,n_real)
			fprintf(fid,' %-7s|',labcv(j,:));
			switch (method)
				case 'count'
					fprintf(fid,'%5i  ',C(j,:)');
					fprintf(fid,'|');
					fprintf(fid,'%5i',sum(C(j,:)));
				case 'disagreement'
					fprintf(fid,' %5.3f ',C(j,:)');
					fprintf(fid,'|');
					fprintf(fid,' %5.3f ',sum(C(j,:)));
			end
			fprintf(fid,'\n ');
		end

		fprintf(fid,repmat('-',1,8+nspace));
		fprintf(fid,'|%s',repmat('-',1,7*n));
		fprintf(fid,'|-------');
		fprintf(fid,['\n  Totals ' cspace '|']);

		switch (method)
			case 'count'
				fprintf(fid,'%5i  ',sum(C));
				fprintf(fid,'|');
				fprintf(fid,'%5i',sum(C(:)));
			case 'disagreement'
				fprintf(fid,' %5.3f ',sum(C));
				fprintf(fid,'|');
				fprintf(fid,' %5.3f ',sum(C(:)));
		end
		fprintf(fid,'\n\n');
	end

	if nargout > 0
		CC = C;
	end

return

author: marianux

repo: ecg-kit

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


text:
function w = arrivals2waveforms(dbpath, arrivals, pretrig, posttrig, taper_seconds, nwaveforms)
% ARRIVALS2WAVEFORMS Load waveform objects corresponding to an arrivals structure
%	w = ARRIVALS2WAVEFORMS(dbpath, arrivals, pretrig, posttrig, taper_seconds, maxNumberOfWaveforms)
	w = waveform(); % Initialize output

	% Check input variables.
	if nargin<4
		disp('Incorrect number of input variables')
		return
	end

	% Create datasource
	ds = datasource('antelope', dbpath);

	% Loop over arrivals structure
	for i=1:min([nwaveforms numel(arrivals.arid)])
        scnl = scnlobject(arrivals.sta{i}, arrivals.chan{i});
		snum = epoch2datenum(arrivals.time(i) - pretrig - taper_seconds);
		enum = epoch2datenum(arrivals.time(i) + posttrig + taper_seconds);
        try
            thisw = waveform(ds, scnl, snum, enum);
        end
        if numel(thisw)==1
            w(i) = thisw;
        else
            w(i)=waveform();
        end
	end
end

author: geoscience-community-codes

repo: GISMO


text:
%% meshCleave
% Below is a demonstration of the features of the |meshCleave| function

%%
clear; close all; clc;

%% Syntax
% |[logicAt]=meshCleave(E,V);|
% |[logicAt]=meshCleave(E,V,P,n);|
% |[logicAt,logicAbove,logicBelow]=meshCleave(E,V,P,n,inclusiveSwitch);|
% |[logicAt,logicAbove,logicBelow]=meshCleave(E,V,P,n,inclusiveSwitch);|

%% Description
% This function creates logic arrays for the mesh components (e.g. elements
% or faces) which are at, above, or below a plane defined by the point P,
% and the normal direction n.
% The optional inclusiveSwitch is a 2-component vector (default [0 0]) and
% sets how "inclusive", the below/above logic is, i.e. they set wether <
% and > is used ([0 0]), or <= and >= are used ([1 1]). A combination may
% also be used e.g. [1 0] results in below checks which features <= and
% above checks using >.

%% Examples
clear; close all; clc;

%% Example: Cleaving a surface mesh

%%
% Create example patch data
[F,V]=stanford_bunny;

%%

n=[0 0 1]; %Normal direction to plane
P=mean(V,1); %Point on plane
[logicAt,logicAbove,logicBelow]=meshCleave(F,V,P,n);

%%
% Visualize
cFigure;
subplot(1,3,1); hold on;
title('At cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F(logicAt,:),V,'bw','k',1);
camlight headlight;
axisGeom;

subplot(1,3,2); hold on;
title('Above cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F(logicAbove,:),V,'bw','k',1);
camlight headlight;
axisGeom;

subplot(1,3,3); hold on;
title('Below cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F(logicBelow,:),V,'bw','k',1);
camlight headlight;
axisGeom;

gdrawnow;

%% Example 2: Cleaving a tetrahedral mesh

%%
% Creating an example tetrahedral mesh
boxDim=[5 5 5]; % Box dimenstions
pointSpacing=0.25;
[meshStruct]=tetMeshBox(boxDim,pointSpacing);
E=meshStruct.elements;
V=meshStruct.nodes;
F=meshStruct.facesBoundary;
VE=patchCentre(E,V);
C=minDist(VE,mean(VE,1));

%%

n=[0 0 1]; %Normal direction to plane
P=mean(V,1); %Point on plane
[logicAt,logicAbove,logicBelow]=meshCleave(E,V,P,n);

% Get faces and matching color data for visualization
[F_cleave,CF_cleave]=element2patch(E(logicAt,:),C(logicAt));
[F_above,CF_above]=element2patch(E(logicAbove,:),C(logicAbove));
[F_below,CF_below]=element2patch(E(logicBelow,:),C(logicBelow));

%%
% Visualize
cFigure;
subplot(1,3,1); hold on;
title('At cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F_cleave,V,CF_cleave,'k',1);
camlight headlight;
axisGeom;

subplot(1,3,2); hold on;
title('Above cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F_above,V,CF_above,'k',1);
camlight headlight;
axisGeom;

subplot(1,3,3); hold on;
title('Below cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F_below,V,CF_below,'k',1);
camlight headlight;
axisGeom;

gdrawnow;

%%
% Visualizing cleaving operation for varying angles

hf=cFigure; hold on;
gpatch(F,V,'w','none',0.2);
hp1=gpatch(F_cleave,V,CF_cleave,'k',1);
axisGeom; axis manual;
camlight headligth;
gdrawnow;

nSteps=50; %Number of animation steps

%Create the time vector
animStruct.Time=linspace(0,1,nSteps);

%The vector lengths
a=linspace(0,2*pi,nSteps);
b=linspace(0,2*pi,nSteps);
for q=1:1:nSteps
    R=euler2DCM([a(q) b(q) 0]);
    nn=n*R;

    logicAt=meshCleave(E,V,P,nn,[1 0]);

    % Get faces and matching color data for cleaves elements
    [F_cleave,CF_cleave]=element2patch(E(logicAt,:),C(logicAt));

    %Set entries in animation structure
    animStruct.Handles{q}=[hp1 hp1]; %Handles of objects to animate
    animStruct.Props{q}={'Faces','CData'}; %Properties of objects to animate
    animStruct.Set{q}={F_cleave,CF_cleave}; %Property values for to set in order to animate
end
anim8(hf,animStruct);

%% Example 3: Cleaving a hexahedral mesh

boxDim=[5 5 5]; % Box dimenstions
boxEl=[20 20 20];
[meshStruct]=hexMeshBox(boxDim,boxEl,2);

E=meshStruct.elements;
V=meshStruct.nodes;
F=meshStruct.facesBoundary;
VE=patchCentre(E,V);
C=minDist(VE,mean(VE,1));

%%

n=[0 0 1]; %Normal direction to plane
P=mean(V,1); %Point on plane
inclusiveSwitch=[1 0];
[logicAt,logicAbove,logicBelow]=meshCleave(E,V,P,n,inclusiveSwitch);

logicPlot=logicAt;

% Get faces and matching color data for visualization
[F_cleave,CF_cleave]=element2patch(E(logicAt,:),C(logicAt));
[F_above,CF_above]=element2patch(E(logicAbove,:),C(logicAbove));
[F_below,CF_below]=element2patch(E(logicBelow,:),C(logicBelow));

%%
% Visualize
cFigure;
subplot(1,3,1); hold on;
title('At cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F_cleave,V,CF_cleave,'k',1);
camlight headlight;
axisGeom;

subplot(1,3,2); hold on;
title('Above cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F_above,V,CF_above,'k',1);
camlight headlight;
axisGeom;

subplot(1,3,3); hold on;
title('Below cleave plane');
gpatch(F,V,'w','none',0.2);
gpatch(F_below,V,CF_below,'k',1);
camlight headlight;
axisGeom;

gdrawnow;

%%
% Visualizing slicing operation for varying angles

hf=cFigure; hold on;
gpatch(F,V,'w','none',0.2);
hp1=gpatch(F_cleave,V,CF_cleave,'k',1);
axisGeom; axis manual;
camlight headligth;
gdrawnow;

nSteps=50; %Number of animation steps

%Create the time vector
animStruct.Time=linspace(0,1,nSteps);

%The vector lengths
a=linspace(0,2*pi,nSteps);
b=linspace(0,2*pi,nSteps);
for q=1:1:nSteps
    R=euler2DCM([a(q) b(q) 0]);
    nn=n*R;

    logicAt=meshCleave(E,V,P,nn,inclusiveSwitch);

    % Get faces and matching color data for cleaves elements
    [F_cleave,CF_cleave]=element2patch(E(logicAt,:),C(logicAt));

    %Set entries in animation structure
    animStruct.Handles{q}=[hp1 hp1]; %Handles of objects to animate
    animStruct.Props{q}={'Faces','CData'}; %Properties of objects to animate
    animStruct.Set{q}={F_cleave,CF_cleave}; %Property values for to set in order to animate
end
anim8(hf,animStruct);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function outlet = onl_lslsend_init(streamname,samplingrate,numchans,datatype)
if nargin<2
    samplingrate = 0;
end
if nargin<3
    numchans = 1;
end
if nargin<4
    datatype = 'cf_string';  % can also be cft_float32, etc...
end

% instantiate the library
lib = lsl_loadlib();
info = lsl_streaminfo(lib,streamname,'MATLAB-serialized',numchans,samplingrate,datatype,'whatevergd456546');
outlet = lsl_outlet(info);

repo: REST


text:
function idx = argmin(x)

    [~,idx] = min(x);

author: gordonberman

repo: MotionMapper


text:
function fx = p00_fun ( problem, n, x )

%*****************************************************************************80
%
%% P00_FUN evaluates the integrand for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 September 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(2,N), the evaluation points.
%
%    Output, real FX(N,1), the integrand values.
%
  fx = zeros ( n, 1 );

  if ( problem == 1 )
    fx = p01_fun ( n, x );
  elseif ( problem == 2 )
    fx = p02_fun ( n, x );
  elseif ( problem == 3 )
    fx = p03_fun ( n, x );
  elseif ( problem == 4 )
    fx = p04_fun ( n, x );
  elseif ( problem == 5 )
    fx = p05_fun ( n, x );
  elseif ( problem == 6 )
    fx = p06_fun ( n, x );
  elseif ( problem == 7 )
    fx = p07_fun ( n, x );
  elseif ( problem == 8 )
    fx = p08_fun ( n, x );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_FUN - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem index = %d\n', problem );
    error ( 'P00_FUN - Fatal error!\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function alpha = circ_ang2rad(alpha)

% alpha = circ_ang2rad(alpha)
%   converts values in degree to radians
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html

alpha = alpha * pi /180;

author: brainstorm-tools


text:
function calpak_test77 ( )

%*****************************************************************************80
%
%% CALPAK_TEST77 tests YMDF_TO_WEEKDAY_ENGLISH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST77\n' );
  fprintf ( 1, '  For the English calendar:\n' );
  fprintf ( 1, '  YMDF_TO_WEEKDAY_ENGLISH returns the day of the week.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED   YMDF    Day of the week\n' );
  fprintf ( 1, '\n' );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0 )
      break
    end

    jed2 = jed_to_next_noon ( jed1 );

    [ y1, m1, d1, f1 ] = jed_to_ymdf_english ( jed2 );
    s1 = ymdf_to_s_english ( y1, m1, d1, f1 );

    w2 = ymdf_to_weekday_english ( y1, m1, d1, f1 );
    s2 = weekday_to_name_common ( w2 );

    fprintf ( 1, '  %11.2f  %20s  %2d  %20s\n', jed2, s1, w2, s2 );

  end

  return
end


text:
function imdb = vocSetup(varargin)
opts.edition = '07' ;
opts.dataDir = fullfile('data','voc07') ;
opts.archiveDir = fullfile('data','archives') ;
opts.includeDetection = false ;
opts.includeSegmentation = false ;
opts.includeTest = false ;
opts = vl_argparse(opts, varargin) ;

% Download data
if ~exist(fullfile(opts.dataDir,'Annotations')), download(opts) ; end

% Source images and classes
imdb.paths.image = esc(fullfile(opts.dataDir, 'JPEGImages', '%s.jpg')) ;
imdb.sets.id = uint8([1 2 3]) ;
imdb.sets.name = {'train', 'val', 'test'} ;
imdb.classes.id = uint8(1:20) ;
imdb.classes.name = {...
  'aeroplane', 'bicycle', 'bird', 'boat', 'bottle', 'bus', 'car', ...
  'cat', 'chair', 'cow', 'diningtable', 'dog', 'horse', 'motorbike', ...
  'person', 'pottedplant', 'sheep', 'sofa', 'train', 'tvmonitor'} ;
imdb.classes.images = cell(1,20) ;
imdb.images.id = [] ;
imdb.images.name = {} ;
imdb.images.set = [] ;
index = containers.Map() ;
[imdb, index] = addImageSet(opts, imdb, index, 'train', 1) ;
[imdb, index] = addImageSet(opts, imdb, index, 'val', 2) ;
if opts.includeTest, [imdb, index] = addImageSet(opts, imdb, index, 'test', 3) ; end

% Source segmentations
if opts.includeSegmentation
  n = numel(imdb.images.id) ;
  imdb.paths.objectSegmentation = esc(fullfile(opts.dataDir, 'SegmentationObject', '%s.png')) ;
  imdb.paths.classSegmentation = esc(fullfile(opts.dataDir, 'SegmentationClass', '%s.png')) ;
  imdb.images.segmentation = false(1, n) ;
  [imdb, index] = addSegmentationSet(opts, imdb, index, 'train', 1) ;
  [imdb, index] = addSegmentationSet(opts, imdb, index, 'val', 2) ;
  if opts.includeTest, [imdb, index] = addSegmentationSet(opts, imdb, index, 'test', 3) ; end
end

% Compress data types
imdb.images.id = uint32(imdb.images.id) ;
imdb.images.set = uint8(imdb.images.set) ;
for i=1:20
  imdb.classes.images{i} = uint32(imdb.classes.images{i}) ;
end

% Source detections
if opts.includeDetection
  imdb.aspects.id = uint8(1:5) ;
  imdb.aspects.name = {'front', 'rear', 'left', 'right', 'misc'} ;
  imdb = addDetections(opts, imdb) ;
end

% Check images on disk and get their size
imdb = getImageSizes(imdb) ;

% -------------------------------------------------------------------------
function download(opts)
% -------------------------------------------------------------------------
switch opts.edition
  case '07',
    url='' ;
  case '08',
    url='' ;
  case '09'
    url='' ;
  case '10'
    url='http://host.robots.ox.ac.uk:8080/pascal_trainval/VOCtrainval_03-May-2010.tar' ;
  case '11'
    url='http://host.robots.ox.ac.uk:8080/pascal_trainval/VOCtrainval_25-May-2011.tar' ;
  case '12'
    url='http://host.robots.ox.ac.uk:8080/pascal_trainval/VOCtrainval_11-May-2012.tar' ;
end
trainvalData=sprintf('VOC20%strainval.tar',opts.edition) ;
testData=sprintf('VOC20%stest.tar', opts.edition) ;

if ~exist(opts.archiveDir), mkdir(opts.archiveDir) ; end
archivePath=fullfile(opts.archiveDir, trainvalData) ;
if ~exist(archivePath)
  fprintf('%s: downloading %s to %s\n', mfilename, url, archivePath) ;
  urlwrite(url, archivePath) ;
end
fprintf('%s: decompressing and rearranging %s\n', mfilename, archivePath) ;
untar(archivePath, opts.dataDir) ;
switch opts.edition
  case '11'
    % Decompresses in TrainVal/VOCdevkit, Test/VOCdevkit (sigh!)
    movefile(fullfile(opts.dataDir, 'TrainVal'), fullfile(opts.dataDir,'Test')) ;
  otherwise
    % Decomporesses directly in VOCdevkit
end

if opts.includeTest
  archivePath=fullfile(opts.archiveDir, testData) ;
  if ~exist(archivePath)
    error('Cannot download the test data automatically. Please download the file %s manually from the PASCAL test server.', archivePath)
  end
  fprintf('%s: decompressing and rearranging %s\n', mfilename, archivePath) ;
  untar(archivePath, opts.dataDir) ;
end
switch opts.edition
  case '11'
    movefile(fullfile(opts.dataDir, 'Test', 'VOCdevkit', sprintf('VOC20%s',opts.edition), '*'), opts.dataDir) ;
    rmdir(fullfile(opts.dataDir, 'Test'),'s') ;
  otherwise
    movefile(fullfile(opts.dataDir, 'VOCdevkit', sprintf('VOC20%s',opts.edition), '*'), opts.dataDir) ;
    rmdir(fullfile(opts.dataDir, 'VOCdevkit'),'s') ;
end

% -------------------------------------------------------------------------
function [imdb, index] = addImageSet(opts, imdb, index, setName, setCode)
% -------------------------------------------------------------------------
j = length(imdb.images.id) ;
for ci = 1:length(imdb.classes.name)
  className = imdb.classes.name{ci} ;
  annoPath = fullfile(opts.dataDir, 'ImageSets', 'Main', ...
    [className '_' setName '.txt']) ;
  fprintf('%s: reading %s\n', mfilename, annoPath) ;
  [names,labels] = textread(annoPath, '%s %f') ;
  for i=1:length(names)
    if ~index.isKey(names{i})
      j = j + 1 ;
      index(names{i}) = j ;
      imdb.images.id(j) = j ;
      imdb.images.set(j) = setCode ;
      imdb.images.name{j} = names{i} ;
      imdb.images.classification(j) = true ;
    else
      j = index(names{i}) ;
    end
    if labels(i) > 0, imdb.classes.images{ci}(end+1) = j ; end
  end
end

% -------------------------------------------------------------------------
function [imdb, index] = addSegmentationSet(opts, imdb, index, setName, setCode)
% -------------------------------------------------------------------------
segAnnoPath = fullfile(opts.dataDir, 'ImageSets', 'Segmentation', [setName '.txt']) ;
fprintf('%s: reading %s\n', mfilename, segAnnoPath) ;
segNames = textread(segAnnoPath, '%s') ;
j = numel(imdb.images.id) ;
for i=1:length(segNames)
  if index.isKey(segNames{i})
    k = index(segNames{i}) ;
    imdb.images.segmentation(k) = true ;
    imdb.images.set(k) = setCode ;
  else
    j = j + 1 ;
    index(segNames{i}) = j ;
    imdb.images.id(j) = j ;
    imdb.images.set(j) = setCode ;
    imdb.images.name{j} = segNames{i} ;
    imdb.images.classification(j) = false ;
    imdb.images.segmentation(j) = true ;
  end
end

% -------------------------------------------------------------------------
function imdb = getImageSizes(imdb)
% -------------------------------------------------------------------------
for j=1:numel(imdb.images.id)
  info = imfinfo(sprintf(imdb.paths.image, imdb.images.name{j})) ;
  imdb.images.size(:,j) = uint16([info.Width ; info.Height]) ;
  fprintf('%s: checked image %s [%d x %d]\n', mfilename, imdb.images.name{j}, info.Height, info.Width) ;
end

% -------------------------------------------------------------------------
function imdb = addDetections(opts, imdb)
% -------------------------------------------------------------------------
rois = {} ;
k = 0 ;
fprintf('%s: getting detections for %d images\n', mfilename, numel(imdb.images.id)) ;
for j=1:numel(imdb.images.id)
  fprintf('.') ; if mod(j,80)==0,fprintf('\n') ; end
  name = imdb.images.name{j} ;
  annoPath = fullfile(opts.dataDir, 'Annotations', [name '.xml']) ;
  if ~exist(annoPath, 'file')
    if imdb.images.classification(j) && imdb.images.set(j) ~= 3
      warning('Could not find detection annotations for image ''%s''. Skipping.', name) ;
    end
    continue ;
  end

  doc = xmlread(annoPath) ;
  x = parseXML(doc, doc.getDocumentElement()) ;

  %   figure(1) ; clf ;
  %   imagesc(imread(sprintf(imdb.paths.image,imdb.images.name{j}))) ;

  for q = 1:numel(x.object)
    xmin = sscanf(x.object(q).bndbox.xmin,'%d') ;
    ymin = sscanf(x.object(q).bndbox.ymin,'%d') ;
    xmax = sscanf(x.object(q).bndbox.xmax,'%d') - 1 ;
    ymax = sscanf(x.object(q).bndbox.ymax,'%d') - 1 ;

    k = k + 1 ;
    roi.id = k ;
    roi.image = imdb.images.id(j) ;
    roi.class = find(strcmp(x.object(q).name, imdb.classes.name)) ;
    roi.box = [xmin;ymin;xmax;ymax] ;
    roi.difficult = logical(sscanf(x.object(q).difficult,'%d')) ;
    roi.truncated = logical(sscanf(x.object(q).truncated,'%d')) ;
    if isfield(x.object(q),'occluded')
      roi.occluded = logical(sscanf(x.object(q).occluded,'%d')) ;
    else
      roi.occluded = false ;
    end
    switch x.object(q).pose
      case 'frontal', roi.aspect = 1 ;
      case 'rear', roi.aspect = 2 ;
      case {'sidefaceleft', 'left'}, roi.aspect = 3 ;
      case {'sidefaceright', 'right'}, roi.aspect = 4 ;
      case {'','unspecified'}, roi.aspect = 5 ;
      otherwise, error(sprintf('Unknown view ''%s''', x.object(q).pose)) ;
    end
    rois{k} = roi ;

    %     hold on ;
    %     label=sprintf('%s %s d:%d t:%d o:%d',...
    %       imdb.classes.name{roi.class}, ...
    %       imdb.aspects{roi.aspect}, ...
    %       roi.difficult, ...
    %       roi.truncated, ...
    %       roi.occluded) ;
    %     plotbox(roi.box,'label',label) ;
    %     drawnow ;
  end
end
fprintf('\n') ;

rois = horzcat(rois{:}) ;
imdb.objects = struct(...
    'id', uint32([rois.id]), ...
    'image', uint32([rois.image]), ...
    'class', uint8([rois.class]), ...
    'box', single([rois.box]), ...
    'difficult', [rois.difficult], ...
    'truncated', [rois.truncated], ...
    'occluded', [rois.occluded], ...
    'aspect', uint8([rois.aspect])) ;

% -------------------------------------------------------------------------
function value = parseXML(doc, x)
% -------------------------------------------------------------------------
text = {''} ;
opts = struct ;
if x.hasChildNodes
  for c = 1:x.getChildNodes().getLength()
    y = x.getChildNodes().item(c-1) ;
    switch y.getNodeType()
      case doc.TEXT_NODE
        text{end+1} = lower(char(y.getData())) ;
      case doc.ELEMENT_NODE
        param = lower(char(y.getNodeName())) ;
        if strcmp(param, 'part'), continue ; end
        value = parseXML(doc, y) ;
        if ~isfield(opts, param)
          opts.(param) = value ;
        else
          opts.(param)(end+1) = value ;
        end
    end
  end
  if numel(fieldnames(opts)) > 0
    value = opts ;
  else
    value = strtrim(horzcat(text{:})) ;
  end
end

% -------------------------------------------------------------------------
function str=esc(str)
% -------------------------------------------------------------------------
str = strrep(str, '\', '\\') ;

author: vlfeat

repo: matconvnet-fcn


text:
% Count the number of multiple edges in the graph.
%
% INPUT: adjacency matrix, nxn
% OUTPUT: integer, number of multiple edges
%
% Examples: multiEdges([0 2; 2 0])=1, and multiEdges([0 0 1; 2 0 0; 0 1 0])=1
%
% Note 1: The definition of number of multi-arcs/edges (node pairs
%         that have multiple edges across them) here is:
%         mA = length(find(adj>1)); (normalized by 2 if the graph is directed).
%
% Note 2: This creates a natural difference in counting for
% undirected and directed graphs.
%
% Other routines used: isSymmetric.m
% GB: last updated, Sep 26 2014

function mE=multiEdges(adj)

if isSymmetric(adj) % here use "is symmetric" as surrogate for "is directed"
  mE=length(find(adj>1))/2;
else
  mE=length(find(adj>1));
end

author: aeolianine

repo: octave-networks-toolbox


text:
function [eu,cEu]=grIsEulerian(E)
% Function eu=grIsEulerian(E) returns 1 for Eulerian graph,
% 0.5 for semi-Eulerian and 0 otherwise.
% Input parameter:
%   E(m,2) - the edges of graph;
%     1st and 2nd elements of each row is numbers of vertexes;
%     m - number of edges.
% Output parameter:
%   eu = 1 for Eulerian graph;
%   eu = 0.5 for semi-Eulerian graph;
%   eu = 0 otherwise.
% The graph is Eulerian, if he is connected and powers of
% all vertexes is even.
% The graph is semi-Eulerian, if he is connected and only
% two vertexes have odd powers, and other powers of vertexes is even.
% [eu,cEu]=grIsEulerian(E) return also
%   cEu(m,1) - the vector-column with numbers of edges
%     included to Eulerian cycle (if eu=1) or Eulerian way (if eu=0.5)
%     and cEu=[] otherwise. The Fleury algorithm is used.
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

if nargin<1,
  error('There are no input data!')
end
eu = 0; % the graph is not Eulerian
cEu=[];
ncV=grComp(E); % the number of components
if max(ncV)>1, % 2 or more components
  return % the graph is not Eulerian
end
n=max(max(E(:,1:2))); % number of vertexes
E1=find([diff(sort(E(:)));1]);
p=[E1(1);diff(E1)]; % powers of vertexes
rp=rem(p,2); % remainder after division to 2
srp=sum(rp); % summa of remainders
switch srp
  case 0, % Eulerian graph
    eu=1;
  case 2, % semi-Eulerian graph
    eu=0.5;
  otherwise, % not Eulerian graph
    return
end
%=========== we find the Eulerian cycle or way =============
if nargout>1,
  if srp==0, % Eulerian graph
    v1=1; % first vertex of Eulerian cycle
  else % semi-Eulerian graph
    v1=find(rp);
    v1=v1(1); % first vertex of Eulerian way
  end
  vc=v1; % the current vertex
  m=size(E,1); % number of edges
  E1=[E(:,1:2), [1:m]']; % all edges with numbers
  while ~isempty(E1), % the Fleury algorithm
    evc=find((E1(:,1)==vc)|(E1(:,2)==vc)); % edges connected with vc
    levc=length(evc); % number of edges connected with vertex vc
    if levc==1, % only one way
      cEu=[cEu;E1(evc,3)]; % we add new edge to Eulerian cycle (way)
      vcold=vc;
      vc=sum(E1(evc,1:2))-vc; % new current vertex
      E1=E1(setdiff([1:size(E1,1)],evc),:); % we delete isolated vertex
      E2=E1(:,1:2);
      E2gv=E2>vcold;
      E2(E2gv)=E2(E2gv)-1;
      E1(:,1:2)=E2;
      if vc>vcold,
        vc=vc-1;
      end
      if v1>vcold,
        v1=v1-1;
      end
    else % several ways from vertex vc
      for k=1:levc,
        E2=E1(setdiff([1:size(E1,1)],evc(k)),:);
        ncv=grComp(E2); % number of components
        nco=max(ncv);
        if (max(ncv)==1), % this edge is not bridge
          cEu=[cEu;E1(evc(k),3)]; % we add new edge to Eulerian cycle (way)
          vc=sum(E1(evc(k),1:2))-vc; % new current vertex
          E1=E2;
          break;
        end
      end
    end
  end
end
return

author: HuangCongQing


text:
function C = blue()
  % BLUE Returns a nice bluecolor
  %
  % Outputs:
  %   [0.2 0.3 0.8]
  C = [0.2 0.3 0.8];
end

author: alecjacobson

repo: gptoolbox


text:
% Visualize all the inception_4c filters in bvlc_googlenet.
% Because of the size, I splitted the visualization map into 32*16 grids,
% with each grid contains visualization of 4 filters.
% This code is messay, I will refactor it in the future.
% Based on paper:
% Feng Wang, Haijun Liu, Jian Cheng,
% Visualizing Deep Neural Network by Alternately Image Blurring and Deblurring
caffe.reset_all();
caffe.set_mode_gpu();
gpu_id = 0;  % we will use the first gpu in this demo
caffe.set_device(gpu_id);
original_prototxt = 'D:\deepLearning\caffe-windows\models\bvlc_googlenet\googlenet.prototxt';
net_weights = 'D:\deepLearning\caffe-windows\models\bvlc_googlenet\thinned.caffemodel';
layer_name = 'inception_4c/output';
pattern_index=400;
channels = 512;
initial_size = [400 400];
[field_size, field_stride] = getReceptiveField(original_prototxt, net_weights, layer_name, pattern_index, initial_size);
field_size = ceil(field_size ./ field_stride) .* field_stride;
border = field_stride(1);
map_border = border ./ field_stride;
if max(field_size) > 224
%     error('field_size:(%d,%d),field_stride:(%d,%d), please increase the initial size.',field_size(1),field_size(2),field_stride(1),field_stride(2));
    field_size = [224 224];
    field_size = ceil(field_size ./ field_stride) .* field_stride;
    map_stride = ceil((field_size + border) ./ field_stride);
else
    map_stride = ceil((field_size + border) ./ field_stride);
end;

vert_num = 32;
hori_num = 16;
total_map = reshape(1:vert_num*hori_num,hori_num,vert_num)';
vert_split = vert_num/16;
hori_split = hori_num/8;
for vert_ind = 1:16
for hori_ind = 1:8
%     disp([layer_name '_' num2str(vert_ind) num2str(hori_ind) '.png']);
%     if exist([layer_name '_' num2str(vert_ind) num2str(hori_ind) '.png'],'file')
%         continue;
%     end;
% vert_ind = 2;
% hori_ind = 2;
caffe.reset_all();
current_map = total_map((vert_ind-1)*vert_split+1:vert_ind*vert_split,(hori_ind-1)*hori_split+1:hori_ind*hori_split)';

height = vert_split * (field_size(1)+border(1)) + border(1);
width = hori_split * (field_size(2)+border(1)) + border(1);

vgg_mean =  [103.939, 116.779, 123.68];
mean_image = permute(repmat(vgg_mean',[1,width,height]),[2,3,1]);
input_data = randn(width, height, 3, 1, 'single')*50;

for need_negative=1:1
original_net_model = fileread(original_prototxt);
original_net_model = strrep(original_net_model,'negative_slope:0#4c','negative_slope:1#4c');

visualize_prototxt = strrep(original_prototxt,'.prototxt','_visualize.prototxt');

fid = fopen(visualize_prototxt,'w');
proto_txt{1} = 'name: "Visualize"';
proto_txt{2} = 'input: "data"';
proto_txt{3} = 'input_dim: 1';
proto_txt{4} = 'input_dim: 3';
proto_txt{5} = ['input_dim: ' num2str(height)];
proto_txt{6} = ['input_dim: ' num2str(width)];
for i=1:6
    fprintf(fid,'%s\r\n',proto_txt{i});
end;

fprintf(fid,'%s\r\n',original_net_model);
fclose(fid);

if need_negative == 0
    input_data = randn(width, height, 3, 1, 'single')*50;
end;
all_grad = zeros(width, height, 3, 1, 'single');

visualize_net = caffe.Net(visualize_prototxt,net_weights,'test');
visualize_net.blobs(visualize_net.inputs{1}).set_data(input_data);
visualize_net.forward_to(layer_name);
target_blob = visualize_net.blob_vec(visualize_net.name2blob_index(layer_name));
output_data = target_blob.get_data();
backward_mask = zeros(size(output_data,1), size(output_data,2),'uint16');
backward_label = zeros(size(output_data,1), size(output_data,2),'uint16');
backward_data = zeros(size(output_data,1), size(output_data,2),channels, 'uint16');
label_mat = repmat(1:channels,size(output_data,1) * size(output_data,2),1)';
if map_stride(1) > 0
    backward_mask(ceil(map_stride(1)/2+map_border(1)):map_stride(1):end, ceil(map_stride(2)/2+map_border(2)):map_stride(2):end) = 1;
    backward_label(backward_mask==1) = current_map(:);
    for i=1:channels
        backward_data(:,:,i) = (backward_label == i);
    end;
else
    backward_mask(floor(size(backward_mask,1)/2),floor(size(backward_mask,2)/2)) = 1;
end;
item_num = sum(backward_mask(:));
target_blob = visualize_net.blob_vec(visualize_net.name2blob_index(layer_name));
data_blob = visualize_net.blob_vec(visualize_net.name2blob_index('data'));
for i=1:channels
    if i>384 && i<=448
    backward_data(:,:,i) =  backward_data(:,:,i) * 2;
    end;
end;

weight_decay = 0;
use_color_prior = false;
num_cluster=6;
color_prior = 0.5;
long_size = 512;
if need_negative==0
    tv_norm = 0;
    use_image_blur = false;
    use_image_deblur = false;
else
    tv_norm = 0;
    use_image_blur = true;
    use_image_deblur = true;
end;

if need_negative==0
    lr = 10;
else
    lr=200;
end;
max_lr = 50;
momentum = 0.8;
momentum2 = 0.99;
lastgrad = zeros(size(mean_image));
lastgrad2 = zeros(size(mean_image));
last_cost = -9999999999999;
blurred = false;
numLast = 0;

for iter = 1:1000
    bak_data = input_data;
    bak_grad = lastgrad;
    visualize_net.blobs(visualize_net.inputs{1}).set_data(input_data);
    visualize_net.forward_to(layer_name);
    output_data = target_blob.get_data();
    output_data(backward_data==0) = -999;
    output_map = max(output_data,[],3);
    min_cost = min(output_map(backward_mask==1));
    if min_cost>300 && need_negative==0
        disp(min_cost);
        break;
    end;
    cost = sum(output_map(backward_mask==1)) / item_num;
    if need_negative == 0
        for i=1:channels
            if output_map(backward_label == i) > 300
                backward_data(:,:,i) = 0;
                backward_mask(backward_label == i) = 0;
            end;
        end;
        fprintf('iter=%d,lr=%f,this_cost=%f,last_cost=%f,min_cost=%f,min_num=%d\n',iter,lr,cost, last_cost,min_cost,sum(sum(backward_mask==1)));
    else
        fprintf('iter=%d,lr=%f,this_cost=%f,last_cost=%f,image_norm=%f,min_cost=%f\n',iter,lr,cost, last_cost,norm(input_data(:)),min_cost);
    end;
%     output_data(:) = 0;
%     output_data(:,:,pattern_index) = backward_mask;
    target_blob.set_diff(backward_data);
    visualize_net.backward_from(layer_name);
    grad = data_blob.get_diff();
    all_grad = all_grad + grad;
    if iter==20
        input_data(all_grad==0) = 0;
    end;

    if tv_norm > 0
        I = input_data(:,:,:,1);
        Gx = (I(2:end-1,:,:) - I(1:end-2,:,:)) - (I(3:end,:,:) - I(2:end-1,:,:));
        Gx = [(I(1,:,:) - I(2,:,:)); Gx; (I(end,:,:) - I(end-1,:,:))];
        Gy = (I(:,2:end-1,:) - I(:,1:end-2,:)) - (I(:,3:end,:) - I(:,2:end-1,:));
        Gy = [(I(:,1,:) - I(:,2,:)) Gy (I(:,end,:) - I(:,end-1,:))];
        grad = grad - tv_norm * (Gx+Gy);
    end;

    if weight_decay > 0
        grad = grad - weight_decay * I;
    end;

    if use_color_prior
        gmm_prior = gaussian_net.forward({input_data});
        sum_gp = zeros(size(mean_image,1),size(mean_image,2));
        sum_prob_gradient = zeros(size(mean_image));
        for i=1:num_cluster
            gp = bsxfun(@minus,input_data(:,:,:,1),reshape(colorObj.mu(i,:),[1 1 3])) .* gmm_prior{1}(:,:,(i-1)*3+1:i*3);
            gp = sum(gp,3);
            gp = colorObj.PComponents(i) * exp(-gp);
            sum_prob_gradient = sum_prob_gradient + bsxfun(@times,gp,gmm_prior{1}(:,:,(i-1)*3+1:i*3));
            sum_gp = sum_gp + gp;
        end;
        sum_prob_gradient = bsxfun(@rdivide,sum_prob_gradient,sum_gp);
        sum_prob_gradient(isnan(sum_prob_gradient)) = 0;
        input_data(:,:,:,1) = input_data(:,:,:,1) - lr * color_prior * sum_prob_gradient;
    end;
    if need_negative==1
        %%%%%%%%%%%%%%%%%%%%%%%%gd linear search
        lastgrad = (1 - momentum) * lr * grad   + momentum * lastgrad;%/ norm(res(:))
        input_data(:,:,:,1) = input_data(:,:,:,1) + lastgrad;
    else
        %%%%%%%%%%%%%%%%%%%%%%%%adam
        lastgrad = (1 - momentum) * grad   + momentum * lastgrad;%/ norm(res(:))
        lastgrad2 = (1 - momentum2) * grad.^2 + momentum2 * lastgrad2;%/ norm(res(:))
        lg_correct = lastgrad ./ (1 - momentum^iter);
        lg2_correct = lastgrad2 ./ (1 - momentum2^iter);
        input_data(:,:,:,1) = input_data(:,:,:,1) + lr * lg_correct ./ (sqrt(lg2_correct) + 1e-8);
    end;
%     input_data(:,:,:,1) = input_data(:,:,:,1) / norm(input_data(:));

    k = mod(iter,10);
    if k==1
        H = fspecial('gaussian',[7 7],1.2);
        if use_image_blur
            input_data(:,:,:,1) = imfilter(input_data(:,:,:,1),H,'same');
            blurred = true;
        end
    end;
    if k==6
        H = fspecial('gaussian',[7 7],1);
        if use_image_deblur
            input_data(:,:,:,1) = deconvlucy(input_data(:,:,:,1), H);
            blurred = false;
        end;
    end;

    if (need_negative==1 && cost>last_cost)
        last_cost = cost;
        numLast = 0;
    else
        numLast = numLast + 1;
    end;
    if (need_negative==0 && min_cost>last_cost)
        last_cost = min_cost;
        numLast = 0;
    else
        numLast = numLast + 1;
    end;
    if numLast>100
        break;
    end;
    if lr<1e-5
        break;
    end;

    %%%%%%%%%%%%%%%%%%%%%%gd

    if mod(iter,100)==0
        output = mean_image + input_data(:,:,:,1);
        output = output(:, :, [3, 2, 1]);
        output = permute(output, [2 1 3]);
        figure(3);
        imshow(uint8(output));
%         title('generated image');
        I = output;
    end;
end;
end;
% imwrite(uint8(output),[layer_name '_' num2str(vert_ind) num2str(hori_ind) '.png']);
end;
end;

author: happynear

repo: DeepVisualization


text:
function n = ndims(t)
%NDIMS Return the number of dimensions of a tensor.
%
%   NDIMS(X) returns the number of dimensions of tensor X.
%
%   Examples
%   A = rand(4,3,1); ndims(A) %<-- Returns 2
%   X = tensor(A); ndims(X) %<-- Returns 2
%   X = tensor(A,[4 3 1]); ndims(X) %<-- Returns 3
%
%   See also TENSOR
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

n = numel(t.size);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function done = mission7(occupiers, target)

for i = 1:6
    if any(occupiers{i} == target)
        done = false;
        return
    end
end

done = true;


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
function amgausst
%AMGAUSST Unit test for the function AMGAUSS.

%	O. Lemoine - February 1996.

N=256; t0=149; T=50;
sig=amgauss(N,t0,T);
if abs(sig(t0)-1)>sqrt(eps),				% sig(t0)=1
 error('amgauss test 1 failed ');
end
[tm,T1]=loctime(sig);
if abs(T-T1)>sqrt(eps),					% width
 error('amgauss test 2 failed ');
end
dist=1:min([N-t0,t0-1]);
if any(abs(sig(t0-dist)-sig(t0+dist))>sqrt(eps))~=0, 	% symmetry
 error('amgauss test 3 failed ');
end

N=121; t0=60; T=27;
sig=amgauss(N,t0,T);
if abs(sig(t0)-1)>sqrt(eps),				% sig(t0)=1
 error('amgauss test 4 failed ');
end
[tm,T1]=loctime(sig);
if abs(T-T1)>sqrt(eps),					% width
 error('amgauss test 5 failed ');
end
dist=1:min([N-t0,t0-1]);
if any(abs(sig(t0-dist)-sig(t0+dist))>sqrt(eps))~=0, 	% symmetry
 error('amgauss test 6 failed ');
end

N=534; t0=333; T=70;
sig=amgauss(N,t0,T);
if abs(sig(t0)-1)>sqrt(eps),				% sig(t0)=1
 error('amgauss test 7 failed ');
end
[tm,T1]=loctime(sig);
if abs(T-T1)>sqrt(eps),					% width
 error('amgauss test 8 failed ');
end
dist=1:min([N-t0,t0-1]);
if any(abs(sig(t0-dist)-sig(t0+dist))>sqrt(eps))~=0, 	% symmetry
 error('amgauss test 9 failed ');
end

author: HeLiangHIT


text: %NOT Logical NOT. Dataset overload


text:
function x = r8sto_sl ( n, a, b )

%*****************************************************************************80
%
%% R8STO_SL solves a R8STO system.
%
%  Discussion:
%
%    The R8STO storage format is used for a symmetric Toeplitz matrix.
%    It stores the N elements of the first row.
%
%    For this routine, the matrix is also required to be positive definite.
%
%    This implementation of the algorithm assumes that the diagonal element
%    is 1.
%
%    The real symmetric Toeplitz matrix can be described by N numbers, which,
%    for convenience, we will label A(0:N-1).
%
%    Note that there is a typographical error in the presentation
%    of this algorithm in the reference, and another in the presentation
%    of a sample problem.  Both involve sign errors.  A minor error
%    makes the algorithm incorrect for the case N = 1.
%
%  Example:
%
%    To solve
%
%     1.0 0.5 0.2    x1    4.0
%     0.5 1.0 0.5 *  x2 = -1.0
%     0.2 0.5 1.0    x3    3.0
%
%    we input:
%
%      N = 3
%      A(0:N-1) = (/ 1.0, 0.5, 0.2 /)
%      B(1:3) = (/ 4.0, -1.0, 3.0 /)
%
%    with output:
%
%      X(1:3) = (/ 355, -376, 285 /) / 56
%             = (/ 6.339, -6.714, 5.089 /)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gene Golub and Charles Van Loan,
%    Section 4.7.3, "The General Right Hand Side Problem",
%    Matrix Computations,
%    Third Edition,
%    Johns Hopkins, 1996.
%
%  Parameters:
%
%    Input, integer N, the order of the system.
%
%    Input, real A(N), the R8STO matrix.
%
%    Input, real B(N), the right hand side of the linear system.
%
%    Output, real X(N), the solution of the linear system.
%
  k = 0;

  beta = 1.0E+00;
  x(k+1) = b(k+1) / beta;

  if ( k < n-1 )
    y(k+1) = -a(k+2) / beta;
  end

  for k = 1 : n-1

    beta = ( 1.0E+00 - y(k) * y(k) ) * beta;

    x(k+1) = ( b(k+1) - a(2:k+1) * x(k:-1:1)' ) / beta;

    x(1:k) = x(1:k) + x(k+1) * y(k:-1:1);

    if ( k < n - 1 )
      y(k+1) = ( -a(k+2) - a(2:k+1) * y(k:-1:1)' ) / beta;
      y(1:k) = y(1:k) + y(k+1) * y(k:-1:1);
    end

  end

  return
end


text:
function root_num = p05_root_num ( )

%*****************************************************************************80
%
%% P05_ROOT_NUM returns the number of known roots for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer ROOT_NUM, the number of known roots.
%
  root_num = 3;

  return
end


text:
function   [A, D, alpha] = LowRankDictionaryUpdate_m(A_prev, D_prev, alpha_prev, lambda)

% min \|E\|_1+ \lambda*(\|D\|_f^2+\|\alpha\|_f^2) s.t. D*\lambda+E=X

% this function is to update the bases D and coefficients \alpha, over
% the constraint D*\alpha = A_prev, by one iteration of ADM

% Xianbiao Shu (xshu2@illinois.edu)
% Updated on Aug-10-2011
% Copyright: Mitsubishi Electric Research Lab

    D = A_prev*alpha_prev'*inv(alpha_prev*alpha_prev'+lambda);
    alpha = inv(D'*D+ lambda)*D'*A_prev;

    A = D*alpha;

end

repo: lrslibrary


text:
function saverd3(fname,A,timewindow,timeinterval,antennasep)
% saverd3(fname,A,timewindow,timeinterval,antennasep)
%
% timewindow: [ns]
% timeinterval: seconds between traces
%
% Aslak Grinsted feb 2004

fname=strrep(lower(fname),'.rd3','');

fid=fopen([fname '.rd3'],'w');
A=(A-mean(mean(A)));
A=A*32767./max(max(abs(A)));
fwrite(fid,A,'short'); % should it be transposed?
fclose(fid);

fid=fopen([fname '.rad'],'w');

freq=size(A,1)*1000/timewindow;
freqsteps=1;
syscal=1/(freqsteps*freq);

stacks=1;
stacktime=size(A,1)*stacks/100000;

fprintf(fid,'SAMPLES:%i\r\n',size(A,1)); % SAMPLES:1024
fprintf(fid,'FREQUENCY:%0.6f\r\n',freq); % FREQUENCY:9410.717437
fprintf(fid,'FREQUENCY STEPS:%i\r\n',freqsteps);% FREQUENCY STEPS:12                    ??????????????
fprintf(fid,'SIGNAL POSITION:%0.6f\r\n',0); % SIGNAL POSITION:-0.156152          ????????????
fprintf(fid,'RAW SIGNAL POSITION:%i\r\n',1); % RAW SIGNAL POSITION:50401      ????????????
fprintf(fid,'DISTANCE FLAG:%i\r\n',0); % DISTANCE FLAG:0
fprintf(fid,'TIME FLAG:%i\r\n',1); % TIME FLAG:1
fprintf(fid,'PROGRAM FLAG:%i\r\n',0); % PROGRAM FLAG:0
fprintf(fid,'EXTERNAL FLAG:%i\r\n',0); % EXTERNAL FLAG:0
fprintf(fid,'TIME INTERVAL:%0.6f\r\n',timeinterval); % TIME INTERVAL:0.100000
fprintf(fid,'DISTANCE INTERVAL:%0.6f\r\n',0.1); % DISTANCE INTERVAL:0.100000
fprintf(fid,'OPERATOR:\r\n'); % OPERATOR:
fprintf(fid,'CUSTOMER:\r\n');% CUSTOMER:
fprintf(fid,'SITE:\r\n');% SITE:
fprintf(fid,'ANTENNAS:unknown\r\n');% ANTENNAS:800 MHz
fprintf(fid,'ANTENNA ORIENTATION:NOT VALID FIELD\r\n');% ANTENNA ORIENTATION:NOT VALID FIELD
fprintf(fid,'ANTENNA SEPARATION:%0.6f\r\n',antennasep);% ANTENNA SEPARATION:0.150000    ?????????????????
fprintf(fid,'COMMENT:Saved from matlab\r\n'); % COMMENT:
fprintf(fid,'TIMEWINDOW:%0.6f\r\n',timewindow); % TIMEWINDOW:108.812108   --nanoseconds for all samples
fprintf(fid,'STACKS:%i\r\n',stacks); % STACKS:1
fprintf(fid,'STACK EXPONENT:%i\r\n',0); % STACK EXPONENT:0
fprintf(fid,'STACKING TIME:%0.6f\r\n',stacktime); % STACKING TIME:0.010240                ??????????
fprintf(fid,'LAST TRACE:%i\r\n',size(A,2)); %LAST TRACE:18461

fprintf(fid,'STOP POSITION:%0.2f\r\n',size(A,2)*timeinterval); % STOP POSITION:53.16                 ??????????
fprintf(fid,'SYSTEM CALIBRATION:%0.10f\r\n',syscal);% SYSTEM CALIBRATION:0.0000088552   ?????? this has some effect on y-axis???
fprintf(fid,'START POSITION:%0.2f\r\n',0);% START POSITION:0.00
fprintf(fid,'SHORT FLAG:%i\r\n',1);% SHORT FLAG:1
fprintf(fid,'INTERMEDIATE FLAG:%i\r\n',0);% INTERMEDIATE FLAG:0
fprintf(fid,'LONG FLAG:%i\r\n',0);% LONG FLAG:0
fprintf(fid,'PREPROCESSING:%i\r\n',0);% PREPROCESSING:0
fprintf(fid,'HIGH:%i\r\n',5);% HIGH:5
fprintf(fid,'LOW:%i\r\n',15);% LOW:15
fprintf(fid,'FIXED INCREMENT:%0.3f\r\n',0.3);% FIXED INCREMENT:0.300
fprintf(fid,'FIXED MOVES UP:%i\r\n',0);% FIXED MOVES UP:0
fprintf(fid,'FIXED MOVES DOWN:%i\r\n',1);% FIXED MOVES DOWN:1
fprintf(fid,'FIXED POSITION:%0.3f\r\n',0); % FIXED POSITION:0.000
fclose(fid);


text:
%roscore=rosmatlab.roscore()

%node = rosmatlab.node('MyNode', '10.0.0.76')
% publisher = node.addPublisher('MyTopic', 'std_msgs/String');
% subscriber = node.addSubscriber('MyTopic', 'std_msgs/String', 1);
% msg = node.newMessage('std_msgs/String');
% msg.setData('hello world!');
% subscriber.setOnNewMessageListeners(@function3)
% publisher.publish(msg)

msg = rosmessage('std_msgs/String');
msg.Data = 'hello world!';

% pub = rospublisher('/MyTopic');
% pub.send(msg);

rospublisher('/MyTopic', msg);

author: petercorke

repo: robotics-toolbox-matlab


text:
function [F,TimeVector] = in_fread_fif(sFile, iEpoch, SamplesBounds, iChannels)
% IN_READ_FIF:  Read a block of recordings from a FIF file
%
% USAGE:  [F,TimeVector] = in_fread_fif(sFile, iEpoch, SamplesBounds, iChannels)
%         [F,TimeVector] = in_fread_fif(sFile, iEpoch, SamplesBounds)            : Read all the channels

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

if (nargin < 4)
    iChannels = [];
end

% Epoched data
% if isempty(SamplesBounds) || ~isfield(sFile.header, 'raw') || isempty(sFile.header.raw)
if ~isfield(sFile.header, 'raw') || isempty(sFile.header.raw)
    % Use data already read
    if isfield(sFile.header, 'epochData') && ~isempty(sFile.header.epochData)
        F = permute(sFile.header.epochData(iEpoch,:,:), [2,3,1]);
        TimeVector = linspace(sFile.epochs(iEpoch).times(1), sFile.epochs(iEpoch).times(2), size(F,2));
    % Read data from file
    else
        sfid = fopen(sFile.filename, 'r', sFile.byteorder);
        [F, TimeVector] = fif_read_evoked(sFile, sfid, iEpoch);
        fclose(sfid);
    end
    % Specific selection of channels
    if ~isempty(iChannels)
        F = F(iChannels, :);
    end
    % Specific time selection
    if ~isempty(SamplesBounds)
        iTime = SamplesBounds - round(sFile.epochs(iEpoch).times(1) .* sFile.prop.sfreq) + 1;
        F = F(:, iTime(1):iTime(2));
        TimeVector = TimeVector(iTime(1):iTime(2));
    end
    % Calibration matrix
    Calibration = diag([sFile.header.info.chs.cal]);
% Raw data
else
    % If time not specified, read the entire file
    if isempty(SamplesBounds)
        % Multiple files
        if isfield(sFile.header, 'fif_headers') && (length(sFile.header.fif_headers) > 1)
            SamplesBounds = [sFile.header.fif_headers{1}.raw.first_samp, sFile.header.fif_headers{end}.raw.last_samp];
        % Single file
        else
            SamplesBounds = [sFile.header.raw.first_samp, sFile.header.raw.last_samp];
        end
    end
    % If there are multiple FIF files: check in which one the data should be read
    if isfield(sFile.header, 'fif_list') && isfield(sFile.header, 'fif_times') && (length(sFile.header.fif_list) >= 2)
        % Check if all the samples are gathered from the same file
        fif_samples = round(sFile.header.fif_times .* sFile.prop.sfreq);
        iFile = find((SamplesBounds(1) >= fif_samples(:,1)) & (SamplesBounds(2) <= fif_samples(:,2)));
        % If this requires reading multiple files, call this function recursively
        if isempty(iFile)
            F = [];
            TimeVector = [];
            % Select the files to read from
            iFileStart = find((SamplesBounds(1) >= fif_samples(:,1)) & (SamplesBounds(1) <= fif_samples(:,2)));
            iFileStop  = find((SamplesBounds(2) >= fif_samples(:,1)) & (SamplesBounds(2) <= fif_samples(:,2)));
            for iFile = iFileStart:iFileStop
                % Create local structures for this file
                sFile_i = sFile;
                sFile_i.header = sFile.header.fif_headers{iFile};
                % Read the samples available in this file
                SamplesBounds_i = [max(min(SamplesBounds(1), fif_samples(iFile,2)), fif_samples(iFile,1)), ...
                                   max(min(SamplesBounds(2), fif_samples(iFile,2)), fif_samples(iFile,1))];
                [F_i,TimeVector_i] = in_fread_fif(sFile, iEpoch, SamplesBounds_i, iChannels);
                % Concatenate with previous files
                F = [F, F_i];
                TimeVector = [TimeVector, TimeVector_i];
            end
            return;
        else
            FifFile = sFile.header.fif_list{iFile(1)};
            sFile.header = sFile.header.fif_headers{iFile};
        end
    else
        FifFile = sFile.filename;
    end
    % Read block of data
    sfid = fopen(FifFile, 'r', sFile.byteorder);
    [F, TimeVector] = fif_read_raw_segment(sFile, sfid, SamplesBounds, iChannels);
    fclose(sfid);
    % Calibration matrix
    Calibration = diag([sFile.header.info.chs.range] .* [sFile.header.info.chs.cal]);
end

% Apply calibration
if ~isempty(iChannels)
    F = Calibration(iChannels,iChannels) * F;
else
    F = Calibration * F;
end


text:
function dtiMakeBrainMaskSafe(dtiDir,  dtiAlignedAvgFile)

%E.g., dtiMakeBrainSafeMask('dti06rt', 'dti_g712_b2000_aligned_trilin_avg.nii.gz')
%Take brainMask.nii.gz and remove the voxels that have a zero in the raw data (at
%least in one gradient directions)
%Input: directory with raw dti data
%Will overwrite your brainMask.nii.gz with a new safer one.

%ER 02/23/2009 pulled code out of mtrNotes to perform this operation
%ER 09/09 removed redundant code and changed input (dtiRawDir->dtiAlignedAvgFile)

dirPrev = pwd;
ni=niftiRead(dtiAlignedAvgFile);

% Update brain mask so that bad voxels are removed
cd([dtiDir filesep 'bin']);
bm = niftiRead('brainMask.nii.gz');
rawm = min(ni.data,[],4);
bm.data(rawm<1) = 0;

%bm.fname='brainMaskRawSafe.nii.gz';
writeFileNifti(bm);

cd(dirPrev);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function convertDicom2Vista_UU_donders(subjectName,SessionName,dicomServer)
% convertDicom2Vista_UU_donders(subjectName,SessionName,dicomServer)
%
% adapted from prepareUCSFParnassusUniversal
%
% e.g. prepareUCSFParnassusUniversal('norcia','LOCMT',2);
% The 1st argument is the subject name (corresponds to that subject's
% anatomy directory under /raid/MRI/anatomy/
% The 2nd argument is a plain text string to tag the exam. Avoid slashes,
% minus signs, spaces etc.
% The third argument indicates which dicom server the data came from. In
% theory we can auto-detec this but for now use
% 1: For NIC data via MPRAGE
% 2: For NIC data via Hurricane
% 3: For Donders' data
%
%
% This is a script to recon the dicom files that come from the UCSF MRI
% scanner at China Basin
% The first part uses the medcon program to convert those files into
% analyze format. After that, we use FSL to do some preprocessing like
% motion correction.
% The final step is to generate a mrLoadRet project from the set of 4d
% Analyze files.
% Motion correction is performed in several stages: We first do within-scan
% MC and store all the resulting transform matrices for each scan. Then we
% do between-scan MC using the mean image from each scan to align everything to the
% first scan in the session. We store the transforms for these alignments
% as well.
%
% Then we go back and combine the within and between scan transforms and
% re-apply the entire transform to the original data set. In this way we
% avoid re-sampling the original data more than we have to.
% EXAMPLE: prepareUCSF('appelbaum','LOC_Localizer');
% 2006.04.20 Mark Schira wrote it (mark@ski.org)
% 2007.06.13 Alex Wade modified it for use with data from Parnassus NIC
% Siemens scanner. Uses some ideas from dicom2analyze_v4.m : SPM and GazLab
% people
% Issue is that Siemens system uses mosaic image format: Each dicom file
% has an entire volume in it sometimes (and sometimes not in the case of
% anatomies..)
% The SPM5 routine dicom2analyze takes care of this.
% The critical step is...
%  % convert dicoms to analyzer
%        hdr = spm_dicom_headers(files);
%        spm_dicom_convert(hdr,opts,root_dir,format);
% Rather than using medcon (great though medcon is!)
% This 'universal' version is supposed to cope with the slightly different
% filename conventions in the two DICOM suervers 'hurricane' and 'MPRAGE'
% 080907: Can now use any directory names inside the RawDicom subdirectory

%%

% we need spm
if(isempty(which('spm_dicom_headers')))
    error('Requires spm2+ tools!');
end

if not(exist('subjectName','var'))
    error('need the lastname of subject as input');
end

if not(exist('SessionName','var'))
    SessionName='Whatever';
end

if not(exist('dicomServer','var'))
    error ('You must  specify a dicom server index: 1 is UCSF NIC / MPRAGE, 2 is UCSF NIC / Hurricane, 3 is Donders Institute');
end

ActionList={'Nothing','Inplanes','Functional','RefFunc','Anatomical','DTI'};

%%  Opening dialogs
% We present the user with a list of options. Not all stages have to be
% performed at once.

initOptions = {'clean previous files (be careful!)',...
    'DB Query',...
    'Do Detect data',...
    'Generate MLR structure',...
    'Generare BV structure',...
    'Motion cor. in scan',...
    'Motion cor. between scan',...
    'Perform blocked analysis',...
    'Perform event related analysis',...
    'Clean up, delete analyze files',...
    'Clean up. Zip the Dicom files',...
    };
initReply=[0 0 1 1 0 0 0 0 0 0 0]; % These are the defaults
initReply = buttondlg('mrInit UU-Data', initOptions,initReply);
if isempty(initReply), return; end

pause(0.2);
cleanfirst=initReply(1);
doDBQuery=initReply(2); %#ok<NASGU>
detectData=initReply(3);
doMLR=initReply(4);
doBV=initReply(5);%#ok<NASGU>
motionInScan=initReply(6);%#ok<NASGU>
motionBetween=initReply(7);%#ok<NASGU>
blockedData=initReply(8);%#ok<NASGU>
eventData=initReply(9);%#ok<NASGU>
doCleanUp=initReply(10);%#ok<NASGU>
doPackDicom=initReply(11);%#ok<NASGU>

baseDir = pwd;

%% we use some unix commands
operatingSystem=computer;
if not(strcmp(operatingSystem(1:4),'GLNX'))
    error('This function  is for Linux machines only');
end

%% remove analyze file from previous runs

if cleanfirst
    ! rm ./*/*/*/*.hdr
    ! rm ./*/*/*/*.img
    ! rm -r ./Inplane
    ! rm -r ./Gray
    ! rm -r ./Volume
    ! rm mrSESSION*.mat
    ! rm prepare*vars.mat
end

%% Automatic detection of the Directories, Inplanes etc
if detectData
    disp('Reading the directory structure...');
    RawDir = './Raw/Functionals/';
    baseDirListing=dir(RawDir);

    % Sort this so that it is in ascending numeric order
    nDirs=length(baseDirListing);
    tempDirName=zeros(nDirs,1);

    for t=1:length(baseDirListing);
        switch dicomServer
            case {1,3}		% MPRAGE
                tempDirNum = str2double( baseDirListing(t).name(find(baseDirListing(t).name=='_',1,'last')+1:end) );

            otherwise
                tempDirNum=str2double(baseDirListing(t).name);
        end

        if (tempDirNum>0) && isfinite(tempDirNum)
            tempDirName(t)=tempDirNum;
        end

    end

    [tempDirName,i]=sort(tempDirName); % Because we don't know how DIR might order the output
    baseDirListing=baseDirListing(i); % Apply the sort

    baseDirListing=baseDirListing(3:end); % Strip off stuff that didn't start with a number - like . and ..
    disp({baseDirListing.name}')

    foundScans(1).DirName='';
    foundScans(1).Files=0;
    foundScans(1).Slices=0;
    foundScans(1).Volumes=0;
    foundScans(1).Sequence='';
    foundScans(1).Action='';
    foundScans(1).Filenames='';
    foundScans(1).SkipVols='';
    foundScans(1).Cycles='';
    foundScans(1).Filename4d='';
    foundScans(1).lastfile='';

    inplanedir=0;

    for thisEntry=1:length(baseDirListing)
        fprintf('\nProcessing %d',thisEntry);

        if baseDirListing(thisEntry).isdir

            thisdir=baseDirListing(thisEntry).name;

            switch dicomServer
                case 1 % Coming from MPRAGE at the NIC
                    tempdir=dir([baseDir,filesep,RawDir,filesep,baseDirListing(thisEntry).name,filesep,'*.dcm']); % Only list the Dicom files.
                case 2 % Coming from Hurricane
                    tempdir=dir([baseDir,filesep,RawDir,filesep,baseDirListing(thisEntry).name,filesep,'*.DCM']); % Only list the Dicom files.
                case 3 % Coming from Donders
                     tempdir=dir([baseDir,filesep,RawDir,filesep,baseDirListing(thisEntry).name,filesep,'*.IMA']); % Only list the Dicom files.
                otherwise
                    error('No valid dicomServer type specified');
            end

            foundScans(thisEntry).RawDir = RawDir;
            foundScans(thisEntry).DicomServer = dicomServer;
            foundScans(thisEntry).Filenames=tempdir(1).name;
            foundScans(thisEntry).Filenames=...
            foundScans(thisEntry).Filenames(1:find(foundScans(thisEntry).Filenames=='I'));
            foundScans(thisEntry).DirName=baseDirListing(thisEntry).name;
            foundScans(thisEntry).Files=length(tempdir);
            header=dicominfo([baseDir,filesep,RawDir,filesep,baseDirListing(thisEntry).name,filesep,tempdir(1).name]);
            if (isfield(header,'Private_0019_100a')) % This is a functional image stored in Mosaic format
                disp(tempdir);
                foundScans(thisEntry).Slices=header.Private_0019_100a;%Private_0021_104f;
                foundScans(thisEntry).Volumes=(length(tempdir));%/foundScans(thisdir).Slices;
            else % Otherwise it's a T1 image stored slice-wise
                foundScans(thisEntry).Slices=(length(tempdir));
                foundScans(thisEntry).Volumes=1;
            end

%             if (sum(findstr(header.SeriesDescription,'ep2d'))) && (foundScans(thisEntry).Volumes>40) % Currently, this is hard-coded. We should be more flexible in case we change the protocol name.
            if (~isempty(strfind(header.SeriesDescription,'ep2d')) || ~isempty(strfind(header.SeriesDescription,'bas_MoCo'))) && (foundScans(thisEntry).Volumes>40) % Currently, this is hard-coded. We should be more flexible in case we change the protocol name.
                foundScans(thisEntry).Action='Functional';
            else
                foundScans(thisEntry).Action='Nothing';
            end
            if strcomp(header.SeriesDescription,'inplanes')
                inplanedir=thisdir;
            end

            foundScans(thisEntry).Sequence=header.SeriesDescription;

        end
    end
    if not(inplanedir==0)
        foundScans(inplanedir).Action='Inplanes';
    else
        disp(sprintf('[%s]:No Inplanes found',mfilename));
    end
    foundScans=selectInitParameters(foundScans);

else %if not scan the directory one would have to read the info from a previous run

    foundScans=selectInitParameters(foundScans);
end

%% The above does seem to work quite well. Now we can pass the dir names
% into the SPM5 dicom converters to generate NIFTI or ANALYZE format
% data...
%
% Recon all the functionals to Analyse

disp(sprintf('[%s]:Reconning all scans...',mfilename));drawnow;
procScans = cell(numel(foundScans),1);
hdrList   = cell(numel(foundScans),1);
for thisDir=1:numel(foundScans)
    if strcmp(foundScans(thisDir).Action,'Functional') ||  strcmp(foundScans(thisDir).Action,'RefFunc')
        tic
        [procScans{thisDir},hdrList{thisDir}]=Dicom2Analyze_SPM5(foundScans(thisDir));
        toc
    end
end
%
%
% for thisDir=1:length(foundScans)
%     if strcmp(foundScans(thisDir).Action,'Inplanes') ||  strcmp(foundScans(thisDir).Action,'Anatomical')
%         outFileName=[foundScans(thisDir).Filenames,'_a4d.hdr']; % We will generate a 3d analyze (Nifti) file from anatomicals
%
%         tic
%         [procScans{thisDir},hdrlist{thisDir}]=Dicom2ThreeDAnalyze_SPM5(foundScans(thisDir),pwd,dicomServer);
%         toc
%
%     end
% end
% datestr(now)
% clear thisDir:

%%

cd (baseDir)
if doMLR
    disp('Making directories');
    !mkdir Inplane
    !mkdir Volume
    !mkdir Gray
end

save prepareUUvars

status=mrInitRetFromSPM(foundScans,subjectName,SessionName);

disp('Result of mrInitRetFromSPM');
disp(status);

if blockedData
    mrvista('inplane');
    inplane=getSelectedInplane;
    inplane=computeMeanMap(inplane,0);
    inplane=computeCorAnal(inplane,0);
    inplane = selectDataType(inplane,2);
    inplane=computeMeanMap(inplane,0);
    inplane=computeCorAnal(inplane,0);
    %    closeInplaneWindow();
    datestr(now)
end