text:
classdef TopOpt_Problem < handle

    properties (GetAccess = public, SetAccess = public)
        designVariable
        dualVariable
        cost
        constraint
        optimizer
        incrementalScheme
        optimizerSettings
    end

    properties (Access = private)
        mesh
        videoMaker
        homogenizedVarComputer
    end

    methods (Access = public)

        function obj = TopOpt_Problem(cParams)
            obj.createMesh(cParams);
            obj.createIncrementalScheme(cParams);
            obj.createDesignVariable(cParams);
            obj.createHomogenizedVarComputer(cParams)
            obj.createCostAndConstraint(cParams);
            obj.createDualVariable();
            obj.createOptimizer(cParams);
            obj.createVideoMaker(cParams);
        end

        function createOptimizer(obj,settings)
            obj.completeOptimizerSettings(settings);
            obj.computeBounds();
            obj.optimizerSettings.outputFunction.type        = 'Topology';
            obj.optimizerSettings.outputFunction.iterDisplay = 'none';
            obj.optimizerSettings.outputFunction.monitoring  = MonitoringManager(obj.optimizerSettings);
            obj.optimizer = Optimizer.create(obj.optimizerSettings);
        end

        function computeBounds(obj)
            switch obj.designVariable.type
                case 'Density'
                obj.optimizerSettings.ub = 1;
                obj.optimizerSettings.lb = 0;
                otherwise

            end
        end

        function completeOptimizerSettings(obj,cParams)
            s = cParams.optimizerSettings;
            s.uncOptimizerSettings.scalarProductSettings = obj.designVariable.scalarProduct;

            s.uncOptimizerSettings.targetParameters = obj.incrementalScheme.targetParams;
            s.uncOptimizerSettings.designVariable   = obj.designVariable;

            s.monitoringDockerSettings.mesh = obj.mesh;

            s.designVar         = obj.designVariable;
            s.targetParameters  = obj.incrementalScheme.targetParams;
            s.cost              = obj.cost;
            s.constraint        = obj.constraint;
            s.incrementalScheme = obj.incrementalScheme;
            s.dualVariable      = obj.dualVariable;

            obj.optimizerSettings = s;
        end

        function computeVariables(obj)
            while obj.incrementalScheme.hasNext()
                obj.incrementalScheme.next();
                obj.optimizer.solveProblem();
            end
%             obj.optimizer.saveMonitoring();
%             obj.optimizer.simulationPrinter.print();
        end

        function postProcess(obj)
            iter = 0:obj.optimizer.nIter;
            obj.videoMaker.iterations = iter;
            obj.videoMaker.makeVideo();
        end

    end

    methods (Access = private)

        function optSet = obtainOptimizersSettings(obj,settings)
            epsilon = obj.incrementalScheme.targetParams.epsilon;
            settings.optimizerSettings.uncOptimizerSettings.lineSearchSettings.epsilon = epsilon;
            settings.optimizerSettings.uncOptimizerSettings.scalarProductSettings.epsilon = epsilon;
            set = settings.clone();
            optSet = set.optimizerSettings;
        end

        function createIncrementalScheme(obj,cParams)
            s = cParams.incrementalSchemeSettings;
            s.mesh = obj.mesh;
        %    s.targetParamsSettings.epsilonPerInitial = 10*s.targetParamsSettings.epsilonPerFinal;
            obj.incrementalScheme = IncrementalScheme(s);
        end

        function createDesignVariable(obj,cParams)
            s = cParams.designVarSettings;
            s.mesh = obj.mesh;
            s.scalarProductSettings.epsilon = obj.incrementalScheme.targetParams.epsilon;
            s.scalarProductSettings.mesh    = obj.mesh;
            obj.designVariable = DesignVariable.create(s);
        end

        function createDualVariable(obj)
            s.nConstraints = obj.constraint.nSF;
            obj.dualVariable = DualVariable(s);
        end

        function createHomogenizedVarComputer(obj,cParams)
            s = cParams.homogenizedVarComputerSettings;
            obj.homogenizedVarComputer = HomogenizedVarComputer.create(s);
        end

        function createCostAndConstraint(obj,cParams)
            obj.createCost(cParams);
            obj.createConstraint(cParams);
        end

        function createCost(obj,cParams)
            s = cParams.costSettings;
            s.designVar = obj.designVariable;
            s.homogenizedVarComputer = obj.homogenizedVarComputer;
            s.targetParameters = obj.incrementalScheme.targetParams;
            obj.cost = Cost(s);
        end

        function createConstraint(obj,cParams)
            s = cParams.constraintSettings;
            s.designVar = obj.designVariable;
            s.homogenizedVarComputer = obj.homogenizedVarComputer;
            s.targetParameters = obj.incrementalScheme.targetParams;
            s.dualVariable = obj.dualVariable;
            obj.constraint = Constraint(s);
        end

        function createVideoMaker(obj,cParams)
            s = cParams.videoMakerSettings;
            obj.videoMaker = VideoMaker.create(s);
        end

        function createMesh(obj,cParams)
            s = cParams.designVarSettings;

            a.coord  = s.femData.mesh.coord;
            a.connec = s.femData.mesh.connec;
            innerMesh = Mesh(a);
            boundMesh = innerMesh.createBoundaryMesh();
            sM.backgroundMesh = innerMesh;
            sM.boundaryMesh   = boundMesh;
%             obj.mesh = UnfittedMesh(sM);
            obj.mesh = innerMesh;
%             sM.coord  = s.femData.mesh.coord;
%             sM.connec = s.femData.mesh.connec;
%             obj.mesh = Mesh_Total(sM);
%             obj.mesh.innerMeshOLD.setMasterSlaveNodes(s.femData.mesh.masterSlaveNodes);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [err,tree] = l2boost(train_feats, train_y, val_feats, val_y, depth, max_iter, alpha)
% Run L2-Boosting and test parameters
    DefaultVal('max_iter',100);
    residual = train_y;
    err_increase = 1;
    pred_train = zeros(numel(train_y),1);
    pred_val = zeros(numel(val_y),1);
    for i=1:max_iter
        tree{i} = classregtree_fuxin(train_feats, residual, 'minleaf',numel(train_y) / 2^(depth-1));
        if exist('alpha','var')
            tree{i} = tree{i}.prune('alpha',alpha);
        end
        pred_train = pred_train + tree{i}.eval(train_feats);
        residual = train_y - pred_train;
        pred_val = pred_val + tree{i}.eval(val_feats);
        residual_val = val_y - pred_val;
        err(i) = mean(abs(residual_val))
        if i > 1 && err(i) >= err(i-1)
            err_increase = err_increase + 1;
        else
            err_increase = 0;
        end
        if err_increase > 3
            break;
        end
    end
    [a,b] = min(err(1:i));
    tree = tree(1:b);
end

author: Cloud-CV

repo: object-proposals


text:
function title = p11_title ( )

%*****************************************************************************80
%
%% P11_TITLE returns the title for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = '1 / ( (1+x) * sqrt(x) )';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function centroids = computeCentroids(X, idx, K)
%COMPUTECENTROIDS returs the new centroids by computing the means of the
%data points assigned to each centroid.
%   centroids = COMPUTECENTROIDS(X, idx, K) returns the new centroids by
%   computing the means of the data points assigned to each centroid. It is
%   given a dataset X where each row is a single data point, a vector
%   idx of centroid assignments (i.e. each entry in range [1..K]) for each
%   example, and K, the number of centroids. You should return a matrix
%   centroids, where each row of centroids is the mean of the data points
%   assigned to it.
%

% Useful variables
[m n] = size(X);

% You need to return the following variables correctly.
centroids = zeros(K, n);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every centroid and compute mean of all points that
%               belong to it. Concretely, the row vector centroids(i, :)
%               should contain the mean of the data points assigned to
%               centroid i.
%
% Note: You can use a for-loop over the centroids to compute this.
%

numberOfElementsHavingCentroid_k = zeros(K,1);
sumOfElementsHavingCentroid_k = zeros(K,n);
for i = 1:size(idx,1)
	z = idx(i);
	numberOfElementsHavingCentroid_k(z) += 1;
	sumOfElementsHavingCentroid_k(z,:) += X(i,:);
end

centroids = sumOfElementsHavingCentroid_k./numberOfElementsHavingCentroid_k;

% =============================================================

end

author: AvaisP


text:
caffe.reset_all();
caffe.set_mode_gpu();
gpu_id = 0;  % we will use the first gpu in this demo
caffe.set_device(gpu_id);

batch_size = 50;
feature_dim = 512;
mean_value = 128;
scale = 0.0078125;
ROIx = 1:96;
ROIy = 1:112;
height = length(ROIx);
width = length(ROIy);

