text:
function [net, stats] = cnn_train_modified(net, imdb, getBatch, varargin)
%CNN_TRAIN  An example implementation of SGD for training CNNs
%    CNN_TRAIN() is an example learner implementing stochastic
%    gradient descent with momentum to train a CNN. It can be used
%    with different datasets and tasks by providing a suitable
%    getBatch function.
%
%    The function automatically restarts after each training epoch by
%    checkpointing.
%
%    The function supports training on CPU or on one or more GPUs
%    (specify the list of GPU IDs in the `gpus` option).

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005;
opts.momentum = 0.9 ;
opts.saveMomentum = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.conserveMemory = true ;
opts.backPropDepth = +Inf ; %+Inf
opts.sync = false ;
opts.cudnn = true ;
opts.errorFunction = 'softmax_regression' ;
opts.errorLabels = {} ;
opts.plotDiagnostics = false ;
opts.plotStatistics = true;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end
if isnan(opts.val), opts.val = [] ; end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

net = vl_simplenn_tidy(net); % fill in some eventually missing values
net.layers{end-1}.precious = 1; % do not remove predictions, used for error
vl_simplenn_display(net, 'batchSize', opts.batchSize) ;

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  for i=1:numel(net.layers)
    J = numel(net.layers{i}.weights) ;
    if ~isfield(net.layers{i}, 'learningRate')
      net.layers{i}.learningRate = ones(1, J) ;
    end
    if ~isfield(net.layers{i}, 'weightDecay')
      net.layers{i}.weightDecay = ones(1, J) ;
    end
  end
end

% setup error calculation function
hasError = true ;
if isstr(opts.errorFunction)
  switch opts.errorFunction
    case 'none'
      opts.errorFunction = @error_none ;
      hasError = false ;
    case 'multiclass'
      opts.errorFunction = @error_multiclass ;
      if isempty(opts.errorLabels), opts.errorLabels = {'top1err', 'top5err'} ; end
    case 'binary'
      opts.errorFunction = @error_binary ;
      if isempty(opts.errorLabels), opts.errorLabels = {'binerr'} ; end
    case 'softmax_regression'
      opts.errorFunction = @error_crossentropy_or_l2norm ;
      if isempty(opts.errorLabels), opts.errorLabels = {'crossentropy_or_l2norm'} ; end
    otherwise
      error('Unknown error function ''%s''.', opts.errorFunction) ;
  end
end

state.getBatch = getBatch ;
stats = [] ;

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
%   prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(params.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;

    if size(imdb.images.labels,1) > 1
        [net.meta.svm_model1,net.meta.svm_err1] = training_SVR(net,params);
    end

    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

%   if params.plotStatistics
%     switchFigure(1) ; clf ;
%     plots = setdiff(...
%       cat(2,...
%       fieldnames(stats.train)', ...
%       fieldnames(stats.val)'), {'num', 'time'}) ;
%     for p = plots
%       p = char(p) ;
%       values = zeros(0, epoch) ;
%       leg = {} ;
%       for f = {'train', 'val'}
%         f = char(f) ;
%         if isfield(stats.(f), p)
%           tmp = [stats.(f).(p)] ;
%           values(end+1,:) = tmp(1,:)' ;
%           leg{end+1} = f ;
%         end
%       end
%       subplot(1,numel(plots),find(strcmp(p,plots))) ;
%       plot(1:epoch, values','o-') ;
%       xlabel('epoch') ;
%       title(p) ;
%       legend(leg{:}) ;
%       grid on ;
%     end
%     drawnow ;
%     print(1, modelFigPath, '-dpdf') ;
%   end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function err = error_multiclass(params, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
[~,predictions] = sort(predictions, 3, 'descend') ;

% be resilient to badly formatted labels
if numel(labels) == size(predictions, 4)
  labels = reshape(labels,1,1,1,[]) ;
end

% skip null labels
mass = single(labels(:,:,1,:) > 0) ;
if size(labels,3) == 2
  % if there is a second channel in labels, used it as weights
  mass = mass .* labels(:,:,2,:) ;
  labels(:,:,2,:) = [] ;
end

m = min(5, size(predictions,3)) ;

error = ~bsxfun(@eq, predictions, labels) ;
err(1,1) = sum(sum(sum(mass .* error(:,:,1,:)))) ;
err(2,1) = sum(sum(sum(mass .* min(error(:,:,1:m,:),[],3)))) ;

% -------------------------------------------------------------------------
function err = error_crossentropy_or_l2norm(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
if size(predictions,3) > 1
    error = sum( - labels(:) .* log(predictions(:)+eps(1)));
else
    error = sum((predictions(:) - labels(:)).^2);
end
err = error;

% -------------------------------------------------------------------------
function err = error_binary(params, labels, res)
% -------------------------------------------------------------------------
% predictions = gather(res(end-1).x) ;
predictions = squeeze(gather(res(end-1).x))' ;
error = bsxfun(@times, predictions, labels) < 0 ;
err = sum(error(:)) ;

% -------------------------------------------------------------------------
function err = error_none(params, labels, res)
% -------------------------------------------------------------------------
err = zeros(0,1) ;

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.momentum)
  for i = 1:numel(net.layers)
    for j = 1:numel(net.layers{i}.weights)
      state.momentum{i}{j} = 0 ;
    end
  end
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net = vl_simplenn_move(net, 'gpu') ;
  for i = 1:numel(state.momentum)
    for j = 1:numel(state.momentum{i})
      state.momentum{i}{j} = gpuArray(state.momentum{i}{j}) ;
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  vl_simplenn_start_parserv(net, parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

subset = params.(mode) ;
num = 0 ;
stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;
adjustTime = 0 ;
res = [] ;
error = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, params.epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    [im, labels] = params.getBatch(params.imdb, batch) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch) ;
    end

    if numGpus >= 1
      im = gpuArray(im) ;

    end

    if strcmp(mode, 'train')
      dzdy = 1 ;
      evalMode = 'normal' ;
    else
      dzdy = [] ;
      evalMode = 'test' ;
    end
    net.layers{end}.class = labels ;
    res = vl_simplenn_modified(net, im, dzdy, res, ...
                      'accumulate', s ~= 1, ...
                      'mode', evalMode, ...
                      'conserveMemory', params.conserveMemory, ...
                      'backPropDepth', params.backPropDepth, ...
                      'sync', params.sync, ...
                      'cudnn', params.cudnn, ...
                      'parameterServer', parserv, ...
                      'holdOn', s < params.numSubBatches) ;

    % accumulate errors
    error = error + sum(double(gather(res(end).x))) ;
  end

  % accumulate gradient
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    [net, res, state] = accumulateGradients(net, res, state, params, batchSize, parserv) ;
  end

  % get statistics
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
%   stats = extractStats(net, params, error / num) ;
  stats.objective = error / num;
  stats.num = num ;
  stats.time = time ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.6f', f, stats.(f)) ;
  end
  fprintf('\n') ;

  % collect diagnostic statistics
  if strcmp(mode, 'train') && params.plotDiagnostics
    switchFigure(2) ; clf ;
    diagn = [res.stats] ;
    diagnvar = horzcat(diagn.variation) ;
    diagnpow = horzcat(diagn.power) ;
    subplot(2,2,1) ; barh(diagnvar) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTick', 1:numel(diagnvar), ...
      'YTickLabel',horzcat(diagn.label), ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1], ...
      'XTick', 10.^(-5:1)) ;
    grid on ;
    subplot(2,2,2) ; barh(sqrt(diagnpow)) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTick', 1:numel(diagnpow), ...
      'YTickLabel',{diagn.powerLabel}, ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1e5], ...
      'XTick', 10.^(-5:5)) ;
    grid on ;
    subplot(2,2,3); plot(squeeze(res(end-1).x)) ;
    drawnow ;
  end
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveMomentum
  state.momentum = [] ;
else
  for i = 1:numel(state.momentum)
    for j = 1:numel(state.momentum{i})
      state.momentum{i}{j} = gather(state.momentum{i}{j}) ;
    end
  end
end

net = vl_simplenn_move(net, 'cpu') ;

% -------------------------------------------------------------------------
function [net, res, state] = accumulateGradients(net, res, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for l=numel(net.layers):-1:1
  for j=numel(res(l).dzdw):-1:1

    if ~isempty(parserv)
      tag = sprintf('l%d_%d',l,j) ;
      parDer = parserv.pull(tag) ;
    else
      parDer = res(l).dzdw{j}  ;
    end

    if j == 3 && strcmp(net.layers{l}.type, 'bnorm')
      % special case for learning bnorm moments
      thisLR = net.layers{l}.learningRate(j) ;
      net.layers{l}.weights{j} = vl_taccum(...
        1 - thisLR, ...
        net.layers{l}.weights{j}, ...
        thisLR / batchSize, ...
        parDer) ;
    else
      % Standard gradient training.
      thisDecay = params.weightDecay * net.layers{l}.weightDecay(j) ;
      thisLR = params.learningRate * net.layers{l}.learningRate(j) ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, net.layers{l}.weights{j}) ;

        % Update momentum.
        state.momentum{l}{j} = vl_taccum(...
          params.momentum, state.momentum{l}{j}, ...
          -1, parDer) ;

        % Nesterov update (aka one step ahead).
        if params.nesterovUpdate
          delta = vl_taccum(...
            params.momentum, state.momentum{l}{j}, ...
            -1, parDer) ;
        else
          delta = state.momentum{l}{j} ;
        end

        % Update parameters.
        net.layers{l}.weights{j} = vl_taccum(...
          1, net.layers{l}.weights{j}, ...
          thisLR, delta) ;
      end
    end

    % if requested, collect some useful stats for debugging
    if params.plotDiagnostics
      variation = [] ;
      label = '' ;
      switch net.layers{l}.type
        case {'conv','convt'}
          variation = thisLR * mean(abs(state.momentum{l}{j}(:))) ;
          power = mean(res(l+1).x(:).^2) ;
          if j == 1 % fiters
            base = mean(net.layers{l}.weights{j}(:).^2) ;
            label = 'filters' ;
          else % biases
            base = sqrt(power) ;%mean(abs(res(l+1).x(:))) ;
            label = 'biases' ;
          end
          variation = variation / base ;
          label = sprintf('%s_%s', net.layers{l}.name, label) ;
      end
      res(l).stats.variation(j) = variation ;
      res(l).stats.power = power ;
      res(l).stats.powerLabel = net.layers{l}.name ;
      res(l).stats.label{j} = label ;
    end
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(net, params, errors)
% -------------------------------------------------------------------------
stats.objective = errors(1) ;
for i = 1:numel(params.errorLabels)
  stats.(params.errorLabels{i}) = errors(i+1) ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net, state)
% -------------------------------------------------------------------------
% save(fileName, 'net', 'state') ;
save(fileName, 'net') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = vl_simplenn_tidy(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
%clear vl_tmove vl_imreadjpeg ;
disp('Clearing mex files') ;
clear mex ;
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(params, cold)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end
end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename) ;
  clearMex() ;
  if numGpus == 1
    disp(gpuDevice(params.gpus)) ;
  else
    spmd
      clearMex() ;
      disp(gpuDevice(params.gpus(labindex))) ;
    end
  end
end

function [svm_model,svm_err] = training_SVR(net,params)

net = vl_simplenn_move(net, 'gpu') ;
net.layers = net.layers(1:end-1);
batch = 256;
sampling = (length(params.imdb.images.score) / 10000);
scores_train = params.imdb.images.score(1:sampling:end);
data = params.imdb.images.data(:,:,:,1:sampling:end);
num_batch = floor(length(scores_train) / batch);
for i = 1:num_batch
    images = gpuArray(data(:,:,:,(i-1)*batch+1:i*batch));
    res = vl_simplenn(net, images, [],[],'Mode','test','ConserveMemory',true);
    feats(:,(i-1)*batch+1:i*batch) = squeeze(gather(res(end).x));
end
images = gpuArray(data(:,:,:,num_batch*batch+1:length(scores_train)));
res = vl_simplenn(net, images, [],[],'Mode','test','ConserveMemory',true);
feats(:,num_batch*batch+1:length(scores_train)) = squeeze(gather(res(end).x));

svm_model = svmtrain(scores_train', double(feats'), '-s 4 -t 0');
[scores_pred, ~,~] = svmpredict(scores_train', double(feats'), svm_model,'-q');
svm_err = (sum(abs(scores_pred - scores_train'))) / length(scores_train);

author: HuiZeng


text:
function plotmainmap(obj)
    % PLOTMAINMAP set up main map window

    % TAG : PURPOSE
    % "active quakes" : selected events
    % "big evens" : selected events, above a threshhold magnitude

    % %
    % this probably belongs somewhere else...
    if ~isempty(obj.shape) && ~isvalid(obj.shape)
        msg.dbdisp('shape had been deleted');
        obj.shape = ShapeGeneral();
    end
    % %
    axm = obj.map_axes;
    axm.Visible = 'off';
    if isempty(axm)||~isvalid(axm)
        error('Somehow lost track of main map');
    end

    % update the active earthquakes
    eq=findobj(axm,'Tag','active quakes');

    mainEventOpts = obj.mainEventProps; % local copy
    szFcn = str2func(mainEventOpts.MarkerSizeFcn);

    if mainEventOpts.UseDotsForTooManyEvents && obj.catalog.Count > mainEventOpts.HowManyAreTooMany
        mainEventOpts.Marker = '.';
    end
    if isempty(eq)
        % CREATE the plot

        axm.NextPlot='add';
        dispname = replace(obj.catalog.Name,'_','\_');

        szFcn = str2func(mainEventOpts.MarkerSizeFcn);
        eq=scatter(axm, obj.catalog.X, obj.catalog.Y, ...
            szFcn(obj.catalog.Magnitude), getLegalColors(),...
            'Tag','active quakes',...
            'HitTest','off',...
            'DisplayName',dispname);
        eq.ZData = obj.catalog.Z;
        axm.NextPlot='replace';
        %obj.do_colorbar(axm);

    else

        % REUSE the plot
        eq.XData = obj.catalog.X;
        eq.YData = obj.catalog.Y;
        eq.ZData = obj.catalog.Z;
        eq.SizeData = szFcn(obj.catalog.Magnitude);
        eq.MarkerEdgeColor='flat';
        eq.CData = getLegalColors();

        % this is a kludge, because if a MarkerEdgeColor is defined that isn't specifically 'flat'
        % then it overrides the CData colors.
        if size(eq.CData(:,1)) > 1
            mainEventOpts=renameStructField(mainEventOpts,'MarkerEdgeColor','Marker_Edge_Color');
        else
            mainEventOpts=renameStructField(mainEventOpts,'Marker_Edge_Color','MarkerEdgeColor');
        end
        dispname = replace(obj.catalog.Name, '_', '\_');
        if ~strcmp(eq.DisplayName, dispname)
            eq.DisplayName = dispname;
        end
    end

    fix_colorbar()

    set_valid_properties(eq, mainEventOpts);

    % update the largest events
    update_large()

    % update the shape
    axm.NextPlot='add';
    if ~isempty(obj.shape)
        obj.shape.plot(axm);
    end
    axm.NextPlot='replace';

    % update the grid
    if ~isempty(obj.Grid)
        if isempty(obj.shape) && all(obj.Grid.ActivePoints(:))
            % do nothing needs to be done.
            obj.Grid.plot(obj.map_axes, 'HitTest', 'off', 'ActiveOnly');
        else
            maskedGrid = obj.Grid.MaskWithShape(obj.shape);
            if ~isequal(maskedGrid.ActivePoints, obj.Grid.ActivePoints)
                obj.Grid.ActivePoints = maskedGrid.ActivePoints;
                obj.Grid.plot(obj.map_axes, 'HitTest', 'off', 'ActiveOnly');
            end
        end
    end
    axm.Visible='on';

    function update_large()
        beq = findobj(axm,'Tag','big events');

        if ~isempty(obj.bigEvents)
            beq.XData = obj.bigEvents.X;
            beq.YData = obj.bigEvents.Y;
            beq.ZData = obj.bigEvents.Z;
            beq.SizeData=mag2dotsize(obj.bigEvents.Magnitude);
        else
            [beq.XData, beq.YData, beq.ZData, beq.SizeData]=deal([]);
        end

    end

    function c = getLegalColors()
        % because datetime isn't allowed
        switch  obj.colorField
            case '-none-'
                c=[0 0 .15];
            case 'Date'
                c=datenum(obj.catalog.Date);
            otherwise
                c=obj.catalog.(obj.colorField);
        end
    end

    function fix_colorbar()
        h = obj.map_axes.Colorbar;
        if isempty(h) && obj.colorField ~= "-none-"
            h = colorbar(obj.map_axes);
        end

        if isempty(h) || h.Label.String == string(obj.colorField)
            return
        end

        h.Label.String = obj.colorField;

        switch obj.colorField
            case '-none-'
                delete(h)

            case 'Depth'
                h.Direction = 'reverse';
                h.TickLabels = h.Ticks;

            case 'Date'
                    h.TickLabels   = datestr(h.Ticks,'yyyy-mm-dd');
                    h.Label.String = 'Date';
                    h.Direction    = 'normal';

            otherwise
                h.Direction = 'normal';
                h.TickLabels = h.Ticks;
        end

    end

end

author: CelsoReyes


text:
function [tmpl, cord] = sample_pos(frame, p, out_size, random, num_scale)
    %% sample positive samples with different scales and translations

    len = (p(3)+p(4)) / 2;
    for iter = 1 : 2
        for i = 1 : num_scale
            context_scale = 0.5 * (i);
            translationX = randi(floor(len / 2 * context_scale)) - 1;
            translationY = randi(floor(len / 2 * context_scale)) - 1;
            if randi(2) == 1
                translationX = -translationX;
            end
            if randi(2) == 1
                translationY = -translationY;
            end
            if random == 0
                translationX = 0;
                translationY = 0;
            end
            p_sample = [p(1) + translationX, p(2) + translationY, p(3) + len * context_scale, p(4) + len * context_scale];
            crop_img = im_crop(frame, p_sample);
            scale_height = size(crop_img,1) / out_size;
            scale_width = size(crop_img,2) / out_size;
            center_x = (1+out_size) / 2 - translationX / scale_width;
            center_y = (1+out_size) / 2 - translationY / scale_height;
            x1 = floor(center_x - (p(3) / 2 ) / scale_width);
            y1 = floor(center_y - (p(4) / 2 ) / scale_height);
            x2 = ceil(center_x + (p(3) / 2 ) / scale_width);
            y2 = ceil(center_y + (p(4) / 2 ) / scale_height);
            %minus 1 for numeric reasons
            cord(iter * num_scale - num_scale + i, :) = [x1, y1, x2, y2] - 1 ;
            crop_img = imresize(crop_img, [out_size, out_size], 'Antialiasing', false);
            tmpl(:,iter * num_scale - num_scale + i) = crop_img(:);
        end
    end
end

author: flyers

repo: drone-tracking


text:
% TD | Tucker-ADAL | Tucker Decomposition solved by ADAL (Goldfarb and Qin, 2013)
% process_video('TD', 'Tucker-ADAL', 'dataset/demo.avi', 'output/demo_Tucker-ADAL.avi');

alg_path_aux = fullfile(lrs_conf.td_path,'RLRT');
addpath(genpath(alg_path_aux));

pdata.T = T;
pdata.X = T;
N = ndims(pdata.T);
r = 1/sqrt(max(size(pdata.T)));
params.E0 = tenzeros(size(pdata.T));
params.X0 = tenzeros(size(pdata.T));
params.V0 = cell(1, N);
for i = 1:N
  params.V0{i} = tenzeros(size(pdata.T));
end
params.mu0 = 1/(N+1);
params.mode = N;
params.IsTC = false; % is tensor completion
params.rRatio = 1/4;
params.opt_tol = 1e-3;
params.eta = 1/(N+1);
params.max_iter = 1000;
params.mu1fac = 10;
params.mu1 = params.mu1fac*std(T(:));
params.mu2 = params.mu1;
params.mu_min = 1e-4;
params.mu_max = 1e2;
params.lambdaS = 1;
params.lambda = params.lambdaS*r*params.rRatio;
params.verbose = 1;
params.use_cont = true;
params.k = [size(T,1) size(T,2) 1];
%%%%%%%%%% for PROPACK %%%%%%%%%%%%
% declare global var 'sv'
global sv;
global tmode;
global use_propack;
global curr_mu;
sv =  ceil(min(size(pdata.T)) * 0.1) * ones( 1, N );
use_propack = true;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
results = tensor_tucker_adal_ncx(pdata, params);
L = double(results.X);
S = double(results.E);
clear sv tmode use_propack curr_mu;

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
%  doinvers
% This file calculates orintation of the stress tensor
% based on Gephard's algorithm.
% stress tensor orientation. The actual calculation is done
% using a call to a fortran program.
%
% Stefan Wiemer 03/96

global mi mif1 mif2  hndl3 a newcat2 mi2
global tmp cumu2
report_this_filefun(mfilename('fullpath'));
think

tmp = [newt2(:,10:14)];
save /home/stefan/ZMAP/invers/data.inp tmp -ascii
infi = ['/home/stefan/ZMAP//invers/data.inp'];
outfi = ['/home/stefan/ZMAP/tmpout.dat'];

cd /home/stefan/ZMAP/invers
com1 =input('Which computer?','s');
tic
comm = ['! rsh ' com1 ' /home/stefan/ZMAP/invers/invshell1 ',...
    num2str(length(tmp(:,1))) ' ' num2str(i)  ' &']
eval(comm)

t = toc/60


text:
%High speed MATLAB codes

%This program use SCL with small list size as a "filter", i.e., under the
%same noise realization, If CRC-SCL with smaller L is correct, then CRC-SCL
%with larger L must be correct. You may not believe in this conjecture, but you can have
%a try, This is true with high probability (almost 1)

%This program has another accelerator. If in lower snr, CRC-SCL is correct, then at higher snr,
%the same CRC-SCL must be correct.  You may not believe in this conjecture, but you can have
%a try, This is true with high probability (almost 1).

%I know above methods sound dangerous, but it is safe under following 4
%conditions

%1. The same code word
%2. The same AWGN noise realization with distribution N(0, 1)
%3. CRC must be used (Only SCL is not permitted)
%4. The same code construction in all SNR range.

%This program satisfies above 4 conditions. You can verify the bler performance by comparisons with existing results.

%Since MATLAB is not good at recursive function (too many parameters to be passed)
%I cancel the well-known recursiveCalcP() and recursiveCalcB() proposed by
%I. Tal
%Instead, 'For' function is used. You may argue that we can use objective
%oriented (OO) style. However, OO in matlab is also slow.

%Besides, the algorithms of following papers are provided.

%How to Construct Polar Codes
%Fast Successive-Cancellation Decoding of Polar Codes: Identification and Decoding of New Nodes
%beta-expansion: A Theoretical Framework for Fast and Recursive Construction of Polar Codes

%Above three algorithms are made by myself so the correctness is not guaranteed.

clear
addpath('GA/')
% addpath('HowToConstructPolarCode/')
addpath('NodeProcess/')
% addpath('BECconstruction/')
% addpath('PolarizaedChannelsPartialOrder/')
%adding above folders will take round 2 seconds

design_epsilon = 0.32;
crc_length = 16;
[gen, det, g] = get_crc_objective(crc_length);
n = 10;
N = 2^n;
K = N/2 + crc_length;
ebno_vec = [2 2.5]; %row vec, you can write it like [1 1.5 2 2.5 3]
list_vec = [1 16];  %row vec, you can write it like [1 4 16 32 ...]. The first element is always 1 for acceleration purpose. The ramaining elements are power of two.
max_runs = 1e7;
max_err = 100;
resolution = 1e4;%the results are shown per max_runs/resolution.
[bler, ber] = simulation(N, K, design_epsilon, max_runs, max_err, resolution,  ebno_vec, list_vec, gen, det, g, crc_length);

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
disp('Fitting AR model to AR-5 process data with known coefficients');

secs=5;
ns=128;
t=[1/ns:1/ns:secs];
N=length(t);
noise_var=1^2;
noise=sqrt(noise_var)*randn(1,N);
a_true=[-1.8517,1.3741,0.1421,-0.6852,0.3506];
y=filter(1,[1,a_true],noise);
y=y(1:N);
figure
plot(t,y);
xlabel('Seconds');
title('Sample of AR-5 process');

total_variance=std(y)^2;
signal_variance=total_variance-noise_var;
snr=sqrt(signal_variance)/sqrt(noise_var);
disp(sprintf('SNR=%1.3f',snr));

disp('For model order p=5');
ar=spm_ar (y,5);

disp(' ');
disp('True coefficients');
disp(a_true);

disp(' ');
disp('Estimated coefficients');
disp(ar.a_mean');

for p=1:10,
    disp(sprintf('Now fitting model with p=%d coefficients',p));
    ar=spm_ar (y,p,0);
    logev(p)=ar.fm;
end
logev=logev-min(logev);
figure
bar(logev);
ylabel('Log Evidence');
xlabel('Model order');

author: spm


text:
%%**********************************************************************
%% HSDNTpred: Compute (dX,dy,dZ) for NT direction.
%%
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'.
%% compute W satisfying W*Z*W = X.
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************************

function [par,dX,dy,dZ,coeff,L,hRd] = ...
    HSDNTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol)

global schurfun schurfun_par
%%
%% compute NT scaling matrix
%%
[par.W,par.G,par.sv,par.gamx,par.gamz,par.dd,par.ee,par.ff] = ...
    NTscaling(blk,X,Z,Zchol,invZchol);
%%
%% compute schur matrix
%%
m = par.m;
schur = sparse(m+2,m+2);
UU = []; EE = [];
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if strcmp(pblk{1},'l')
        [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,par.dd);
    elseif strcmp(pblk{1},'q');
        [schur,UU,EE] = schurmat_qblk(blk,At,par,schur,UU,EE,p,par.dd,par.ee);
    elseif strcmp(pblk{1},'s')
        if isempty(schurfun{p})
            schur = schurmat_sblk(blk,At,par,schur,p,par.W);
        elseif ischar(schurfun{p})
            if ~isempty(par.permZ{p})
                Wp = par.W{p}(par.permZ{p},par.permZ{p});
            else
                Wp = par.W{p};
            end
            schurtmp = feval(schurfun{p},Wp,Wp,schurfun_par(p,:));
            schur = schur + schurtmp;
        end
    end
end
%%
%% compute rhs
%%
[rhs,EinvRc,hRd] = HSDNTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu);
%%
%% solve linear system
%%
par.addschur   = par.kap/par.tau;
schur(m+1,m+1) = schur(m+1,m+1) + par.kap/par.tau;
schur(m+2,m+2) = schur(m+2,m+2) + par.addschur;
[xx,coeff,L] = HSDlinsysolve(par,schur,UU,EE,par.Umat,rhs);
%%
%% compute (dX,dZ)
%%
[par,dX,dy,dZ] = HSDNTdirfun(blk,At,par,Rd,EinvRc,xx);
%%**********************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [M_low]=voxelate(M,reduction_factor)

% function [M_low]=voxelate(M,reduction_factor)
% ------------------------------------------------------------------------
% This function lowers the resulution of a 3D image based on
% 'reduction_factor'.
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 25/05/2008
% ------------------------------------------------------------------------

%%
num_dims=ndims(M);

switch num_dims

    case 2
        %Setting up cell shape
        cell_shape_1=reduction_factor*ones(size(M,1)/reduction_factor,1);
        cell_shape_2=reduction_factor*ones(size(M,2)/reduction_factor,1);

        %Converting matrix to cell
        M_cell=mat2cell(M, cell_shape_1, cell_shape_2);

        %Calculating average for each cell
        M_cell=cellfun(@mean, cellfun(@mean, M_cell, 'UniformOutput',0), 'UniformOutput',0);

    case 3
        %Setting up cell shape
        cell_shape_1=reduction_factor*ones(size(M,1)/reduction_factor,1);
        cell_shape_2=reduction_factor*ones(size(M,2)/reduction_factor,1);
        cell_shape_3=reduction_factor*ones(size(M,3)/reduction_factor,1);

        %Converting matrix to cell
        M_cell=mat2cell(M, cell_shape_1, cell_shape_2, cell_shape_3);

        %Calculating average for each cell
        M_cell=cellfun(@mean, cellfun(@mean, cellfun(@mean, M_cell, 'UniformOutput',0), 'UniformOutput',0), 'UniformOutput',0);

end

%Converting cell to matrix
M_low=cell2mat(M_cell);
%% OLD VERSION
%
% % Setting field of view FOV size
% FOV=[size(M,1)*voxeldim_high, size(M,2)*voxeldim_high, size(M,3)*voxeldim_high];
% no_elements=FOV(1)/voxeldim_low;
%
% reduction_factor=size(M,1)/no_elements;
% cell_shape=reduction_factor*ones(size(M,1)/reduction_factor,1);
%
% for k=1:size(M,3)
%     M_slice=M(:,:,k);
%     cell_slice=mat2cell(M_slice, cell_shape, cell_shape);
%     column_means = cellfun(@mean, cell_slice, 'UniformOutput',0);
%     total_means = cellfun(@mean, column_means, 'UniformOutput',0);
%     M_means(:,:,k)=cell2mat(total_means);
% end
%
% i=1;
% for k=1:size(M,3)/reduction_factor
%     M_low(:,:,k)=sum(M_means(:,:,(i:(i+reduction_factor-1))),3) /reduction_factor;
%     i=i+reduction_factor;
% end
%
% %% END

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function rimgs = slresizeimg(imgs, newsiz, interpker)
%SLRESIZEIMG Resizes the images by interpolation
%
% $ Syntax $
%   - rimgs = slresizeimg(imgs, newsiz, interpker)
%
% $ Arguments $
%   - imgs:         The set of images
%   - newsiz:       The new image size to be resized to
%                   It can be in two forms:
%                   - [new_height, new_width]
%                   - ratio to the original size
%   - interpker:    The interpolation kernel (default = 'linear')
%   - rimgs:        The resized images
%
% $ Description $
%   - rimgs = slresizeimg(imgs, newsiz, interpker) resizes the image
%     to new size by interpolation.
%
% $ Remarks $
%   - The implementation is based on slimginterp.
%
% $ History $
%   - Created by Dahua Lin, on Sep 3, 2006
%

%% parse and verify input

if nargin < 2
    raise_lackinput('slresizeimg', 2);
end
h0 = size(imgs, 1);
w0 = size(imgs, 2);

if isnumeric(newsiz)
    if length(newsiz) == 1
        h = newsiz * h0;
        w = newsiz * w0;
    elseif length(newsiz) == 2
        [h, w] = sltakeval(newsiz);
    else
        error('sltoolbox:invalidarg', 'The newsiz is invalid');
    end
else
    error('sltoolbox:invalidarg', 'The newsiz is invalid');
end

if nargin < 3 || isempty(interpker)
    interpker = 'linear';
end

%% generate coordinate map

I = linspace(1, h0, h)';
J = linspace(1, w0, w);
I = I(:, ones(1, w));
J = J(ones(h, 1), :);

%% Do interpolation

rimgs = slimginterp(imgs, I, J, interpker);

author: lmthang


text:
clear
        % Incarcarea datelor din tabel
I=0:0.5:5;
x=linspace(0,pi,length(I));
for k=1:length(I)
    yy=3*I(k)*sin(x);
    y(:,k)=yy';
end
        % Incarcarea tabelului de cautare
tab=zeros(12);
tab(2:12,1)=x';
tab(1,2:12)=I;
tab(2:12,2:12)=y;
        % Cautarea in tabel
a=table2(tab,pi/10,1)
b=table2(tab,0.5,0.75)

author: Sable

repo: mcbench-benchmarks


text:
% Numerical solution for wave propagation

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory in Section 9.3 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% Equation (9.34) is solved with the kappa scheme
% dy/dt + c dy/dx = 0
% Known exact solution f(x-ct)
% Boundary conditions using upwind schemes or extrapolation
% Uniform vertex-centred grid
% This program generates Figs. 9.6 and 9.7 in the book

% Functions called: exact_solution

		% ...............Input.........................................
kappa = 0;	% Parameter of kappa-scheme
c = 1;		% Velocity
sigma = 0.7;	% CFL number
J = 31;		% Number of grid cells
bc = 2;		% Governs choice of numerical boundary conditions; see page 350
		% Enter 1 for boundary conditions using upwind schemes
		%    or 2 for extrapolation using  equation (9.36)
T = 0.4;	% Final time
		% ..............End of input...................................

h = 1.0/(J-1);		% Size of grid cells
dt = sigma*h/c;		% Time step
n = floor(T/dt);	% Number of time steps

% 	Definition of grid numbering

%     x=0    			  x=1
% grid |------|------|------|------|
%      1      2      3             J

x = [0:J-1]'*h;		% Location of grid nodes

sig2 = sigma^2; sig4 = sigma/4;

ynew = zeros(J,1); yright = zeros(n+1,1);
t = 0;	yright(1) = exact_solution(t,1,c);  yold = exact_solution(t,x,c);
for i = 1:n,  t = t + dt;
  a = 2*yold(1) - yold(2);	% Extrapolated inflow value
  ynew(1) = exact_solution(t,0,c);
  if bc == 1
    ynew(2) =  yold(2) - sigma*(yold(2) - yold(1));
  else
    ynew(2) =  yold(2) +sig4*(-1+kappa + (1-3*kappa)*sigma +2*kappa*sig2)*a +...
    sig4*(5-3*kappa+(9*kappa-1)*sigma-6*kappa*sig2)*yold(1)...
    +sig4*(-3+3*kappa -(1+9*kappa)*sigma+6*kappa*sig2)*yold(2)+...
    sig4*(-1-kappa+(1+3*kappa)*sigma-2*kappa*sig2)*yold(3);
  end
  for j = 3:J-1
    ynew(j) = yold(j) +sig4*(-1+kappa + (1-3*kappa)*sigma +2*kappa*sig2)*...
    yold(j-2)+sig4*(5-3*kappa+(9*kappa-1)*sigma-6*kappa*sig2)*yold(j-1)...
    +sig4*(-3+3*kappa -(1+9*kappa)*sigma+6*kappa*sig2)*yold(j)+...
    sig4*(-1-kappa+(1+3*kappa)*sigma-2*kappa*sig2)*yold(j+1);
  end
  if bc == 1		% Fully upwind scheme at outflow
    ynew(J) = (1-1.5*sigma+0.5*sig2)*yold(J) + sigma*(2-sigma)*yold(J-1)...
    +sigma*0.5*(sigma-1)*yold(J-2);
  else			% Extrapolation at outflow
      ynew(J) = yold(J) +sig4*(-1+kappa + (1-3*kappa)*sigma +2*kappa*sig2)*...
    yold(J-2)+sig4*(5-3*kappa+(9*kappa-1)*sigma-6*kappa*sig2)*yold(J-1)...
    +sig4*(-3+3*kappa -(1+9*kappa)*sigma+6*kappa*sig2)*yold(J)+...
      sig4*(-1-kappa+(1+3*kappa)*sigma-2*kappa*sig2)*(2*yold(J) - yold(J-1));
  end
  yright(i+1) = ynew(J);  yold = ynew;
end

yexact =  exact_solution(t,x,c);	error = ynew - yexact;
norm1 = norm(error,1)/J			% Compute error norms
norm2 = norm(error,2)/sqrt(J)
norminf = norm(error,inf)

figure(1), clf, hold on
xx = 0:0.005:1;	plot (xx,exact_solution(t,xx,c),'-'); plot (x,ynew,'o');
s1 = ['kappa=',num2str(kappa),' bc =', num2str(bc),' sigma=',...
	num2str(sigma),' cells=',num2str(J-1),' t=',num2str(t)];
title(s1,'fontsize',18);

repo: Aero-matlab


text:
%
% Original source
% http://www.mat.univie.ac.at/~neum/glopt/coconut/Benchmark/Library1/alkyl.mod
%
function test_alkyl

  addpath('../lib') ;

  auxdata = {} ;
  options.lb = [ 0, 0,   0,   0, 0, 0.85, 0.9,  3,  1.2, 1.45, 0.99, 0.99, 0.9, 0.99 ] ;
  options.ub = [ 2, 1.6, 1.2, 5, 2, 0.93, 0.95, 12, 4,   1.62, 1.01010101010101, 1.01010101010101, 1.11111111111111, 1.01010101010101 ] ;

  % The constraint functions are bounded to zero
  options.cl = zeros(1,7); %  constraints
  options.cu = zeros(1,7);

  % Set up the auxiliary data.
  options.auxdata = auxdata ;

  % Set the IPOPT options.
  options.ipopt.jac_d_constant   = 'no';
  options.ipopt.hessian_constant = 'no';
  options.ipopt.mu_strategy      = 'adaptive';
  options.ipopt.max_iter         = 400;
  options.ipopt.tol              = 1e-10;
  options.ipopt.linear_solver    = 'ma57';
  %options.ipopt.linear_solver    = 'mumps';
  %options.ipopt.linear_solver    = 'pardiso';

  % The callback functions.
  funcs.objective         = @objective;
  funcs.constraints       = @constraints;
  funcs.gradient          = @gradient;
  funcs.jacobian          = @jacobian;
  funcs.jacobianstructure = @jacobianstructure;
  if true
    funcs.hessian           = @hessian;
    funcs.hessianstructure  = @hessianstructure;
    options.ipopt.derivative_test = 'second-order';
  else
    options.ipopt.hessian_approximation      = 'limited-memory';
    options.ipopt.limited_memory_update_type = 'bfgs' ; % {bfgs}, sr1 = 6; % {6}
    options.ipopt.derivative_test            = 'first-order';
  end

  % Run IPOPT.
  x0 = [ 1.745,   ...
         1.2,     ...
         1.1,     ...
         3.048,   ...
         1.974,   ...
         0.893,   ...
         0.928,   ...
         8,       ...
         3.6,     ...
         1.45,    ...
         1,       ...
         1,       ...
         1,       ...
         1 ] ;

  tic
  [x, info] = ipopt_auxdata(x0,funcs,options);
  elapsed = toc ;

  info;

  x

end

%%
% map the indices with the corresponding index in the spase matrix
function f = objective( x, auxdata )
  x2  = x(1) ;
  x3  = x(2) ;
  x4  = x(3) ;
  x5  = x(4) ;
  x6  = x(5) ;
  x7  = x(6) ;
  x8  = x(7) ;
  x9  = x(8) ;
  x10 = x(9) ;
  x11 = x(10) ;
  x12 = x(11) ;
  x13 = x(12) ;
  x14 = x(13) ;
  x15 = x(14) ;
  f   = - 6.3*x5*x8 + 5.04*x2 + 0.35*x3 + x4 + 3.36*x6 ;
end

%%
% map the indices with the corresponding index in the spase matrix
function g = gradient(x,auxdata)
  x2  = x(1) ;
  x3  = x(2) ;
  x4  = x(3) ;
  x5  = x(4) ;
  x6  = x(5) ;
  x7  = x(6) ;
  x8  = x(7) ;
  x9  = x(8) ;
  x10 = x(9) ;
  x11 = x(10) ;
  x12 = x(11) ;
  x13 = x(12) ;
  x14 = x(13) ;
  x15 = x(14) ;
  g   = [ 5.04, 0.35, 1, -6.3*x8, 3.36, 0, -6.3*x5, 0, 0, 0, 0, 0, 0, 0 ] ;
end

function c = constraints(x,auxdata)
  x2  = x(1) ;
  x3  = x(2) ;
  x4  = x(3) ;
  x5  = x(4) ;
  x6  = x(5) ;
  x7  = x(6) ;
  x8  = x(7) ;
  x9  = x(8) ;
  x10 = x(9) ;
  x11 = x(10) ;
  x12 = x(11) ;
  x13 = x(12) ;
  x14 = x(13) ;
  x15 = x(14) ;
  e2  =  - 0.819672131147541*x2 + x5 - 0.819672131147541*x6 ;
  e3  = 0.98*x4 - x7*(0.01*x5*x10 + x4) ;
  e4  =  - x2*x9 + 10*x3 + x6 ;
  e5  = x5*x12 - x2*(1.12 + 0.13167*x9 - 0.0067*x9*x9) ;
  e6  = x8*x13 - 0.01*(1.098*x9 - 0.038*x9*x9) - 0.325*x7 - 0.57425;
  e7  = x10*x14 + 22.2*x11 - 35.82;
  e8  = x11*x15 - 3*x8 + 1.33 ;
  c   = [e2;e3;e4;e5;e6;e7;e8] ;
end

function jac = jacobian(x,auxdata)
  x2  = x(1) ;
  x3  = x(2) ;
  x4  = x(3) ;
  x5  = x(4) ;
  x6  = x(5) ;
  x7  = x(6) ;
  x8  = x(7) ;
  x9  = x(8) ;
  x10 = x(9) ;
  x11 = x(10) ;
  x12 = x(11) ;
  x13 = x(12) ;
  x14 = x(13) ;
  x15 = x(14) ;

  jac = sparse(7,14) ;

  jac(1,1) = -0.819672131147541 ;
  jac(1,4) = 1 ;
  jac(1,5) = -0.819672131147541 ;

  jac(2,3) = 0.98 - x7 ;
  jac(2,4) = -x7*(0.01*x10);
  jac(2,6) = -(0.01*x5*x10 + x4);
  jac(2,9) = - x7*(0.01*x5);

  jac(3,1) = -x9 ;
  jac(3,2) = 10 ;
  jac(3,5) = 1 ;
  jac(3,8) = -x2 ;

  jac(4,1)  = -(1.12 + 0.13167*x9 - 0.0067*x9*x9) ;
  jac(4,4)  = x12 ;
  jac(4,8)  = -x2*(0.13167 - 2*0.0067*x9) ;
  jac(4,11) = x5 ;

  jac(5,6)  = -0.325;
  jac(5,7)  = x13 ;
  jac(5,8)  = - 0.01*(1.098 - 2*0.038*x9) ;
  jac(5,12) = x8 ;

  jac(6,9)  = x14;
  jac(6,10) = 22.2;
  jac(6,13) = x10;

  jac(7,7)  = -3 ;
  jac(7,10) = x15 ;
  jac(7,14) = x11 ;

end

function jac = jacobianstructure(auxdata)
  jac = sparse(7,14) ;

  jac(1,1) = 1 ;
  jac(1,4) = 1 ;
  jac(1,5) = 1 ;

  jac(2,3) = 1 ;
  jac(2,4) = 1 ;
  jac(2,6) = 1 ;
  jac(2,9) = 1 ;

  jac(3,1) = 1 ;
  jac(3,2) = 1 ;
  jac(3,5) = 1 ;
  jac(3,8) = 1 ;

  jac(4,1)  = 1 ;
  jac(4,4)  = 1 ;
  jac(4,8)  = 1 ;
  jac(4,11) = 1 ;

  jac(5,6)  = 1 ;
  jac(5,7)  = 1 ;
  jac(5,8)  = 1 ;
  jac(5,12) = 1 ;

  jac(6,9)  = 1 ;
  jac(6,10) = 1 ;
  jac(6,13) = 1 ;

  jac(7,7)  = 1 ;
  jac(7,10) = 1 ;
  jac(7,14) = 1 ;
end

function H = hessian(x, sigma, L, auxdata)

  x2  = x(1) ;
  x3  = x(2) ;
  x4  = x(3) ;
  x5  = x(4) ;
  x6  = x(5) ;
  x7  = x(6) ;
  x8  = x(7) ;
  x9  = x(8) ;
  x10 = x(9) ;
  x11 = x(10) ;
  x12 = x(11) ;
  x13 = x(12) ;
  x14 = x(13) ;
  x15 = x(14) ;

  Hf = sparse(14,14) ;
  Hf(4,7) = -6.3 ;
  Hf(7,4) = -6.3 ;

  H2 = sparse(14,14) ;
  H2(3,6) = -1 ;
  H2(6,3) = -1;
  H2(4,6) = -(0.01*x10) ;
  H2(6,4) = -(0.01*x10);
  H2(4,9) = -x7*0.01 ;
  H2(9,4) = -x7*0.01 ;
  H2(6,9) = -(0.01*x5);
  H2(9,6) = -(0.01*x5);

  H3 = sparse(14,14) ;
  H3(1,8)  = -1 ;
  H3(8,1)  = -1 ;

  H4 = sparse(14,14) ;
  H4(1,8)  = -(0.13167 - 2*0.0067*x9) ;
  H4(8,1)  = -(0.13167 - 2*0.0067*x9) ;
  H4(4,11) = 1 ;
  H4(11,4) = 1 ;
  H4(8,8)  = x2*2*0.0067 ;

  H5 = sparse(14,14) ;
  H5(7,12) = 1 ;
  H5(12,7) = 1 ;
  H5(8,8)  = 0.01*(2*0.038) ;

  H6 = sparse(14,14) ;
  H6(9,13) = 1;
  H6(13,9) = 1;

  H7 = sparse(14,14) ;
  H7(10,14) = 1 ;
  H7(14,10) = 1 ;

  H  = tril(sparse(sigma * Hf + L(2)*H2 + L(3)*H3 + L(4)*H4 + L(5)*H5 + L(6)*H6 + L(7)*H7)) ;
end

function H = hessianstructure(auxdata)

  Hf = sparse(14,14) ;
  Hf(4,7) = 1;
  Hf(7,4) = 1;

  H2 = sparse(14,14) ;
  H2(3,6) = 1;
  H2(6,3) = 1;
  H2(4,6) = 1;
  H2(6,4) = 1;
  H2(4,9) = 1;
  H2(9,4) = 1;
  H2(6,9) = 1;
  H2(9,6) = 1;

  H3 = sparse(14,14) ;
  H3(1,8) = 1;
  H3(8,1) = 1;

  H4 = sparse(14,14) ;
  H4(1,8)  = 1;
  H4(8,1)  = 1;
  H4(4,11) = 1;
  H4(11,4) = 1;
  H4(8,8)  = 1;

  H5 = sparse(14,14) ;
  H5(7,12) = 1;
  H5(12,7) = 1;
  H5(8,8)  = 1;

  H6 = sparse(14,14) ;
  H6(9,13) = 1;
  H6(13,9) = 1;

  H7 = sparse(14,14) ;
  H7(10,14) = 1 ;
  H7(14,10) = 1 ;

  H = tril(sparse(Hf+H2+H3+H4+H5+H6+H7)) ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [timelock] = ft_datatype_timelock(timelock, varargin)

% FT_DATATYPE_TIMELOCK describes the FieldTrip MATLAB structure for timelock data
%
% The timelock data structure represents averaged or non-averaged event-releted
% potentials (ERPs, in case of EEG) or ERFs (in case of MEG). This data structure is
% usually generated with the FT_TIMELOCKANALYSIS or FT_TIMELOCKGRANDAVERAGE function.
%
% An example of a timelock structure containing the ERF for 151 channels MEG data is
%
%     dimord: 'chan_time'       defines how the numeric data should be interpreted
%        avg: [151x600 double]  the average values of the activity for 151 channels x 600 timepoints
%        var: [151x600 double]  the variance of the activity for 151 channels x 600 timepoints
%      label: {151x1 cell}      the channel labels (e.g. 'MRC13')
%       time: [1x600 double]    the timepoints in seconds
%       grad: [1x1 struct]      information about the sensor array (for EEG data it is called elec)
%        cfg: [1x1 struct]      the configuration used by the function that generated this data structure
%
% Required fields:
%   - label, dimord, time
%
% Optional fields:
%   - avg, var, dof, cov, trial, trialinfo, sampleinfo, grad, elec, opto, cfg
%
% Deprecated fields:
%   - <none>
%
% Obsoleted fields:
%   - fsample
%
% Revision history:
%
% (2017/latest) The data structure cannot contain an average and simultaneously single
% trial information any more, i.e. avg/var/dof and trial/individual are mutually exclusive.
%
% (2011v2) The description of the sensors has changed, see FT_DATATYPE_SENS
% for further information.
%
% (2011) The field 'fsample' was removed, as it was redundant.
%
% (2003) The initial version was defined.
%
% See also FT_DATATYPE, FT_DATATYPE_COMP, FT_DATATYPE_FREQ, FT_DATATYPE_RAW

% Copyright (C) 2011, 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$

% get the optional input arguments, which should be specified as key-value pairs
version       = ft_getopt(varargin, 'version', 'latest');
hassampleinfo = ft_getopt(varargin, 'hassampleinfo', 'ifmakessense'); % can be yes/no/ifmakessense
hastrialinfo  = ft_getopt(varargin, 'hastrialinfo',  'ifmakessense'); % can be yes/no/ifmakessense

% convert it into true/false
if isequal(hassampleinfo, 'ifmakessense')
  hassampleinfo = makessense(timelock, 'sampleinfo');
else
  hassampleinfo = istrue(hassampleinfo);
end
if isequal(hastrialinfo, 'ifmakessense')
  hastrialinfo = makessense(timelock, 'trialinfo');
else
  hastrialinfo = istrue(hastrialinfo);
end

if strcmp(version, 'latest')
  version = '2017';
end

if isempty(timelock)
  return;
end

% do some sanity checks
assert(isfield(timelock, 'time') && isfield(timelock, 'label'), 'inconsistent timelock data structure, some field is missing');
assert(length(unique(timelock.label))==length(timelock.label), 'channel labels must be unique');

% ensure consistency between the dimord string and the axes that describe the data dimensions
timelock = fixdimord(timelock);

% remove these very obsolete fields, it is unclear when they were precisely used
if isfield(timelock, 'numsamples'),       timelock = rmfield(timelock, 'numsamples');       end
if isfield(timelock, 'numcovsamples'),    timelock = rmfield(timelock, 'numcovsamples');    end
if isfield(timelock, 'numblcovsamples'),  timelock = rmfield(timelock, 'numblcovsamples');  end

if ~iscolumn(timelock.label)
  timelock.label = timelock.label';
end
if ~isrow(timelock.time)
  timelock.time = timelock.time';
end
if ~isfield(timelock, 'label')
  ft_warning('data structure is incorrect since it has no channel labels');
end

switch version
  case '2017'
    % ensure that the sensor structures are up to date
    if isfield(timelock, 'grad')
      timelock.grad = ft_datatype_sens(timelock.grad);
    end
    if isfield(timelock, 'elec')
      timelock.elec = ft_datatype_sens(timelock.elec);
    end
    if isfield(timelock, 'opto')
      timelock.opto = ft_datatype_sens(timelock.opto);
    end

    fn = fieldnames(timelock);
    fn = setdiff(fn, ignorefields('appendtimelock'));
    fn = fn(~endsWith(fn, 'dimord'));
    dimord = cell(size(fn));
    hasrpt = false(size(fn));
    for i=1:numel(fn)
      dimord{i} = getdimord(timelock, fn{i});
      hasrpt(i) = ~isempty(strfind(dimord{i}, 'rpt')) || ~isempty(strfind(dimord{i}, 'subj'));
    end
    if any(hasrpt) && ~all(hasrpt)
      ft_warning('timelock structure contains field with and without repetitions');
      str = sprintf('%s, ', fn{hasrpt});
      str = str(1:end-2);
      ft_notice('selecting these fields that have repetitions: %s', str);
      str = sprintf('%s, ', fn{~hasrpt});
      str = str(1:end-2);
      ft_notice('removing these fields that do not have repetitions: %s', str);
      timelock = removefields(timelock, fn(~hasrpt));
      if isfield(timelock, 'dimord') && ~ismember(timelock.dimord, dimord(hasrpt))
        % the dimord does not apply to any of the existing fields any more
        timelock = rmfield(timelock, 'dimord');
        timelock = fixdimord(timelock);
      end
    end

    if (hassampleinfo && ~isfield(timelock, 'sampleinfo')) || (hastrialinfo && ~isfield(timelock, 'trialinfo'))
      % try to reconstruct the sampleinfo and trialinfo
      timelock = fixsampleinfo(timelock);
    end

    if ~hassampleinfo && isfield(timelock, 'sampleinfo')
      timelock = rmfield(timelock, 'sampleinfo');
    end

    if ~hastrialinfo && isfield(timelock, 'trialinfo')
      timelock = rmfield(timelock, 'trialinfo');
    end

    % this field can be present in raw data, but is not desired in timelock data
    timelock = removefields(timelock, {'fsample'});

    % ensure that the structure has all required fields
    % note that dimord is listed as required field, but it might also be xxxdimord, or dynamically determined with GETDIMORD
    for required={'label' 'time'}
      assert(isfield(timelock, required), 'required field "%s" is missing', required{:});
    end

  case '2011v2'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ensure that the sensor structures are up to date
    if isfield(timelock, 'grad')
      timelock.grad = ft_datatype_sens(timelock.grad);
    end
    if isfield(timelock, 'elec')
      timelock.elec = ft_datatype_sens(timelock.elec);
    end
    if isfield(timelock, 'opto')
      timelock.opto = ft_datatype_sens(timelock.opto);
    end

    % these fields can be present in raw data, but are not desired in timelock data
    timelock = removefields(timelock, {'sampleinfo', 'fsample'});

  case '2003'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % there are no known conversions for backward or forward compatibility support

  otherwise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_error('unsupported version "%s" for timelock datatype', version);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Channels,FileInfo]=importfamos(FullRawData)
%__________________________________________________________________________
% The sequnce importfamos.m was produced to convert imc raw data(*.raw;
% *.dat) to matlab data. Here, struct is used to manage the channel
% information.
%
% For more information of FAMOS file format, please see the
% manufacturer's website: http://www.imc-berlin.de
%
% Corresponding to Data Structure in Matlab, the channels are stored
% as struct struct: Channel + name (channel name)
%                           + comment
%                           + data (channel Value)
%                           + length
%                           + yUnit
%                           + t0
%                           + dt(sampling period)
%                           + xUnit
% Version history:
% Version 1 (2011.1.19); Current version only can deal with analog rawdata
% from imc devices. The digital and group function is ongoning and will be
% released in version 2.
%%-------------------------------------------------------------------------
% Author: Liang
% Danke.Liang@gmail.com
% Started on Dec.14, 2010
%__________________________________________________________________________
if nargin == 0;
    %[FamosName, FamosPath] = uigetfile({'*.raw';'*.dat'},'Select Famos Raw data');
   FullRawData = uigetfile({'*.raw';'*.dat'},'Select Famos Raw data');
    %pwd current path
    %filesep e.g "/" under windows os
%elseif nargin == 1;
    %FamosPath = [pwd filesep];
end

fid = fopen(FullRawData,'r');
KeyFlag={'|CF';'|Ck';'|NO';'|CT';'|CB';'|CI';'|CG';'|CD';'|NT';'|CZ';'|CC';'|CP';'|Cb';'|CR';'|ND';'|CN';'|CS';'|NU'};

disp(['Info: Read data @' datestr(now) ' from' ])
% disp([ ' ' FamosPath FamosName '...'])
disp([ ' ' FullRawData '...'])
%initiation
BinaryStart=0;
Str3='';
%while
BinaryStart=FindBinaryStart(fid);
ChannelID=1;
Channel=struct('name','','comment','','data',[],'length',0,'yUnit','','t0','','dt','','xUnit','');
FileInfo=struct('Origin','','Company','','Comment','','ChannelNums',[]);

%Struct TriggerTime
FlagCS='|CS';
fseek(fid, 1, 'bof');
 while strcmp(Str3,FlagCS)==0
     lastlocation=ftell(fid);
     tline=fgetl(fid);
     Str3=tline(1:3);
     switch (Str3)
         case KeyFlag(1)
         CFfini=ProcessCFandCK(tline);
         if CFfini==0
         Str3='|CS'; %force to exit while loop
         end
         case KeyFlag(3);
        [NOKeyOrigin,NOKeyCompany,NOKeyComment]=ProcessNO(tline);
            FileInfo.Origin=NOKeyOrigin;
            FileInfo.Company=NOKeyCompany;
            FileInfo.Comment=NOKeyComment;
         case KeyFlag(4);
         ProcessCT(tline);
          case KeyFlag(5);
         ProcessCBbig(tline);
         case KeyFlag(6);
         ProcessCI(tline);
         case KeyFlag(7);
         [CGNumberComponents,CGFieldType,CGDimension]=ProcessCG(tline);
         case KeyFlag(8);
         [CDSample, CDUnit]=ProcessCD(tline);
         case KeyFlag(9);
         TriggerTime = ProcessNT(tline);
         case KeyFlag(10);
         ProcessCZ(tline);
         case KeyFlag(11);
         KeyComponentIndex= ProcessCC(tline);
         case KeyFlag(12);
         [KyBufferRef,KeyBytes,KeyNumberFormat,KeySignBits]=ProcessCP(tline);
         case KeyFlag(13);
         [KeyBufferRefIndex,KeyBufferRefCb,KeyOffsetBufferInSamplesKey,KeyBufferFilledBytes]=ProcessCblittle(tline);
         case KeyFlag(14);
         [KeyTransformation,KeyCRfactor,KeyCRoffset,KeyUnit]=ProcessCR(tline);
         case KeyFlag(15);
         ProcessND();
         case KeyFlag(16);
             CurrentLocation=ftell(fid);
             [ChannelName,ChannelComment]=ProcessCN(tline);
             %Read Data from Buffer
             BinaryRead= BinaryStart+KeyOffsetBufferInSamplesKey;
             ChannelLength=KeyBufferFilledBytes*8/KeySignBits;
             TempChannel=ReadChannel(fid,BinaryRead,ChannelLength,KeyNumberFormat,KeyCRfactor,KeyCRoffset);
             Channel(ChannelID).name=ChannelName;
             Channel(ChannelID).comment=ChannelComment;
             Channel(ChannelID).length=ChannelLength;
             Channel(ChannelID).data=TempChannel;
             Channel(ChannelID).yUnit=KeyUnit;
             Channel(ChannelID).t0=TriggerTime;
             Channel(ChannelID).dt=CDSample;
             Channel(ChannelID).xUnit=CDUnit;
%              plot(TempChannel)
             clear TempChannel;
             ChannelID=ChannelID+1;
         fseek(fid,CurrentLocation,'bof');
         case KeyFlag(17);
         ProcessCS();
         case KeyFlag(18);
         ProcessNU();
     end
 end
fclose(fid);
FileInfo.ChannelNums=ChannelID-1;
Channels=Channel;
end

function BinaryStart=FindBinaryStart(fileid)
Str3='';
while not(strcmp(Str3,'|CS'))
LastLocation=ftell(fileid);
  tline=fgetl(fileid);
  Str3=tline(1:3);
end
CommaLocation=find(tline==',');
BinaryStart=CommaLocation(4)+LastLocation;

end
%--------------------------------------------------------------------------
%%
% CF,2,1,1;
%  Dateiformat,Prozessor
%  Dateiformat = 2, Keyl?nge =1, Prozessor = 1
%--------------------------------------------------------------------------
function  Finished=ProcessCFandCK(TxtString)
disp('Info: Processing key CF and CK...');
FileFormat=TxtString(5);
KeyLength=TxtString(7);
temp=TxtString(21);
Finished=str2double(temp);
if Finished==0
    msgbox('the rawdata is not finished storage','incorrect data','error');
end
disp('Info: Finished Processing key CF and CK!')
end

%-------------------------------------------------------------------------------
% CK,1,3,1,Abgeschlossen;
%-------------------------------------------------------------------------------
function  [keyVersion,keyLength,keyData]=ProcessCK(keyVersion,keyLength,keyData)
disp('Info: Processing key CK...')
if (keyLength ~= 3)
    msgbox('ERROR: CKs key length should always be 3.');
end
if (dataClosed == 0)
msgbox('WARNING: CK key indicates not all data was written');
msgbox('done processing key');
end
disp('Info: Finished Processing key CK!')
end
%-------------------------------------------------------------------------------
% NO,1,KeyLang,Ursprung,NameLang,Name, KommLang,Kommentar;
%-------------------------------------------------------------------------------
function [KeyOrigin,KeyCompany,KeyComment] = ProcessNO(TxtString)
%disp('Info: Processing key NO...');

CommaLocation=find(TxtString==',');
KeyOrigin=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
KeyCompany=TxtString(CommaLocation(5)+1:CommaLocation(6)-1);
KeyCommLength=TxtString(CommaLocation(6)+1:CommaLocation(7)-1);
if KeyCommLength=='0';
KeyComment='';
else
    temp=str2double(KeyCommLength);
    KeyComment=TxtString(CommaLocation(7)+1,CommaLocation(7)+temp);
end

%disp('Info: Finished Processing key NO!')
end

% '-------------------------------------------------------------------------------
% ' CB,1,KeyLang,IndexGruppe,NameLang,Name,KommLang,Kommentar;
% '-------------------------------------------------------------------------------
function tempCBbig= ProcessCBbig(TxtString)
% This function is relatvie to Group function
tempCBbig=0;
end

% '-------------------------------------------------------------------------------
% ' CT,1,KeyLang,IndexGruppe, NameLang,Name,TextLang,Text, CommLeng,Comment;
% '
% '  NOTE: The use of the word 'name' here is misleading.  This key contains extra
% '  information about a given channel group, not the channel group's actual name.
% '  'name', as used by this key, refers to the name of the attribute which
% '  will appear in the channel group's properties.
% '-------------------------------------------------------------------------------
function ProcessCT(TxtString)
% function is relative to Group
end

% '-------------------------------------------------------------------------------
% ' CG,1,KeyLang,NumberComponent,Fieldtype,Dimension;
% '-------------------------------------------------------------------------------
function [KeyNumnerCommponents,KeyFieldType,KeyDimension]=ProcessCG(TxtString)
% 	TracePrint(" Processing key CG...")
%disp('Info: Processing key CG...');
global CGCount;
CGCount=CGCount+1;
global CGFlag;
CGFlag=1;
% 	Dim paramList : paramList = Split(keyData,",")
CommaLocation=find(TxtString==',');

    Txtemp=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
    KeyNumnerCommponents=str2num(Txtemp);
    Txtemp=TxtString(CommaLocation(4)+1:CommaLocation(5)-1);
    KeyFieldType=str2num(Txtemp);
    Txtemp=TxtString(CommaLocation(5)+1);
    KeyDimension=str2num(Txtemp);
%disp('Info: Finished Processing key CG!')
end

% '-------------------------------------------------------------------------------
% ' CD,1,KeyLength,dx,kalibriert,EinheitLang, Einheit, 0,0,0;
% ' CD,2,KeyLang,dx,kalibriert,EinheitLang, Einheit, Reduktion, IsMultiEvents,
% '  SortiereBuffer, x0, PretriggerVerwendung;
% '-------------------------------------------------------------------------------
function [KeyDx, KeyUnit]= ProcessCD(TxtString)
% disp('Info: Processing key CD...');

CommaLocation=find(TxtString==',');
Txtemp=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
KeyDx=str2double(Txtemp);
KeyUnit=TxtString(CommaLocation(6)+1:CommaLocation(7)-1);
% disp('Info: Finished Process key CD!');
end

% '-------------------------------------------------------------------------------
% ' NT,1,KeyLang,Tag,Monat,Jahr,Stunden,Minuten,Sekunden;
% '-------------------------------------------------------------------------------
function TimeStart = ProcessNT(TxtString)

    CommaLocation=find(TxtString==',');
    Txtemp=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
    KeyDay=str2num(Txtemp);
    Txtemp=TxtString(CommaLocation(4)+1:CommaLocation(5)-1);
    KeyMonth=str2num(Txtemp);
    Txtemp=TxtString(CommaLocation(5)+1:CommaLocation(6)-1);
    KeyYear=str2num(Txtemp);

    Txtemp=TxtString(CommaLocation(6)+1:CommaLocation(7)-1);
    KeyHours=str2num(Txtemp);
    Txtemp=TxtString(CommaLocation(7)+1:CommaLocation(8)-1);
    KeyMinutes=str2num(Txtemp);
    Txtemp=TxtString(CommaLocation(8)+1:length(TxtString));
    KeySeconds=str2num(Txtemp);

    TimeStart=datestr(datenum([KeyYear, KeyMonth,KeyDay,KeyHours,KeyMinutes,KeySeconds]),'yyyy-mm-dd HH:MM:SS');
% disp('Info: Finished Processing key NT!');
end

% '-------------------------------------------------------------------------------
% ' CZ,1,KeyLang,dz,dzkali,z0,z0kali,EinheitLang,Einheit,SegmentLang;
% '-------------------------------------------------------------------------------
function ProcessCZ(TxtString)
% 	TracePrint(" Processing key CZ...")
%
%   TracePrint(" ignoring key")
end

% '-------------------------------------------------------------------------------
% ' CC,1,KeyLang,KomponentenIndex,AnalogDigital;
% '
% ' Indicates a new channel.
% '-------------------------------------------------------------------------------
function KeyComponentIndex= ProcessCC(TxtString)
% 	TracePrint(" Processing key CC...")
% disp('Info: Processing key CC...');
% 	Dim paramList : paramList = Split(keyData,",")
%   Dim chan : Set chan = currentChanGroup.AddChan()
CommaLocation=find(TxtString==',');
Txtemp=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
KeyComponentIndex=str2num(Txtemp);
end

% '-------------------------------------------------------------------------------
% ' CP,1,KeyLang,BufferReferenz,Bytes,Zahlenformat,SignBits,Maske,Offset,
% '  DirekteFolgeAnzahl,AbstandBytes;
% '-------------------------------------------------------------------------------
function [KeyBufferRef,KeyBytes,KeyNumerFormat,KeySignBits]= ProcessCP(TxtString)
% disp('Info: Processing key CP...');
CommaLocation=find(TxtString==',');
Txtemp=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
KeyBufferRef=str2num(Txtemp);
Txtemp=TxtString(CommaLocation(4)+1:CommaLocation(5)-1);
KeyBytes=str2num(Txtemp);
Txtemp=TxtString(CommaLocation(5)+1:CommaLocation(6)-1);
NumberFormat=str2num(Txtemp);
switch (NumberFormat)
    case 1
        KeyNumerFormat='*uint';
    case 2
        KeyNumerFormat='*int';
    case 3
        KeyNumerFormat='*ushort';
    case 4
        KeyNumerFormat='*short';
    case 5
        KeyNumerFormat='*ulong';
    case 6
        KeyNumerFormat='*long';
    case 7
        KeyNumerFormat='*float';
    case 8
        KeyNumerFormat='*float32';
    case 9
        KeyNumerFormat='*'; % imc Device Transitional Recording
    case 10
        KeyNumerFormat='*TimeStampASII' ;% TimeStamp is famos type
    case 11
        KeyNumberFormat='*bit16'; %2-byte-word digital
    case 13
         KeyNumberFormat='*bit48';
end
Txtemp=TxtString(CommaLocation(6)+1:CommaLocation(7)-1);
KeySignBits=str2num(Txtemp);
end
% '-------------------------------------------------------------------------------
% ' Cb,1,KeyLang,AnzahlBufferInKey,BytesInUserInfo,BufferReferenz, IndexSamplesKey,
% '  OffsetBufferInSamplesKey, BufferLangBytes,OffsetFirstSampleInBuffer,
% '  BufferFilledBytes, 0, X0, AddZeit, UserInfo, [BufferReferenz, ... ];
% '
% ' NOTE: This description above is taken from the specs, but all examples use a
% '  format as given below, with a zero after AddZeit.
% ' Cb,1,KeyLang,AnzahlBufferInKey,BytesInUserInfo,BufferReferenz, IndexSamplesKey,
% '  OffsetBufferInSamplesKey, BufferLangBytes,OffsetFirstSampleInBuffer,
% '  BufferFilledBytes, X0, AddZeit, 0, UserInfo, [BufferReferenz, ... ];
% '-------------------------------------------------------------------------------
function [KeyBufferRefIndex,KeyBufferRefCb,KeyOffsetBufferInSamplesKey,KeyBufferFilledBytes] = ProcessCblittle(TxtString)
% disp('Info: Processing key Cb...');
CommaLocation=find(TxtString==',');
Txtemp=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
KeyBufferRefIndex=str2num(Txtemp);
Txtemp=TxtString(CommaLocation(5)+1:CommaLocation(6)-1);
KeyBufferRefCb=str2double(Txtemp);
Txtemp=TxtString(CommaLocation(7)+1:CommaLocation(8)-1);
KeyOffsetBufferInSamplesKey=str2double(Txtemp);
Txtemp=TxtString(CommaLocation(10)+1:CommaLocation(11)-1);
KeyBufferFilledBytes=str2double(Txtemp);
%disp('Info: Finished Processing key Cb!');
end

% '-------------------------------------------------------------------------------
% ' CR,1,KeyLang,Transformieren,Faktor,Offset,Kalibriert,EinheitLang, Einheit;
% '-------------------------------------------------------------------------------
function [KeyTransformation,KeyCRfactor,KeyCRoffset,KeyUnit]= ProcessCR(TxtString)
%   disp('Info: Processing key CR...')
%
CommaLocation=find(TxtString==',');
Txtemp=TxtString(CommaLocation(3)+1:CommaLocation(4)-1);
KeyTransformation=str2num(Txtemp);
Txtemp=TxtString(CommaLocation(4)+1:CommaLocation(5)-1);
KeyCRfactor=str2double(Txtemp);
Txtemp=TxtString(CommaLocation(5)+1:CommaLocation(6)-1);
KeyCRoffset=str2double(Txtemp);
Txtemp=TxtString(CommaLocation(7)+1:CommaLocation(8)-1);
KeyUnitLength=str2double(Txtemp);
KeyUnit=TxtString(CommaLocation(8)+1:CommaLocation(8)+KeyUnitLength);

%   disp('Info: Finished Processing key CR!');
end

% '-------------------------------------------------------------------------------
% ' ND,1,KeyLang,FarbeR,FarbeG,FarbeB,yMin,yMax;
% '-------------------------------------------------------------------------------
function  ProcessND(TxtString)
% this function is relative to color and MaxMin value of channel
end

% '-------------------------------------------------------------------------------
% ' CN,1,KeyLang,IndexGruppe,0,IndexBit,NameLang,Name,KommLang,Kommentar;
% '-------------------------------------------------------------------------------
function [ChannelName,ChannelComment]= ProcessCN(TxtString)
% 	TracePrint(" Processing key CN...")
%disp('Info: Processing key CR...');
CommaLocation=find(TxtString==',');
ChannelName=TxtString(CommaLocation(7)+1:CommaLocation(8)-1);
%   NOTE: It is unclear why this key has an index group, as it seems that the
%   index group in the CB and CT keys are channel group numbers and this key
%   relates only to a single channel
ChannelCommLength=TxtString(CommaLocation(8)+1:CommaLocation(9)-1);
if ChannelCommLength=='0';
ChannelComment='';
else
    temp=str2double(ChannelCommLength);
    ChnannelComment=TxtString(CommaLocation(9)+1,CommaLocation(9)+temp);
end
%disp('Info: Finished Processing key CN!');
end

function tempCS= ProcessCS(TxtString)
% 	TracePrint(" Processing key CS...")
%
%   globalFile.Formatter.Delimiters  = ","
%   Dim dataSampleIndex : dataSampleIndex = globalFile.GetNextStringValue(defaultIntType)
%   Dim filePositionKeyStart : filePositionKeyStart = keyData
%   Dim filePositionDataStart : filePositionDataStart = globalFile.Position
%   'file position should now be at beginning of data
%
%   Dim currentChannelGroup
%   Dim currentChannel
%   Dim chunkLength,chanLength
%   Dim chan, chanGroup
%   Dim wordType, bytesPerWord
%   'If the channels are intermixed (channelsMixedTogether = true), one BinaryBlock
%   ' must be created to contain all the direct access channels
%   Dim binBlock : Set binBlock = globalFile.GetBinaryBlock()
%
%   Call DetermineExactFileFormat(filePositionKeyStart,filePositionDataStart,keyLength,avlFormat,singleBinaryBlock,samplesOddlyClumped,channelsMixedTogether)
%
%   For Each chanGroup In chanGroups
%     TracePrint(" looking at channelGroup: " & chanGroup.CBname)
%
%     'processing channel group...
%     If (Not chanGroup.processed) Then
%       TracePrint("  this channel group has never been processed; adding it")
%       Set currentChannelGroup = Root.ChannelGroups.Add(chanGroup.CBname)
%       Set chanGroup.channelGroup = currentChannelGroup
%       chanGroup.processed = true
%
%       'Fill in the channel group's meta-data
%       If (chanGroup.existsCT) Then
%         Call currentChannelGroup.Properties.Add(chanGroup.CTname,chanGroup.CTtext)
%         Call currentChannelGroup.Properties.Add(chanGroup.CTname & " comment",chanGroup.CTcomment)
%       End If
%       Call currentChannelGroup.Properties.Add("description",chanGroup.CBcomment)
%     Else
%       TracePrint("  this channel group has ALREADY been processed")
%       Set currentChannelGroup = chanGroup.channelGroup
%     End If
%     'currentChannelGroup has now been correctly set
%
%     For Each chan In chanGroup.chans
%       TracePrint(" looking at channel: " & chan.CNname)
%       If (chan.Cbindex = dataSampleIndex) Then
%         'this channel has data in this CS key
%         TracePrint("  this channel shares dataSampleIndex: " & dataSampleIndex)
%
%         'Add an implicit channel for each channel, if possible.  Sizes will
%         ' be updated at the end.
%         If ( (Not chan.processed) And chan.subGroup.existsCD) Then
%           Dim dummySize : dummySize = 10
%           Set currentChannel = currentChannelGroup.Channels.AddImplicitChannel(chan.CNname & " x-axis",chan.CbtimeOffset,chan.subGroup.CDdx,dummySize,defaultRealType)
%           Set chan.associatedImplicitChannel = currentChannel
%           Call currentChannel.Properties.Add("unit_string",chan.subGroup.CDunit)
%
%           TracePrint("  adding implicit channel, offset: " & chan.CbtimeOffset & " incr: " & chan.subGroup.CDdx & " length: " & chanLength)
%         End If
%
%         'Determine if direct access channels can be used.
%         If ( (Not forceStableDataRead) And singleBinaryBlock And (Not samplesOddlyClumped) ) Then
%           TracePrint("  using direct access channels")
%
%           If (chan.processed) Then
%             TracePrint("ERROR: direct access channels cannot process channels twice")
%             ErrorPrint("ERROR: The file format appears to be incorrect.  See log file for more details.")
%           End If
%           chan.processed = true
%
%           'chunkLength is measured in bytes.
%           chunkLength = (chan.CPnumberSamplesTogether * chan.CPbytesPerMeasurement) + chan.CPspacerBytes
%
%           'chanLength is measured in number of samples
%           chanLength = (chan.Cblength \ chunkLength) * chan.CPnumberSamplesTogether
%           TracePrint("  chunkLength: " & chunkLength & "   and chanLength: " & chanLength)
%
%           If (channelsMixedTogether) Then
%             'Typical channel structure, e.g. A1 B1 A2 B2 A3 B3...
%             'All direct channels will come from one BinaryBlock.
%
%             'NOTE: Currently, when the channels are mixed together in the binary
%             ' block, this plugin assumes that the first channel described in the
%             ' FAMOS file (i.e. with a CC key) is the first channel that appears
%             ' in the binary block.  If it is the case that a FAMOS file can
%             ' describe channel Y then channel X, but in the binary block the data
%             ' looks like X1 Y1 X2 Y2... (i.e. channel X starts first), then
%             ' the channels must be sorted before being added to the BinaryBlock
%             TracePrint("  channels are mixed together, using one BinaryBlock")
%           Else
%             'Channels follow each other, e.g. A1 A2 A3... B1 B2 B3...
%             'Each channel gets its own BinaryBlock.
%             Set binBlock = globalFile.GetBinaryBlock()
%             'NOTE: the interaction among all the offsets is not entirely clear
%             binBlock.Position = filePositionDataStart + chan.CPinitialOffset + chan.CbbufferOffset + chan.CbsampleOffset
%             TracePrint("  channels are NOT mixed together, using many BinBlocks")
%           End If
%
%           binBlock.BlockLength = chanLength
%
%           Call FAMOStoUSI(chan.CPdataType,wordType,bytesPerWord)
%
%           If (wordType = eNoType) Then
%             TracePrint("ERROR: invalid type for reading data.  FAMOS type # " & chan.CPdataType)
%             ErrorPrint("ERROR: Unable to process part of the source file.  See log file for more details.")
%           End If
%
%           Set currentChannel = binBlock.Channels.Add(chan.CNname,wordType)
%           Set chan.channel = currentChannel
%
%           If (chan.CRtransformationType = 1) Then
%             TracePrint("  scaled, factor: " & chan.CRfactor & "   offset: " & chan.CRoffset)
%             currentChannel.Factor = chan.CRfactor
%             currentChannel.Offset = chan.CRoffset
%           Else
%             TracePrint("  not scaled")
%           End If
%
%           TracePrint("  adding direct access channel to channel group")
%           currentChannelGroup.Channels.AddDirectAccessChannel(currentChannel)
%
%         Else
%           TracePrint("  reading binary data serially with normal (not direct access) channels")
%
%           'processing channel...
%           If (Not chan.processed) Then
%             TracePrint("  this channel has never been processed; adding it")
%             chan.processed = true
%
%             'Before reading binary data, fill in the channel's meta-data
%             Set currentChannel = currentChannelGroup.Channels.Add(chan.CNname,defaultRealType)'chan.CPdataType
%             Set chan.channel = currentChannel
%
%           Else
%             TracePrint("  this channel has ALREADY been processed")
%             Set currentChannel = chan.channel
%           End If
%           'currentChannel has now been correctly set
%
%           Call ReadBinaryDataSerial(chan,currentChannel,filePositionDataStart)
%
%         End If ' If (direct access channels can be used or not)
%
%         'update meta-data
%         Call currentChannel.Properties.Add("unit_string",chan.CRunit)
%         Call currentChannel.Properties.Add("description",chan.CNcomment)
%         Call currentChannel.Properties.Add("datetime",chan.subGroup.NTdateAndTime.variantdate)
%         If (chan.existsND) Then
%           Call currentChannel.Properties.Add("colorRed",chan.NDcolorRed)
%           Call currentChannel.Properties.Add("colorGreen",chan.NDcolorGreen)
%           Call currentChannel.Properties.Add("colorBlue",chan.NDcolorBlue)
%           Call currentChannel.Properties.Add("yMin",chan.NDyMin)
%           Call currentChannel.Properties.Add("yMax",chan.NDyMax)
%         End If
%
%       End If   ' If (channel has data in this CS key or not)
%       TracePrint(" processing next channel...")
%     Next
%     TracePrint(" processing next chanGroup...")
%   Next
%
%   'adjust the file position to the end of the CS key
%   TracePrint(" done looping; adjusting file position from " & globalFile.Position & " to " & filePositionKeyStart + keyLength)
%   globalFile.Position = filePositionKeyStart + keyLength
%
%   'consume the trailing semicolon unless this is an AVL format file
%   If (Not avlFormat) Then
%     If (Not ByteStream(1) = ";") Then
%       TracePrint("ERROR: CS key improperly formatted, no ; found at end (and file is not an AVL file)")
%       ErrorPrint("ERROR: The file is incorrectly formatted.  See log file for more details.")
%     End If
%   End If
%
%   TracePrint(" done processing CS key")
tempCS=0;
end

% '-------------------------------------------------------------------------------
% ' NU,1,KeyLang,KennLang,Kennwort,Daten;
% ' This is a binary key so keydata is the file position directly after
% '  the comma following the key length
% '-------------------------------------------------------------------------------
function tempNU= ProcessNU(TxtString)
%
tempNU=0;
end

% '-------------------------------------------------------------------------------
% ' CI,1,KeyLang, IndexBlockKey, Zahlenformat, NameLang, Name, Wert, EinheitLang,
% '  Einheit, KommentarLang, Kommentar, Zeit;
% '-------------------------------------------------------------------------------
function tempCI= ProcessCI(TxtString)
%This function is relative to single value
tempCI=0;
end

% '-------------------------------------------------------------------------------
% ' MakeNewChanGroup: Makes a new chanGroup, adds the group to the array of
% '  chanGroups, and sets the parameter 'currentChanGroup' to the new
% '  chanGroup.  currentChanGroup is also a global variable, but it made most of
% '  the program more easily understandable to pass a variable to MakeNewChanGroup.
% ' name: the name of the new chanGroup
% ' currentChanGroup: output - refers to the newly made ChanGroup
% '-------------------------------------------------------------------------------
function MakeNewChanGroup(name,currentChanGroup)
%The funcion is relative to Group
end

function tempChannel=ReadChannel(FileID, ReadStart,ChannelLength, Datatype,factor,offset)
fseek(FileID,ReadStart,'bof');

tempChannel=double(fread(FileID,ChannelLength,Datatype))*factor+offset;
%disp('Info: a Channle was imported.... ');
end


text:
function cvx_optpnt = complex_lorentz( sx, dim )

%COMPLEX_LORENTZ   Complex second-order cone.
%   COMPLEX_LORENTZ(N), where N is a positive integer, creates a column
%   variable of length N and a scalar variable, and constrains them
%   to lie in a second-order cone. That is, given the declaration
%       variable x(n) complex
%       variable y
%   the constraint
%       {x,y} == complex_lorentz(n)
%   is equivalent to
%       norm(x,2) <= y
%   The inequality form is more natural, and preferred in most cases. But
%   in fact, the COMPLEX_LORENTZ set form is used by CVX itself to convert
%   complex NORM()-based constraints to solvable form.
%
%   COMPLEX_LORENTZ(SX,DIM), where SX is a valid size vector and DIM is a
%   positive integer, creates an array variable of size SX and an array
%   variable of size SY (see below) and applies the second-order cone
%   constraint along dimension DIM. That is, given the declarations
%       sy = sx; sy(min(dim,length(sx)+1))=1;
%       variable x(sx) complex
%       variable y(sy)
%   the constraint
%       {x,y} == complex_lorentz(sx,dim)
%   is equivalent to
%       norms(x,2,dim) <= y
%   Again, the inequality form is preferred, but CVX uses the set form
%   internally. DIM is optional; if it is omitted, the first non-singleton
%   dimension is used.
%
%   LORENTZ(SX,DIM,CPLX) creates real second-order cones if CPLX is FALSE,
%   and complex second-order cones if CPLX is TRUE. The latter case is
%   equivalent to COMPLEX_LORENTZ(SX,DIM).
%
%   Disciplined convex programming information:
%       LORENTZ is a cvx set specification. See the user guide for
%       details on how to use sets.

narginchk(1,2);
if nargin == 1,
    cvx_optpnt = lorentz( sx, [], true );
else
    cvx_optpnt = lorentz( sx, dim, true );
end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX