text:
classdef TraceData
   %TraceData   Handles the data associated with timeseries.
   %
   %   TraceData contains the functionality associated with data
   %   manipulation, including dealing with units.
   %
   %   TraceData might be considered a "light" version of MATLAB's timeseries
   %   class. Whereas the timeseries class has lots of functionality, it
   %   suffers (at the time of this writing) from much slower execution times.
   %
   % About TraceData vs waveform
   %
   % TraceData Properties:
   %       data       - column of evenly sampled data
   %       samplerate - number of samples per second
   %       units      - units for the data, ex. 'counts', 'm / sec'
   %       duration   - time between first and last samples
   %       hasnan     - true if any value is nan
   %
   % TraceData Methods:
   %
   %    Derived properties:
   %       nyquist - Nyquist frequency (samplerate / 2)
   %       period - Period  (1/samplerate)
   %       formattedduration - retrieve duration as formatted text
   %
   %    additional Methods:
   %       nsamples  - return number of data samples
   %
   %    Mathamatical operations:
   %       plus - (+)addition
   %       minus - (-)subtraction
   %       times - (.*) element multiplication
   %       mtimes - (*) matrix multiplication
   %       rdivide - (./) element division
   %
   %       uminus - (-A) unary minus
   %       sign - signum of data (returns array of +1, 0, or -1)
   %       abs - Absolute value
   %
   %       forEach - apply a function element-wise Tnew(n) = f(T(n), values(n))
   %
   %    Binary Operators:
   %       eq - (A==B) true if data, samplerate and units match.
   %       ne - (A~=B) false if data, samplerate and units match.
   %
   %    Statistics (returns a single value for each trace):
   %       min - Minimum value of data
   %       max - Maximum  value of data
   %       mean - Average value of data
   %       median - Median value of data
   %       var - Variance of data
   %       std - Standard Deviation of data
   %
   %    Advance Mathamatical operations:
   %       diff - Difference and approximate derivative for traces
   %       integrate - cumulative sums of data
   %
   %       hilbert - Hilbert envelope (real only)
   %       taper - Apply a taper to the data
   %       resample - Resample the data
   %       demean - Remove the mean
   %       detrend - Remove trend from data
   %       clip - Clip the data
   %
   %       extract - Retrieve a subset of the data
   %
   %    Conversion operations:
   %       double - retrieve data as 1xN array of double
   %
   %       stack - stack all waveforms
   %       binStack - stack N waveforms at a time (with optional overlap)
   %
   %       compatibleWith - Compare units, samplerate and datalength
   %       assertCompatibleWith - Error if units, samplerate and datalength do not match
   %
   %       fillgaps - replace nan values (can replace other values)
   %       zero2nan - Replace values close to zero with nans
   %       amplitude_spectrum
   %       setlength - adjust length of data to allow batch processing
   %
   % See also Seismictrace, waveform, timeseries

   % Error recovery and trapping with traces
   % trying new tactic. I won't try to anticipate all the various ways
   % someone can provide incompatible data. Instead, I'm going to provide a
   % comment with my expectations which will show up automatically in the
   % displayed error.

   % Waveform functions included in Tracedata:
   %  - functions that manipulate the data itself
   %
   % Waveform functions NOT added to Tracedata:
   %  - functions that require knowledge of time
   %  - functions that require knowledge of location
   %  - function that access user defined fields
   %  - functions that access history
   %  - additinal not-so-sure-we-need-these functions
   %    : getpeaks

   properties
      data = []% time-series data, kept in a column
      samplerate  = NaN % in samples/sec
      units = 'none' % text description of data units
   end

   properties(Dependent)
      duration % duration in seconds. To get matlab duration, divide by 86400
      hasnan;
   end

   properties(Hidden=true)
      % when trust_assignments is false, then data
      trust_assignments = false; % if not trusted, then slows down computations somewhat, but is useful for debugging.
      debug_level
   end

   methods
      function obj = TraceData(varargin)
         %TraceData construct a TraceData object
         switch nargin
            case 1
               if isa(varargin{1}, 'waveform')
                  obj.samplerate = get(varargin{1},'freq');
                  obj.data = get(varargin{1},'data');
                  obj.units = get(varargin{1}, 'units');

               elseif isa(varargin{1},'TraceData') && ~strcmp(class(varargin{1}),'TraceData')
                  if numel(varargin{1})==1;
                  obj.data = varargin{1}.data;
                  obj.samplerate = varargin{1}.samplerate;
                  obj.units = varargin{1}.units;
                  else
                     for n=1:numel(varargin{1})
                        obj(n) = TraceData(varargin{1}(n));
                     end
                     reshape(obj,size(varargin{1}));
                  end
               else
                  error('Unknown conversion');
               end
            case 3 % TraceData(data, samplerate, units);
               obj.data = varargin{1};
               obj.samplerate = varargin{2};
               obj.units = varargin{3};
         end %switch
      end

      function tf = get.hasnan(obj)
         %get.hasnan   functional hasnan allows for
         % NOTE: technically, hasnan could be set as data is entered.
         %   any([T.hasnan])
         tf = any(isnan(obj.data));
      end
      function obj = set.data(obj, values)
         %set.data   Assign values to data as a column
         obj.data = values(:);
      end

      function N = nyquist(T)
         %nyquist   the nyquist frequency calculated as (samplerate / 2)
         N=[T.samplerate] ./ 2;
         reshape(N,size(T));
      end

      function p = period(T)
         %period   Period, calculated as (1 / samplerate)
         p = 1 ./ [T.samplerate];
         % NOT reshaped
      end

      function val = sampletimes(obj)
         %sampletimes  MATLAB time offset of each sample.
         %
         % equivalent to waveform's get(w,'timevector');
         assert(numel(obj) == 1, 'only works on one TraceData at a time');
         val = (0:(numel(obj.data)-1)) .* obj.period / 86400;
         val = val(:);
      end

      function secondsOfData = get.duration(obj)
         %secondsOfData   Duration of samples, in seconds
         if isempty(obj.data) || isempty(obj.samplerate)
            secondsOfData = 0;
         else
            secondsOfData = numel(obj.data) / obj.samplerate;
         end
      end

      function n = nsamples(obj)
         %nsamples   number of samples in the trace(s)
         %   n = nsamples(traces) will return the number of samples for
         %   each trace in an array of the same size as traces
         %
         %   See also numel
         n = zeros(size(obj));
         for m=1:numel(obj)
            n(m) = numel(obj(m).data);
         end
      end

      function s = formattedduration(obj, fmt)
         %formattedduration   Duration as a formatted string
         %  s = trace.formattedduration() retrieves the duration in the
         %  default format as 'dd:hh:mm:ss.SSS'
         %  s = trace.formattedduration(fmt) retrieves duration in format
         %  specified.
         %
         %  depends upon the duration class, introduced in r2014b
         %
         %  See also duration
         secsOfData = obj.duration;
         if exist('duration','class') %available in recent r2014b and later of matlab
            if exist('fmt','var') && ~isempty(fmt)
               s = char(duration(0,0,secsOfData,'Format',fmt));
            else
               s = char(duration(0,0,secsOfData,'Format','dd:hh:mm:ss.SSS'));
            end
         else
            error('advanced duration functions supported in r2014b and later');
         end
      end

      %% Mathamatical - BASIC OPERATIONS
      function A = plus(A, B)
         %+   Plus.
         %   C=A+B add something to the TraceData's data.
         %
         %   valid combinations
         %      TraceData + NumericVector; % same length as TraceData.data
         %      TraceData + Scalar;
         %
         %  To avoid ambiguity with the metadata, two traces are added
         %  together by explicitly adding the data from one to the other:
         %
         %      % assertCompatiblewith(TraceData1, TraceData2); % for debug
         %      TraceData1 + TraceData2.data
         %
         % See also minus, arrayApply, compatiblewith, assertCompatiblewith

         if ~isa(A, 'TraceData')
            [A, B] = deal(B, A); % swap values
         end

         % A is guaranteed to be TraceData
         if isnumeric(B)
            for n = 1: numel(A)
               A(n).data = A(n).data + B; % add to either a scalar or a COLUMN of numbers (same length as TraceData's data)
            end
         elseif isa(B,'TraceData')
            error('TraceData:plus:ambiguousOperation',...
               ['Adding two TraceData objects results in ambiguous metadata.'...
               '\nInstead, add the data explicitly to the Trace whose'...
               ' metadata you wish to keep.\nEx. T = T1 + T2.data']);
         else
            error('TraceData:plus:unknownClass','do not know how to add a %s to a TraceData object', class(B));
         end
      end

      function A = minus(A, B)
         %-  Subtract something from the Trace's data,
         %   C=A-B This will return a TraceData object
         %
         %   valid combinations
         %      TraceData - NumericVector; % same length as TraceData.data
         %      TraceData - Scalar;
         %
         %   To avoid ambiguity with the metadata, two traces are subtracted
         %   by explicitly subtracting the one's data from the other:
         %
         %      % assertCompatiblewith(TraceData1, TraceData2); % for debug
         %      TraceData1 - TraceData2.data
         %
         %   See also plus, arrayApply, compatiblewith, assertCompatiblewith

         if isnumeric(B)
            % A is guaranteed to be a TraceData
            for n = 1:numel(A)
               A(n).data = A(n).data - B; % subtract either a scalar or a COLUMN of numbers (same length as TraceData's data)
            end
         elseif isa(B,'TraceData')
            error('TraceData:minus:ambiguousOperation',...
               ['Subtracting a Trace from a constant is not supported.\n'...
               'For equivalent functionality, add the negative.\nEx. ans = 5 + (-T)']);
         else
            error('TraceData:minus:unknownClass','do not know how to subtract a %s from a %s', class(B), class(A));
         end
      end

      function A = times(A,B)
         %.*   Elementwise Trace data multiplication
         %   C=A.*B or C=B.*A when A is a TraceData object
         %   Either A or B can be a scalar, or a vector of same size as the
         %   data elements of the TraceData object.
         %
         %   To avoid ambiguity with the metadata, two traces are multiplied
         %   by explicitly multiplying one's data with the other:
         %
         %      % assertCompatiblewith(TraceData1, TraceData2); % for debug
         %      TraceData1 .* TraceData2.data
         %
         %  See also arrayApply, mtimes, compatiblewith, assertCompatiblewith
         if isnumeric(B)
            for n=1:numel(A)
               A(n).data = A(n).data .* B; % B should be either scalar or same size as obj.data
            end
         elseif isnumeric(A)
            [A, B] = deal(B, A); % swap values
            for n=1:numel(A)
               A(n).data = A(n).data .* B; % B should be either scalar or same size as obj.data
            end
         else
            error('TraceData:times:unknownClass','do not know how to multiply a %s with a TraceData object', class(B));
         end
      end

      function C = mtimes(A, B)
         %*   Matrix multiplication against data within a trace
         %   C=A*B matrix multiplication against data within a trace
         %   result is a matrix, vector, or scalar. (NOT a TraceDataObject)
         %
         %   See also times, arrayApply
         if isa(A,'TraceData')
            C = A.data * B;
         else
            C = A * B.data;
         end
      end

      function A = rdivide(A, B)
         %./  Divide data elements of A by B
         %   C=A./B performs elementwise division on trace data
         %   A must be a TraceData object
         %   B can be either a scalar or vector of numbers the same size as A.data
         %
         %   To avoid ambiguity with the metadata, two traces are multiplied
         %   by explicitly dividing one's data from the other:
         %
         %   % assertCompatiblewith(TraceData1, TraceData2); % for debug
         %   TraceData1 ./ TraceData2.data
         %
         %   See also times, mtimes, arrayApply, compatiblewith, assertCompatiblewith

         if isnumeric(B)
            % A is guaranteed to be a TraceData
            for n=1:numel(A)
               A(n).data = A(n).data ./ B; % subtract either a scalar or a COLUMN of numbers (same length as TraceData's data)
            end
         else
            error('TraceData:rdivide:unknownClass','do not know how to divide a %s from a TraceData object', class(B));
         end
      end

      function A = power(A, B)
         %.^   Array power for TraceData
         %   C=A.^B raises each data element of A to the power B and then
         %   returns the resulting trace(s)
         %
         %   See also power, arrayApply

         assert(isa(A,'TraceData'),'TraceData:power:invalidType',...
            'for A .^ B, B cannot be a TraceData object');
         assert(isnumeric(B),'TraceData:power:invalidType',...
            'for A .^ B, B must be numeric');
         for n=1:numel(A)
            A(n).data = A(n).data .^ B; % B should be scalar or same length as A
         end
      end

      function A = uminus(A)
         %-   Unary minus.
         %   -A negates the data in A and returns the resulting trace(s)
         %
         %   see also uminus
         for n=1:numel(A)
            A(n).data = -A(n).data;
         end
      end

      function trace = abs(trace)
         %Absolute value of trace data
         %  T = abs(trace) returns traces containing the absolute values of
         %  the data.
         %
         %See also abs, sign
         for n=1:numel(trace)
            trace(n).data = abs(trace(n).data);
         end
      end

      function trace = sign(trace)
         %sign   Convert each data point to its signum (+1, 0, +1)
         % T = sign(trace) returns a a trace containing the signs of the
         % data, instead of the data.
         %
         % See also sign
         for n=1:numel(trace)
            trace(n).data = sign(trace(n).data);
            %trace(n).units = ['sign(', trace(n).units, ')'];
         end
      end

      %% more complicated
      function T = diff(T, varargin)
         %diff   Difference and approximate derivative for traces
         %  A = diff(trace)
         %  A = diff(trace, options) see the builtin diff for details on
         %  available options.
         %
         %  trace must have samplerate and data assigned, otherwise you may
         %  get "a Matrix dimensions must agree" error
         %
         %  units are automatically changed. Assuming the sample rate is
         %  samples/sec, then the new trace is in currentunits / sec.
         %
         %  See also diff
         if isempty(varargin)
            for I=1:numel(T)
               T(I).data = diff(T(I).data) .* T(I).samplerate; % must have data and sample rate
            end
         else
            error('not implemented yet');
         end
         for I=1:numel(T)
            tempUnits = T(I).units;
            whereInUnits = strfind(tempUnits,' * sec');
            if isempty(whereInUnits)
               T(I).units = [tempUnits, ' / sec'];
            else
               tempUnits(whereInUnits(1) :whereInUnits(1)+5) = [];
               T(I).units = tempUnits;
            end
         end
      end

      function T = integrate (T,method)
         %integrate   Integrate a tracedata signal
         %   trace = trace.integrate([method])
         %   goes from Acceleration -> Velocity, and from Velocity -> displacement
         %
         %   tr = trace.integrate()
         %   tr = integrate(trace,'cumsum') performs integration by summing the
         %   data points with the cumsum function, taking into account time interval
         %   and updating the units as appropriate.
         %
         %   tr = trace.integrate('trapz') as above, but uses matlab's
         %   cumtrapz function to perform the integration.
         %
         %   Input Arguments
         %       trace: tracedata   N-DIMENSIONAL
         %       METHOD: either 'cumtrapz' or 'cumsum'  [default is cumsum]
         %
         %   Actual implementation  merely does a cumulative sum of the trace's
         %   samples, and updates the units accordingly. These units may be a
         %   little kludgey.
         %
         %   See also cumsum, cumtrapz, diff

         Nmax = numel(T);
         allfreq = [T.samplerate];

         if ~exist('method','var')
            method = 'cumsum';
         end

         switch lower(method)
            case 'cumsum'
               integratefn = str2func('cumsum');
            case 'trapz'
               integratefn = str2func('cumtrapz');
            otherwise
               error('TraceData:integrate:unknownMethod',...
                  'Unknown integration method. Valid methods are ''cumsum'' and ''trap''');
         end

         for I = 1 : Nmax
            T(I).data = integratefn(T(I).data) ./ allfreq(I);
            tempUnits = T(I).units;
            whereInUnits = strfind(tempUnits,' / sec');
            if isempty(whereInUnits)
               T(I).units = [tempUnits, ' * sec'];
            else
               tempUnits(whereInUnits(1) :whereInUnits(1)+5) = [];
               T(I).units = tempUnits;
            end
         end
      end

      function A = demean(A)
         %demean   Remove the average from TraceData
         %   A = demean(A) will subtract the mean of a trace from its data.
         %
         %   See also demean
         for n=1:numel(A);
            A(n) = A(n) - mean(A(n).data(~isnan(A(n).data)));
         end
      end

      function A = detrend(A, varargin)
         %detrend   Remove the trend from TraceData
         %   A = detrend(A) will subtract the trend of a trace from its
         %   data.
         %   A = detrand(A, ...) allows for additional options.  For an
         %   explanation of these, see MATLAB's builtin detrend.
         %
         %   See also detrend
         for n=1:numel(A);
            A(n).data = detrend(A(n).data,varargin{:});
         end
      end

      function val = bulkCalculate(T, F)
         %bulkCalculate   Run a function against the data
         %  vals = T.bulkCalculate(funcHandle) will against the data field
         %  of each element in T shape is preserved, and all empty
         %  traces return nan FH must return a single value
         %
         %  basically, this will run:
         %    funcHandle(T.data)
         %
         %  for each trace in T
         %
         %  See also max, min, mean, median

         val = nan(size(T));
         for n=numel(T) : -1 : 1
            hasdata(n) = ~isempty(T(n).data);
         end
         val(hasdata) = arrayfun(F,T(hasdata));
      end

      function T = forEach(T, f, values, varargin)
         %forEach  applies a function element-wise to an array B(n) = f(T(n), values(n))
         %   VALUES must be a numeric array of the same size (shape) as T.
         %   Either T or VALUES could also be individual values.
         %
         %   When numel(T)==1:
         %   the output will be an array of T the same size as VALUES.
         %
         %   When numel(V)==1,
         %   you might as well use the standard operator instead.
         %       eg.  T = T + VALUE;
         %
         %  the function f can either be a function handle, string name of
         %  a function  (such as 'plus'), or the valid representation of a
         %  simple mathamatical function, such as: '.*'  './'  '+'  '-'  '.^'
         %
         %  This means that
         %       Tsq = forEach(T, @power, valueArray);
         %  is equivelent to
         %       Tsq = T.forEach('power', valueArray);
         %  which is equivelent to
         %       Tsq = T.forEach('.^', valueArray);
         %
         %
         %  apply(...,'loose') ignores the size restriction, but
         %  enforces the number of elements.  Most useful in getting around
         %  applying something 1xN to something Nx1
         %
         %  example:
         %    % add offset to each trace
         %    offsets = 1:numel(T);
         %    Toffset = T.apply('+', offsets, 'loose');
         %
         %  example:
         %    % you can also create your own function of the form f(T,x):
         %    myfn = @(T, x) T .* x + x;
         %    newT = T.apply(T, myfn, values)
         %
         % See also times, rdivide, plus, minus, power
         TisScalar = numel(T)==1;
         VisScalar = numel(values)==1;
         eitherIsScalar = TisScalar || VisScalar;
         ignoreShape = eitherIsScalar || (~isempty(varargin) && strcmpi(varargin{end},'loose'));
         if ignoreShape
            assert(eitherIsScalar || numel(values)==numel(T), ...
               'Both the traces and values must have the same number of elements T:[%d] vs values:[%d',...
               numel(T), numel(values) );
         else
            assert(eitherIsScalar || all(size(values)==size(T)),...
               ['Both the traces and the values must be the same size.'...
               'T:[%s] vs values:[%s]\n',...
               'Use arrayApply(...,''loose'') if only the number of elements matters'],...
               TraceData.size2str(T), TraceData.size2str(values));
         end
         if ischar(f)
            f=str2func(f);
         end

         if VisScalar
            T = f(T, values); % not necessary to use arrayApply when either is a scalar
            return
         elseif TisScalar
            T = repmat(T,size(values));
            for n=1:numel(values)
               T(n) = f(T(n), values(n));
            end
            return
         end

         for n=1:numel(T)
            T(n) = f(T(n), values(n));
         end
      end

      function val = max(T)
         %max   Maximum value for trace data
         %   maxVals = max(traces);
         %
         %   See also max
         F = @(X) max(X.data);
         val = T.bulkCalculate(F);
      end
      function val = min(T)
         %min   Minimum value for trace data
         %   minVals = min(traces);
         %
         %   See also min
         F = @(X) min(X.data);
         val = T.bulkCalculate(F);
      end
      function val = mean(T)
         %mean   Average or mean value for trace data
         %   meanVals = mean(traces);
         %
         %   See also mean
         F = @(X) mean(X.data(~isnan(X.data)));
         val = T.bulkCalculate(F);
      end
      function val = median(T)
         %median   Median value for trace data
         %   medianVals = median(traces);
         %
         %   See also median
         F = @(X) median(X.data(~isnan(X.data)));
         val = T.bulkCalculate(F);
      end
      function val = std(T, varargin)
         %std   standard deviation for traces
         %  vals = std(traces) retrieves the std from each trace, returning
         %  a traces array of the same shape as T.
         %  vals = std(traces, options); lets you declare options as per the
         %  builtin version of std
         %  useful option:
         %
         %     vals = T.std('omitnan') % or 'includenan'
         %
         %  See also std, bulkCalculate
         if exist('varargin','var')
            F = @(X) std(X.data,varargin{:});
         else
            F = @(X) std(X.data);
         end
         val = T.bulkCalculate(F);
      end
      function val = var(T, varargin)
         %var   calulate variance for traces
         %  vals = var(traces);
         %  vals = var(traces,options); lets you declare options as per the
         %  builtin version of var.
         %  useful option:
         %    vals = traces.var(taces,'omitnan') % or 'includenan'
         %
         %  See also var, bulkCalculate
         if exist('varargin','var')
            F = @(X) var(X.data,varargin{:});
         else
            F = @(X) var(X.data);
         end
         val = T.bulkCalculate(F);
      end
      %% extended functionality
      [A, phi, f] = amplitude_spectrum(td)

      function obj = clip(obj, vals)
         %clip   clips a trace's data at a particular max/min value range
         %   clippedtraces = clip(traces, values)
         %
         %   Input Arguments
         %       TRACES: an N-dimensional TraceData object
         %       VALUES: a number. If a scalar, then amplitudes will be
         %           clipped at +/- value. If a pair, eg. [Max, Min] then
         %           traces will be clipped between these two values
         %
         %   All values beyond maximum ranges will be set to the maximum range

         if nargin < 2,
            vals = [];
         end

         switch numel(vals)
            case 0
               warning('TraceData:clip:noClipValue',...
                  'No values given, TraceData remains unclipped');
               return

            case 1 % clip at +/- input value
               if ~isnumeric(vals)
                  error('TraceData:clip:invalidClipValue','Non-numeric clip value');
               end
               vals = abs(vals);
               disp(['Clipping at +/- ' num2str(vals)]);
               vals = [-vals vals];

            case 2
               if ~isnumeric(vals)
                  error('TraceData:clip:invalidClipValue','Non-numeric clip value');
               end
               if vals(1) > vals(2)
                  error('TraceData:clip:reversedValues',...
                     'Clipping values are  reversed: Value1 %f > Value2 %f', vals(1),vals(2));
               end
               %everything is AOK
            otherwise
               warning('TraceData:clip:invalidClipValue',...
                  'Invalid clip values. TraceData remains unclipped. specify [min max] values');
               return
         end

         for N = 1: numel(obj)
            d = obj(N).data;
            d(d < vals(1)) = vals(1);
            d(d > vals(2)) = vals(2);
            obj(N).data = d;
         end
      end

      function d = double(obj,option)
         %double   TraceData data as double
         %   d = double(T) returns the data within T as a column. If T is
         %   an array with N elements, then data is returned in array with
         %   size N x maximum_data_length_in_T.
         %
         %   d = double(T, fillfunction) allows you to specify how the array
         %   is filled. FILLFUNCTION is a handle or name of any function
         %   that can accept (row, col) arguments, and return an array of some
         %   value. By default, the array is filled with zeros.
         %
         %   Example. Retrieve array using nan as filler
         %   d = double(T, @nan);
         %
         %   See also zeros, nan

         if nargin == 1
            createDefaultArray = @zeros;
         else
            if ischar(option)
               createDefaultArray = str2func(option);
            elseif isa(option,'function_handle')
               createDefaultArray = option;
            else
               error('TraceData:double:invalidOption',...
                  '''option'' must be either a function handle or function name.');
            end
         end
         dlens = obj.nsamples();
         maxlen = max(dlens(:));
         d = createDefaultArray(maxlen, numel(obj));
         for n = 1:numel(obj)
            d(1:dlens(n),n) = double(obj(n).data);
         end
      end

      function tf = eq(A, B)
         %== Equal for TraceData
         %  eq(A,B) is called for A==B when A or B is a TraceData object.
         %  will return equal if the data , samplerate, and units match.
         %
         %  See also ne
         tf = all(A.data == B.data) && A.samplerate==B.samplerate && strcmp(A.units, B.units);
      end

      function tf = ne(A,B)
         %~=   Not equal for TraceData.
         %  A~=B compares the data, samplerate, and units
         %
         %  C = ne(A,B) is called for the syntax 'A ~= B' when A or B is a
         %  TraceData object.
         %  See also eq
         tf = ~eq(A,B);
      end;

      function outT = extract(Tr, method, startPos, endPos)
         %extract   Retrieve subset of TraceData
         %
         %    subTraces = extract(traces, 'index', startPos, endPos)
         %
         %
         %   waveform = T.extract('INDEX', startPos, endPos)
         %       returns traces with the subset of data from startIndex to
         %       endIndex.  Essentially, this is equivelent to
         %       T.data = T.data(startPos, endPos)
         %
         %       data into an array, as in D = get(W,'data'), then returning a
         %       waveform with the subset of data,
         %       ie. waveform = set(waveform,'data', D(startIndex:endIndex));
         %
         %    See also SeismicTrace/extract
         switch lower(method)
            case 'index'
               assert(numel(startPos) == numel(endPos), 'number of start and end indices must match');
               outT = repmat(TraceData,numel(Tr),numel(startPos));
               for n=1:numel(T);
                  for t=1:numel(startPos)
                     outT(n,t) = Tr(n);
                     outT(n,t).data = Tr(n).data(startPos(t), endPos(t));
                  end
               end
         end
      end

      function td = fillgaps(td,value, gapvalue)
         %fillgaps   fill missing data with values of your choice
         %  T = T.fillgaps(value) fills data with the number of your choice
         %  VALUE can also be nan or inf or -inf
         %
         %  T = T.fillgaps([]) removes missing data from trace. Warning, the
         %  resulting timing issues are NOT corrected for!
         %
         %  T = T.fillgaps('method') replaces data using an interpolation method of
         %  your choice. Valid methods are:
         %
         %    'meanall' - replace missing values with the mean of the whole
         %    dataset
         %
         %    'number' - replace missing value with a numeric value of your
         %    choice (can also be Inf, -Inf, NaN)
         %
         %    'interp' - assuming missing values are marked by NaN, this will use
         %    cubic interpolation to estimate the missing values
         %
         %  FILLGAPS is designed to replace NaN values. However, if if you
         %  use T = T.fillgaps(number, gapvalue), then ALL data points with
         %  the value GAPVALUE will be replaced by NUMBER.

         if ~(isnumeric(value) && (isscalar(value) || isempty(value)) || ischar(value))
            warning('TraceData:fillgaps:invalidGapValue',...
               'Value to replace data with must be string or scalar or []');
         end

         if ~exist('gapvalue','var') || isnan(gapvalue)
            getgaps = @(x) isnan(x.data);
         else
            getgaps = @(x) x.data == gapvalue;
         end;

         if ischar(value)
            fillmethod = lower(value);
         else
            fillmethod = 'number';
         end

         switch fillmethod
            case 'meanall'
               for N = 1:numel(td);
                  allgaps = getgaps(td(N));
                  % do not include the values to be replaced
                  meanVal = mean(td(N).data(~allgaps));
                  if isnan(meanVal)
                     meanVal = 0;
                  end
                  td(N).data(allgaps) = meanVal;
               end
            case 'number'
               for N = 1:numel(td);
                  td(N).data(getgaps(td(N))) = value;
               end
            case 'interp' % blame Glenn, inspired by http://www.mathworks.com/matlabcentral/fileexchange/8225-naninterp by E Rodriguez
               for N = 1:numel(td);
                  X = td(N).data;
                  X(isnan(X)) = interp1(find(~isnan(X)), X(~isnan(X)), find(isnan(X)),'cubic');
                  td(N).data = X;
               end
            otherwise
               error('TraceData:fillgaps:unimplementedMethod',...
                  'Unimplemented fillgaps method [%s]', fillmethod);
         end
      end

      function T = zero2nan(T,mgl)
         %zero2nan   replaces zeros with nan values.
         % Should be replaced with T.fillgaps(nan, 0)
         % however, it was designed to have a minimum gap length Perhaps this
         % should be added to fillgaps
         %
         %This function replaces gaps that have been filled with zeros and
         %converts them to NaN values. This is the inverse of T =
         %T.fillgaps(0). An input mgl defines the minimum gap length to be
         %converted to NaN gaps, i.e. if only 5 consecutive zero values
         %exist in a given small gap, they will be converted to NaN values
         %if mgl <= 5 and left as zero values if mgl > 5
         %
         %USAGE: T = T.zero2nan(mgl)
         %
         %REQUIRED INPUTS:
         %   mgl - minimum gap length (datapoints) to convert to NaN values
         %
         %OUTPUTS: w - trace with gaps converted to NaN
         %
         %See also nan, fillgaps

         % Author: Dane Ketner, Alaska Volcano Observatory
         % Modified: Celso Reyes: rewrote algorithm to elimenate looping (2x faster)
         %                        results differ because old method converted
         %                        5 zeros only if mgl <5 (not <=5)

         for nw = 1:numel(T)
            closeToZero = abs(T(nw).data) < 0.1; % 0.1 was already chosen -CR

            % --- the logic below should be pulled into a new function, since it is
            % shared across a couple different areas, such as trace/clean -- %
            firstZeros = find(diff([false; closeToZero(:)]) == 1);
            lastZeros = find(diff([closeToZero(:); false]) == -1);
            assert(numel(firstZeros) == numel(lastZeros));
            nContiguousZeros = lastZeros - firstZeros + 1;
            firstZeros(nContiguousZeros < mgl) = [];
            lastZeros(nContiguousZeros < mgl) = [];

            for c=1:numel(firstZeros)
               T(nw).data(firstZeros(c) : lastZeros(c)) = NaN;
            end
         end
      end

      function t = setlength(t, maxlen)
         %setlength   adjust length of trace data to allow batch processing
         %   trace = traces.setlength()
         %       adjusts all traces to the length of the largest, while
         %       zero-padding all shorter traces
         %
         %   trace = traces.setlength(maxlength)
         %       sets all data lengths to maxlength, padding with zero or
         %       truncating as necessary.
         %
         %  examples
         %       % let traces be a 1x2 TraceData object
         %       % 3000 samples in traces(1)
         %       % 10025 samples in traces(2)
         %
         %       % set both waves' data to a length of to 10025 while padding the
         %       % smaller of the two with zeroes.
         %       outTraces = traces.setlength
         %
         %       % set both sample lengths to 500 truncating both of them...
         %       outTraces = traces.setlength(500)
         datalengths = arrayfun(@(x) numel(x.data), t);
         if ~exist('maxlen','var')
            maxlen = max(datalengths);
         end

         for n=find(datalengths < maxlen)
            t(n).data(maxlen) = 0;
         end

         for n=find(datalengths > maxlen);
            t(n).data = t(n).data(1:maxlen);
         end
      end

      function T = hilbert(T, n)
         %hilbert   Discrete-time analytic Hilbert transform for traces.
         %   trace = trace.hilbert()
         %   trace = trace.hilbert(N);
         %
         % THIS version only returns the abs value in the trace. If you
         % want to keep the imaginary values, then you should use the
         % built-in hilbert transform. ie. Don't feed it a trace, feed it
         % a vector... - CR
         %
         % See also FFT, IFFT, HILBERT

         if nargin==2
            T = arrayfun(@myHilbertN,T);
         else
            T = arrayfun(@myHilbert,T);
         end

         function Tr = myHilbert(Tr)
            Tr.data = abs(hilbert(Tr.data));
            Tr.units = 'abs(hilbert)';
         end

         function Tr = myHilbertN(Tr)
            Tr.data = abs(hilbert(Tr.data, n));
            Tr.units = 'abs(hilbert)';
         end
      end

      function T = resample(T, method, crunchFactor)
         %resample   Resample a trace over a specified interval
         %   T = trace.resample(method, crunchfactor)
         %
         %   Input Arguments
         %       Trace: TraceData or Trace       N-dimensional
         %
         %       METHOD: which method of sampling to perform within each sample
         %                window
         %           'max' : maximum value
         %           'min' : minimum value
         %           'mean': average value
         %           'median' : median value
         %           'rms' : rms value (added 2011/06/01)
         %           'absmax': absolute maximum value (greatest deviation from zero)
         %           'absmin': absolute minimum value (smallest deviation from zero)
         %           'absmean' : mean deviation from zero (added 2011/06/01)
         %           'absmedian' : median deviation from zero (added 2011/06/01)
         %           'builtin': Use MATLAB's built in resample routine
         %
         %       CRUNCHFACTOR : the number of samples making up the sample window
         %
         % For example, T.resample('max',5) would grab the max value of every 5
         % samples and return that in a trace with an adjusted frequency. as a
         % result, the trace will have 1/5 of the samples
         %
         %
         % To use matlab's built-in RESAMPLE method...
         %       % assume T is an existing trace
         %       D = double(T);
         %       ResampleD = D.resample(P,Q);  % see matlab's RESAMPLE for specifics
         %
         %       %put back into trace, but don't forget to update the frequency
         %       T.data = ResampleD;
         %       T.samplerate = NewFrequency;
         %
         % See also RESAMPLE, MIN, MAX, MEAN, MEDIAN.

         % AUTHOR: Celso Reyes, Glenn Thompson

         persistent STATS_INSTALLED;
         if isempty(STATS_INSTALLED)
            STATS_INSTALLED = ~isempty(ver('stats'));
         end

         if ~(round(crunchFactor) == crunchFactor)
            disp ('val needs to be an integer');
            return;
         end;
         method = lower(method);

         % determine which function to use on the data. This will
         % automatically determine whether it can use the NaN-safe version
         if ~strcmp(method,'builtin')
            if STATS_INSTALLED
               if numel(method) > 3 && strcmpi(method(1:3),'ABS')
                  methodfn = str2func(['nan' method(4:end)]);
               else
                  methodfn = str2func(['nan', method]);
               end
            else
               if numel(method) > 3 && strcmpi(method(1:3),'ABS')
                  methodfn = str2func(method(4:end));
               else
                  methodfn = str2func(method);
               end
            end
         end

         for i=1:numel(T)
            nRows = ceil(length(T(i).data) / crunchFactor);
            totVals = nRows * crunchFactor; % total number of values that can be accomodated
            if length(T(i).data) < totVals
               T(i).data(end+1:totVals) = mean(T(i).data((nRows-1)*totVals : end)); %pad it with the avg value
            end;

            d = reshape(T(i).data,crunchFactor,nRows); % produces ( val x rowcount) matrix
            switch upper(method)
               case {'MAX', 'MIN', 'RMS'}
                  T(i).data = methodfn(d, [], 1);
               case {'MEAN', 'MEDIAN'}
                  T(i).data = methodfn(d, 1);
               case {'ABSMAX', 'ABSMIN'}
                  T(i).data = methodfn(abs(d),[],1);
               case {'ABSMEAN', 'ABSMEDIAN'}
                  T(i).data = methodfn(abs(d), 1);
               case 'BUILTIN'
                  % assume T is an existing trace
                  ResampleD = resample(T(i).data,1,crunchFactor);  % see matlab's RESAMPLE for specifics
                  T(i).data = ResampleD(:);
               otherwise
                  error('TraceData:resample:UnknownSampleMethod',...
                     'Don''t know what you mean by resample via %s', method);
            end;
            T(i).samplerate = T(i).samplerate ./ crunchFactor;
         end
      end

      function A = smooth(A, varargin)
         %smooth   Smooth response data
         %  A = smooth(T) smooths the Trace's data using a moving average
         %  filter.
         %
         %  smooth requires the signal fitting toolbox.  For additional
         %  options, see that smooth
         % See also smooth
         for n = 1:numel(A)
            A(n).data = smooth(A(n).data,varargin{:});
         end
      end

      function T = taper(T, style, R)
         %taper   apply a taper to a trace
         % trace = trace.taper() applies a cosine (tukey) taper to the ends
         % of a trace with a default taper to the first and last 10% of the
         % trace. This is same as trace.taper('tukey', 0.2)
         %
         % trace = trace.taper(style) applies the tapering window function
         % STYLE to the trace. STYLE can be any valid windowing function
         % some possible taper styles include hanning, tukey, and gaussian.
         % see help for WINDOW for a more complete list.
         %
         % trace = trace.taper(style, R)
         % R varies in meaning according to the window style. For the
         % default cosine (tukey) taper, R is the ratio of tapered to
         % constant sections and is between 0 and 1. For example, if R =
         % 0.1 then the taper at each end of the trace is 5% of the total
         % trace length. R can be either a scalar or the same size as
         % TRACE. If R is a scalar, it is applied uniformly to each
         % trace.
         %
         % for a tukey taper, setting R ito 1 results in a hanning taper
         %
         % All tapers require the signal processing toolbox.
         %
         % See also window
         % AUTHOR: Michael West
         % Modified: Celso Reyes
         a=ver;
         HAVE_SIGBOX = ismember('Signal Processing Toolbox',{a.Name});
         if ~HAVE_SIGBOX
            error('TraceData:taper:signalToolboxNotInstalled',...
               'Using a taper requires access to the Signal Processing Toolbox');
         end
         if ~exist('style','var') || isempty(style) || strcmpi(style,'cosine')
            style = 'tukeywin';
         end
         if exist(lower(style))
            taperfun = str2func(lower(style));
         else
            error('TraceData:taper:invalidTaperType',...
               'This style of taper is not recognized.');
         end

         %% massage R into place
         if strcmpi(style,'tukeywin') && (~exist('R','var') || isempty(R))
            R = 0.2; %assign default taper
         end
         if exist('R','var') &&  ~isnumeric(R)
            error('TraceData:taper:InvalidRValue',...
               'R, if specified, must be numeric');
         end

         if isscalar(R)
            R = repmat(R,size(T));
         end

         if (isvector(R) && isvector(T)) && numel(R) == numel(T)
            if all(size(T)) ~= size(R)
               % same number of elements, but R is Nx1 and w is 1xN or vice-versa
               warning('TraceData:taper:columnsVsRows',...
                  ['One input (either R or the wavform) is arranged in '...
                  'columns while the other is arranged in Rows. While they '...
                  'should be the same shape, taper is continuing with R''']);
               R = R';
            end
         end

         if ~all(size(T) == size(R))
            error('TraceData:taper:InvalidRSize',...
               'R must either be a scalar value, or must be the same size as the input traces');
         end
         nsamples = T.nsamples();
         %% Do the window processing
         if exist('R','var')
            for N=1:numel(T)
               T(N) = T(N) .* window(taperfun, nsamples(N), R(N));
            end
         else
            for N=1:numel(T)
               T(N) = T(N) .* window(taperfun, nsamples(N));
            end
         end

      end

      function [tf, msg] = compatiblewith(A, B)
         %compatiblewith   confirm units, samplerate, and data length match
         %  tf = A.compatiblewith(B) will make sure that A and B have
         %  matching units, sample frequencies, and data length.
         %  If units are empty for either A or B, then units will pass
         %  If samplerate is empty for either A or B, then samplerate will
         %  pass. samplerate is compared within a tolerance of 10e-2
         %
         %  [tf, msg] = A.compatiblewith(B) will also return a message
         %  describing how the test failed.
         %
         %  sample usage:
         %    % assume A & B are TraceData
         %    if A.compatiblewith(B)
         %       dosomething(A, B)
         %    end
         %
         % See also assertCompatiblewith
         TOL = 10e-2;
         if ~(numel(A.data) == numel(B.data))
            tf = false;
            msg = 'Data lengths do not match';
         elseif ~(isempty(B.units) ||isempty(A.units) || strcmp(A.units, B.units))
            tf = false;
            msg = sprintf('Units do not match [%s] vs [%s]', A.units, B.units);
         elseif ~(isempty(A.samplerate) || isempty(B.samplerate) || isnan(A.samplerate) || isnan(B.samplerate) || ismembertol(A.samplerate,B.samplerate, TOL))
            tf = false;
            msg = sprintf('Frequencies do not match [%f] vs [%f]',A.samplerate, B.samplerate);
         else
            tf = true;
            msg = '';
         end
      end

      function assertCompatiblewith(A, B)
         %assertCompatiblewith  asserts units, samplerate, and data length match
         %  assertCompatiblewith(A, B) will error unless A and B have
         %  matching units, sample frequencies, and data length. The error
         %  message will describe what is wrong
         %
         %  use this within a try-catch, or when debugging
         %
         % sample debugging usage:
         %   % assume A & B are TraceData
         %   A.assertCompatiblewith(B)
         %   C = A + B.data;
         %
         % sample prduction usage:
         %   try
         %     A.assertCompatiblewith(B)
         %     C = A + B.data;
         %   catch er
         %     % handle this failure somehow
         %   end
         %
         % See also compatiblewith
         [tf, msg] = compatiblewith(A, B);
         assert(tf, 'TraceData:compatabilityCheckFailed', msg);
      end

      %% stacking functions
      function out = stack(T)
         %stack  stacks data from array of traces
         %   StackedTraces = stack(traces)
         %   ASSUMES frequencies are the same. data does not need to be the same
         %   length, but shorter traces will be padded with zeros at the end
         %   prior to stacking.
         %
         %   Stacks all traces, regardless of dimension
         %
         %   Data is summed, but the average is not taken, nor is it normalized. You
         %   may wish to change the station and/or channel names to reflect the
         %   properties of this trace. Possibly change the units, also, if that
         %   makes sense.
         %
         %   Output retains the same info as the very first trace (minus history)
         %   the station name becomes the original name with "- stack (N)" tacked onto
         %   the end.
         %
         %   To ensure frequency and time matching, use ALIGN
         %
         %   See also trace.align

         out = T(1);
         out.station = [out.station ' - stack (' num2str(numel(T)) ')'];
         out.data = sum(double(T),2);
      end
      function stackedTraces = binStack(T,nBins,binOverlap)
         %binStack   Stack traces with specified number of traces per stack.
         %
         %USAGE: stk = binStack(T,nBins,binOverlap)
         %    If input T contains 100 traces then:
         %    stk = binStack(T,20,0) % stk has 5 stacks of 20 traces
         %    stk = binStack(T,5,0)  % stk has 20 stacks of 5 traces
         %    stk = binStack(T,20,10)  % stk has 9 stacks of 20 traces
         %
         %OUTPUTS: stk - Array of stacked traces

         % Author: Dane Ketner, Alaska Volcano Observatory

         nw = numel(T);
         inc = nBins - binOverlap;
         n=1;
         while 1
            if (n-1)*inc+nBins > nw
               return
            else
               stackedTraces(n) = stack(T((n-1)*inc+1:(n-1)*inc+nBins));
               n=n+1;
            end
         end
      end
      function info(T)
            fprintf('[%s] TraceData array with:\n', TraceData.size2str(T));
            fprintf('  nSaples  sampleRate      Units Duration(sec)   min      max   hasnan?\n');
         for n=1:numel(T)
            if T(n).hasnan; hnt = 'Y'; else hnt = 'N';end
            fprintf('%8d    %f %10s   %8.3f  %8.3f  %8.3f     %s \n',...
               T(n).nsamples, T(n).samplerate, T(n).units, T(n).duration, min(T(n)), max(T(n)), hnt);
         end
      end
       %}
   end

   %% protected methods
   methods(Static, Access=protected)
      function s = size2str(A)
         sz = size(A);
         s = num2str(sz(1));
         for n=2:numel(sz)
            s = [s, 'x', num2str(sz(n))];
         end
      end
   end

   methods(Static)
      function setParameter(name, val)
         %setParameter  controlstate behavior for the TraceData objects
         %
         %     TRUST_ASSIGNMENTS
         %     set_parameter('trust_assignments', true) will verify the
         %     shapes and types of each calculation before storing the data
         %     in the object
         %     set_parameter('trust_assignments', false) will not perform
         %     size/class checks. Use this only when you are looking for
         %     extra speed and tightly control the inputs to this class.
         %
         %     DEBUG_LEVEL
         disp(val)
         disp('not active yet');
         switch lower(name)
            case 'trust_assignments'
               % set the parameter trust_assignments to logical(val);
            case 'debug_level'
               % set the parameter debug_level to numerical val
            otherwise
         end
      end

      newunit = autoscale(axishandle, oldunit);

   end
end

author: geoscience-community-codes

repo: GISMO


text:
function T = multiply_pots(T1, T2)
% MULTIPLY_POTS Multiply a pair of dpots together pointwise (cgpot)
% T = multiply_pots(pots)

ddom = myunion(T1.ddom, T2.ddom);
cdom = myunion(T1.cdom, T2.cdom);
dom = myunion(ddom, cdom);
ns = zeros(1, max(dom));
ns(T1.ddom) = T1.dsizes;
ns(T2.ddom) = T2.dsizes;
ns(T1.cdom) = T1.csizes;
ns(T2.cdom) = T2.csizes;

T = cgpot(ddom, cdom, ns);
T = multiply_by_pot(T, T1);
T = multiply_by_pot(T, T2);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


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

%%
clear; close all; clc;

%% Syntax
% |[Nr,Nr_units]=MRimageUnitScale(N,N_info);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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:
function [jac,err] = jacobianest(fun,x0)
% gradest: estimate of the Jacobian matrix of a vector valued function of n variables
% usage: [jac,err] = jacobianest(fun,x0)
%
%
% arguments: (input)
%  fun - (vector valued) analytical function to differentiate.
%        fun must be a function of the vector or array x0.
%
%  x0  - vector location at which to differentiate fun
%        If x0 is an nxm array, then fun is assumed to be
%        a function of n*m variables.
%
%
% arguments: (output)
%  jac - array of first partial derivatives of fun.
%        Assuming that x0 is a vector of length p
%        and fun returns a vector of length n, then
%        jac will be an array of size (n,p)
%
%  err - vector of error estimates corresponding to
%        each partial derivative in jac.
%
%
% Example: (nonlinear least squares)
%  xdata = (0:.1:1)';
%  ydata = 1+2*exp(0.75*xdata);
%  fun = @(c) ((c(1)+c(2)*exp(c(3)*xdata)) - ydata).^2;
%
%  [jac,err] = jacobianest(fun,[1 1 1])
%
%  jac =
%           -2           -2            0
%      -2.1012      -2.3222     -0.23222
%      -2.2045      -2.6926     -0.53852
%      -2.3096      -3.1176     -0.93528
%      -2.4158      -3.6039      -1.4416
%      -2.5225      -4.1589      -2.0795
%       -2.629      -4.7904      -2.8742
%      -2.7343      -5.5063      -3.8544
%      -2.8374      -6.3147      -5.0518
%      -2.9369      -7.2237      -6.5013
%      -3.0314      -8.2403      -8.2403
%
%  err =
%   5.0134e-15   5.0134e-15            0
%   5.0134e-15            0   2.8211e-14
%   5.0134e-15   8.6834e-15   1.5804e-14
%            0     7.09e-15   3.8227e-13
%   5.0134e-15   5.0134e-15   7.5201e-15
%   5.0134e-15   1.0027e-14   2.9233e-14
%   5.0134e-15            0   6.0585e-13
%   5.0134e-15   1.0027e-14   7.2673e-13
%   5.0134e-15   1.0027e-14   3.0495e-13
%   5.0134e-15   1.0027e-14   3.1707e-14
%   5.0134e-15   2.0053e-14   1.4013e-12
%
%  (At [1 2 0.75], jac should be numerically zero)
%
%
% See also: derivest, gradient, gradest
%
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 3/6/2007

% get the length of x0 for the size of jac
nx = numel(x0);

MaxStep = 100;
StepRatio = 2.0000001;

% was a string supplied?
if ischar(fun)
    fun = str2func(fun);
end

% get fun at the center point
f0 = fun(x0);
f0 = f0(:);
n = length(f0);
if n==0
    % empty begets empty
    jac = zeros(0,nx);
    err = jac;
    return
end

relativedelta = MaxStep*StepRatio .^(0:-1:-25);
nsteps = length(relativedelta);

% total number of derivatives we will need to take
jac = zeros(n,nx);
err = jac;
for i = 1:nx
    x0_i = x0(i);
    if x0_i ~= 0
        delta = x0_i*relativedelta;
    else
        delta = relativedelta;
    end

    % evaluate at each step, centered around x0_i
    % difference to give a second order estimate
    fdel = zeros(n,nsteps);
    for j = 1:nsteps
        fdif = fun(swapelement(x0,i,x0_i + delta(j))) - ...
            fun(swapelement(x0,i,x0_i - delta(j)));

        fdel(:,j) = fdif(:);
    end

    % these are pure second order estimates of the
    % first derivative, for each trial delta.
    derest = fdel.*repmat(0.5 ./ delta,n,1);

    % The error term on these estimates has a second order
    % component, but also some 4th and 6th order terms in it.
    % Use Romberg exrapolation to improve the estimates to
    % 6th order, as well as to provide the error estimate.

    % loop here, as rombextrap coupled with the trimming
    % will get complicated otherwise.
    for j = 1:n
        [der_romb,errest] = rombextrap(StepRatio,derest(j,:),[2 4]);

        % trim off 3 estimates at each end of the scale
        nest = length(der_romb);
        trim = [1:3, nest+(-2:0)];
        [der_romb,tags] = sort(der_romb);
        der_romb(trim) = [];
        tags(trim) = [];

        errest = errest(tags);

        % now pick the estimate with the lowest predicted error
        [err(j,i),ind] = min(errest);
        jac(j,i) = der_romb(ind);
    end
end

end % mainline function end

% =======================================
%      sub-functions
% =======================================
function vec = swapelement(vec,ind,val)
% swaps val as element ind, into the vector vec
vec(ind) = val;

end % sub-function end

% ============================================
% subfunction - romberg extrapolation
% ============================================
function [der_romb,errest] = rombextrap(StepRatio,der_init,rombexpon)
% do romberg extrapolation for each estimate
%
%  StepRatio - Ratio decrease in step
%  der_init - initial derivative estimates
%  rombexpon - higher order terms to cancel using the romberg step
%
%  der_romb - derivative estimates returned
%  errest - error estimates
%  amp - noise amplification factor due to the romberg step

srinv = 1/StepRatio;

% do nothing if no romberg terms
nexpon = length(rombexpon);
rmat = ones(nexpon+2,nexpon+1);
% two romberg terms
rmat(2,2:3) = srinv.^rombexpon;
rmat(3,2:3) = srinv.^(2*rombexpon);
rmat(4,2:3) = srinv.^(3*rombexpon);

% qr factorization used for the extrapolation as well
% as the uncertainty estimates
[qromb,rromb] = qr(rmat,0);

% the noise amplification is further amplified by the Romberg step.
% amp = cond(rromb);

% this does the extrapolation to a zero step size.
ne = length(der_init);
rhs = vec2mat(der_init,nexpon+2,ne - (nexpon+2));
rombcoefs = rromb\(qromb'*rhs);
der_romb = rombcoefs(1,:)';

% uncertainty estimate of derivative prediction
s = sqrt(sum((rhs - rmat*rombcoefs).^2,1));
rinv = rromb\eye(nexpon+1);
cov1 = sum(rinv.^2,2); % 1 spare dof
errest = s'*12.7062047361747*sqrt(cov1(1));

end % rombextrap

% ============================================
% subfunction - vec2mat
% ============================================
function mat = vec2mat(vec,n,m)
% forms the matrix M, such that M(i,j) = vec(i+j-1)
[i,j] = ndgrid(1:n,0:m-1);
ind = i+j;
mat = vec(ind);
if n==1
    mat = mat';
end

end % vec2mat

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
randinit
clf
hold on
g = PGraph(3);

sim_time = 2;

xinit = [0, 0, 0];
plot(xinit(1), xinit(2), 'go');
rho = 1;

g.add_node(xinit);
t =2;

for k=1:500
    % Step 3
    % find random state x,y in [-10, 10], th in [0 2pi]
    x = (rand-0.5)*20;
    y = (rand-0.5)*20;
    theta = rand*2*pi;
    xrand = [x, y, theta]';
    xrand'

    % Step 4
    % find the existing node closest in state space
    dmin = Inf;
    for v=1:g.n
        xv = g.coord(v);
        d = sum( (xv(1:2)-xrand(1:2)).^2 ) + t*angdiff(xv(3)-xrand(3))^2;
        if d < dmin
            xnear = xv;
            vnear = v;
            dmin = d;
        end
    end

    % Step 5
    % figure how to drive the robot from xnear to xrand
    x0 = xnear';
    xg = xrand';
    %set_param('sl_movepoint2/Bicycle', 'x0', sprintf('[%f,%f,%f]', xnear) );
    %set_param('sl_movepoint2/Goal pose', 'Value', sprintf('trans2([%f,%f,0])', xrand(1:2)) );

    r = sim('sl_drivepose', sim_time);
    t = r.find('tout');
    y = r.find('yout');
    plot2(y)
    drawnow
    xnew = y(end,:)';
    plot2(xnew', 'go');

    % ensure that the path is collision free

    % Step 7,8
    % add xnew to the graph, with an edge to xnear
    g.add_node(xnew, vnear);
end

grid on
xyzlabel
zlabel('\theta');
axis([-10 10 -10 10])
iprint('rrt_paths');
view(0, 0);
iprint('rrt_paths_xt');

author: petercorke

repo: robotics-toolbox-matlab


text:
function fea = NormalizeFea(fea,row)
% if row == 1, normalize each row of fea to have unit norm;
% if row == 0, normalize each column of fea to have unit norm;

if ~exist('row','var')
    row = 1;
end

% if row
%     nSmp = size(fea,1);
%     feaNorm = max(1e-14,full(sum(fea.^2,2)));
%     fea = spdiags(feaNorm.^-.5,0,nSmp,nSmp)*fea;
% else
%     nSmp = size(fea,2);
%     feaNorm = max(1e-14,full(sum(fea.^2,1))');
%     fea = fea*spdiags(feaNorm.^-.5,0,nSmp,nSmp);
% end
%
% return;

if row
    [nSmp, mFea] = size(fea);
    if issparse(fea)
        fea2 = fea';
        feaNorm = mynorm(fea2,1);
        for i = 1:nSmp
            fea2(:,i) = fea2(:,i) ./ max(1e-10,feaNorm(i));
        end
        fea = fea2';
    else
        feaNorm = sum(fea.^2,2).^.5;
        fea = fea./feaNorm(:,ones(1,mFea));
    end
else
    [mFea, nSmp] = size(fea);
    if issparse(fea)
        feaNorm = mynorm(fea,1);
        for i = 1:nSmp
            fea(:,i) = fea(:,i) ./ max(1e-10,feaNorm(i));
        end
    else
        feaNorm = sum(fea.^2,1).^.5;
        fea = fea./feaNorm(ones(1,mFea),:);
    end
end

author: Elyorcv

repo: SAE

save_path: github-repos/MATLAB/Elyorcv-SAE


text:
function [dlnZ_dmu, dlnZ_dvs] = cmpndNoiseGradVals(noise, mu, varsigma, y)

% CMPNDNOISEGRADVALS Gradient of CMPND noise log Z with respect to input mean and variance.
% FORMAT
% DESC computes the gradient of the compound
% noise with respect to the input mean and the input variance.
% ARG noise : noise structure for which gradients are being
% computed.
% ARG mu : mean input locations with respect to which gradients are
% being computed.
% ARG varSigma : variance input locations with respect to which
% gradients are being computed.
% ARG y : noise model output observed values associated with the given points.
% RETURN dlnZ_dmu : the gradient of log Z with respect to the input mean.
% RETURN dlnZ_dvs : the gradient of log Z with respect to the input variance.
%
% SEEALSO cmpndNoiseParamInit, cmpndNoiseGradientParam, noiseGradVals,
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

startVal = 1;
endVal = 0;
dlnZ_dmu = zeros(size(mu));
dlnZ_dvs = zeros(size(varsigma));
for i = 1:length(noise.comp)
  fhandle = str2func([noise.comp{i}.type 'NoiseGradVals']);
  [dlnZ_dmu(:, i), dlnZ_dvs(:, i)]  = fhandle(noise.comp{i}, ...
                                              mu(:, i), ...
                                              varsigma(:, i), ...
                                              y(:, i));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [BNEW,U,EXITFLAG,OUTPUT] = genNpReg2(A,B,options,defaultopt,regDim,varargin);
%GENNPREG2 solves general nonparametric image registration problems
%using a fixed point iteration
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% check whether the combination of options is supported
[SimMeasFlag,RegularizerFlag,BoundCondFlag] = checkOptions(options,defaultopt);

% initialize similarity measure value and type of comparison
[SimMeas,simMeasComparison] = initializeSimilarityMeasure(SimMeasFlag);
SimMeasNew = SimMeas;

% get appropriate function handle for computing body force
computeBodyForce = getBodyForceFunction(SimMeasFlag);

% get appropriate function handle for solving core PDE based on regularizer
% and boundary conditions
solveCorePDE = getCorePDESolverFunction(RegularizerFlag,BoundCondFlag,regDim);

% get various fields from options structure
DisplayFlag = npRegGet(options,'Display',defaultopt,'fast');
VoxSize = [npRegGet(options,'VoxSizeX',defaultopt,'fast') ...
    npRegGet(options,'VoxSizeY',defaultopt,'fast') ...
    npRegGet(options,'VoxSizeZ',defaultopt,'fast')];
MaxIter = npRegGet(options,'MaxIter',defaultopt,'fast');
UDiffTol = npRegGet(options,'UDiffTol',defaultopt,'fast');
VDiffTol = npRegGet(options,'VDiffTol',defaultopt,'fast');
BodyForceTol = npRegGet(options,'BodyForceTol',defaultopt,'fast');
SimMeasTol = npRegGet(options,'SimMeasTol',defaultopt,'fast');
BodyForceDiffTol = npRegGet(options,'BodyForceDiffTol',defaultopt,'fast');
SimMeasDiffTol = npRegGet(options,'SimMeasDiffTol',defaultopt,'fast');
SimMeasPercentDiffTol = npRegGet(options,'SimMeasPercentDiffTol',defaultopt,'fast');
FixedPointMaxFlowDistance = npRegGet(options,'FixedPointMaxFlowDistance',defaultopt,'fast');
RegridTol = npRegGet(options,'RegridTol',defaultopt,'fast');
Mu = npRegGet(options,'Mu',defaultopt,'fast');
Lambda = npRegGet(options,'Lambda',defaultopt,'fast');
ForceFactor = npRegGet(options,'ForceFactor',defaultopt,'fast');
RegularizerFactor = npRegGet(options,'RegularizerFactor',defaultopt,'fast');
StabilityConstant = npRegGet(options,'StabilityConstant',defaultopt,'fast');

% get size of images
[NumRows,NumCols,NumPages] = size(A);

% initialize various variables
if regDim == 2 % 2-D case
    [U,V,BodyForce,BodyForceNew,DU,X,BGrad] = deal(zeros(NumRows,NumCols,2));
    DefGrad = zeros(NumRows,NumCols,2,2);
    DefJac = zeros(NumRows,NumCols);
else % 3-D case
    [U,V,BodyForce,BodyForceNew,DU,X,BGrad] = deal(zeros(NumRows,NumCols,NumPages,3));
    DefGrad = zeros(NumRows,NumCols,NumPages,3,3);
    DefJac = zeros(NumRows,NumCols,NumPages);
end
MaxDUNorm = inf;
MeanDUNorm = inf;
MedDUNorm = inf;

% initialize vectors of grid positions
if regDim == 2 % 2-D case
    [X(:,:,1),X(:,:,2)] = ndgrid((0:(NumRows-1))*VoxSize(1),...
        (0:(NumCols-1))*VoxSize(2));
else % 3-D case
    [X(:,:,:,1),X(:,:,:,2),X(:,:,:,3)] = ndgrid((0:(NumRows-1))*VoxSize(1),...
        (0:(NumCols-1))*VoxSize(2),(0:(NumPages-1))*VoxSize(3));
end

% compute gradient of floating image, as this will be needed repeatedly
% within incremental updating
[HX,HY,HZ] = constructGradientFilters(VoxSize);
if regDim == 2 % 2-D case
    %HX = imfilter(HX,fspecial('gaussian',7,2),'full');
    %HY = imfilter(HY,fspecial('gaussian',7,2),'full');
    BGrad(:,:,1) = imfilter(B,HX,'replicate','same');
    BGrad(:,:,2) = imfilter(B,HY,'replicate','same');
else % 3-D case
    BGrad(:,:,:,1) = imfilter(B,HX,'replicate','same');
    BGrad(:,:,:,2) = imfilter(B,HY,'replicate','same');
    BGrad(:,:,:,3) = imfilter(B,HZ,'replicate','same');
end

% follow general framework of Bro-Nielsen and Gramkow, but allowing for
% different body force and regularizers

if isequal(DisplayFlag,'iter')
    displayString = sprintf('Iteration\tSimilarity\tBodyForceMax\tDisplaceMax\tMean\tMedian');
    disp(displayString);
end

% fixed point iteration loop
i = 1;
while i <= MaxIter

    % compute body force
    [BodyForceNew(:),SimMeasNew(:)] = computeBodyForce(A,B,BGrad,U,X,NumRows,NumCols,NumPages,VoxSize,regDim,HX,HY,HZ);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%% try to make robust
    %if SimMeasNew > SimMeas
    %    ForceFactor = ForceFactor./((1+sqrt(5))/2);
    %    disp('Force Factor reduced.');
    %end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    BodyForceNew(:) = BodyForceNew*ForceFactor;

    % solve PDE at current iteration
    %DU(:) = solveCorePDE(U,BodyForceNew,StabilityConstant,VoxSize,NumRows,NumCols,NumPages,Mu,Lambda,RegularizerFactor,DU,HX,HY,HZ);  % why not -U?
    DU(2:end-1,2:end-1,:) = RegularizerFactor.*imfilter(U(2:end-1,2:end-1,:),[0 1 0;1 -4 1;0 1 0],0) - BodyForceNew(2:end-1,2:end-1,:);
    %keyboard
    % compute time step based on maximum flow distance
    [TimeStep,MaxDUNorm(:),MeanDUNorm(:),MedDUNorm(:)] = ...
        computeTimeStep(DU,FixedPointMaxFlowDistance);

    % now add back to displacement field
    U(:) = U + TimeStep.*DU;

    % display information
    if isequal(DisplayFlag,'iter')
       displayString = sprintf('%5d\t\t%g\t\t%g\t\t%g\t\t%g\t\t%g',i,SimMeasNew,max(BodyForceNew(:)),MaxDUNorm,MeanDUNorm,MedDUNorm);
       disp(displayString);
    end

    % check tolerances for convergence
    [BreakCond,BreakMsg] = checkThresholds(BodyForce,BodyForceNew,SimMeas,...
        SimMeasNew,MaxDUNorm,UDiffTol,VDiffTol,BodyForceTol,SimMeasTol,...
        BodyForceDiffTol,SimMeasDiffTol,SimMeasPercentDiffTol,simMeasComparison);
    if BreakCond
        disp(BreakMsg);
        break;
    end

    % update similarity measure and body force
    SimMeas(:) = SimMeasNew;
    BodyForce(:) = BodyForceNew;

    % finally increment counter
    i = i+1;

end

% update image B to final deformation grid
if regDim == 2 % 2-D case
    BNEW = interp2(X(:,:,2),X(:,:,1),B,X(:,:,2)-U(:,:,2),...
        X(:,:,1)-U(:,:,1),'*linear',0);
else % 3-D case
    BNEW = interp3(X(:,:,:,2),X(:,:,:,1),X(:,:,:,3),B,...
        X(:,:,:,2)-U(:,:,:,2),X(:,:,:,1)-U(:,:,:,1),...
        X(:,:,:,3)-U(:,:,:,3),'*linear',0);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% need to add some code to return collection of deformation fields
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set EXITFLAG
if i==(MaxIter+1) % maximum number of iterations reached
    EXITFLAG = 0;
else % solution found in fewer than maximum number of iterations
    EXITFLAG = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% need to add some code to catch errors thrown by any called function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set OUTPUT
OUTPUT.fixedPointIterations = i;
switch EXITFLAG
    case 1
        OUTPUT.message = 'NPREG2D converged to a solution.';
        if isequal(DisplayFlag,'iter') | isequal(DisplayFlag,'final')
            displayString = sprintf(OUTPUT.message);
            disp(displayString);
        end
    case 0
        OUTPUT.message = 'Maximum number of iterations reached.';
        if isequal(DisplayFlag,'iter') | isequal(DisplayFlag,'final')
            displayString = sprintf('%s%s',OUTPUT.message,' Terminating.');
            disp(displayString);
        end
    otherwise
        OUTPUT.message = 'If you see this message, you''re fucked.';
end

%-----------------------------------------------------------------------------------------
function [BreakCond,BreakMsg] = checkThresholds(BodyForce,BodyForceNew,SimMeas,SimMeasNew,MaxDUNorm,UDiffTol,VDiffTol,BodyForceTol,SimMeasTol,BodyForceDiffTol,SimMeasDiffTol,SimMeasPercentDiffTol,simMeasComparison);

% initialize break condition to zero
BreakCond = false;
BreakMsg = '';

% if body force is below a threshold for all points, then stop
MaxBodyForce = max(abs(BodyForceNew(:)));
if MaxBodyForce < BodyForceTol
    BreakMsg = 'Body force threshold met.';
    BreakCond = true;
end

% if body force difference is below a threshold for all points, then stop
MaxBodyForceDiff = max(abs(BodyForceNew(:)-BodyForce(:)));
if MaxBodyForceDiff < BodyForceDiffTol
    BreakMsg = 'Body force difference threshold met.';
    BreakCond = true;
end

% if similarity measure is beneath its threshold, then stop
if simMeasComparison(SimMeasNew,SimMeasTol)
    BreakMsg = 'Similarity Measure threshold met.';
    BreakCond = true;
end

% if similarity measure difference is beneath its threshold, then stop
SimMeasDiff = SimMeasNew - SimMeas;
if simMeasComparison(SimMeasNew,SimMeas)
    if abs(SimMeasDiff) < SimMeasDiffTol
        BreakMsg = 'Similarity measure difference threshold met.';
        BreakCond = true;
    end
end

% if similarity measure percentage decrease is beneath its threshold,
% then stop
if abs(SimMeas) < eps
    SimMeasPercentDiff = inf;
else
    SimMeasPercentDiff = SimMeasDiff./SimMeas;
end
if simMeasComparison(SimMeasNew,SimMeas)
    if abs(SimMeasPercentDiff) < SimMeasPercentDiffTol
        BreakMsg = 'Similarity measure percentage decrease threshold met.';
        BreakCond = true;
    end
end

% if maximum flow (velocity or displacement) is less than tolerance, then
% stop
if MaxDUNorm < UDiffTol
    BreakMsg = 'Maximum displacement threshold met.';
    BreakCond = true;
elseif MaxDUNorm < VDiffTol % this is wrong - needs to be changed
    BreakMsg = 'Maximum velocity threshold met.';
    BreakCond = true;
end

%-----------------------------------------------------------------------------------------
function [TimeStep,MaxDUNorm,MeanDUNorm,MedDUNorm] = computeTimeStep(DU,MaxFlow);

DUNorm = sqrt(sum(DU.^2,3));
MaxDUNorm = max(DUNorm(:));
MeanDUNorm = mean(DUNorm(:));
MedDUNorm = median(DUNorm(:));

if MaxDUNorm < MaxFlow
    TimeStep = 1;
else
    TimeStep = MaxFlow/MaxDUNorm;
end

%-----------------------------------------------------------------------------------------
function [HX,HY,HZ] = constructGradientFilters(VoxSize);

HX = [-1;0;1]/(2*VoxSize(1));
HY = [-1 0 1]/(2*VoxSize(2));
HZ = cat(3,-1,0,1)/(2*VoxSize(3));

%-----------------------------------------------------------------------------------------
function Y = getCorePDESolverFunction(Regularizer,BoundCond,regDim);

RegName = lower(Regularizer);
switch lower(BoundCond)
    case 'dirichlet'
        BoundCondName = 'Dirichlet';
    case 'periodic'
        BoundCondName = 'Periodic';
    case 'neumann'
        BoundCondName = 'Neumann';
end

Y = eval(['@',RegName,BoundCondName,int2str(regDim),'D']);

%-----------------------------------------------------------------------------------------
function Y = getBodyForceFunction(SimMeasFlag);

switch SimMeasFlag
    case 'SSD'
        Y = @bodyForceSSD;
    case 'NCC'
        Y = @bodyForceNCC;
    case 'CR'
        Y = @bodyForceCR;
    case 'MI'
        Y = @bodyForceMI;
    case 'NMI'
        Y = @bodyForceNMI;
end

%-----------------------------------------------------------------------------------------
function [S,comp] = initializeSimilarityMeasure(Flag);

switch Flag
    case 'SSD'
        S = inf;
        comp = @lt; % set comparison function to less than, as we are minimizing
    case {'NCC','CR','MI','NMI'}
        S = -inf;
        comp = @gt; % set comparison function to greater than, as we are maximizing
end

%-----------------------------------------------------------------------------------------
function [simMeasFlag,regularizerFlag,boundaryCondFlag] = checkOptions(options,defaultopt);

simMeasFlag = npRegGet(options,'SimilarityMeasure',defaultopt,'fast');
regularizerFlag = npRegGet(options,'Regularizer',defaultopt,'fast');
boundaryCondFlag = npRegGet(options,'BoundaryCond',defaultopt,'fast');

% test similarity measure support
if isequal(simMeasFlag,'NCC')
    errid = 'npRegLib:genNpReg:checkOptions:NCCNotSupported';
    errmsg = 'Normalized Cross Correlation not currently supported as similarity measure.';
elseif isequal(simMeasFlag,'CR')
    errid = 'npRegLib:genNpReg:checkOptions:CRNotSupported';
    errmsg = 'Correlation ratio not currently supported as similarity measure.';
elseif isequal(simMeasFlag,'MI')
    errid = 'npRegLib:genNpReg:checkOptions:MINotSupported';
    errmsg = 'Mutual information not currently supported as similarity measure.';
elseif isequal(simMeasFlag,'NMI')
    errid = 'npRegLib:genNpReg:checkOptions:NMINotSupported';
    errmsg = 'Normalized mututal information not currently supported as similarity measure.';
else
    errid = '';
    errmsg = '';
end

% report error
error(errid,errmsg);

author: tomdoel

repo: pulmonarytoolkit


text:
function garchplot(data, ht, resids)
%{
-----------------------------------------------------------------------
 PURPOSE:
 Plots the data, volatility and residuals
-----------------------------------------------------------------------
 USAGE:
 garch(data, ht, resids)

 INPUTS:
 data:      (T x 1) vector of data
 ht:        (T x 1) vector of conditional variance
 resids:    (T x 1) vector of residuals

-----------------------------------------------------------------------
Author:
 Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
 Date: 08/2011
-----------------------------------------------------------------------
%}
if nargin == 0
    error('Data, Conditional Variance, Residuals')
end

if size(data,2) > 1 | size(resids,2) > 1 | size(resids,2) > 2
     error('Data, variance and residual vectors should be column vectors')
end

figure
subplot(3,1,1), plot(data);
title('Returns');
subplot(3,1,2), plot(resids);
title('Residuals');
subplot(3,1,3), plot( ht);
title('Conditional Variance');

end

author: Sable

repo: mcbench-benchmarks


text:
function [tri, x] = remove_vertex_from_tri(tri, x, idx)
% REMOVE_VERTEX_FROM_TRI  Remove vertices from a triangular mesh
%
% This function removes vertices from a mesh and all triangles they are
% part of. It also relabels vertex indices accordingly.
%
% [TRI2, X2] = remove_vertex_from_tri(TRI, X, IDX)
%
%   TRI is a 3-column matrix. Each row represents the indices of the three
%   vertices that form a triangle. TRI as a whole represents the closed
%   surface.
%
%   X is a 3-column matrix. Each row represents the Cartesian coordinates
%   of a vertex on the surface, indexed by TRI values.
%
%   IDX is a vector of vertex indices.
%
%   TRI2 and X2 are the matrices that describe the mesh after the vertices
%   have been removed.
%
% See also: addpoint2tri.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(3, 3);
nargoutchk(0, 2);

% map between node indices map(10) = 6 means that vertex 10 will be renamed
% as 6
map = 1:size(x, 1);

% remove duplicates
idx = unique(idx);

% number of indices
N = length(idx);

% loop every index to remove
for I = 1:N

    % if we remove one vertex, all vertices after it decrease their index
    map(idx(I):end) = map(idx(I):end) - 1;

end

% find triangles connected to the vertex that are going to be removed
[I, J] = find(ismember(tri, idx));

% remove triangles
tri(I, :) = [];

% remove vertices
x(idx, :) = [];

% rename vertices in triangles
tri = map(tri);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% Functions for data management.
%
% Files
%   inInterval     - Check apartenance to interval.
%   inSquare       - Points inside square.
%   split          - Split vectors into scalars, or matrices into row vectors.
%   closeLogFile   - Close log file.
%   diginto        - Pretty print entire structure trees.
%   neesAnalysis   - NEES analysis for slamtb.
%   neesPlots      - Plot the results of NEES analysis for slamtb.
%   openLogFile    - Open log file.
%   writeLogData   - Write log data in file.
%   writeLogHeader - Write log file header.
%   errorAnalysis  - Error analysis for slamtb.
%   rmseAnalysis   - RMSE analysis for slamtb.
%   rmsePlots      - Plot the results of RMSE analysis for slamtb.

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% Test split for gesture labels, for video chunks data
% run in directory of dtfv_learning
addpath('other');
addpath('../../../toolbox/libsvm/matlab');
load('../dataset_trial/saved_anno.mat');
load('../dataset_trial/frame_rates.mat');

% Build feature Map of easy load of FV feats
FV_path = '../dataset_trial/FVs_chunks';
CNN_path = '../dataset_trial/CNNfeats';

FVfeaMap = FV_feaMap(FV_path);
CNNfeaMap = CNN_feaMap(CNN_path);

gesture_score = cell(1,10);
gesture_acc = cell(1,10);
gesture_pred = cell(1,10);
train_lab = cell(1, 10);
test_lab = cell(1, 10);
gesture_scorei_train = cell(1,10);
num_chunks_train = cell(1,10);
num_chunks_test = cell(1,10);

num_ges = 5;

use_CNN = 1;
for cv = 0:9
    gesture_score{cv+1} = cell(1, num_ges);
    gesture_acc{cv+1} = cell(1, num_ges);
    gesture_pred{cv+1} = cell(1, num_ges);
    train_lab{cv+1} = cell(1, num_ges);
    test_lab{cv+1} = cell(1, num_ges);

    total_posnum_tr = zeros(1,5); % number of total pos for training
    total_negnum_tr = zeros(1,5); % number of total neg for training
    train_au_lab = []; % train facial unit labels
    test_au_lab = [];
    current_i = 1;  % used for accumulating id of each sample
    fprintf('CV: %d\n', cv);

    %% Loading data
    % Train data
    trainfile = ['../Scripts_by_ID/trainVideo',num2str(cv),'.txt'];
    testfile = ['../Scripts_by_ID/testVideo',num2str(cv),'.txt'];
    fid = fopen(trainfile);
    C = textscan(fid, '%s');
    fclose(fid);
    num_v = length(C{1})/2;
    feat_train_CNN = [];
    feat_train_FV = [];
    num_chunks_train{cv+1} = zeros(num_v,1);
    for i = 1:num_v
        [pathstr,name,ext] = fileparts(C{1}{2*i-1});
        disp(name)
        chunk_names = get_chunkname(name);
        num_chunks_train{cv+1}(i) = length(chunk_names);
        for chunk_id = 1: length(chunk_names)
            % check whether is a bad chunk
            if ~bad_chunk(chunk_names{chunk_id})
                eval(sprintf('chunk_label=fea_dict.%s;',chunk_names{chunk_id}));
                chunk_name = chunk_names{chunk_id};
                chunk_label = str2num(chunk_label);
                train_au_lab = [train_au_lab; chunk_label'];  % #sample x #AU (5)
                % get features CNN
                CNN_feats = get_CNNfeats(chunk_name, CNNfeaMap, frame_rates);
                %CNN_feats = mean(CNN_feats, 1);  % 1-by-4096
                % get features FV
                FV_feats = get_FVfeats(chunk_name, FVfeaMap);
                % update step
                feat_train_CNN(current_i,:) = CNN_feats;
                feat_train_FV(current_i,:) = FV_feats;
                current_i = current_i +1;
            end
        end
    end
    disp('train#')
    sum(train_au_lab, 1)

    % Test data
    fid = fopen(testfile);
    C = textscan(fid, '%s');
    fclose(fid);
    num_v = length(C{1})/2;
    feat_test_CNN = [];
    feat_test_FV = [];
    current_i = 1;
    num_chunks_test{cv+1} = zeros(num_v,1);
    for i = 1:num_v
        [pathstr,name,ext] = fileparts(C{1}{2*i-1});
        chunk_names = get_chunkname(name);
        num_chunks_test{cv+1}(i) = length(chunk_names);
        for chunk_id = 1: length(chunk_names)
            % check whether is a bad chunk
            if ~bad_chunk(chunk_names{chunk_id})
                eval(sprintf('chunk_label=fea_dict.%s;',chunk_names{chunk_id}));
                chunk_name = chunk_names{chunk_id};
                chunk_label = str2num(chunk_label);
                test_au_lab = [test_au_lab; chunk_label'];  % #sample x #AU (5)
                % get features CNN
                CNN_feats = get_CNNfeats(chunk_name, CNNfeaMap, frame_rates);
                %CNN_feats = mean(CNN_feats, 1);  % 1-by-4096
                % get features FV
                FV_feats = get_FVfeats(chunk_name, FVfeaMap);
                % update step
                feat_test_CNN(current_i,:) = CNN_feats;
                feat_test_FV(current_i,:) = FV_feats;
                current_i = current_i +1;
            end
        end
    end
    disp('test#')
    sum(test_au_lab, 1)

    %% Training classifiers
    for ges_i = 1:num_ges
        fprintf('CV: %d, ges_id: %d\n', cv, ges_i);
        lab_vec = train_au_lab(:, ges_i);
        if use_CNN == 1
            fea_mat = feat_train_CNN;
        else
            fea_mat = feat_train_FV;
        end
        train_lab{cv+1}{ges_i} = lab_vec;
        model = svmtrain(lab_vec, fea_mat, '-t 0 -q -b 1');
        [pred, acc, prob] = svmpredict(lab_vec,fea_mat, model, '-b 1');
        gesture_score_train{cv+1}{ges_i} = prob(:,2);
        % Test
        lab_vec = test_au_lab(:, ges_i);
        if use_CNN == 1
            fea_mat = feat_test_CNN;
        else
            fea_mat = feat_test_FV;
        end
        [pred, acc, prob] = svmpredict(lab_vec,fea_mat, model, '-b 1');
        pred
        gesture_score{cv+1}{ges_i} = prob(:,2);
        gesture_acc{cv+1}{ges_i} = acc;
        gesture_pred{cv+1}{ges_i} = pred;
        test_lab{cv+1}{ges_i} = lab_vec;
    end
    cv_gesture_score = gesture_score{cv+1};
    %save(['CNN_gesture_score',num2str(cv),'.mat'],'cv_gesture_score');

end
if use_CNN == 1
    appname = 'CNN';
else
    appname = 'FV';
end
save(['total_',appname,'_gesture_score.mat'], 'gesture_score');
save(['total_',appname,'_gesture_pred.mat'], 'gesture_pred');
save(['total_',appname,'ges_results.mat'], 'gesture_score', 'gesture_pred', 'train_lab', 'test_lab','gesture_score_train','num_chunks_train','num_chunks_test');

author: Doubaibai

repo: DARE

save_path: github-repos/MATLAB/Doubaibai-DARE


text:
function ccfi_1_tabulate ( )

%*****************************************************************************80
%
%% CCFI_1_TABULATE tabulates CCFI_1 quadrature rules for the Laguerre integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CCFI_1_TABULATE\n' );
  fprintf ( 1, '  Tabulate CCFI_1 quadrature rules for the Laguerre integral.\n' );
  fprintf ( 1, '  Density function rho(x) = 1.\n' );
  fprintf ( 1, '  Region: 0 <= x < +oo.\n' );
  fprintf ( 1, '  Exactness: NONE.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Print the first 5 rules:\n' );
  fprintf ( 1, '   I          X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  ell = 1.0;

  for n = 1 : 5

    [ x, w ] = ccfi_1 ( n, ell );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %14.6f  %14.6f\n', i, x(i), w(i) );
    end
    fprintf ( 1, ' Sum                  %14.6f\n', sum ( w(1:n) ) );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Print the first 5 nested rules:\n' );
  fprintf ( 1, '   I          X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  ell = 1.0;

  for j = 1 : 5

    n = 2^j - 1;

    [ x, w ] = ccfi_1 ( n, ell );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %14.6f  %14.6f\n', i, x(i), w(i) );
    end
    fprintf ( 1, ' Sum                  %14.6f\n', sum ( w(1:n) ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% RMSE - Evaluates root mean squared error.
%
% E = RMSE(Y)
%
% E = RMSE(Y1,Y2)

% Last modified 2010-10-21
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function e = rmse(Y, Yh)

if nargin < 2
  e = sqrt(Y(:)'*Y(:)/numel(Y));
else
  z = Y(:)-Yh(:);
  e = sqrt(z'*z/numel(z));
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function h = compute_histogram_rbf(f, sigma, x, options)

% compute_histogram_rbf - parzen windows density estimation
%
%   h = compute_histogram_rbf(f, sigma, x);
%
%   f is the signal, h is an estimate of the histogram,
%   where h(i) is the density of the estimation around value x(i).
%
%   sigma is the bandwidth used to estimate the histogram (approx. size of
%       the bins).
%
%   If f is (n,2) matrix, then a joint histogram is estimated.
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;

if size(f,1)<size(f,2)
    f = f';
end
if size(x,1)<size(x,2)
    x = x';
end

nb_samples = Inf;
if isfield(options, 'nb_samples')
    nb_samples = options.nb_samples;
end
nb_samples = min(nb_samples,size(f,1));
sel = randperm(size(f,1)); sel = sel(1:nb_samples);
f = f(sel,:);

if size(f,2)==2
    h = compute_histogram_rbf_2d(f, sigma, x, options);
    return;
end

if isfield(options, 'renormalize')
    renormalize = options.renormalize;
else
    renormalize = 0;
end

n = length(x);
p = length(f);  % number of samples
f = f(:); x = x(:);

% h is of size p x n
h = repmat(f, [1 n]) - repmat( x(:)', [p 1] );
h = exp( -h.^2/(2*sigma^2) );
if renormalize
    d = repmat( sum(h,2), [1 n] );
    d(d<eps) = 1;
end
h = h./d;
h = sum(h,1);
h = h(:)/sum(h(:));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = compute_histogram_rbf_2d(f, sigma, x, options)

options.null = 0;

if isfield(options, 'nb_max')
    nb_max = options.nb_max;
else
    nb_max = 0;
end

if isfield(options, 'renormalize')
    renormalize = options.renormalize;
else
    renormalize = 0;
end

n = size(x,1);
p = size(f,1);

if p>nb_max
    h = zeros(n);
    niter = ceil(p/nb_max);
    for i=1:niter
        progressbar(i,niter);
        sel = (i-1)*nb_max+1:min(i*nb_max,p);
        h = h + compute_histogram_rbf_2d( f(sel,:) , sigma, x, options);
    end
    h = h/sum(h(:));
    return;
end

% hx is of size p x n x n
hx = repmat( f(:,1), [1 n n]) - repmat( reshape(x(:,1),[1 n 1]),[p 1 n] );
hy = repmat( f(:,2), [1 n n]) - repmat( reshape(x(:,2),[1 1 n]),[p n 1] );
h = exp( -(hx.^2/(2*sigma(1)^2)+hy.^2/(2*sigma(2)^2)) );
if renormalize
    d = repmat( sum(sum(h,2),3), [1 n n] );
    d(d<eps) = 1;
    h = h./d;
end
h = sum(h,1);
h = reshape(h,n,n);
h = h./sum(sum(h));

author: gpeyre

repo: matlab-toolboxes


text:
% Fig. 9.30   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

h=0.1;
N=1;
nn=1;
nb=nn*nn+1;
ai=0.1:0.001:nn;
NA=(4*N./(pi*ai)).*exp(-sqrt(-1)*asin(h./ai));
subplot(2,1,1)
plot(ai,abs(NA));
grid;
title('Describing function for hysteresis nonlinearity')
xlabel('a');
ylabel('Magnitude, |K_{eq}|');
pause;
ff=180/pi;
subplot(2,1,2)
plot(ai,ff*angle(NA));
grid;
xlabel('a')
ylabel('Phase, \angle K_{eq}, deg')
hold off


text:
% This file should NOT change when run through MBeautify
% If you find anything that is difficult or failed in some MBeautify
% version, please add it here.

% unary operator testcases

+2
+2.
+.2
'   string  with  lots  of  spaces   '
1 + 2
f(+1) + 1

x = y
x + 1. + 2
x + 1. + +.1
x + 1 + 2
x = 1
x = -1
x = +1
x = +.1
+(-[-.1])
z = [1, 2, 3, 4]

if 1 > +2
    return
end; % comment +-+-+- +++ 123 ***
if 1 > -2
end
% different meanings of 'end'
if any(z == -[-1, -2, -3, -4])
    ifmyvariablenamecontainsif = z(1:end);
end

% old-style function calls
disp +end+ this is not any keyword if else endif while +1
% bracket handling
while (1)
    a = [0, 1];
    a(1) = 2 * [a(0)];
    break
end;

% transpose
-x' + +1 + x'' + 2 * x''' * 1
a = eye(27)
a(3, 4:5) = [3, 1]

% norm notation
1.e-6

% #36
a(1, 1:2) = [3, 1]
a(3, :) = [3, 1]
b = zeros(3, 3, 3);
b(:, :, 2) = rand(2, 2);

 %if AddCommasToMatrices=1
a=[@(x) minus(x,1),]

%if AddCommasToCellArrays=1
a={@(x) minus(x,1),  @(x,y) minus(x,y)}

% #34
if -1 > -2
end
if +1 > -2
end

% #59
a = [1, 2, 3];
a = [1, a...
.* 2]

% #58
a = [1, a ... % 111
    .* 2, ... % 222
    123, 4- ...
    5] % 333

a = {'', ... % hello
    1, 4+ ... %hello2
     2} % hello3

if true || ... % aaa
        true || ... asd
        false % //
end

%80
a + sprintf("%d", b)%comment
a + sprintf("'%d'", b) %comment
a + sprintf("""%d""", b) % comment
a + sprintf('"%d"', b)
a.' + sprintf('''%d''', b)
a' + sprintf('%d', b)

% Remove extra space after @
f = @ (x) a
% Remove extra space before unary operator
f = @(x) - a
num@MySuper(obj) - a
% remove spaces around @
num @ MySuper(obj) - a

% if AddCommasToMatrices=1
% add comma after b
[a, - b c + d]
% add comma before b
[a -b]
% one comma added after b
[a *b (c -d)]

% treat whitespace as delimiter regardless of AddCommasToMatrices
[1 (2) {3}]
% same for cells
{1 (2) {3}}

author: davidvarga

repo: MBeautifier


text:
function t=gpst2time(week,sec)

gpst0=[1980,1, 6,0,0,0]; %gps time reference

t0=epoch2time(gpst0);

if sec<-1e9||sec>1e9,sec=0;end

t.time=t0.time+week*7*86400+fix(sec);
t.sec=sec-fix(sec);

return;

repo: GINav


text:
function [ outputTexture ] = Remap( inputTexture, mapX, mapY )
%REMAP Summary of this function goes here
%   Detailed explanation goes here

    outputTexture = zeros(size(mapX));

    [X,Y] = meshgrid(0:size(inputTexture,2)-1,0:size(inputTexture,1)-1);

    inds = find(mapX ~= -1);

    xSources = mapX(inds);
    ySources = mapY(inds);

    Z = interp2(X, Y, double(inputTexture), xSources, ySources);
    outputTexture(inds) = Z;
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function himmelblau_test ( )

%*****************************************************************************80
%
%% HIMMELBLAU_TEST works with the Himmelblau function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HIMMELBLAU_TEST:\n' );
  fprintf ( 1, '  Test COMPASS_SEARCH with the Himmelblau function.\n' );
  m = 2;
  delta_tol = 0.00001;
  delta = 0.3;
  k_max = 20000;

  x = [ 1.0, 1.0 ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', himmelblau ( m, x ) );

  [ x, fx, k ] = compass_search ( @himmelblau, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );

  x = [ -1.0, 1.0 ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', himmelblau ( m, x ) );

  [ x, fx, k ] = compass_search ( @himmelblau, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );

  x = [ -1.0, -1.0 ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', himmelblau ( m, x ) );

  [ x, fx, k ] = compass_search ( @himmelblau, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );

  x = [ 1.0, -1.0 ];
  r8vec_print ( m, x, '  Initial point X0:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X0) = %g\n', himmelblau ( m, x ) );

  [ x, fx, k ] = compass_search ( @himmelblau, m, x, delta_tol, delta, k_max );
  r8vec_print ( m, x, '  Estimated minimizer X1:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X1) = %g, number of steps = %d\n', fx, k );
%
%  Demonstrate Himmelblau minimizers.
%
  x = [ 3.0, 2.0 ];
  r8vec_print ( m, x, '  Correct Himmelblau minimizer X1*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', himmelblau ( m, x ) );

  x = [ 3.58439, -1.84813 ];
  r8vec_print ( m, x, '  Correct Himmelblau minimizer X2*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', himmelblau ( m, x ) );

  x = [ -3.77934, -3.28317 ];
  r8vec_print ( m, x, '  Correct Himmelblau minimizer X3*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', himmelblau ( m, x ) );

  x = [ -2.80512,  3.13134 ];
  r8vec_print ( m, x, '  Correct Himmelblau minimizer X4*:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X*) = %g\n', himmelblau ( m, x ) );

  return
end


text:
function A = transformMatrix (varargin)
%
%  CONSTRUCTOR FOR transfromMatrix OBJECT
%
%  The transformMatrix class is based on a structure with the
%  two fields:
%    transform  - character string indicating the transform:
%                 'fft'  or 'dct'
%    transpose  - indicates if the matrix has been transposed.
%
%   Calling Syntax:
%       A = transformMatrix           (returns object with empty fields)
%       A = transformMatrix(transformMatrixObj)
%       A = transformMatrix(transform)
%
%    where
%       * transformMatrixObj is an already existing transformMatrix object
%       * transform is a character string indicating the transform:
%         'fft' or 'dct'
%

%  J. Nagy  6/2/02

switch nargin
case 0
% default
   A.transform = '';
   A.transpose = 0;
   A = class(A, 'transformMatrix');
case 1
%
   if isa (varargin{1}, 'kronMatrix')
      A = varargin{1};
   elseif ischar(varargin{1})
      A.transform = varargin{1};
      A.transpose = 0;
      A = class(A, 'transformMatrix');
   end
otherwise
   error ('Incorrect number of arguments');
end

repo: IRtools