feature_suffix = '_hard64nofd';
aligned_image_folder = 'F:\datasets\MegaFace\FaceScrub\aligned\';
feature_folder = ['F:\datasets\MegaFace\FaceScrub\feature' feature_suffix '\'];
if exist(feature_folder, 'dir')==0
    mkdir(feature_folder);
end;

test_list_file = 'E:\datasets\MegaFace\devkit\templatelists\facescrub_uncropped_features_list.json';
json_string = fileread(test_list_file);
json_string = json_string(strfind(json_string,'path')+8:end);
test_list = regexp(json_string(8:end), '"(.*?)"','tokens');
test_list_raw = cell(length(test_list),1);
for i=1:length(test_list)
    test_list_raw{i} = [aligned_image_folder test_list{i}{1}];
    dot_pos = strfind(test_list{i}{1},'.');
    if isempty(dot_pos)||(dot_pos(end) ~= length(test_list{i}{1}) - 3&&dot_pos(end) ~= length(test_list{i}{1}) - 4)
%         disp(test_list{i}{1});
        test_list{i}{1} = [test_list{i}{1} '.jpg'];
    end;
    if ~isempty(dot_pos) && strcmp(test_list{i}{1}(dot_pos(end):end), '.gif')
        test_list{i}{1} = [test_list{i}{1}(1:end-3) 'jpg'];
    end;
    test_list{i} = [aligned_image_folder test_list{i}{1}];
end;

total_image = length(test_list);
total_iter = ceil(total_image / batch_size);

net = caffe.Net('D:\face project\experiment\96_112_l2_distance\face_deploy_64.prototxt','D:\face project\norm_face_model\hard-margin-64-nofd.caffemodel', 'test');%face_model_my

image_p = 1;
for i=1:total_iter
    fprintf('%d/%d\n',i, total_iter);
    J = zeros(height,width,3,batch_size,'single');
    feature_names = cell(batch_size,1);
    for j = 1 : batch_size
        if image_p <= total_image
            I = imread(test_list{image_p});
            I = permute(I,[2 1 3]);
            I = I(:,:,[3 2 1]);
            I = I(ROIx,ROIy,:);
            I = single(I) - mean_value;
            J(:,:,:,j) = I*scale;
            feature_names{j} = [strrep(test_list_raw{image_p},aligned_image_folder, feature_folder) feature_suffix '.bin'];
            image_p = image_p + 1;
        end;
    end;
    f1 = net.forward({J});
    feature = squeeze(f1{1});
    for j = 1 : batch_size
        if ~isempty(feature_names{j})
            [file_folder, file_name, file_ext] = fileparts(feature_names{j});
            if exist(file_folder,'dir')==0
                mkdir(file_folder);
            end;
            fp = fopen(feature_names{j},'wb');
            fwrite(fp, [feature_dim 1 4 5], 'int32');
            fwrite(fp, feature(:,j), 'float32');
            fclose(fp);
        end;
    end;
end;

author: happynear

repo: FaceVerification


text:
function value = r8vec_in_ab ( n, x, a, b )

%*****************************************************************************80
%
%% R8VEC_IN_AB is TRUE if the entries of an R8VEC are in the range [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real X(N), the vector.
%
%    Input, real A, B, the limits, with A <= B.
%
%    Output, logical VALUE, is TRUE if every entry is
%    between A and B.
%
  if ( any ( x(1:n) < a | b < x(1:n) ) )
    value = 0;
  else
    value = 1;
  end

  return
end


text:
function X=stft(x,wl,olr,window,alpha)
% STFT Short Time Fourier Transform
%
% X = stfth(x,wl,olr,window,alpha) returns the short time Fourier Transform of
% a signal "x" with a window length "wl", overlap ratio "olr" and window
% type "window". Valid window types are as follows:
%
%   'rectwin'   - Rectangular Window
%   'hann'      - Hann window
%   'hamming'   - Hamming window
%   'gausswin'   - Gaussian window
%
% If a Gaussian window is used, alpha is the reciprocal of the standard
% deivation of the Gaussian window (if used), as per the GAUSSWIN function
%
% Modified by Al-Hafeez Dhalla, PhD. Contact email: hafeez.dhalla@gmail.com
%
% Original version produced by Suraj Kamya, available at:
%
% www.mathworks.com/matlabcentral/fileexchange/38035-stft-short-time-fourier-transform

if nargin < 5
    alpha = 2.5;
end

L=length(x);
if L<wl
    z=wl-L;
    x=[x,zeros(1,z)];
end
switch window  % Window functions are (all =0 outside the interval 0<=n<=wl)
    case 'rectwin'
        win=ones(1,wl); % Rectangular Window
    case 'hann'
        win=hamming(wl)'; % Hann window w[n]=0.5-0.5*cos(2*pi*n/M), for 0<=n<=M
    case 'hamming'
        win=hann(wl)';  % Hamming window w[n]=0.54-.46*cos(2*pi*n/M) for 0<=n<=M
    case 'gausswin'
        win=gausswin(wl,3.5)';  % Gaussian window
    otherwise
        win=ones(1,window_length);
        disp('Invalid window type, defaulting to rectangular window.');
end

L=length(x);
hop=ceil(wl * (1-olr));  % Hop size of window

if hop<1
    fprintf('\n Warning: Overlap ratio too large. Using maximum overlap ratio of (wl-1)/wl = %f. \n',(wl-1)/wl);
    fprintf(' This may be very slow. \n')
    hop=1;
end

i=1;str=1; len=wl; X=[];
while (len<=L || i<2)
    if i==1
    if  len>L   % If window size excceds the L of signal for 1st time
        z=len-L;
        x=[x,zeros(1,z)]; % padding zeros
        i=1+1;
    end
    x1=x(str:len);
    X=[X;fft(x1.*win)];  % Matrix mul
    str=str+hop; len=str+wl-1; % to make window overlapping
        end
end

X = fftshift(X,2)';

author: Sable

repo: mcbench-benchmarks


text:
% importevent() - Import experimental events from data file or Matlab
%                 array into a structure.
%
% Usage: >> eventstruct = importevent( event, oldevent, srate);
%        >> eventstruct = importevent( event, oldevent, srate, 'key1', 'value1', ...);
%
% Input:
%   event     - [ 'filename'|array ] Filename of a text file, or name of
%               Matlab array in the global workspace containing an
%               array of events in the folowing format: The first column of
%               the cell array is the type of the event, the second the latency.
%               The others are user-defined. The function can read
%               either numeric or text entries in ascii files.
%   oldevent  - Old event structure. Used for aligning new events with old
%               ones. Enter [] is no such structure exist.
%   srate     - Sampling rate of the underlying data. Event latencies are
%               expressed in terms of latency in sample point with respect
%               to the data.
%
% Optional file or array input:
%  'fields'   - [Cell array] List of the name of each user-defined column,
%               optionally followed by a description. Ex: { 'type', 'latency' }
%  'skipline' - [Interger] Number of header rows to skip in the text file
%  'timeunit' - [ latency unit rel. to seconds ]. Default unit is 1 = seconds.
%               NaN indicates that the latencies are given in time points.
%  'delim'    - [string] String of delimiting characters in the input file.
%               Default is tab|space.
%
% Optional oldevent input:
%  'append'   - ['yes'|'no'] 'yes' = Append events to the current events in
%               the EEG dataset {default}: 'no' = Erase the previous events.
%  'indices'  - {integer vector] Vector indicating the indices of the events to
%               modify.
%  'align'    - [num] Align the first event latency to the latency of existing
%               event number (num), and check latency consistency.
%               A value of 0 indicates that the first events of the pre-defined
%               events and imported events will be aligned. A positive value (num)
%               aligns the first event to the num-th pre-existing event.
%               A negative value can also be used; then event number (-num)
%               is aligned to the first pre-existing event.  Default is 0.
%               (NaN-> no alignment).
%  'optimmeas' - ['median'|'mean'] Uses either the median of the mean
%               distance to align events. Default is 'mean'. Median is
%               preferable if events are missing in the event file.
%  'optimalign' - ['on'|'off'] Optimize the sampling rate of the new events so
%               they best align with old events. Default is 'on'.
%
% Outputs:
%   eventstruct - Event structure containing imported events
%
% Example: >> eventstruct = importevent( 'event_values.txt', [], 256, ...
%         'fields', {'type', 'latency','condition' }, ...
%         'append', 'no', 'align', 0, 'timeunit', 1E-3 );
%
%         This loads the ascii file 'event_values.txt' containing 3 columns
%         (event_type, latency, and condition). Latencies in the file are
%         in ms (1E-3). The first event latency is re-aligned with the
%         beginning of the dataset ('align', 0). Any previous events
%         are erased ('append', 'no').
%
% Author: Arnaud Delorme & Scott Makeig, CNL / Salk Institute, 2004
%
% See also: pop_importevent()

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 2004, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%% REVISION HISTORY
%% function definition, input check
function event = importevent(event, oldevent, srate, varargin)

if nargin < 1
   help importevent;
   return;
end;

I = [];

% remove the event field
% ----------------------
if ~isempty(oldevent), allfields = fieldnames(oldevent);
else                   allfields = {}; end;

g = finputcheck( varargin, { 'fields'    'cell'     []         {};
                         'skipline'  'integer'  [0 Inf]    0;
                         'indices'   'integer'  [1 Inf]    [];
                         'append'    'string'   {'yes';'no';'''yes''';'''no''' }         'yes';
                         'timeunit'  'real'     []         1;
                         'event'     { 'cell';'real';'string' }     []    [];
                         'align'     'integer'  []         NaN;
                         'optimalign' 'string'  { 'on';'off' }         'on';
                         'optimmeas'  'string'  { 'median';'mean' }         'mean';
                         'delim'     {'integer';'string'}   []         char([9 32 44])}, 'importevent');
if isstr(g), error(g); end;
if ~isempty(g.indices), g.append = 'yes'; end;
g.delim = char(g.delim);

% call from pop_importevent
% -------------------------
if ~isempty(g.event)
    event = g.event;
    if ~isstr(event)
        if size(event,2) > size(event,1), event = event'; end;
        if iscell(event)
            eventStr = cellfun(@ischar, event);
            eventNum = cellfun(@isnumeric, event);
            if all(eventStr(1,:)) && ~all(eventStr(:,1))
                event = event';
            end;
            if all(eventNum(1,:)) && ~all(eventNum(:,1))
                event = event';
            end;
        end;
    end;
end;
g.event = event;

% determine latency for old event alignment
% -----------------------------------------
tmpalign    = g.align;
g.align     = [];
g.align.val = tmpalign;
if ~isnan(g.align.val)
    if isempty(oldevent)
        error('Setevent: no pre-existing event, cannot perform alignment');
    end;
    if ~isfield(oldevent, 'latency')
        error('Setevent: pre-existing events do not have a latency field for re-alignment');
    end;
    switch g.append
        case {'yes' '''yes'''}, disp('Setevent warning: cannot align and append events at the same time; disabling event alignment');
    end;
    if g.align.val < 0
        g.align.event = oldevent(1).latency;
    else
        g.align.event = oldevent(g.align.val+1).latency;
    end
    g.align.nbevent = length(oldevent);
    g.oldevents = oldevent;
    g.align.txt = sprintf([ 'Check alignment between pre-existing (old) and loaded event' ...
                          ' latencies:\nOld event latencies (10 first): %s ...\n' ], ...
                          int2str([ oldevent(1:min(10, length(oldevent))).latency ]));
else
    g.oldevents = [];
end;

tmpfields = fieldnames(g);
event = oldevent;
%% scan all the fields of g
% ------------------------
for curfield = tmpfields'
    if ~isempty(event), allfields = fieldnames(event);
    else                    allfields = {}; end;
    switch lower(curfield{1})
        case {'append', 'fields', 'skipline', 'indices', 'timeunit', 'align', 'delim' }, ; % do nothing now
        case 'event', % load an ascii file
            switch g.append
                case { '''no''' 'no' } % ''no'' for backward compatibility
                      if isstr(g.event) && ~exist(g.event), g.event = evalin('caller', g.event); end;
                      event = load_file_or_array( g.event, g.skipline, g.delim );
                      allfields = g.fields(1:min(length(g.fields), size(event,2)));
                      event = eeg_eventformat(event, 'struct', allfields);
					  % generate ori fields
					  % -------------------
                      if ~isnan(g.timeunit)
                          for index = 1:length(event)
                              event(index).init_index = index;
                              event(index).init_time  = event(index).latency*g.timeunit;
                          end;
                      end;
                      event = recomputelatency( event, 1:length(event), srate, ...
                                                    g.timeunit, g.align, g.oldevents, g.optimalign, g.optimmeas);
                case { '''yes''' 'yes' }
                      % match existing fields
                      % ---------------------
                      if isstr(g.event) && ~exist(g.event), g.event = evalin('caller', g.event); end;
                      tmparray = load_file_or_array( g.event, g.skipline, g.delim );
                      if isempty(g.indices) g.indices = [1:size(tmparray,1)] + length(event); end;
                      if length(g.indices) ~= size(tmparray,1)
                            error('Set error: number of row in file does not match the number of event given as input');
                      end;

                      % add field
                      % ---------
                      g.fields = getnewfields( g.fields, size(tmparray,2)-length(g.fields));

                      % add new values
                      % ---------------------
                      for eventfield = 1:size(tmparray,2)
                          event = setstruct( event, g.fields{eventfield}, g.indices, { tmparray{:,eventfield} });
                      end;
					  % generate ori fields
					  % -------------------
                      offset = length(event)-size(tmparray,1);
                      for index = 1:size(tmparray,1)
                          event(index+offset).init_index = index;
                          event(index+offset).init_time  = event(index+offset).latency*g.timeunit;
                      end;
                      event = recomputelatency( event, g.indices, srate, g.timeunit, ...
                                                    g.align, g.oldevents, g.optimalign, g.optimmeas);
            end;
      end;
end;

if isempty(event) % usefull 0xNB empty structure
    event = [];
end;

%% remove the events wit out-of-bound latencies
% --------------------------------------------
if isfield(event, 'latency')
    try
        res = cellfun('isempty', { event.latency });
        res = find(res);
        if ~isempty(res)
            fprintf( 'importevent warning: %d/%d event(s) have invalid latencies and were removed\n', ...
                     length(res), length(event));
            event( res ) = [];
        end;
    end;
end;

%% interpret the variable name
% ---------------------------
function array = load_file_or_array( varname, skipline, delim );
    if isstr(varname) & exist(varname) == 2  % mean that it is a filename
                                             % --------------------------
        array = loadtxt( varname, 'skipline', skipline, 'delim', delim, 'blankcell', 'off' );

    else
         if ~iscell(varname)
             array = mattocell(varname);
         else
             array = varname;
         end;
    end;
return;

%% update latency values
% ---------------------
function event = recomputelatency( event, indices, srate, timeunit, align, oldevents, optimalign, optimmeas)

    % update time unit
    % ----------------
    if ~isfield(event, 'latency'),
        if isfield(event, 'duration')
            error('A duration field cannot be defined if a latency field has not been defined');
        end;
        return;
    end;
    if ~isnan(timeunit)
        for index = indices
            event(index).latency  = event(index).latency*srate*timeunit;
            if isfield(event, 'duration')
                event(index).duration = event(index).duration*srate*timeunit;
            end;
        end;
    end;

    % alignment with old events
    % -------------------------
    if ~isnan( align.val )
        if align.val >= 0, alignlatency = event(1).latency;
        else               alignlatency = event(-align.val+1).latency;
        end;
        for index = indices
             event(index).latency = event(index).latency-alignlatency+align.event;
        end;
        if length(event) ~= align.nbevent
            disp('Setevent warning: the number of pre-existing events do not correspond to the number of events');
            disp('                  that were read, so their latencies may have been wrongly re-aligned');
        end;
        fprintf(align.txt);
        fprintf('New event latencies (10 first): %s ...\n', int2str(round([ event(1:min(10, length(event))).latency ])));
    end;
    if strcmpi(optimalign, 'on') & ~isempty(oldevents)
        newlat = [ event.latency     ];
        oldlat = [ oldevents.latency ];

        newlat = repmat(newlat, [length(oldlat) 1]);
        oldlat = repmat(oldlat', [1 size(newlat,2)]);
        if align.val >= 0
            newlat = newlat-newlat(1);
            oldlat = oldlat-oldlat(1+align.val);
        else
            newlat = newlat-newlat(1-align.val);
            oldlat = oldlat-oldlat(1);
        end;

        try
            newfactor = fminsearch('eventalign',1,[],newlat, oldlat, optimmeas);
        catch
            newfactor = fminsearch('eventalign',1,[],[],newlat, oldlat, optimmeas); % Octave
        end;
        fprintf('Best sampling rate ratio found is %1.7f. Below latencies after adjustment\n', newfactor);
        if newfactor > 1.01 | newfactor < 0.99
            disp('Difference is more than 1%, something is wrong; ignoring ratio');
            newfactor = 1;
        else
            difference1 = eventalign( 1        , newlat, oldlat, optimmeas);
            difference2 = eventalign( newfactor, newlat, oldlat, optimmeas);
            fprintf('The average difference before correction was %f sample points\n', difference1);
            fprintf('The average difference after correction is %f sample points\n', difference2);
        end;

        %diffarray = abs(newfactor*newlat-oldlat)';
        %[allmins poss] = min(diffarray);
        %figure; hist(allmins);
    else
        newfactor = 1;
    end;
    if ~isnan( align.val ) & newfactor ~= 1
        if align.val >= 0
            latfirstevent = event(1).latency;
        else
            latfirstevent = event(-align.val+1).latency;
        end;
        for index = setdiff_bc(indices, 1)
            event(index).latency = round(event(index).latency-latfirstevent)*newfactor+latfirstevent;
        end;
        if ~isempty(oldevents)
            fprintf('Old event latencies (10 first): %s ...\n', int2str(round([ oldevents(1:min(10, length(oldevents))).latency ])));
            fprintf('New event latencies (10 first): %s ...\n', int2str(round([ event(1:min(10, length(event))).latency ])));
        end;
    else
        % must add one (because first sample point has latency 0
        % ------------------------------------------------------
        if ~isnan(timeunit)
            for index = indices
                event(index).latency = round((event(index).latency+1)*1000*newfactor)/1000;
            end;
        end;
    end;

%% create new field names
% ----------------------
function epochfield = getnewfields( epochfield, nbfields )
   count = 1;
   while nbfields > 0
       if isempty( strmatch([ 'var' int2str(count) ], epochfield ) )
               epochfield =  { epochfield{:} [ 'var' int2str(count) ] };
               nbfields = nbfields-1;
       else    count = count+1;
       end;
   end
return;

%%
% ----------------------
function var = setstruct( var, fieldname, indices, values )
    if exist('indices') ~= 1, indices = 1:length(var); end
    if ~isempty(values)
        for index = 1:length(indices)
            var = setfield(var, {indices(index)}, fieldname, values{index});
        end
    else
        for index = 1:length(indices)
            var = setfield(var, {indices(index)}, fieldname, '');
        end
    end
return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [slice] = spm_vb_taylor_R(Y,slice)
% Get Taylor series approximation to posterior correlation matrices
% FORMAT [slice] = spm_vb_taylor_R(Y,slice)
%
% Y        - data
% slice    - VB-GLMAR data structure
%
% See paper VB3.
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_vb_taylor_R.m 6079 2014-06-30 18:25:37Z spm $

% Get mean hyperparameter values
h0 = [];
if slice.p > 0
    if size(slice.ap_mean,2)==1
        % Single voxel in slice
        a      = slice.ap_mean';
        a_cov  = slice.a_cov{1};
    else
        a      = mean(slice.ap_mean');
        a_covs = cat(3,slice.a_cov{:});
        a_cov  = mean(a_covs,3);
    end
    slice.mean.a     = a;
    slice.mean.a_cov = a_cov;
    h0 = a';
end
slice.mean.b = mean(slice.b');

lambda = mean(slice.mean_lambda);
slice.mean.lambda = lambda;
h0 = [h0;lambda];

R = spm_vb_get_R(slice,h0);
slice.mean.R = R;

delta = 0.0001;
% Get first order Taylor terms about slice
% mean values of a and lambda
for i=1:length(h0)
    h    = h0;
    h(i) = h(i)-delta;

    [R1] = spm_vb_get_R(slice,h);

    h    = h0;
    h(i) = h(i) + delta;
    % Loop over hyperparameters
    R2   = spm_vb_get_R(slice,h);

    dR(:,:,i) = (R2-R1) / (2*delta);
end
slice.mean.h0 = h0;
slice.mean.dR = dR;

author: spm


text:
function [A,E,Y] = singular_value_rpca( D, lambda, tau, delta, svdMethod, A0)

%   Solves the Robust PCA relaxation
%
%      min  \tau ( |A|_* + \lambda |E|_1 ) + 1/2 |(A,E)|_F^2
%      subj  A+E = D
%
%   by iterative thresholding.
%
%   Inputs:
%      D      -- the data matrix, m x n.
%      lambda -- relative weight of sparsity of error
%
%   Optional:
%      tau    -- magnitude of L2 relaxation of the pure robust PCA SDP,
%                 higher value is desirable.
%      delta  -- stepsize, should be in (0,1).
%      svdMethod -- SVD routine to be used in each iteration, must be one
%                    of 'svdlibc', 'propack', or 'svds'. If not any of
%                    these, default option is MATLAB's svd command. (May
%                    require additional library dependencies if custom routine is used.)
%      A0 -- true low-rank solution, if known, to enable better display of
%            progress in each iteration.
%
%   Outputs:
%      A      -- estimate of the low-rank generating matrix
%      E      -- estimate of the error or corruption
%
%   Winter '08, John Wright, Shankar Rao. Questions? jnwright@uiuc.edu
%
%   Copyright: Perception and Decision Laboratory
%			University of Illinois, Urbana-Champaign

VERBOSE                         = 2;
EPSILON_PRIMAL                  = 5e-4;

if nargin < 5,  svdMethod = 'svd'; end

if nargin < 4,  delta = 0.9; end;

if nargin < 3, tau = 1e4; end;

MAX_ITER                        = 25000;
DISPLAY_EVERY                   = 100;

[m,n] = size(D);

Y = zeros(m,n);  % Lagrange multiplier
A = zeros(m,n);  % Structure
E = zeros(m,n);  % Error

rankA = 0;

iter      = 0;
converged = false;

while ~converged
    iter = iter + 1;

    switch lower(svdMethod)
        case 'svdlibc'
            [U,diagS,V] = svdlibc(Y, rankA+1);
        case 'propack'
            [U,S,V] = lansvd(Y,rankA+1,'L');
            diagS = diag(S);
        case 'svds'
            [U,S,V] = svds(Y, rankA+1, 'L');
            diagS = diag(S);
        otherwise
            [U,S,V] = svd(Y,0);
            diagS = diag(S);
    end

    A = U * diag(pos(diagS-tau)) * V';
    E = sign(Y) .* pos( abs(Y) - lambda*tau );
    M = D - A - E;

    rankA  = sum(diagS>tau);
    cardE = sum(sum(double(abs(E)>0)));

    Y = Y + delta * M;

%     if VERBOSE > 1 && mod(iter, DISPLAY_EVERY)==0 && nargin>=6,
%         disp(['    Iteration '    num2str(iter)               ...
%             ' |A|_F '         num2str(norm(A,'fro'))      ...
%             ' rank(A) '         num2str(rankA)              ...
%             ' |E|_F '         num2str(norm(E,'fro'))      ...
%             ' |E|_0 '         num2str(cardE) ...
%             ' |D-A-E|_F ' num2str(norm(M,'fro')) ...
%             ' |A-A0|_F / |A0|_F ' num2str(norm(A-A0,'fro')/norm(A0,'fro')) ...
%             ' |D-A-E|_1,inf ' num2str(max(max(abs(M)))) ]);
%     elseif VERBOSE > 0 && mod(iter, DISPLAY_EVERY)==0,
%         disp(['    Iteration '    num2str(iter)               ...
%             ' rank(A) '         num2str(rankA) ...
%             ' ||E||_0 ' num2str(cardE) ]);
%     end

    if ( norm(D-A-E,'fro')/norm(D,'fro') < EPSILON_PRIMAL || iter >= MAX_ITER )
        converged = true;
    end

%     if ( iter >= MAX_ITER )
%         disp('Maximum number of iterations reached.') ;
%     end

end

author: andrewssobral

repo: lrslibrary


text:
function X=image_show(img,grlvls,scale,title)
%Displays image, with additional parameters
%X = image_show(img,grlvls,scale,title)
%
%Input:
% img - image array; if the values fall outside [0,255] they will be
%       scaled
% grlvls - if it is a single value it represents the number of gray levels
%          to be displayed;
%          if grlvls is an array then it represents the colormap
% scale - scale for displaying the image
% title - name of the image figure, if it starts with 'save_' the image will be
%         saved instead. Don't forget the image format(e.g. 'save_test.png')!!
%
%Output:
% X - contains the image (useful if the image values have been scaled to
%     fit in the [0, 255] range
%
%Note:
% Works also for truecolor non-indexed images.
% If the range between the maximum and the minimum value in the image is larger
% than 2^16 then the image values will be rescaled logarithmically.
%
%Example:
% X = image_show(A,256,2,'A'); %zoom 2x

%for wavelet coefficients:
%Z2= 128+32*log2(abs(Y)+1); //high pass
%Z2(1:36,1:44) = Y(1:36,1:44)/8; //LL

LogTreshold = 10;

if isstr(img)
    img=imread(img);
end;
if numel(grlvls)>1
    mapa=grlvls;
else
    mapa=gray(grlvls);
end;
% maximg=max(img(:));
% minimg=min(img(:));
% if (maximg ~= minimg)
%     if maximg>255 | minimg<0
%         if log2(maximg-minimg)>LogTreshold
%             logimg=log2(abs(img)+1);
%             img=255*logimg/max(max(logimg));
%         else
%             img=round(255*(img-minimg)./(maximg-minimg));
%         end;
%     end;
%     if nargin<4
%         title=''
%         if nargin <3
%             scale=1;
%         end;
%     end;
% end;

set(0,'Units','pixels');
rd=get(0,'ScreenSize');
[rws,cls,ndms]=size(img);
rws=scale*rws;
cls=scale*cls;
fwd=round(cls+cls*0.1);
fhg=round(rws+rws*0.1);
if ~isempty(strmatch('save_',title))
    if ndims(img)==3
        imwrite(img,title);
    else
        imwrite(img,mapa,title);
    end;
else
    figure('Name',title,'Menubar','none','Units','pixels','Position',[(rd(3)-fwd)/2 (rd(4)-fhg)/2 fwd fhg]);
    image(img);
    if ndims(img)<3
        colormap(mapa);
    end;
    fh=gca;
    set(gca,'DataAspectRatio',[1 1 1],'TickLength',[0 0],'XTick',[],'YTick',[],'Units','pixels','Position',[cls*0.05 rws*0.05 cls rws]);
end;
X=img;


text:
function [mesh, edge_to_vertex, edge_to_face] = geodesic_new_mesh(points, tri)

global geodesic_library
if ~libisloaded(geodesic_library)
    hfile = 'geodesic_matlab_api.h';
    loadlibrary([geodesic_library '.so'], hfile);
end

dim = find(size(points) == 3);
if dim == 1
    p = points(:);
else
    p = points';
    p = p(:);
end;

dim = find(size(tri) == 3);
if dim == 1
    t = tri(:) - 1;
else
    t = tri';
    t = t(:) - 1;
end;

tmp = libpointer('doublePtrPtr');
[id, tmp1, tmp2, num_edges, edges] = calllib(geodesic_library, 'new_mesh', length(p)/3, p, length(t)/3, t, 1, tmp);
setdatatype(edges, 'doublePtr', 4, num_edges);

mesh.id = id;
mesh.object_type = 'mesh';
edge_to_vertex = (edges.Value(1:2,:) + 1)';
edge_to_face = (edges.Value(3:4,:) + 1)';

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [s, cfg] = ft_statfun_depsamplesFmultivariate(cfg, dat, design)

% FT_STATFUN_DEPSAMPLESFMULTIVARIATE calculates the MANOVA dependent samples
% F-statistic on the biological data in dat (the dependent variable), using the
% information on the independent variable (ivar) in design.
%
% Use this function by calling one of the high-level statistics functions as
%   [stat] = ft_timelockstatistics(cfg, timelock1, timelock2, ...)
%   [stat] = ft_freqstatistics(cfg, freq1, freq2, ...)
%   [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% with the following configuration option
%   cfg.statistic = 'ft_statfun_depsamplesFmultivariate'
%
% Configuration options
%   cfg.contrastcoefs  = matrix of contrast coefficients determining the
%                        effect being tested. The number of columns of this
%                        matrix has to be equal to the number of conditions.
%                        The default is a matrix that specifies the
%                        main effect of the independent variable. This matrix
%                        has size [(ncond-1),ncond].
%   cfg.computestat    = 'yes' or 'no', calculate the statistic (default='yes')
%   cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no')
%   cfg.computeprob    = 'yes' or 'no', calculate the p-values (default='no')
%
% The following options are relevant if cfg.computecritval='yes' and/or
% cfg.computeprob='yes'.
%   cfg.alpha = critical alpha-level of the statistical test (default=0.05)
%   cfg.tail  = -1, 0, or 1, left, two-sided, or right (default=1)
%               cfg.tail in combination with cfg.computecritval='yes'
%               determines whether the critical value is computed at
%               quantile cfg.alpha (with cfg.tail=-1), at quantiles
%               cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at
%               quantile (1-cfg.alpha) (with cfg.tail=1).
%
% Design specification
%   cfg.ivar  = row number of the design that contains the labels of the conditions that must be
%               compared (default=1). The labels range from 1 to the number of conditions.
%   cfg.uvar  = row number of design that contains the labels of the units-of-observation (subjects or trials)
%               (default=2). The labels are assumed to be integers ranging from 1 to
%               the number of units-of-observation.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS

% Copyright (C) 2006, Eric Maris
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% set defaults
if ~isfield(cfg, 'computestat'),       cfg.computestat='yes';     end
if ~isfield(cfg, 'computecritval'),    cfg.computecritval='no';   end
if ~isfield(cfg, 'computeprob'),       cfg.computeprob='no';      end
if ~isfield(cfg, 'alpha'),             cfg.alpha=0.05;            end
if ~isfield(cfg, 'tail'),              cfg.tail=1;                end

nconds=length(unique(design(cfg.ivar,:)));
if ~isfield(cfg,'contrastcoefs')
    % specify the default contrast coefficient matrix.
    ncontrasts = nconds-1;
    cfg.contrastcoefs = zeros(ncontrasts,nconds);
    cfg.contrastcoefs(:,1) = 1;
    for contrastindx=1:ncontrasts
        cfg.contrastcoefs(contrastindx,contrastindx+1)=-1;
    end
else
    ncontrasts = size(cfg.contrastcoefs,1);
end

% perform some checks on the configuration
if strcmp(cfg.computeprob,'yes') && strcmp(cfg.computestat,'no')
    ft_error('P-values can only be calculated if the test statistics are calculated.');
end
if ~isfield(cfg,'uvar') || isempty(cfg.uvar)
    ft_error('uvar must be specified for dependent samples statistics');
end

% perform some checks on the design
nuospercond=zeros(nconds,1);
for condindx=1:nconds
    nuospercond(condindx)=length(find(design(cfg.ivar,:)==condindx));
end
if sum(nuospercond)<size(design,2) || any(nuospercond~=nuospercond(1))
  ft_error('Invalid specification of the design array.');
end
nunits = max(design(cfg.uvar,:));
dfdenom = nunits - ncontrasts;
if dfdenom<1
    ft_error('The data must contain more units-of-observation (usually subjects) than the number of contrasts.')
end
nrepl=nunits*nconds;
if (nrepl~=sum(nuospercond)) || (nrepl~=size(dat,2))
  ft_error('Invalid specification of the design array.');
end
nsmpls = size(dat,1);

if strcmp(cfg.computestat,'yes')
% compute the statistic
    % store the positions of the condition labels nunits-by-nconds array
    poslabelsperunit=zeros(nunits,nconds);
    for condindx=1:nconds
        poslabel=find(design(cfg.ivar,:)==condindx);
        [dum,i]=sort(design(cfg.uvar,poslabel),'ascend');
        poslabelsperunit(:,condindx)=poslabel(i);
    end
    % reshape poslabelsperunit into a row vector that contains the
    % replications of the first condition on the first nunits positions,
    % the replications of the second condition on the second nunits
    % positions, etc.
    poslabelsperunit=reshape(poslabelsperunit,1,nrepl);
    s.stat=zeros(nsmpls,1);
    for smplindx=1:nsmpls
        datonesmpl=reshape(dat(smplindx,poslabelsperunit),nunits,nconds);
        contrasts=datonesmpl*cfg.contrastcoefs';
        contrastavg=mean(contrasts,1);
        dev=contrasts-repmat(contrastavg,nunits,1);
        covmat=(dev'*dev)/(nunits-1);
        s.stat(smplindx)=nunits*contrastavg*inv(covmat)*contrastavg';
    end
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  s.dfnum = ncontrasts;
  s.dfdenom = nunits - ncontrasts;
  if cfg.tail==-1
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a left tail critical value.');
  end
  if cfg.tail==0
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a two-sided critical value.');
  end
  if cfg.tail==1
    s.critval = ((nunits-1).*ncontrasts./(nunits-ncontrasts)).*finv(1-cfg.alpha,s.dfnum,s.dfdenom);
  end
end

if strcmp(cfg.computeprob,'yes')
  % also compute the p-values
  s.dfnum = ncontrasts;
  s.dfdenom = nunits - ncontrasts;
  if cfg.tail==-1
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a left tail p-value.');
  end
  if cfg.tail==0
      ft_error('For a dependent samples F-statistic, it does not make sense to calculate a two-sided p-value.');
  end
  if cfg.tail==1
    scaledstat = ((nunits-ncontrasts)./((nunits-1).*ncontrasts)).*s.stat;
    s.prob = 1-fcdf(scaledstat,s.dfnum,s.dfdenom);
  end
end


text:
clear;
% some integration experiments

% Example 1: E[p(x)], x~N(0,1), p(x) polynomial
degree = 10;
p = rand(1,degree+1); % make some polynomial
fprintf('adaptive quad, tol %.1g = %.10g\n', ...
        1e-10, quad(@(x) polyval(p,x).*normpdf(x), ...
                    -30,30,1e-10));

fprintf(' gauss-hermite quadrature\n');
for n=1:((degree+1)/2+4)
  % use hermite -- will be exact for n>=(degree+1)/2
  int = gaussHermite(n);
  % notice the change of variables!
  ep = polyval(p,int.x*sqrt(2))'*int.w/sqrt(pi);
  if (n==round((degree+1)/2))
    fprintf('--- the rest should be exact ----\n');
  end
  fprintf('n=%2d  E[p(x)] = %.10g\n',n,ep);
end

fprintf('\n monte carlo integration\n')
for n=1:6
  fprintf('n=10^%d  E[p(x)] = %.10g\n',n, ...
          mean(polyval(p,randn(10^n,1))));
end

pause();

% Example 2: E[p(x)|x>-0.1], x~N(0,1)
fprintf('adaptive quad, tol %.1g = %.10g\n', ...
        1e-10, quad(@(x) polyval(p,x).*normpdf(x), ...
                    -0.1,30,1e-10) / ...
        quad(@(x) normpdf(x),-0.1,30,1e-10));

fprintf('\n gauss-hermite quadrature (not the right rule)\n');
for n=1:2:60
  int = gaussHermite(n);
  m = int.x>-0.1;
  ep = polyval(p,int.x(m)*sqrt(2))'* ...
       int.w(m)/sqrt(pi) ...
       * sum(int.w)/sum(int.w(m));
  fprintf('n=%2d  E[p(x)] = %.10g\n',sum(m),ep);
end

fprintf('\n monte carlo integration\n')
for n=1:6
  x = randn(10^n,1);
  x = x(x>-0.1);
  fprintf('n=10^%d  E[p(x)] = %.10g\n',numel(x), ...
          mean(polyval(p,x)));
end
pause();

% Example 2: E[p(x)|x>-0.1], x~N(0,1)
% Example 2: E[p(x)|x>-0.1], x~N(0,1)
fprintf('adaptive quad, tol %.1g = %.10g\n', ...
        1e-10, quad(@(x) exp(x).*normpdf(x), ...
                    -30,30,1e-10));

fprintf('\n gauss-hermite quadrature (not the right rule)\n');
for n=1:20
  int = gaussHermite(n);
  ep = exp(int.x*sqrt(2))'*int.w/sqrt(pi);
  fprintf('n=%2d  E[p(x)] = %.10g\n',n,ep);
end

fprintf('\n monte carlo integration\n')
for n=1:6
  x = randn(10^n,1);
  fprintf('n=10^%d  E[p(x)] = %.10g\n',numel(x), ...
          mean(exp(x)));
end
pause();

repo: Aero-matlab


text:
function [E,J]=SynthMeasWatsonSHCylSingleRadGPD_PGSE(x, grad_dirs, G, delta, smalldel, fibredir, roots)
% Substrate: Impermeable cylinders with one radius in a homogeneous background.
% Orientation distribution: Watson's distribution with SH approximation
% Pulse sequence: Pulsed gradient spin echo
% Signal approximation: Gaussian phase distribution.
%
% [E,J]=SynthMeasWatsonSHCylSingleRadGPD_PGSE(x, grad_dirs, G, delta, smalldel, fibredir, roots)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameters.  The Jacobian does not
% include derivates with respect to the fibre direction.
%
% x is the list of model parameters in SI units:
% x(1) is the volume fraction of the intracellular space.
% x(2) is the free diffusivity of the material inside and outside the cylinders.
% x(3) is the hindered diffusivity outside the cylinders in perpendicular directions.
% x(4) is the radius of the cylinders.
% x(5) is the concentration parameter of the Watson's distribution
%
% grad_dirs is the gradient direction for each measurement.  It has size [N
% 3] where N is the number of measurements.
%
% G, delta and smalldel are the gradient strength, pulse separation and
% pulse length of each measurement in the protocol.  Each has
% size [N 1].
%
% fibredir is a unit vector along the symmetry axis of the Watson's
% distribution.  It must be in Cartesian coordinates [x y z]' with size [3 1].
%
% roots contains solutions to the Bessel function equation from function
% BesselJ_RootsCyl.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

% Duplication with SynthMeasDistributedRadVG remains, because of the
% derivative computation.

f=x(1);
dPar=x(2);
dPerp=x(3);
R=x(4);
kappa=x(5);

% build the input x vector for hindered compartment
x_h = [dPar dPerp kappa];

% build the input x vector for restricted diffusion in Neuman cylinder model
% set diffusion coeff in restricted compartment same as parallel one in
% hindered.
x_r = [dPar R kappa];

% Synthesize measurements from model
if (nargout>1)
	[E_h, J_h] = SynthMeasWatsonHinderedDiffusion_PGSE(x_h, grad_dirs, G, delta, smalldel, fibredir);
	[E_r, J_r] = SynthMeasWatsonSHCylNeuman_PGSE(x_r, grad_dirs, G, delta, smalldel, fibredir, roots);
else
	E_h = SynthMeasWatsonHinderedDiffusion_PGSE(x_h, grad_dirs, G, delta, smalldel, fibredir);
	E_r = SynthMeasWatsonSHCylNeuman_PGSE(x_r, grad_dirs, G, delta, smalldel, fibredir, roots);
end

E=(1-f)*E_h+f*E_r;

% Compute the Jacobian matrix
if(nargout>1)

    % dE_tot/df = E_r - E_h
    dEtdf = E_r - E_h;

    % dE_tot/ddPar
    dEtddPar = (1-f)*J_h(:,1) + f*J_r(:,1);

    % dE_tot/ddPerp
    dEtddPerp = (1-f)*J_h(:,2);

    % dE_tot/dR
    dEtdr = f*J_r(:,2);

    % dE_tot/dk
    dEtdk = (1-f)*J_h(:,3) + f*J_r(:,3);

    % Construct the jacobian matrix.
    J = zeros(length(E), 5);
    J(:,1) = dEtdf;
    J(:,2) = dEtddPar;
    J(:,3) = dEtddPerp;
    J(:,4) = dEtdr;
    J(:,5) = dEtdk;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% [INPUT]
% data = A float t-by-n matrix containing the time series to be smooted.
% s = An integer [1,t] representing the span of the smoothing filter (optional, default=21).
%
% [OUTPUT]
% data = A float t-by-n matrix containing the smoothed time series.

function data = smooth_data(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' '2d' 'nonempty'}));
        ip.addOptional('s',21,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 1 'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    [data,s] = validate_input(ipr.data,ipr.s);

    nargoutchk(1,1);

    data = smooth_data_internal(data,s);

end

function data = smooth_data_internal(data,s)

    w = s + mod(s,2) - 1;

    for i = 1:size(data,2)
        data(:,i) = smoothing_function(data(:,i),s,w);
    end

end

function ys = smoothing_function(y,s,w)

    nan_indices = isnan(y);
    nan_found = any(nan_indices);

    if (~nan_found)
        if (w == 1)
            ys = y;
            return;
        end

        n = numel(y);
        z = filter(ones(w,1) ./ w,1,y);

        ys_begin = cumsum(y(1:w-2));
        ys_begin = ys_begin(1:2:end) ./ (1:2:w-2).';

        ys_end = cumsum(y(n:-1:n-w+3));
        ys_end = ys_end(end:-2:1)./(w-2:-2:1)';

        ys = [ys_begin; z(w:end); ys_end];
    else
        z1 = y;
        z1(nan_indices) = 0;

        z2 = double(~nan_indices);

        ys = smoothing_function(z1,s,w) ./ smoothing_function(z2,s,w);
    end

    ys(nan_indices) = NaN;

end

function [data,s] = validate_input(data,s)

    if (isscalar(data))
        error('The value of ''data'' is invalid. Expected input to be a vector or a 2d matrix.');
    end

    [t,n] = size(data);

    if ((any([t n] == 1)) && (n > t))
        error('The value of ''data'' is invalid. Expected input to be a column vector when a single time series is defined.');
    end

    if (s > t)
        error(['The value of ''s'' is invalid. Expected input to be less than or equal to ' num2str(t) '.']);
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [params] = backwardCompatible(params, fieldNames, defaultValue)
  for ii=1:length(fieldNames)
    field = fieldNames{ii};
    if ~isfield(params, field)
      params.(field) = defaultValue;
    end
  end
end

author: lmthang


text:
function cc_levels_minmax_display ( )

%*****************************************************************************80
%
%% CC_LEVELS_MINMAX_DISPLAY displays grids generated by CC_LEVELS_MINMAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVELS_MINMAX_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display the nested 2D Clenshaw-Curtis grids\n' );
  fprintf ( 1, '  generated by CC_LEVELS_MINMAX.\n' );

  dim_num = 2;

  while ( 1 )
%
%  Get user input.
%
    level_min = input ( 'Enter LEVEL_MIN or RETURN to exit;' );

    if ( isempty ( level_min ) )
      break
    end

    if ( level_min < 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  We require 0 <= LEVEL_MIN\n' );
      continue
    end

    level_max = input ( 'Enter LEVEL_MAX or RETURN to exit;' );

    if ( isempty ( level_max ) )
      break
    end

    if ( level_max < level_min )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  We require LEVEL_MIN <= LEVEL_MAX!\n' );
      continue
    end
%
%  Compute data.
%
    [ grid_num, point_num ] = cc_levels_minmax_size ( dim_num, ...
      level_min, level_max );

    [ grid_level, grid_order, grid_point ] = cc_levels_minmax ( dim_num, ...
      level_min, level_max,  grid_num, point_num );

    clf
%
%  We have to name the axes in order to control the grid.
%
    axes_handle = axes;
%
%  Plot the points.
%
    handle = scatter ( grid_point(1,:), grid_point(2,:), 'filled' );
%
%  Force the plotting region to be square, not rectangular.
%
    axis square
%
%  Request grid lines.
%
    grid on
%
%  Specify the location of the grid lines, and suppress labeling.
%
    set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'xticklabel', [] );
    set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'yticklabel', [] );
%
%  Make the plotting region slightly bigger than the data.
%
    axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
%
%  Title
%
    s = sprintf ( '%d <= LEVEL <= %d', level_min, level_max );
    title ( s );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVELS_MINMAX_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end


text:
function s = read_curry(filename)

% READ_CURRY reads and parses Curry V2 and V4 ascii files and returns the
% content in a structure that is similar to the block-structured layout of
% the file. This function does not interpret the content of the file, but
% is intended as a helper function for READ_CURRY_XXX functions (where XXX
% is the extension of the file).
%
% Use as
%   s = read_curry(filename)

% Copyright (C) 2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen_or_error(filename, 'rt');

s    = [];
line = [];

while ~feof(fid)
  line = fgetl(fid);
  blockname = '';

  % search for the beginning of a block
  if ~isempty(strfind(line, 'START'))
    [blockname, blocktype] = strtok(line);
    blockname = strtrim(blockname);
    blocktype = strtrim(strtok(blocktype));
    val = {};

    % fprintf('reading block %s\n', blockname);
    while isempty(strfind(line, 'END'))

      % read the next line
      line = fgetl(fid);
      % process the data in the line
      if strcmp(blocktype, 'START') && ~isempty(strfind(line, '='))
        % the line looks like "lhs = rhs"
        [lhs, rhs] = strtok(line, '=');
        lhs = strtrim(lhs);
        rhs(1) = ' ';
        rhs = strtrim(rhs);
        val = str2double(rhs);
        if isnan(val)
          val = rhs;
        end
        s = setfield(s, blockname, lhs, val);
      elseif strcmp(blocktype, 'START_LIST')
        % this contains unstructured data, process this line as being part of a list
        val{end+1} = line;
      else
        % process this line differently
        % FIXME
      end

      if isnumeric(line) || feof(fid)
        ft_error('unexpected end of file');
      end

    end % while ~END
  end % if START of block

  % store the collected cell-array if it was part of an unstructured list
  if strcmp(blocktype, 'START_LIST')
    val = val(1:(end-1));   % remove the last item, which is the end of the block
    val = val(:);           % make it a column array
    s = setfield(s, blockname, val);
  end

  % clear them to prevent confusion on the next loop
  blockname = '';
  blocktype = '';

end  % while ~feof

fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function W_RF  = yuweiA2(V_D,V_RF)

global Vn H Nrf Nr Nk;
W_RF = ones(Nr,Nrf);

for k = 1:Nk
    F(:,:,k) = H(:,:,k)*V_RF*V_D(:,:,k)*V_D(:,:,k)'*V_RF'*H(:,:,k)';
end
F = sum(F,3)/Nk;
g = 1/Nr;
a = g/Vn;

for Nloop = 1:10
    for j = 1:Nrf
        VRF = W_RF;
        VRF(:,j)=[];
        C = eye(Nrf-1)+a*VRF'*F*VRF;
        G = a*F-a^2*F*VRF*C^(-1)*VRF'*F;
        for i = 1:Nr
            for l = 1:Nr
                if i~=l
                    x(l)=G(i,l)*W_RF(l,j);
                end
            end
            n = sum(x);
            if n ==0
                W_RF(i,j)=1;
            else
                W_RF(i,j)=n/abs(n);
            end
        end
    end
end

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
 function x = ir_read_mat(file, pick)
%function x = ir_read_mat(file, pick)
%|
%| Read a matlab matrix from file, e.g., 'data.mat'
%| If pick is an integer, it indexes which field in the file is desired.
%| Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if ~isvar('pick') || isempty(pick), pick = 0; end

s = load(file);
names = fieldnames(s);
name = names{1};
if length(names) > 1
	if ~isa(pick, 'numeric')
		error 'only numeric picking done now'
	end
	if pick < 1
		disp(names)
		warn('using default first array: %s', name)
	elseif pick > length(names)
		error 'pick too large'
	else
		name = names{pick};
		printm('Using array %s', name)
	end
end
x = s.(name);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function result=L1median(x,tol);

%L1MEDIAN is an orthogonally equivariant location estimator,
% also known as the spatial median. It is defined as the point which
% minimizes the sum of the Euclidean distances to all observations in the
% data matrix x. It can resist 50% outliers.
%
% Reference (for the algorithm):
%    Hossjer, O. and Croux, C. (1995)
%    "Generalizing univariate signed rank statistics for testing and estimating
%    a multivariate location parameter", Nonparametric Statistics, 4, 293-308.
%
% Required input argument:
%    x : either a data matrix with n observations in rows, p variables in columns
%        or a vector of length n.
%
% Optional input argument:
%  tol : convergence criterium; the iterative process stops when the norm between two solutions < tol.
%        (default = 1.e-08).
%
% I/O: result=L1median(x,tol);
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Original Gauss code by C. Croux, translated to MATLAB by Sabine Verboven
% Last updated 18/02/2004

if nargin <2
   tol=1.e-08;
end;
[n,p]=size(x);
maxstep=200;
%initializing starting value for m
m=median(x);
k=1;
while (k<=maxstep)
   mold=m;
   Xext=sortrows([norme(x-repmat(m,n,1)) x],1);
   dx=Xext(:,1);
   x=Xext(:,2:p+1);
   if all(dx)
      w=1./dx;
   else
      ww=dx(all(dx,2));
      w=1./ww;
      w=[zeros(length(dx)-length(w),1);w];
   end
   delta=sum((x-repmat(m,n,1)).*repmat(w,1,p),1)./sum(w);
   nd=norme(delta);
   if all(nd<tol)
      maxhalf=0;
   else
      maxhalf=log2(nd/tol);
   end
   m=mold+delta;   %computation of a new estimate
   nstep=0;
   while all(mrobj(x,m)>=mrobj(x,mold))&(nstep<=maxhalf)
      nstep=nstep+1;
      m=mold+delta./(2^nstep);
   end
   if (nstep>maxhalf)
      mX=mold;
      break
   end
   k=k+1;
end
if k>maxstep
   display('Iteration failed')
end
result=m;

%------------------------------------------------------------------
function n=norme(x);

%NORME calculates the Euclidian norm of a matrix x
% the output is a column vector containing the norm of each row
% I/O: n=norme(X);

for i=1:size(x,1)
   n(i)=norm(x(i,:));
end;
n=n';

%------------------------------------------------------------------
function s=mrobj(x,m)

%MROBJ computes the objective function in m based on x and a

xm=norme(x-repmat(m,size(x,1),1));
s=sum(xm,1)';

author: marianux

repo: ecg-kit

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


text:
% Modelling Sensor Directivity in 3D Example
%
% This example demonstrates how the sensitivity of a large single element
% detector varies with the angular position of a point-like source. It is a
% 3D version of the Modelling Sensor Directivity in 2D example.
%
% author: Ben Cox
% date: 29th October 2010
% last update: 23rd February 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 64;            % number of grid points in the x direction
Ny = 64;            % number of grid points in the y direction
Nz = 64;            % number of grid points in the z direction
dx = 100e-3/Nx;     % grid point spacing in the x direction [m]
dy = dx;            % grid point spacing in the y direction [m]
dz = dx;            % grid point spacing in the z direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

% define the properties of the propagation medium
medium.sound_speed = 1500;	% [m/s]

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
Nt = length(kgrid.t_array);

% define a large area detector
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(Nx/2+1, (Ny/2-9):(Ny/2+11),(Nz/2-9):(Nz/2+11)) = 1;

% define equally spaced point sources lying on a circle centred at the
% centre of the detector face
Nangles = 11;
circle = makeCartCircle(25*dx, Nangles, [0,0], pi);
circle = [circle; zeros(1,Nangles)];

% find the binary sensor mask most closely corresponding to the cartesian
% points coordinates from makeCartCircle
circle3D = cart2grid(kgrid,circle);

% find the indices of the sources in the binary source mask
source_positions = find(circle3D == 1);

% define a time varying sinusoidal source
source_freq = 0.25e6;
source_mag = 1;
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);

% filter the source to remove high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);

% pre-allocate array for storing the output time series
single_element_data = zeros(Nt,length(source_positions));

% run a simulation for each of these sources to see the effect that the
% angle from the detector has on the measured signal
for source_loop = 1:length(source_positions)

    % select a point source
    source.p_mask = zeros(Nx,Ny,Nz);
    source.p_mask(source_positions(source_loop)) = 1;

    % create a display mask to display the transducer
    display_mask = source.p_mask + sensor.mask;

    % run the simulation
    input_args = {'PMLSize', 10, 'DisplayMask', display_mask, 'PlotScale', [-0.2 0.2], 'PlotFreq', 50,  'DataCast', 'single'};
    sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

    % average the data recorded for each grid point to simulate the
    % measured signal from a large aperture, single element, detector
    single_element_data(:,source_loop) = sum(sum(sensor_data,1),1);

end

% =========================================================================
% VISUALISATION
% =========================================================================

% plot source points and sensor mask
voxelPlot(circle3D + sensor.mask);
view([-34 34])

% plot the time series recorded for each of the sources
figure;
plot(kgrid.t_array, single_element_data);
colormap(getColorMap);
xlabel('time [s]');
ylabel('pressure');
title('time series from each direction');

% calculate angle between source and centre of detector face
angles = atan((kgrid.y(source_positions))./kgrid.x(source_positions));

% plot the maximum amplitudes for each of the sources, showing that the
% detector sensitivity falls off at low angles as expected.
figure;
plot(angles,max(single_element_data),'o')
colormap(getColorMap);
xlabel('angle between source and centre of detector face');
ylabel('maximum detected pressure from each direction')


text:
function g = matRad_objectiveGradient(optiProb,apertureInfoVec,dij,cst)
% matRad IPOPT callback: gradient function for direct aperture optimization
%
% call
%   g = matRad_objectiveGradient(optiProb,apertureInfoVec,dij,cst)
%
% input
%   optiProb:           option struct defining the type of optimization
%   apertureInfoVect:   aperture info in form of vector
%   dij:                matRad dij struct as generated by bixel-based dose calculation
%   cst:                matRad cst struct
%
% output
%   g: gradient
%
% References
%   [1] http://dx.doi.org/10.1118/1.4914863
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% update apertureInfo, bixel weight vector an mapping of leafes to bixels
if ~isequal(apertureInfoVec,optiProb.apertureInfo.apertureVector)
    optiProb.apertureInfo = optiProb.matRad_daoVec2ApertureInfo(optiProb.apertureInfo,apertureInfoVec);
end
apertureInfo = optiProb.apertureInfo;

% bixel based gradient calculation
bixelG = matRad_objectiveGradient@matRad_OptimizationProblem(optiProb,apertureInfo.bixelWeights,dij,cst);

% allocate gradient vector for aperture weights and leaf positions
g = NaN * ones(size(apertureInfoVec,1),1);

% 1. calculate aperatureGrad
% loop over all beams
offset = 0;
for i = 1:numel(apertureInfo.beam)

    % get used bixels in beam
    ix = ~isnan(apertureInfo.beam(i).bixelIndMap);

    % loop over all shapes and add up the gradients x openingFrac for this shape
    for j = 1:apertureInfo.beam(i).numOfShapes
        g(j+offset) = apertureInfo.beam(i).shape(j).shapeMap(ix)' ...
                        * bixelG(apertureInfo.beam(i).bixelIndMap(ix));
    end

    % increment offset
    offset = offset + apertureInfo.beam(i).numOfShapes;

end

% 2. find corresponding bixel to the leaf Positions and aperture
% weights to calculate the gradient
g(apertureInfo.totalNumOfShapes+1:end) = ...
        apertureInfoVec(apertureInfo.mappingMx(apertureInfo.totalNumOfShapes+1:end,2)) ...
     .* bixelG(apertureInfo.bixelIndices(apertureInfo.totalNumOfShapes+1:end)) / apertureInfo.bixelWidth;

% correct the sign for the left leaf positions
g(apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs) = ...
    -g(apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs);

repo: matRad


text:
function sigm = sigmoid(x)

    sigm = 1 ./ (1 + exp(-x));
end

author: singaxiong

repo: SignalGraph


text:
function obj = t_mp_table(quiet)
%T_MP_TABLE  Tests for MP_TABLE (and TABLE).

%   MATPOWER
%   Copyright (c) 2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

% define_constants;
if quiet
    verbose = 0;
else
    verbose = 1;
end

nt = 166;
skip_tests_for_tablicious = 1;
table_classes = {@mp_table};
class_names = {'mp_table'};
if have_feature('table')
    table_classes = {@table, table_classes{:}};
    class_names = {'table', class_names{:}};
end
nc = length(table_classes);

t_begin(nc*nt, quiet);

%% set up example data
var1 = [1:6]';
var2 = var1 + var1/10;
var3 = {'one'; 'two'; 'three'; 'four'; 'five'; 'six'};
var4 = 1./var1;
var5 = [3;-5;0;1;-2;4] <= 0;
var6 = [10 15; 20 25; 30 35; 40 45; 50 55; 60 65];
var_names = {'igr', 'flt', 'str', 'dbl', 'boo', 'mat'};
var_values = {var1, var2, var3, var4, var5, var6};
row_names = {'row1'; 'second row'; 'row 3'; 'fourth row'; '5'; '6th row'};
dim_names = {'records', 'fields'};

%% run tests
for k = 1:nc
    table_class = table_classes{k};
    cls = upper(class_names{k});

    t = sprintf('%s : constructor - no args : ', cls);
    T = table_class();
    skip_oct_tab = isa(T, 'table') && have_feature('octave') && ...
            skip_tests_for_tablicious;
    skip_ml_tab = isa(T, 'table') && have_feature('matlab') && ...
            have_feature('matlab', 'vnum') < 9.010;
    t_ok(isa(T, class_names{k}), [t 'class']);
    t_ok(isempty(T), [t 'isempty']);

    t = sprintf('%s : constructor - ind vars : ', cls);
    T = table_class(var1, var2, var3, var4, var5(:), var6);
    t_ok(isa(T, class_names{k}), [t 'class']);
    t_ok(~isempty(T), [t 'not isempty']);
    t_is(size(T), [6 6], 12, [t 'sz = size(T)']);
    t_is(size(T, 1), 6, 12, [t 'sz = size(T, 1)']);
    t_is(size(T, 2), 6, 12, [t 'sz = size(T, 2)']);
    [nr, nz] = size(T);
    t_is([nr, nz], [6 6], 12, [t '[nr, nz] = size(T)']);
    t_ok(isequal(T.Properties.VariableNames, ...
        {'var1', 'var2', 'var3', 'var4', 'Var5', 'var6'}), [t 'VariableNames'] );
    t_ok(isempty(T.Properties.RowNames), [t 'RowNames'] );
    t_ok(isequal(T.Properties.DimensionNames, {'Row', 'Variables'}), ...
        [t 'DimensionNames'] );
% show_me(T)

    %% subsref () (w/o RowNames)
    t = sprintf('%s : subsref () w/o RowNames : ', cls);
    t_is(size(T([2;4],[1;3;4])), [2 3], 12, [t 'T(ii,jj) : size(T(ii,jj))']);
% show_me(T([2;4],[1;3;4]));
    t_is(size(T(:,1:2:5)), [6 3], 12, [t 'T(:,j1:jN) : size']);
% show_me(T(:,1:2:5));
    t_is(size(T(1:2:5,:)), [3 6], 12, [t 'T(i1:iN,:) : size']);
% show_me(T(1:2:5,:));
    t_is(size(T(5,4)), [1 1], 12, [t 'T(i,j) : size']);
    t_is(size(T(4,6)), [1 1], 12, [t 'T(i,j) : size']);
% show_me(T(5,4));

    t = sprintf('%s : constructor - ind vars, w/names : ', cls);
    if skip_oct_tab || skip_ml_tab
        T = table_class(var1, var2, var3, var4, var5, var6, ...
            'VariableNames', var_names, 'RowNames', row_names);
    else
        T = table_class(var1, var2, var3, var4, var5, var6, ...
            'VariableNames', var_names, 'RowNames', row_names, ...
            'DimensionNames', dim_names);
    end
    t_ok(isa(T, class_names{k}), [t 'class'])
    t_ok(~isempty(T), [t 'not isempty']);
    t_is(size(T), [6 6], 12, [t 'sz = size(T)']);
    t_is(size(T, 1), 6, 12, [t 'sz = size(T, 1)']);
    t_is(size(T, 2), 6, 12, [t 'sz = size(T, 2)']);
    [nr, nz] = size(T);
    t_is([nr, nz], [6 6], 12, [t '[nr, nz] = size(T)']);
    t_ok(isequal(T.Properties.VariableNames, var_names), [t 'VariableNames'] );
    t_ok(isequal(T.Properties.RowNames, row_names), [t 'RowNames'] );
    if skip_oct_tab || skip_ml_tab
        t_skip(1, [t 'DimensionNames not yet supported'] );
    else
        t_ok(isequal(T.Properties.DimensionNames, dim_names), [t 'DimensionNames'] );
    end

    %% . indexing
    %% get full variables
    t = sprintf('%s : subsref . : ', cls);
    t_is(T.igr, var1, 12, [t 'T.igr']);
    t_is(T.flt, var2, 12, [t 'T.flt']);
    t_ok(isequal(T.str, var3), [t 'T.str']);
    t_is(T.dbl, var4, 12, [t 'T.dbl']);
    t_is(T.boo, var5, 12, [t 'T.boo']);
    t_is(T.mat, var6, 12, [t 'T.mat']);

    %% get indexed variables
    t_is(T.igr(2), var1(2), 12, [t 'T.igr(i)']);
    t_is(T.flt(4:6), var2(4:6), 12, [t 'T.flt(i1:iN)']);
    t_ok(isequal(T.str{3}, var3{3}), [t 'T.str{i}']);
    t_ok(isequal(T.str(6:-1:4), var3(6:-1:4)), [t 'T.str(iN:-1:i1)']);
    t_is(T.dbl([5;3]), var4([5;3]), 12, [t 'T.dbl(ii)']);
    t_is(T.boo(var5 == 1), var5(var5 == 1), 12, [t 'T.boo(<logical>)']);
    t_is(T.mat([6;2;4], :), var6([6;2;4], :), 12, [t 'T.mat(ii, :)']);
    t_is(T.mat([6;2;4], 2), var6([6;2;4], 2), 12, [t 'T.mat(ii, j)']);
    t_is(T.mat([6;2;4], [2;1]), var6([6;2;4], [2;1]), 12, [t 'T.mat(ii, jj)']);
    t_is(T.mat(5, :), var6(5, :), 12, [t 'T.mat(i, :)']);
    t_is(T.mat(5, 2), var6(5, 2), 12, [t 'T.mat(i, j)']);
    t_is(T.mat(8:10), var6(8:10), 12, [t 'T.mat(jj)']);

    %% set full variables
    t = sprintf('%s : subsasgn . : ', cls);
    T.igr = var1(end:-1:1);
    T.flt = var2(end:-1:1);
    T.str = var3(end:-1:1);
    T.dbl = var4(end:-1:1);
    T.boo = var5(end:-1:1);
    T.mat = var6(end:-1:1, end:-1:1);
    t_is(T.igr, var1(end:-1:1), 12, [t 'T.igr']);
    t_is(T.flt, var2(end:-1:1), 12, [t 'T.flt']);
    t_ok(isequal(T.str, var3(end:-1:1)), [t 'T.str']);
    t_is(T.dbl, var4(end:-1:1), 12, [t 'T.dbl']);
    t_is(T.boo, var5(end:-1:1), 12, [t 'T.boo']);
    t_is(T.mat, var6(end:-1:1, end:-1:1), 12, [t 'T.mat']);

    %% set indexed variables
    T.igr(2) = 55;
    T.flt(4:6) = [0.4 0.5 0.6];
    T.str{3} = 'tres';
    T.str(6:-1:4) = {'seis'; 'cinco'; 'cuatro'};
    T.dbl([5;3]) = [pi; exp(1)];
    T.boo(var5 == 1) = false;
    T.mat([5;3],:) = [5 0.5; 3 0.3];

    v1 = var1(end:-1:1); v1(2) = 55;
    v2 = var2(end:-1:1); v2(4:6) = [0.4 0.5 0.6];
    v3 = var3(end:-1:1);
    v3(6:-1:3) = {'seis'; 'cinco'; 'cuatro'; 'tres'};
    v4 = var4(end:-1:1); v4([5;3]) = [pi; exp(1)];
    v5 = var5(end:-1:1); v5(var5 == 1) = false;
    v6 = var6(end:-1:1, end:-1:1); v6([5;3],:) = [5 0.5; 3 0.3];
    t_is(T.igr, v1, 12, [t 'T.igr(i)']);
    t_is(T.flt, v2, 12, [t 'T.flt(i1:iN)']);
    t_ok(isequal(T.str, v3), [t 'T.str{i}']);
    t_is(T.dbl, v4, 12, [t 'T.dbl(ii)']);
    t_is(T.boo, v5, 12, [t 'T.boo(<logical>)']);
    t_is(T.mat, v6, 12, [t 'T.mat(ii, :)']);
% show_me(T);

    %% {} indexing
    t = sprintf('%s : subsref {} : ', cls);
    if skip_oct_tab
        t_skip(6, [t 'T{:, j} syntax not yet supported'])
    else
        t_ok(isequal(T{:, 1}, v1), [t 'T{:, 1} == v1']);
        t_ok(isequal(T{:, 2}, v2), [t 'T{:, 2} == v2']);
        t_ok(isequal(T{:, 3}, v3), [t 'T{:, 3} == v3']);
        t_ok(isequal(T{:, 4}, v4), [t 'T{:, 4} == v4']);
        t_ok(isequal(T{:, 5}, v5), [t 'T{:, 5} == v5']);
        t_ok(isequal(T{:, 6}, v6), [t 'T{:, 6} == v6']);
    end

    t_ok(isequal(T{2, 1}, v1(2)), [t 'T{i, j} == v<j>(i)']);
    if skip_oct_tab
        t_skip(1, [t 'T{i, j} syntax not yet supported for matrices'])
    else
        t_ok(isequal(T{4, 6}, v6(4, :)), [t 'T{i, j} == v<j>(i)']);
    end

    %% remove skip_oct_tab when
    %% https://github.com/apjanke/octave-tablicious/pull/89 is merged
    if skip_oct_tab
        t_skip(1, [t 'T{end, end} syntax not yet supported'])
    else
        t_ok(isequal(T{end, end}, v6(end,:)), [t 'T{end, end} == v<end>(end,:)']);
    end
    if skip_oct_tab
        t_skip(8, [t 'T{ii, j} syntax not yet supported'])
    else
        t_ok(isequal(T{1:3, 2}, v2(1:3)), [t 'T{i1:iN, j} == v<j>(i1:iN)']);
        t_ok(isequal(T{1:3, 6}, v6(1:3, :)), [t 'T{i1:iN, j} == v<j>(i1:iN, :)']);
        t_ok(isequal(T{[6;3], 5}, v5([6;3])), [t 'T{ii, j} == v<j>(ii)']);
        t_ok(isequal(T{[6;3], 6}, v6([6;3], :)), [t 'T{ii, j} == v<j>(ii, :)']);
        t_ok(isequal(T{6:-1:3, [2;4;5]}, [v2(6:-1:3) v4(6:-1:3) v5(6:-1:3)]), [t 'T{iN:-1:i1, jj}']);
        t_ok(isequal(T{6:-1:3, [2;4;6]}, [v2(6:-1:3) v4(6:-1:3) v6(6:-1:3, :)]), [t 'T{iN:-1:i1, jj}']);
        t_ok(isequal(T{:, [4;1]}, [v4 v1]), [t 'T{:, jj} == [v<j1> v<j2> ...]']);
        t_ok(isequal(T{:, [4;6;1]}, [v4 v6 v1]), [t 'T{:, jj} == [v<j1> v<j2> ...]']);
    end

    t = sprintf('%s : subsasgn {} : ', cls);
    T{:, 1} = var1;
    t_is(T.igr, var1, 12, [t 'T{:, 1} = var1']);
    T{:, 2} = var2;
    t_is(T.flt, var2, 12, [t 'T{:, 2} = var2']);
    T{:, 3} = var3;
    t_ok(isequal(T.str, var3), [t 'T{:, 3} = var3']);
    T{:, 4} = var4;
    t_is(T.dbl, var4, 12, [t 'T{:, 4} = var4']);
    T{:, 5} = var5;
    t_is(T.boo, var5, 12, [t 'T{:, 5} = var5']);
    if skip_oct_tab
        t_skip(1, 'T{:, j} = M syntax not yet supported for matrices');
    else
        T{:, 6} = var6;
        t_is(T.mat, var6, 12, [t 'T{:, 6} = var6']);
    end

    T{2, 1} = v1(2);
    t_ok(isequal(T{2, 1}, v1(2)), [t 'T{i, j} = v<j>(i)']);
    if skip_oct_tab
        t_skip(1, 'T{i, j} = M syntax not yet supported for matrices');
    else
        T{4, 6} = v6(4, :);
        t_ok(isequal(T{4, 6}, v6(4, :)), [t 'T{i, j} = v<j>(i, :)']);
    end
    T{1:3, 2} = v2(1:3);
    t_ok(isequal(T.flt(1:3), v2(1:3)), [t 'T{i1:iN, j} = v<j>(i1:iN)']);
    if skip_oct_tab
        t_skip(2, 'T{i1:iN, j} = M syntax not yet supported for matrices');
    else
        T{1:3, 6} = v6(1:3, :);
        t_ok(isequal(T.mat(1:3, :), v6(1:3, :)), [t 'T{i1:iN, j} = v<j>(i1:iN, :)']);
        T{1:3, 6} = v6(1:3, [2; 1]);
        t_ok(isequal(T.mat(1:3, :), v6(1:3, [2; 1])), [t 'T{i1:iN, j} = v<j>(i1:iN, jj)']);
    end
    T{[6;3], 5} = v5([6;3]);
    t_ok(isequal(T.boo([6;3]), v5([6;3])), [t 'T{ii, j} = v<j>(ii)']);
    if skip_oct_tab
        t_skip(2, 'T{ii, j} = M syntax not yet supported for matrices');
    else
        T{[6;3], 6} = v6([6;3], :);
        t_ok(isequal(T.mat([6;3], :), v6([6;3], :)), [t 'T{ii, j} = v<j>(ii, :)']);
        T{[6;3], 6} = v6([6;3], [2;1]);
        t_ok(isequal(T.mat([6;3], :), v6([6;3], [2;1])), [t 'T{ii, j} = v<j>(ii, jj)']);
    end
    if skip_oct_tab
        t_skip(2, [t 'T{ii, jj} syntax not yet supported'])
    else
        T{6:-1:3, [2;4;6;5]} = [v2(6:-1:3) v4(6:-1:3) v6(6:-1:3, :) v5(6:-1:3)];
        t_ok(isequal(T{6:-1:3, [2;4;6;5]}, [v2(6:-1:3) v4(6:-1:3) v6(6:-1:3, :) v5(6:-1:3)]), [t 'T{iN:-1:i1, jj}']);
        T{:, [6;4;1]} = [v6 v4 v1];
        t_ok(isequal(T{:, [6;4;1]}, [v6 v4 v1]), [t 'T{:, jj} = [v<j1> v<j2> ...]']);
    end
    T{:, 1} = v1;
    T{:, 2} = v2;
    T{:, 3} = v3;
    T{:, 4} = v4;
    T{:, 5} = v5;
    if ~skip_oct_tab
        T{:, 6} = v6;
    end

    %% () indexing
    t = sprintf('%s : subsref () : T(ii,jj) : ', cls);
    ii = [2;4]; jj = [1;3;4;6];
    T2 = T(ii,jj);
% show_me(T2);
    t_is(size(T2), [2 4], 12, [t 'size']);
    t_ok(isequal(T2.Properties.VariableNames, var_names(jj)), [t 'VariableNames']);
    t_ok(isequal(table_values(T2), {v1(ii), v3(ii), v4(ii), v6(ii, :)}), [t 'VariableValues']);
    t_ok(isequal(T2.Properties.RowNames, row_names(ii)), [t 'RowNames']);
    if skip_oct_tab || skip_ml_tab
        t_skip(1, [t 'DimensionNames not yet supported'])
    else
        t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
    end

    t = sprintf('%s : subsref () : T(:,j1:s:jN) : ', cls);
    jj = [2:2:6];
    T2 = T(:,2:2:6);
% show_me(T2);
    t_is(size(T2), [6 3], 12, [t 'size']);
    t_ok(isequal(T2.Properties.VariableNames, var_names(jj)), [t 'VariableNames']);
    t_ok(isequal(table_values(T2), {v2, v4, v6}), [t 'VariableValues']);
    t_ok(isequal(T2.Properties.RowNames, row_names), [t 'RowNames']);
    if skip_oct_tab || skip_ml_tab
        t_skip(1, [t 'DimensionNames not yet supported'])
    else
        t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
    end

    t = sprintf('%s : subsref () : T(i1:s:iN,:)) : ', cls);
    ii = [1:2:5];
    T2 = T(1:2:5,:);
% show_me(T2);
    t_is(size(T2), [3 6], 12, [t 'size']);
    t_ok(isequal(T2.Properties.VariableNames, var_names), [t 'VariableNames']);
    t_ok(isequal(table_values(T2), {v1(ii), v2(ii), v3(ii), v4(ii), v5(ii), v6(ii, :)}), [t 'VariableValues']);
    t_ok(isequal(T2.Properties.RowNames, row_names(ii)), [t 'RowNames']);
    if skip_oct_tab || skip_ml_tab
        t_skip(1, [t 'DimensionNames not yet supported'])
    else
        t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
    end

    t = sprintf('%s : subsref () : T(i,j) : ', cls);
    T2 = T(5,4);
% show_me(T2);
    t_is(size(T2), [1 1], 12, [t 'size']);
    t_ok(isequal(T2.Properties.VariableNames, var_names(4)), [t 'VariableNames']);
    t_ok(isequal(table_values(T2), {v4(5)}), [t 'VariableValues']);
    t_ok(isequal(T2.Properties.RowNames, row_names(5)), [t 'RowNames']);
    if skip_oct_tab || skip_ml_tab
        t_skip(1, [t 'DimensionNames not yet supported'])
    else
        t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
    end
    T2 = T(3,6);
% show_me(T2);
    t_is(size(T2), [1 1], 12, [t 'size']);
    t_ok(isequal(T2.Properties.VariableNames, var_names(6)), [t 'VariableNames']);
    t_ok(isequal(table_values(T2), {v6(3, :)}), [t 'VariableValues']);
    t_ok(isequal(T2.Properties.RowNames, row_names(3)), [t 'RowNames']);
    if skip_oct_tab || skip_ml_tab
        t_skip(1, [t 'DimensionNames not yet supported'])
    else
        t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
    end

    %% remove skip_oct_tab
    %% when https://github.com/apjanke/octave-tablicious/pull/89 is merged
    if skip_oct_tab
        t_skip(5, sprintf('%s : T(end, end) not yet supported.', cls));
    else
        t = sprintf('%s : subsref () : T(end, end) : ', cls);
        T2 = T(end, end);
    % show_me(T2);
        t_is(size(T2), [1 1], 12, [t 'size']);
        t_ok(isequal(T2.Properties.VariableNames, var_names(6)), [t 'VariableNames']);
        t_ok(isequal(table_values(T2), {v6(end, :)}), [t 'VariableValues']);
        t_ok(isequal(T2.Properties.RowNames, row_names(end)), [t 'RowNames']);
        if skip_oct_tab || skip_ml_tab
            t_skip(1, [t 'DimensionNames not yet supported'])
        else
            t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
        end
    end

    if skip_oct_tab
        t_skip(20, sprintf('%s : subsasgn () not yet supported.', cls));
    else
        t = sprintf('%s : subsasgn () : T(ii,jj) : ', cls);
        ii0 = [2;4]; jj = [1;3;4;6];
        T3 = T(ii0,jj);
        ii = [4;2];
        T(ii,jj) = T3;
        T2 = T(ii,jj);
        t_ok(isequal(T2.Properties.VariableNames, var_names(jj)), [t 'VariableNames']);
        t_ok(isequal(table_values(T2), {v1(ii0), v3(ii0), v4(ii0), v6(ii0, :)}), [t 'VariableValues']);
        t_ok(isequal(T2.Properties.RowNames, row_names(ii)), [t 'RowNames']);
        if skip_oct_tab || skip_ml_tab
            t_skip(1, [t 'DimensionNames not yet supported'])
        else
            t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
        end
        T(ii0, jj) = T2;        %% restore

        t = sprintf('%s : subsasgn () : T(:,j1:s:jN) : ', cls);
        jj = [2:2:6];
        T3 = T(end:-1:1,2:2:6);
        T(:,2:2:6) = T3;
        T2 = T(:,2:2:6);
        t_ok(isequal(T2.Properties.VariableNames, var_names(jj)), [t 'VariableNames']);
        t_ok(isequal(table_values(T2), {v2(end:-1:1), v4(end:-1:1), v6(end:-1:1, :)}), [t 'VariableValues']);
        t_ok(isequal(T2.Properties.RowNames, row_names), [t 'RowNames']);
        if skip_oct_tab || skip_ml_tab
            t_skip(1, [t 'DimensionNames not yet supported'])
        else
            t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
        end
        T(end:-1:1,2:2:6) = T3;     %% restore

        t = sprintf('%s : subsasgn () : T(i1:s:iN,:)) : ', cls);
        ii0 = [1:2:5];
        T3 = T(1:2:5,:);
        ii = [5:-2:1];
        T(5:-2:1, :) = T3;
        T2 = T(5:-2:1, :);
        t_ok(isequal(T2.Properties.VariableNames, var_names), [t 'VariableNames']);
        t_ok(isequal(table_values(T2), {v1(ii0), v2(ii0), v3(ii0), v4(ii0), v5(ii0), v6(ii0, :)}), [t 'VariableValues']);
        t_ok(isequal(T2.Properties.RowNames, row_names(ii)), [t 'RowNames']);
        if skip_oct_tab || skip_ml_tab
            t_skip(1, [t 'DimensionNames not yet supported'])
        else
            t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
        end
        T(1:2:5, :) = T3;   %% restore

        t = sprintf('%s : subsref () : T(i,j) : ', cls);
        T3 = T(5,4);
        T3{1,1} = exp(1);
        T(5,4) = T3;
        T2 = T(5,4);
        t_ok(isequal(T2.Properties.VariableNames, var_names(4)), [t 'VariableNames']);
        t_ok(isequal(table_values(T2), {exp(1)}), [t 'VariableValues']);
        t_ok(isequal(T2.Properties.RowNames, row_names(5)), [t 'RowNames']);
        if skip_oct_tab || skip_ml_tab
            t_skip(1, [t 'DimensionNames not yet supported'])
        else
            t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
        end
        T{5,4} = pi;

        T3 = T(3,6);
        T3{1,1} = [111 222];
        T(3,6) = T3;
        T2 = T(3,6);
        t_ok(isequal(T2.Properties.VariableNames, var_names(6)), [t 'VariableNames']);
        t_ok(isequal(table_values(T2), {[111 222]}), [t 'VariableValues']);
        t_ok(isequal(T2.Properties.RowNames, row_names(3)), [t 'RowNames']);
        if skip_oct_tab || skip_ml_tab
            t_skip(1, [t 'DimensionNames not yet supported'])
        else
            t_ok(isequal(T2.Properties.DimensionNames, dim_names), [t 'DimensionNames']);
        end
        T{3,6} = [3 0.3];
    end

    %% vertical concatenation
    t = sprintf('%s : vertical concatenation : ', cls);
    T1 = T(1:3, :);
    T2 = T(4:6, :);
    T3 = T([4:6 1:3], :);
    T4 = [T2; T1];
    t_ok(isequal(T4, T3), [t '[T1;T2]']);
    T5 = table_class([7;8],[7.7;8.8],{'seven';'eight'},1./[7;8],[-1;2]<=0, [70 75; 80 85], ...
        'VariableNames', var_names);
    if skip_oct_tab
        t_skip(1, [t 'RowNames auto-generation not yet supported']);
    else
        T6 = [T5; T2; T1];
        t_ok(isequal(T6, [T5;T3]), [t '[T1;T2;T3]']);
    end

    %% horizontal concatenation
    t = sprintf('%s : horizontal concatenation : ', cls);
    T1 = T(:, 1:3);
    T2 = T(:, 4:6);
    T3 = T(:, [4:6 1:3]);
    T4 = [T2 T1];
    t_ok(isequal(T4, T3), [t '[T1 T2]']);
    T5 = table_class(var3, var2);
    T6 = [T2 T1 T5];
    t_ok(isequal(T6, [T4 T5]), [t '[T1 T2 T3]']);

    %% deleting variables
    t = sprintf('%s : delete variables : ', cls);
    if skip_oct_tab
        t_skip(1, [t 'not yet supported']);
    else
        T7 = T;
        T7(:, [5 6 4]) = [];
        t_ok(isequal(T7, T1), [t 'T(:, j) = []']);
    end

    %% adding variables
    t = sprintf('%s : add variables : ', cls);
    T7 = T1;
    T7.dbl = T2.dbl;
    T7.boo = T2.boo;
    T7.mat = T2.mat;
    t_ok(isequal(T7, T), [t 'T.new_var = val']);

    %% more {} indexing
    t = sprintf('%s : more subsref {} : ', cls);
    T2 = T(:, [1;2;4;6;5]);
    T3 = T6(:, 6:7);
    ii = [5;3;1];
    jj = [6:7];
    if skip_oct_tab
        t_skip(6, [t 'T{ii, :} syntax not yet supported'])
    else
        ex = horzcat(v3, var3);
        t_is(T2{ii, :}, [v1(ii) v2(ii) v4(ii) v6(ii, :) v5(ii)], 12, [t 'T{ii,:} (double)']);
        t_is(T2{:, :}, [v1 v2 v4 v6 v5], 12, [t 'T{:,:} (double)']);
        t_ok(isequal(T6{ii, 6:7}, ex(ii, :)), [t 'T{ii,j1:j2} (cell)']);
        t_ok(isequal(T6{:, 6:7}, ex), [t 'T{:,j1:j2} (cell)']);
        t_ok(isequal(T3{ii, :}, ex(ii, :)), [t 'T{ii,:} (cell)']);
        t_ok(isequal(T3{:, :}, ex), [t 'T{:,:} (cell)']);
    end

    t = sprintf('%s : more subsasgn {} : ', cls);
    if skip_oct_tab
        t_skip(6, [t 'T{ii, :} syntax not yet supported'])
    else
        T2{ii, :} = [var1(ii) var2(ii) var4(ii) var6(ii, :) var5(ii)];
        t_is(T2{ii, :}, [var1(ii) var2(ii) var4(ii) var6(ii, :) var5(ii)], 12, [t 'T{ii,:} (double)']);
        T2{:, :} = [v1 v2 v4 v6 v5];
        t_is(T2{:, :}, [v1 v2 v4 v6 v5], 12, [t 'T{:,:} (cell)']);
        v1 = horzcat(var3, v3);
        v2 = horzcat(v3, var3);
        T6{ii, 6:7} = v1(ii, :);
        t_ok(isequal(T6{ii, 6:7}, v1(ii, :)), [t 'T{ii,j1:j2} (cell)']);
        T6{:, 6:7} = v2;
        t_ok(isequal(T6{:, 6:7}, v2), [t 'T{:,j1:j2} (cell)']);
        T3{ii, :} = v1(ii, :);
        t_ok(isequal(T3{ii, :}, v1(ii, :)), [t 'T{ii,:} (cell)']);
        T3{:, :} = v2;
        t_ok(isequal(T3{:, :}, v2), [t 'T{:,:} (cell)']);
    end

    %% value class, not handle class
    t = sprintf('%s : value class, not handle', cls);
    T2 = T;
    T.igr(2) = 2;
    t_is(T.igr(2), 2, 12, t);
    t_is(T2.igr(2), 55, 12, t);

    %% nested tables
    t = sprintf('%s : nested tables : ', cls);
    T1 = table_class([1;2;3],[4;5;6]);
    T2 = table_class({'one';'two';'three'},{'four';'five';'six'});
    T3 = table_class([10;20;30],{'uno';'dos';'tres'}, T1, T2);
    t_ok(isequal(T3.T1, T1), [t 'T.T1 == T1']);
    t_ok(isequal(T3.T2, T2), [t 'T.T2 == T2']);
    t_ok(isequal(T3.T1([1;3], :), T1([1;3], :)), [t 'T.T1(ii, :) == T1(ii, :)']);
    t_ok(isequal(T3.T2([1;3], :), T2([1;3], :)), [t 'T.T2(ii, :) == T2(ii, :)']);
    t_ok(isequal(T3.T1.Var1, T1.Var1), [t 'T.T1.Var1 == T1.Var1']);
    t_ok(isequal(T3.T2.Var2([1;3], :), T2.Var2([1;3], :)), [t 'T.T2.Var2(ii) == T.T2.Var2(ii)']);
    %% check for Octave trouble in subsref
    T4 = T3(:, :);
    T5 = T3([1;2;3], :);
    t_ok(isequal(T3, T4), [t 'T == T(:, :))']);
    t_ok(isequal(T3, T5), [t 'T == T(ii, :))']);
end

if nargout
    obj = T;
end

t_end;

function Tv = table_values(T)
if isa(T, 'table') && have_feature('matlab')
    Tv = cellfun(@(c)T{:, c}, num2cell(1:size(T, 2)), 'UniformOutput', false);
else
    Tv = T.Properties.VariableValues;
end

function show_me(T)
if isa(T, 'table') && have_feature('octave')
    prettyprint(T)
else
    T
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower