text:
clear;
%% load file

path_to_package = '';   % path to the folder that contains the package
addpath(genpath(path_to_package));

filename = '';      % path to stack tiff file
%foldername = '';    % path to folder that contains a sequence of tiff files

if exist([filename(1:end-3),'mat'],'file')
    data = matfile([filename(1:end-3),'mat'],'Writable',true);
else
    sframe=1;						% user input: first frame to read (optional, default 1)
    num2read=[];					% user input: how many frames to read   (optional, default until the end)
    chunksize=1000;                 % user input: read and map input in chunks (optional, default read all at once)
    data = memmap_file(filename,sframe,num2read,chunksize);
    %data = memmap_file_sequence(foldername);
end

%% Set parameters
sizY = size(data,'Y');                  % size of data matrix
patch_size = [32,32];                   % size of each patch along each dimension (optional, default: [32,32])
overlap = [4,4];                        % amount of overlap in each dimension (optional, default: [4,4])

patches = construct_patches(sizY(1:end-1),patch_size,overlap);
K = 8;                                            % number of components to be found
tau = 8;                                          % std of gaussian kernel (size of neuron)
p = 0;                                            % order of autoregressive system (p = 0 no dynamics, p=1 just decay, p = 2, both rise and decay)
merge_thr = 0.8;                                  % merging threshold
sizY = data.sizY;

options = CNMFSetParms(...
    'd1',sizY(1),'d2',sizY(2),...
    'search_method','ellipse','dist',3,...      % search locations when updating spatial components
    'deconv_method','constrained_foopsi',...    % activity deconvolution method
    'temporal_iter',2,...                       % number of block-coordinate descent steps
    'ssub',2,...
    'tsub',4,...
    'fudge_factor',0.98,...                     % bias correction for AR coefficients
    'merge_thr',merge_thr,...                    % merging threshold
    'gSig',tau...
    );

%% Run on patches

[A,b,C,f,S,P,RESULTS,YrA] = run_CNMF_patches(data,K,patches,tau,p,options);

%% order and plot

[A_or,C_or,S_or,P] = order_ROIs(A,C,S,P); % order components

contour_threshold = 0.95;                 % amount of energy used for each component to construct contour plot
figure;
[Coor,json_file] = plot_contours(A_or,reshape(P.sn,sizY(1),sizY(2)),contour_threshold,1); % contour plot of spatial footprints
%savejson('jmesh',json_file,'filename');        % optional save json file with component coordinates (requires matlab json library)

author: zhoupc


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Eigen_Vector,Eigen_Value]=Find_K_Max_Gen_Eigen(Matrix1,Matrix2,Eigen_NUM)

[NN,NN]=size(Matrix1);
[V,S]=eig(Matrix1,Matrix2); %Note this is equivalent to; [V,S]=eig(St,SL); also equivalent to [V,S]=eig(Sn,St); %

S=diag(S);
[S,index]=sort(S);

Eigen_Vector=zeros(NN,Eigen_NUM);
Eigen_Value=zeros(1,Eigen_NUM);

p=NN;
for t=1:Eigen_NUM
    Eigen_Vector(:,t)=V(:,index(p));
    Eigen_Value(t)=S(p);
    p=p-1;
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [ gx, gy ] = g00_xy ( gi, gn )

%*****************************************************************************80
%
%% G00_XY returns the grid points for any grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer GI, the index of the grid.
%
%    Input, integer GN, the grid size.
%
%    Output, real GX(GN,1), GY(GN,1), the grid coordinates.
%
  if ( gi == 1 )
    [ gx, gy ] = g01_xy ( gn );
  elseif ( gi == 2 )
    [ gx, gy ] = g02_xy ( gn );
  elseif ( gi == 3 )
    [ gx, gy ] = g03_xy ( gn );
  elseif ( gi == 4 )
    [ gx, gy ] = g04_xy ( gn );
  elseif ( gi == 5 )
    [ gx, gy ] = g05_xy ( gn );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'G00_XY - Fatal error!\n' );
    fprintf ( 1, '  Illegal grid index GI = %d\n', gi );
    error ( 'G00_XY - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef Max < dagnn.ElementWise
  %SUM DagNN max layer
  %   The MAX layer takes the max of all its inputs and store the result
  %   as its only output.

  properties (Transient)
    numInputs
  end

  methods
    function outputs = forward(obj, inputs, params)
      obj.numInputs = numel(inputs) ;
      outputs{1} = vl_nnmax(obj.numInputs, inputs{:}) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs = vl_nnmax(obj.numInputs, inputs{:}, derOutputs{1}) ;
      keyboard
      derParams = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes{1} = inputSizes{1} ;
      for k = 2:numel(inputSizes)
        if all(~isnan(inputSizes{k})) && all(~isnan(outputSizes{1}))
          if ~isequal(inputSizes{k}, outputSizes{1})
            warning('Max layer: the dimensions of the input variables is not the same.') ;
          end
        end
      end
    end

    function rfs = getReceptiveFields(obj)
      numInputs = numel(obj.net.layers(obj.layerIndex).inputs) ;
      rfs.size = [1 1] ;
      rfs.stride = [1 1] ;
      rfs.offset = [1 1] ;
      rfs = repmat(rfs, numInputs, 1) ;
    end

    function obj = Max(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MFT


text:
function phow_caltech101()
% PHOW_CALTECH101 Image classification in the Caltech-101 dataset
%   This program demonstrates how to use VLFeat to construct an image
%   classifier on the Caltech-101 data. The classifier uses PHOW
%   features (dense SIFT), spatial histograms of visual words, and a
%   Chi2 SVM. To speedup computation it uses VLFeat fast dense SIFT,
%   kd-trees, and homogeneous kernel map. The program also
%   demonstrates VLFeat PEGASOS SVM solver, although for this small
%   dataset other solvers such as LIBLINEAR can be more efficient.
%
%   By default 15 training images are used, which should result in
%   about 64% performance (a good performance considering that only a
%   single feature type is being used).
%
%   Call PHOW_CALTECH101 to train and test a classifier on a small
%   subset of the Caltech-101 data. Note that the program
%   automatically downloads a copy of the Caltech-101 data from the
%   Internet if it cannot find a local copy.
%
%   Edit the PHOW_CALTECH101 file to change the program configuration.
%
%   To run on the entire dataset change CONF.TINYPROBLEM to FALSE.
%
%   The Caltech-101 data is saved into CONF.CALDIR, which defaults to
%   'data/caltech-101'. Change this path to the desired location, for
%   instance to point to an existing copy of the Caltech-101 data.
%
%   The program can also be used to train a model on custom data by
%   pointing CONF.CALDIR to it. Just create a subdirectory for each
%   class and put the training images there. Make sure to adjust
%   CONF.NUMTRAIN accordingly.
%
%   Intermediate files are stored in the directory CONF.DATADIR. All
%   such files begin with the prefix CONF.PREFIX, which can be changed
%   to test different parameter settings without overriding previous
%   results.
%
%   The program saves the trained model in
%   <CONF.DATADIR>/<CONF.PREFIX>-model.mat. This model can be used to
%   test novel images independently of the Caltech data.
%
%     load('data/baseline-model.mat') ; # change to the model path
%     label = model.classify(model, im) ;
%

% Author: Andrea Vedaldi

% Copyright (C) 2011-2013 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).

conf.calDir = 'data/caltech-101' ;
conf.dataDir = 'data/' ;
conf.autoDownloadData = true ;
conf.numTrain = 15 ;
conf.numTest = 15 ;
conf.numClasses = 102 ;
conf.numWords = 600 ;
conf.numSpatialX = [2 4] ;
conf.numSpatialY = [2 4] ;
conf.quantizer = 'kdtree' ;
conf.svm.C = 10 ;

conf.svm.solver = 'sdca' ;
%conf.svm.solver = 'sgd' ;
%conf.svm.solver = 'liblinear' ;

conf.svm.biasMultiplier = 1 ;
conf.phowOpts = {'Step', 3} ;
conf.clobber = false ;
conf.tinyProblem = true ;
conf.prefix = 'baseline' ;
conf.randSeed = 1 ;

if conf.tinyProblem
  conf.prefix = 'tiny' ;
  conf.numClasses = 5 ;
  conf.numSpatialX = 2 ;
  conf.numSpatialY = 2 ;
  conf.numWords = 300 ;
  conf.phowOpts = {'Verbose', 2, 'Sizes', 7, 'Step', 5} ;
end

conf.vocabPath = fullfile(conf.dataDir, [conf.prefix '-vocab.mat']) ;
conf.histPath = fullfile(conf.dataDir, [conf.prefix '-hists.mat']) ;
conf.modelPath = fullfile(conf.dataDir, [conf.prefix '-model.mat']) ;
conf.resultPath = fullfile(conf.dataDir, [conf.prefix '-result']) ;

randn('state',conf.randSeed) ;
rand('state',conf.randSeed) ;
vl_twister('state',conf.randSeed) ;

% --------------------------------------------------------------------
%                                            Download Caltech-101 data
% --------------------------------------------------------------------

if ~exist(conf.calDir, 'dir') || ...
   (~exist(fullfile(conf.calDir, 'airplanes'),'dir') && ...
    ~exist(fullfile(conf.calDir, '101_ObjectCategories', 'airplanes')))
  if ~conf.autoDownloadData
    error(...
      ['Caltech-101 data not found. ' ...
       'Set conf.autoDownloadData=true to download the required data.']) ;
  end
  vl_xmkdir(conf.calDir) ;
  calUrl = ['http://www.vision.caltech.edu/Image_Datasets/' ...
    'Caltech101/101_ObjectCategories.tar.gz'] ;
  fprintf('Downloading Caltech-101 data to ''%s''. This will take a while.', conf.calDir) ;
  untar(calUrl, conf.calDir) ;
end

if ~exist(fullfile(conf.calDir, 'airplanes'),'dir')
  conf.calDir = fullfile(conf.calDir, '101_ObjectCategories') ;
end

% --------------------------------------------------------------------
%                                                           Setup data
% --------------------------------------------------------------------
classes = dir(conf.calDir) ;
classes = classes([classes.isdir]) ;
classes = {classes(3:conf.numClasses+2).name} ;

images = {} ;
imageClass = {} ;
for ci = 1:length(classes)
  ims = dir(fullfile(conf.calDir, classes{ci}, '*.jpg'))' ;
  ims = vl_colsubset(ims, conf.numTrain + conf.numTest) ;
  ims = cellfun(@(x)fullfile(classes{ci},x),{ims.name},'UniformOutput',false) ;
  images = {images{:}, ims{:}} ;
  imageClass{end+1} = ci * ones(1,length(ims)) ;
end
selTrain = find(mod(0:length(images)-1, conf.numTrain+conf.numTest) < conf.numTrain) ;
selTest = setdiff(1:length(images), selTrain) ;
imageClass = cat(2, imageClass{:}) ;

model.classes = classes ;
model.phowOpts = conf.phowOpts ;
model.numSpatialX = conf.numSpatialX ;
model.numSpatialY = conf.numSpatialY ;
model.quantizer = conf.quantizer ;
model.vocab = [] ;
model.w = [] ;
model.b = [] ;
model.classify = @classify ;

% --------------------------------------------------------------------
%                                                     Train vocabulary
% --------------------------------------------------------------------

if ~exist(conf.vocabPath) || conf.clobber

  % Get some PHOW descriptors to train the dictionary
  selTrainFeats = vl_colsubset(selTrain, 30) ;
  descrs = {} ;
  %for ii = 1:length(selTrainFeats)
  parfor ii = 1:length(selTrainFeats)
    im = imread(fullfile(conf.calDir, images{selTrainFeats(ii)})) ;
    im = standarizeImage(im) ;
    [drop, descrs{ii}] = vl_phow(im, model.phowOpts{:}) ;
  end

  descrs = vl_colsubset(cat(2, descrs{:}), 10e4) ;
  descrs = single(descrs) ;

  % Quantize the descriptors to get the visual words
  vocab = vl_kmeans(descrs, conf.numWords, 'verbose', 'algorithm', 'elkan', 'MaxNumIterations', 50) ;
  save(conf.vocabPath, 'vocab') ;
else
  load(conf.vocabPath) ;
end

model.vocab = vocab ;

if strcmp(model.quantizer, 'kdtree')
  model.kdtree = vl_kdtreebuild(vocab) ;
end

% --------------------------------------------------------------------
%                                           Compute spatial histograms
% --------------------------------------------------------------------

if ~exist(conf.histPath) || conf.clobber
  hists = {} ;
  parfor ii = 1:length(images)
  % for ii = 1:length(images)
    fprintf('Processing %s (%.2f %%)\n', images{ii}, 100 * ii / length(images)) ;
    im = imread(fullfile(conf.calDir, images{ii})) ;
    hists{ii} = getImageDescriptor(model, im);
  end

  hists = cat(2, hists{:}) ;
  save(conf.histPath, 'hists') ;
else
  load(conf.histPath) ;
end

% --------------------------------------------------------------------
%                                                  Compute feature map
% --------------------------------------------------------------------

psix = vl_homkermap(hists, 1, 'kchi2', 'gamma', .5) ;

% --------------------------------------------------------------------
%                                                            Train SVM
% --------------------------------------------------------------------

if ~exist(conf.modelPath) || conf.clobber
  switch conf.svm.solver
    case {'sgd', 'sdca'}
      lambda = 1 / (conf.svm.C *  length(selTrain)) ;
      w = [] ;
      parfor ci = 1:length(classes)
        perm = randperm(length(selTrain)) ;
        fprintf('Training model for class %s\n', classes{ci}) ;
        y = 2 * (imageClass(selTrain) == ci) - 1 ;
        [w(:,ci) b(ci) info] = vl_svmtrain(psix(:, selTrain(perm)), y(perm), lambda, ...
          'Solver', conf.svm.solver, ...
          'MaxNumIterations', 50/lambda, ...
          'BiasMultiplier', conf.svm.biasMultiplier, ...
          'Epsilon', 1e-3);
      end

    case 'liblinear'
      svm = train(imageClass(selTrain)', ...
                  sparse(double(psix(:,selTrain))),  ...
                  sprintf(' -s 3 -B %f -c %f', ...
                          conf.svm.biasMultiplier, conf.svm.C), ...
                  'col') ;
      w = svm.w(:,1:end-1)' ;
      b =  svm.w(:,end)' ;
  end

  model.b = conf.svm.biasMultiplier * b ;
  model.w = w ;

  save(conf.modelPath, 'model') ;
else
  load(conf.modelPath) ;
end

% --------------------------------------------------------------------
%                                                Test SVM and evaluate
% --------------------------------------------------------------------

% Estimate the class of the test images
scores = model.w' * psix + model.b' * ones(1,size(psix,2)) ;
[drop, imageEstClass] = max(scores, [], 1) ;

% Compute the confusion matrix
idx = sub2ind([length(classes), length(classes)], ...
              imageClass(selTest), imageEstClass(selTest)) ;
confus = zeros(length(classes)) ;
confus = vl_binsum(confus, ones(size(idx)), idx) ;

% Plots
figure(1) ; clf;
subplot(1,2,1) ;
imagesc(scores(:,[selTrain selTest])) ; title('Scores') ;
set(gca, 'ytick', 1:length(classes), 'yticklabel', classes) ;
subplot(1,2,2) ;
imagesc(confus) ;
title(sprintf('Confusion matrix (%.2f %% accuracy)', ...
              100 * mean(diag(confus)/conf.numTest) )) ;
print('-depsc2', [conf.resultPath '.ps']) ;
save([conf.resultPath '.mat'], 'confus', 'conf') ;

% -------------------------------------------------------------------------
function im = standarizeImage(im)
% -------------------------------------------------------------------------

im = im2single(im) ;
if size(im,1) > 480, im = imresize(im, [480 NaN]) ; end

% -------------------------------------------------------------------------
function hist = getImageDescriptor(model, im)
% -------------------------------------------------------------------------

im = standarizeImage(im) ;
width = size(im,2) ;
height = size(im,1) ;
numWords = size(model.vocab, 2) ;

% get PHOW features
[frames, descrs] = vl_phow(im, model.phowOpts{:}) ;

% quantize local descriptors into visual words
switch model.quantizer
  case 'vq'
    [drop, binsa] = min(vl_alldist(model.vocab, single(descrs)), [], 1) ;
  case 'kdtree'
    binsa = double(vl_kdtreequery(model.kdtree, model.vocab, ...
                                  single(descrs), ...
                                  'MaxComparisons', 50)) ;
end

for i = 1:length(model.numSpatialX)
  binsx = vl_binsearch(linspace(1,width,model.numSpatialX(i)+1), frames(1,:)) ;
  binsy = vl_binsearch(linspace(1,height,model.numSpatialY(i)+1), frames(2,:)) ;

  % combined quantization
  bins = sub2ind([model.numSpatialY(i), model.numSpatialX(i), numWords], ...
                 binsy,binsx,binsa) ;
  hist = zeros(model.numSpatialY(i) * model.numSpatialX(i) * numWords, 1) ;
  hist = vl_binsum(hist, ones(size(bins)), bins) ;
  hists{i} = single(hist / sum(hist)) ;
end
hist = cat(1,hists{:}) ;
hist = hist / sum(hist) ;

% -------------------------------------------------------------------------
function [className, score] = classify(model, im)
% -------------------------------------------------------------------------

hist = getImageDescriptor(model, im) ;
psix = vl_homkermap(hist, 1, 'kchi2', 'period', .7) ;
scores = model.w' * psix + model.b' ;
[score, best] = max(scores) ;
className = model.classes{best} ;

author: jianxiongxiao

repo: ProfXkit


text:
clear all; close all; clc; rng('default');

% Ambient space dimensions
D = 50;
% subspace dimension
K = 8;
% Angle between subspaces A-B and B-C.
thetas = [4:2:40]; % in degree
nt  = numel(thetas);
% Number of points per subspace
ng_values = [4:2:60];
nng = numel(ng_values);
SNRs = [15, 25 35];
nsnrs = numel(SNRs);
Trials = 400;

average_fmeasures = zeros(nt, nsnrs, nng);
average_precisions = zeros(nt, nsnrs, nng);
average_recalls = zeros(nt, nsnrs, nng);
average_clustering_ratios = zeros(nt, nsnrs, nng);

for t=1:nt
    theta = thetas(t);
    for sn=1:nsnrs
        SNR = SNRs(sn);
        for i=1:nng
            Ng = ng_values(i);
            fmeasures = zeros(1, Trials);
            precisions = zeros(1, Trials);
            recalls = zeros(1, Trials);
            clustering_ratios = zeros(1, Trials);
            for s=1:Trials
                if mod(s, 1) == 0
                    fprintf('Theta: %d, SNR: %d dB,  Signals per subspace: %d, Trial: %d\n', theta, SNR, Ng, s);
                end
                result = SimulateSSCOMP_3Spaces(D, K, Ng, theta, SNR);
                cmpr = result.comparison;
                fmeasures(s) = cmpr.fMeasure;
                precisions(s) = cmpr.precision;
                recalls(s) = cmpr.recall;
                clustering_ratios(s) = cmpr.clusteringRatio;
            end
            average_fmeasures(t, sn, i) = mean(fmeasures);
            average_precisions(t, sn, i) = mean(precisions);
            average_recalls(t, sn, i) = mean(recalls);
            average_clustering_ratios(t, sn, i) = mean(clustering_ratios);
        end
    end
end
save(sprintf('bin/phase_transition_theta_signals_K=%d.mat',  K));

author: indigits

repo: sparse-plex


text:
function [Alp,blp] = coordinate_constraint(V,C)
  % COORDINATE_CONSTRAINT  Build a "coordinate constraint" matrix for enforcing
  % linear precision in generalized barycentric coordinates or skinning
  % weights.
  %
  % [Alp,blp] = coordinate_constraint(V,C)
  %
  % Inputs:
  %   V  #V by dim list of sample locations
  %   C  #C by dim list of cage/control vertex positions
  % Outputs:
  %   Alp  #V*dim by #V*c sparse matrix of constraint coefficients
  %   blp  #V*dim list of right hand sides
  %
  % Example:
  %   % fine mesh in (V,F) cage in (C,CF), boundary conditions in (b,bc)
  %   [Alp,blp] = coordinate_constraint(V,C)
  %   W = min_quad_with_fixed(Q,[],b,bc,Alp,blp)
  %   max(abs(Alp*W(:) - blp))
  %

  % dimensions
  n = size(V,1);
  d = size(C,2);
  assert(d == size(V,2));
  c = size(C,1);
  AlpI = reshape(repmat(1:n,c,1),n*c,1);
  AlpJ = reshape(bsxfun(@plus,repmat((1:c)-1,n,1)*n,(1:n)')',n*c,1);
  %% Handle each dimension
  AlpI = reshape(bsxfun(@plus,AlpI,((1:d)-1)*n),n*c*d,1);
  AlpJ = repmat(AlpJ(:),d,1);
  AlpV = reshape(repmat(C,n,1),n*c*d,1);
  Alp = sparse(AlpI(:),AlpJ(:),AlpV(:),n*d,c*n);
  blp = V(:);

end

author: alecjacobson

repo: gptoolbox


text:
function c8vec_uniform_01_test ( )

%*****************************************************************************80
%
%% C8VEC_UNIFORM_01_TEST tests C8VEC_UNIFORM_01;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2015
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8VEC_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  C8VEC_UNIFORM_01_TEST returns a random C8VEC\n' );

  seed = 123456789;

  [ a, seed ] = c8vec_uniform_01 ( n, seed );

  c8vec_print ( n, a, '  The random C8VEC:' );

  return
end


text:
%%                   runSignalReconstructionDemo.m
%
% This script will create phaseless measurements from a 1d test signal, and
% then recover the image using phase retrieval methods.  We now describe
% the details of the simple recovery problem that this script implements.
%
%                         Recovery Problem
% This script creates a complex-valued random Gaussian signal. Measurements
% of the signal are then obtained by applying a linear operator to the
% signal, and computing the magnitude (i.e., removing the phase) of
% the results.
%
%                       Measurement Operator
% Measurement are obtained using a linear operator, called 'A', that
% contains random Gaussian entries.
%
%                      The Recovery Algorithm
% The image is recovered by calling the method 'solvePhaseRetrieval', and
% handing the measurement operator and linear measurements in as arguments.
% A struct containing options is also handed to 'solvePhaseRetrieval'.
% The entries in this struct specify which recovery algorithm is used.
%
% For more details, see the Phasepack user guide.
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

function runSignalReconstructionDemo()

%% Specify the signal length, and number of measurements
n = 100;  % The signal length
m = 8*n;  % The number of measurements

%% Build the target signal
x_true = randn(n,1)+1i*randn(n,1);

%% Create the measurement operator
% Note: we use a dense matrix in this example, but PhasePack also supports
% function handles.  See the more complex 'runImageReconstructionDemo.m'
% script for an example using the fast Fourier transform.
A = randn(m,n)+1i*randn(m,n);

%% Compute phaseless measurements
b = abs(A*x_true);

%% Set options for PhasePack - this is where we choose the recovery algorithm
opts = struct;                  % Create an empty struct to store options
opts.algorithm = 'Fienup';      % Use the Fienup method to solve the retrieval problem.  Try changing this to 'twf' for truncated Wirtinger flow.
opts.initMethod = 'optimal';    % Use the optimal spectral initializer method to generate an initial starting point for the solver
opts.tol = 1e-3;                % The tolerance - make this smaller for more accurate solutions, or larger for faster runtimes
opts.verbose = 2;               % Print out lots of information as the solver runs (set this to 1 or 0 for less output)

%% Run the Phase retrieval Algorithm
fprintf('Running %s algorithm\n',opts.algorithm);
% Call the solver using the measurement operator 'A', the
% measurements 'b', the length of the signal to be recovered, and the
% options.  Note, the measurement operator can be either a function handle
% or a matrix.   Here, we use a matrix.  In this case, we have omitted the
% second argument. If 'A' had been a function handle, we would have
% handed the transpose of 'A' in as the second argument.
[x, outs] = solvePhaseRetrieval(A, [], b, n, opts);
% Note: 'outs' is a struct containing convergene information.

%% Remove phase ambiguity
% Phase retrieval can only recover images up to a phase ambiguity.
% Let's apply a phase rotation to align the recovered signal with the
% original so they look the same when we plot them.
rotation = sign(x'*x_true(:));
x = x*rotation;

% Print some useful info to the console
fprintf('Signal recovery required %d iterations (%f secs)\n',outs.iterationCount, outs.solveTimes(end));

%% Plot results
figure;
% Plot the true vs recovered signal.  Ideally, this scatter plot should be
% clustered around the 45-degree line.
subplot(1,2,1);
scatter(real(x_true),real(x));
xlabel('Original signal value');
ylabel('Recovered signal value');
title('Original vs recovered signal');

% Plot a convergence curve
subplot(1,3,3);
convergedCurve = semilogy(outs.solveTimes, outs.residuals);
set(convergedCurve, 'linewidth',1.75);
grid on;
xlabel('Time (sec)');
ylabel('Error');
title('Convergence Curve');
set(gcf,'units','points','position',[0,0,1200,300]);

end

author: tomgoldstein

repo: phasepack-matlab


text:
function pa_wavplay(varargin)
%    pa_wavplay: playback a matrix to multichannel sound hardware
%
%    Usage: pa_wavplay([devicetype])
%           pa_wavplay(buffer, [samplerate], [deviceid], [devicetype])
%
%    pa_wavplay is a tool for playing multi-channel audio through your ASIO
%    device. Arguments in [] are optional
%
%    pa_wavplay([devicetype]) will list all your audio devices for that
%    device type then exit.
%
%    - buffer       is the matrix to play
%    - samplerate	is the sampling frequency. Default: 44100
%    - deviceid 	is the device id to use for output. Default: 0
%    - devicetype   determines which sound driver to use
%        'win'      Windows Multimedia Device
%        'dx'       DirectX DirectSound driver
%        'asio'     ASIO Driver (default)
%
%    For stereo playback, buffer should be an N-by-2 matrix. The number of
%    audio channels supported is hardware dependent.
%
%    samplerate should be a valid integer sample rate, eg. 22050 44100 etc.
%
%    As of pa_wavplay 0.2, buffer should be either of type double or
%    single. pa_wavplay uses 32-bit floats (single) internally.
%
%    SEE ALSO: pa_wavrecord, pa_wavplayrecord

% check right num of args
error(nargchk(0,4,nargin));

% defaults
device_opt = 1;
device = 0;
fs = 44100;

% no args, print devices for asio
if (nargin==0)
    pawavplaya;
    return;
end

% if devtype specified
if ischar(varargin{end}),
   s=varargin{end};
   varargin(end)=[];

   device_opt = strmatch(s,{'asio','win', 'dx'});
   if isempty(device_opt),
       error(['Unrecognized DEVICETYPE: ' s]);
   end

   if (nargin==1)
       if device_opt==1, % asio
           pawavplaya;
       end
       if device_opt==2, % win
           pawavplayw;
       end
       if device_opt==3, % dx
           pawavplayx;
       end
       return;
    end
end

% data buffer
y = varargin{1};

% sample rate
if length(varargin)>=2,
   fs=varargin{2};
end

% device id
if length(varargin) >= 3,
    if (~(ischar(varargin{3})))
        device = varargin{3};
    end
end

if device_opt==1, % asio
    fprintf('Using ASIO driver\n');
    pawavplaya(y, device, fs, 0, 0, 0, -1);

elseif device_opt==2, % win
    fprintf('Using WMME driver\n');
    pawavplayw(y, device, fs, 0, 0, 0, -1);

elseif device_opt==3, % dx
    fprintf('Using DirectX driver\n');
    pawavplayx(y, device, fs, 0, 0, 0, -1);
end

% [EOF] pa_wavplay.m

author: Sable

repo: mcbench-benchmarks


text:
function [DCM] = spm_dcm_erp_results(DCM,Action,fig)
% Results for ERP Dynamic Causal Modeling (DCM)
% FORMAT spm_dcm_erp_results(DCM,'ERPs (mode)');
% FORMAT spm_dcm_erp_results(DCM,'ERPs (sources)');
% FORMAT spm_dcm_erp_results(DCM,'Coupling (A)');
% FORMAT spm_dcm_erp_results(DCM,'Coupling (B)');
% FORMAT spm_dcm_erp_results(DCM,'Coupling (C)');
% FORMAT spm_dcm_erp_results(DCM,'trial-specific effects');
% FORMAT spm_dcm_erp_results(DCM,'Input');
% FORMAT spm_dcm_erp_results(DCM,'Response');
% FORMAT spm_dcm_erp_results(DCM,'Response (image)');
% FORMAT spm_dcm_erp_results(DCM,'Scalp maps');
% FORMAT spm_dcm_erp_results(DCM,'Data');
%
%__________________________________________________________________________
%
% DCM is a causal modelling procedure for dynamical systems in which
% causality is inherent in the differential equations that specify the model.
% The basic idea is to treat the system of interest, in this case the brain,
% as an input-state-output system.  By perturbing the system with known
% inputs, measured responses are used to estimate various parameters that
% govern the evolution of brain states.  Although there are no restrictions
% on the parameterisation of the model, a bilinear approximation affords a
% simple re-parameterisation in terms of effective connectivity.  This
% effective connectivity can be latent or intrinsic or, through bilinear
% terms, model input-dependent changes in effective connectivity.  Parameter
% estimation proceeds using fairly standard approaches to system
% identification that rest upon Bayesian inference.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_erp_results.m 6644 2015-12-12 14:53:37Z karl $

% get Action if necessary
%--------------------------------------------------------------------------
if nargin < 2

    str{1}  = 'ERPs (mode)';
    str{2}  = 'ERPs (sources)';
    str{3}  = 'Coupling (A)';
    str{4}  = 'Coupling (B)';
    str{5}  = 'Coupling (C)';
    str{6}  = 'trial-specific effects';
    str{7}  = 'Input';
    str{8}  = 'Response';
    str{9}  = 'Response (image)';
    str{10} = 'Scalp maps';
    str{11} = 'Data';

    s       = listdlg('PromptString','Select an option:',...
        'SelectionMode','single',...
        'ListString',str);

    Action = str{s};

end

% get figure
%--------------------------------------------------------------------------
if nargin < 3, spm_figure('GetWin','Graphics'); end
colormap(gray)
clf

% trial data
%--------------------------------------------------------------------------
xY  = DCM.xY;                   % data
nt  = length(xY.y);             % Nr trial types
ne  = size(xY.y{1},2);          % Nr electrodes
nb  = size(xY.y{1},1);          % Nr time bins
t   = xY.pst;                   % PST

% plot data
%--------------------------------------------------------------------------
switch(lower(Action))
    case{lower('Data')}
        try
            A = [0 0];
            for i = 1:nt

                % confounds if specified
                %----------------------------------------------------------
                try
                    X0 = spm_orth(xY.X0(1:nb,:),'norm');
                    R  = speye(nb,nb) - X0*X0';
                catch
                    R  = speye(nb,nb);
                end

                % plot data
                %----------------------------------------------------------
                subplot(nt,2,(i - 1)*2 + 1)
                plot(t,R*xY.y{i})
                xlabel('time (ms)')
                try
                    title(sprintf('Observed response (code:%i)',xY.code(i)))
                catch
                    title(sprintf('Observed response %i',i))
                end
                axis square
                a    = axis;
                A(1) = min(A(1),a(3));
                A(2) = max(A(2),a(4));

                % image data
                %----------------------------------------------------------
                subplot(nt,2,(i - 1)*2 + 2)
                imagesc([1:ne],t,R*xY.y{i})
                xlabel('channels');ylabel('peri-stimulus time (ms)')
                axis square
                try
                    title(sprintf('Observed response (code:%i)',xY.code(i)))
                catch
                    title(sprintf('Observed response %i',i))
                end
            end

            % set axis
            %--------------------------------------------------------------
            for i = 1:nt
                subplot(nt,2,(i - 1)*2 + 1)
                set(gca,'YLim',A)
            end
        end
        return
end

% post inversion parameters
%--------------------------------------------------------------------------
nu  = length(DCM.B);          % Nr inputs
nc  = size(DCM.H{1},2);       % Nr modes
ns  = size(DCM.A{1},1);       % Nr of sources
np  = size(DCM.K{1},2)/ns;    % Nr of population per source

% switch
%--------------------------------------------------------------------------
switch(lower(Action))

    case{lower('ERPs (mode)')}

        % spm_dcm_erp_results(DCM,'ERPs (mode)');
        %------------------------------------------------------------------
        co = {'b', 'r', 'g', 'm', 'y', 'k'};
        lo = {'-', '--'}; A = [0 0];

        for i = 1:nc
            subplot(ceil(nc/2),2,i), hold on
            str   = {};
            for j = 1:nt
                plot(t,DCM.H{j}(:,i), lo{1},...
                    'Color', co{j},...
                    'LineWidth',2);
                str{end + 1} = sprintf('trial %i (predicted)',j);
                plot(t,DCM.H{j}(:,i) + DCM.R{j}(:,i), lo{2},...
                    'Color',co{j});
                str{end + 1} = sprintf('trial %i (observed)',j);
                set(gca, 'XLim', [t(1) t(end)]);

            end
            hold off
            title(sprintf('mode %i',i))
            grid on
            axis square
            a    = axis;
            A(1) = min(A(1),a(3));
            A(2) = max(A(2),a(4));
        end
        xlabel('time (ms)')
        legend(str)

        % set axis
        %------------------------------------------------------------------
        for i = 1:nc
            subplot(ceil(nc/2),2,i)
            set(gca,'YLim',A)
        end

    case{lower('ERPs (sources)')}

        % spm_dcm_erp_results(DCM,'ERPs (sources)');
        %------------------------------------------------------------------
        col   = {'b','r','g','m','y','c'}; A = [0 0];
        sty   = {':','-.','-','--','-'};
        for i = 1:ns
            str   = {};
            subplot(ceil(ns/2),2,i), hold on

            % if J maps from states to sources, use J (a matrix)
            %--------------------------------------------------------------
            if strcmpi(DCM.options.model,'DEM')

                for j = find(DCM.Eg.J(i,:))
                    for k = 1:nt
                        plot(t, DCM.K{k}(:,j) , ...
                            'Color',col{k}, ...
                            'LineWidth',1);
                        str{end + 1} = sprintf('trial %i (pop. %i)',k,j);
                    end
                end

            else

                % otherwise assume normal form for states (source x states)
                %----------------------------------------------------------
                for j = 1:np
                    for k = 1:nt

                            plot(t, DCM.K{k}(:,i + ns*(j - 1)), ...
                                'Color',col{k}, ...
                                'LineStyle',sty{j}, ...
                                'LineWidth',2);

                        str{end + 1} = sprintf('trial %i (pop. %i)',k,j);
                    end
                end
            end

            hold off
            title(DCM.Sname{i},'FontSize',16)
            grid on
            axis square
            a    = axis;
            A(1) = min(A(1),a(3));
            A(2) = max(A(2),a(4));
        end
        xlabel('time (ms)','FontSize',14)
        legend(str)

        % set axis
        %------------------------------------------------------------------
        for i = 1:ns
            subplot(ceil(ns/2),2,i)
            axis([t(1) t(end) A(1) A(2)]);

            % or
            %--------------------------------------------------------------
            spm_axis tight

        end

    case{lower('Coupling (A)')}

        % spm_dcm_erp_results(DCM,'coupling (A)');
        %------------------------------------------------------------------
        if ~isfield(DCM.Ep,'A'), return, end
        n    = length(DCM.Ep.A);
        if n == 2
            str = {'Forward','Backward'};
        elseif n == 3
            str = {'Forward','Backward','Lateral'};
        else
            str = {'Forward (i)','Forward (ii)','Backward (i)','Backward (ii)'};
        end

        for i = 1:n

            % images
            %--------------------------------------------------------------
            subplot(4,n,i)
            imagesc(exp(DCM.Ep.A{i}))
            title(str{i},'FontSize',10)
            set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname,'FontSize',8)
            set(gca,'XTick',[])
            xlabel('from','FontSize',8)
            ylabel('to','FontSize',8)
            axis square

            % table
            %--------------------------------------------------------------
            subplot(4,n,i + n)
            text(0,1/2,num2str(full(exp(DCM.Ep.A{i})),' %.2f'),'FontSize',8)
            axis off,axis square

            % PPM
            %--------------------------------------------------------------
            subplot(4,n,i + n + n)
            image(64*DCM.Pp.A{i})
            set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname,'FontSize',8)
            set(gca,'XTick',[])
            title('PPM')
            axis square

            % table
            %--------------------------------------------------------------
            subplot(4,n,i + n + n + n)
            text(0,1/2,num2str(DCM.Pp.A{i},' %.2f'),'FontSize',8)
            axis off, axis square

        end

    case{lower('Coupling (C)')}

        % spm_dcm_erp_results(DCM,'coupling (C)');
        %------------------------------------------------------------------
        if ~isfield(DCM.Ep,'C'), return, end

        % images
        %------------------------------------------------------------------
        subplot(2,4,1)
        imagesc(exp(DCM.Ep.C))
        title('Factors','FontSize',10)
        set(gca,'XTick',[1:nu],'XTickLabel','Input','FontSize',8)
        set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname, 'FontSize',8)
        axis square

        % PPM
        %------------------------------------------------------------------
        subplot(2,4,3)
        image(64*DCM.Pp.C)
        title('Factors','FontSize',10)
        set(gca,'XTick',[1:nu],'XTickLabel','Input','FontSize',8)
        set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname, 'FontSize',8)
        axis square
        title('PPM')

        % table
        %------------------------------------------------------------------
        subplot(2,4,2)
        text(0,1/2,num2str(full(exp(DCM.Ep.C)),' %.2f'),'FontSize',8)
        axis off

        % table
        %------------------------------------------------------------------
        subplot(2,4,4)
        text(0,1/2,num2str(DCM.Pp.C,' %.2f'),'FontSize',8)
        axis off

    case{lower('Coupling (B)')}

        % spm_dcm_erp_results(DCM,'coupling (B)');
        %------------------------------------------------------------------
        if ~isfield(DCM.Ep,'B'), return, end
        if  isfield(DCM.Ep,'N')
            if any(spm_vec(DCM.Ep.N))
                xU.name = [DCM.xU.name DCM.xU.name];
                Ep.B    = [DCM.Ep.B DCM.Ep.N];
                Pp.B    = [DCM.Pp.B DCM.Pp.N];
            else
                xU.name = DCM.xU.name;
                Ep.B    = DCM.Ep.B;
                Pp.B    = DCM.Pp.B;
            end
        else
            xU.name = DCM.xU.name;
            Ep.B    = DCM.Ep.B;
            Pp.B    = DCM.Pp.B;
        end

        nb    = length(Ep.B);
        for i = 1:nb

            % images
            %--------------------------------------------------------------
            subplot(4,nb,i)
            imagesc(exp(Ep.B{i}))
            title(xU.name{i},'FontSize',10)
            set(gca,'YTick',1:ns,'YTickLabel',DCM.Sname,'FontSize',8)
            set(gca,'XTick',[])
            xlabel('from','FontSize',8)
            ylabel('to','FontSize',8)
            axis square

            % tables
            %--------------------------------------------------------------
            subplot(4,nb,i + nb)
            text(0,1/2,num2str(full(exp(Ep.B{i})),' %.2f'),'FontSize',8)
            axis off
            axis square

            % PPM
            %--------------------------------------------------------------
            subplot(4,nb,i + 2*nb)
            image(64*Pp.B{i})
            set(gca,'YTick',[1:ns],'YTickLabel',DCM.Sname,'FontSize',8)
            set(gca,'XTick',[])
            title('PPM')
            axis square

            % tables
            %--------------------------------------------------------------
            subplot(4,nb,i + 3*nb)
            text(0,1/2,num2str(Pp.B{i},' %.2f'),'FontSize',8)
            axis off
            axis square

        end

    case{lower('trial-specific effects')}

        % spm_dcm_erp_results(DCM,'trial-specific effects');
        %------------------------------------------------------------------
        if ~any(isfield(DCM.Ep,{'B','N'})), return, end

        for i = 1:ns
            for j = 1:ns

                % ensure connection is enabled
                %----------------------------------------------------------
                q     = 0;
                for k = 1:nu
                    try, q = q | DCM.Ep.B{k}(i,j); end
                    try, q = q | DCM.Ep.N{k}(i,j); end
                end

                % plot trial-specific effects
                %----------------------------------------------------------
                if q
                    D     = zeros(nt,0);
                    B     = zeros(nt,1);
                    N     = zeros(nt,1);
                    for k = 1:nu
                        try, B = B + DCM.xU.X(:,k)*DCM.Ep.B{k}(i,j); end
                        try, N = N + DCM.xU.X(:,k)*DCM.Ep.N{k}(i,j); end
                    end

                    if any(B(:)), D = B;     end
                    if any(N(:)), D = [D N]; end

                    subplot(ns,ns,(i - 1)*ns + j)
                    bar(exp(D)*100)
                    title([DCM.Sname{j}, ' to ' DCM.Sname{i}],'FontSize',10)
                    xlabel('trial',  'FontSize',8)
                    ylabel('strength (%)','FontSize',8)
                    set(gca,'XLim',[0 nt + 1])
                end
            end
        end

    case{lower('Input')}

        % plot data
        % -----------------------------------------------------------------
        tU    = t - t(1);
        U     = spm_erp_u(tU/1000,DCM.Ep,DCM.M);

        subplot(2,1,1)
        plot(t,U)
        xlabel('time (ms)')
        title('input')
        axis square, grid on
        for i = 1:length(DCM.M.ons)
            str{i} = sprintf('input (%i)',i);
        end
        legend(str)

    case{lower('Response')}

        % get spatial projector
        % -----------------------------------------------------------------
        try
            U = DCM.M.U';
        catch
            U = 1;
        end

        % plot data
        % -----------------------------------------------------------------
        try
            A     = [];
            for i = 1:nt
                subplot(nt,2,2*i - 1)
                plot(t,(DCM.H{i} + DCM.R{i})*U)
                xlabel('time (ms)')
                try
                    title(sprintf('Observed (adjusted-code:%i)',xY.code(i)))
                catch
                    title(sprintf('Observed (adjusted) %i',i))
                end
                A(end + 1,:) = axis;

                subplot(nt,2,2*i - 0)
                plot(t,DCM.H{i}*U)
                xlabel('channels')
                ylabel('time (ms)')
                title('Predicted')
                A(end + 1,:) = axis;
            end
            a(1)  = min(A(:,1));
            a(2)  = max(A(:,2));
            a(3)  = min(A(:,3));
            a(4)  = max(A(:,4));
            for i = 1:nt
                subplot(nt,2,2*i - 1)
                axis(a); axis square, grid on
                subplot(nt,2,2*i - 0)
                axis(a); axis square, grid on
            end

        end

    case{lower('Response (image)')}

        % get spatial projector
        % -----------------------------------------------------------------
        try
            U = DCM.M.U';
        catch
            U = 1;
        end

        % plot data
        % -----------------------------------------------------------------
        try
            for i = 1:nt
                subplot(nt,2,2*i - 1)
                imagesc([1:ne],t,(DCM.H{i} + DCM.R{i})*U)
                xlabel('channels')
                ylabel('time (ms)')
                try
                    title(sprintf('Observed (adjusted-code:%i)',xY.code(i)))
                catch
                    title(sprintf('Observed (adjusted) %i',i))
                end
                axis square, grid on, A = axis;

                subplot(nt,2,2*i - 0)
                imagesc(1:ne,t,DCM.H{i}*U)
                xlabel('channels')
                ylabel('time (ms)')
                title('Predicted')
                axis(A); axis square, grid on
            end
        end

    case{lower('Scalp maps')}

        % get spatial projector
        % -----------------------------------------------------------------
        try
            U = DCM.M.U';
        catch
            U = 1;
        end

        try
            pos = DCM.xY.coor2D;
        catch
            [xy, label]  = spm_eeg_project3D(DCM.M.dipfit.sens, DCM.xY.modality);
            [sel1, sel2] = spm_match_str(DCM.xY.name, label);
            pos = xy(:, sel2);
        end

        ns = 5; %number of time frames

        in           = [];
        in.type      = DCM.xY.modality;
        in.f         = gcf;
        in.noButtons = 1;
        in.cbar      = 0;
        in.plotpos   = 0;

        % plot data
        % -----------------------------------------------------------------

        for i = 1:nt
            Yo  = (DCM.H{i} + DCM.R{i})*U;
            Yp  = DCM.H{i}*U;

            for j = 1:ns
                ind    = ((j-1)*floor(nb/ns)+1):j*floor(nb/ns);

                in.max = max(abs(mean(Yo(ind, :))));
                in.min = -in.max;

                in.ParentAxes = subplot(nt*2,ns,(i - 1)*2*ns + j);
                spm_eeg_plotScalpData(mean(Yo(ind, :))', pos , DCM.xY.name, in);
                title(sprintf('observed\n%s\n%.0f ms', DCM.xY.code{i}, mean(t(ind))));

                in.ParentAxes = subplot(nt*2,ns,(i - 1)*2*ns + ns + j);
                spm_eeg_plotScalpData(mean(Yp(ind, :))', pos, DCM.xY.name, in);
                title(sprintf('predicted\n%s\n%.0f ms', DCM.xY.code{i}, mean(t(ind))));
            end
        end

    case{lower('Dipoles')}

        % return if LFP
        % -----------------------------------------------------------------
        if strcmpi(DCM.xY.modality,'lfp')
            warndlg('There are no ECDs for these LFP data')
            return
        end

        % plot dipoles
        % -----------------------------------------------------------------
        switch DCM.options.spatial

            case{'ECD'}
                P            = DCM.Eg;
                P.L          = spm_cat(P.L);
                np           = size(P.L,2)/size(P.Lpos,2);
                sdip.n_seeds = 1;
                sdip.n_dip   = np*ns;
                sdip.Mtb     = 1;
                sdip.j{1}    = full(P.L);
                sdip.j{1}    = sdip.j{1}./repmat(sqrt(sum(sdip.j{1}.^2)), 3, 1);
                sdip.j{1}    = sdip.j{1}(:);
                sdip.loc{1}  = kron(ones(1,np),full(P.Lpos));
                spm_eeg_inv_ecd_DrawDip('Init', sdip)

            case{'LFP'}
                warndlg('This is a LFP spatial model')

            case{'IMG'}
                warndlg('use the render API button to see reconstruction')

            otherwise
                return
        end

end

author: spm


text:
function [nx, cx, done, rollback, evnts, cb_data, results] = cpf_nose_event_cb(...
        k, nx, cx, px, done, rollback, evnts, cb_data, cb_args, results)
%CPF_NOSE_EVENT_CB  Callback to handle NOSE events
%   [NX, CX, DONE, ROLLBACK, EVNTS, CB_DATA, RESULTS] =
%       CPF_NOSE_EVENT_CB(K, NX, CX, PX, DONE, ROLLBACK, EVNTS, ...
%                               CB_DATA, CB_ARGS, RESULTS)
%
%   Callback to handle NOSE events, triggered by event function
%   CPF_NOSE_EVENT to indicate the nose point of the continuation curve.
%
%   This function sets the msg field of the event when the nose point has
%   been found, raises the DONE.flag and sets the DONE.msg.
%
%   See CPF_DEFAULT_CALLBACK for details of the input and output arguments.

%   MATPOWER
%   Copyright (c) 2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Shrirang Abhyankar, Argonne National Laboratory
%
%   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.

%% skip if initialize, finalize or done
if k <= 0 || done.flag
    return;
end

%% initialize return value
stop_at = cb_data.mpopt.cpf.stop_at;

%% handle event
if ~rollback || nx.step == 0
    for i = 1:length(evnts)
        if strcmp(evnts(i).name, 'NOSE') && evnts(i).zero
            if nx.step == 0
                evnts(i).msg = ...
                    sprintf('Nose point eliminated by limit induced bifurcation at %d continuation steps, lambda = %.4g.', k, nx.lam);
            else
                evnts(i).msg = ...
                    sprintf('Reached steady state loading limit in %d continuation steps, lambda = %.4g.', k, nx.lam);
            end

            %% the following conditional is only necessary if we also allow
            %% finding the location of the nose-point without terminating
            if ischar(stop_at) && strcmp(stop_at, 'NOSE');
                done.flag = 1;
                done.msg = evnts(i).msg;
            end
            break;
        end
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function F= GetMeyerBasedFilter(level,dBand,dataClass)
% Generates windowing filters
%Input:
%        level      : In multilevel decomposition number of level
%        dataClass  : 'single' or 'double'
%Output:
%        F          : Windowing filter at different level
%                   :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function F= GetMeyerBasedFilter(level,dBand,filterSize,dataClass)

F=cell(3,level);

for l=1:level

cubeSize=filterSize(l);

numDir=dBand{level}{l}(1,1);

[P,PF]=GeneratePyramidSection(cubeSize);
SF=PF{1}+PF{2}+PF{3};
PF{1}=SF;
PF{2}=SF;
PF{3}=SF;
shift=floor(cubeSize/numDir);
F{3,l}=cell(numDir,numDir);

A=zeros(cubeSize,cubeSize,cubeSize,dataClass);

for c=1:3
	mRadial=ones(shift+1,cubeSize,shift+1,dataClass);
  mRadial( [1 shift+1],1:cubeSize,1: shift+1)=.5;
  mRadial(1: shift+1,1: cubeSize,[1 shift+1])=.5;
  mRadial(1,1: cubeSize,1)=1/3;

  l2=1;
  l1=1;
  mRadialIdx=[(l1-1)*shift+1  l1*shift+1 ;1 cubeSize ; (l2-1)*shift+1  l2*shift+1  ];

  F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
  A=A+F{c,l}{l2,l1};

  mRadial=ones(shift+1,cubeSize,shift+1,dataClass);
  mRadial( [1 shift+1],1:cubeSize,1: shift+1)=.5;
  mRadial(1: shift+1,1: cubeSize,[1 shift+1])=.5;

   for l2=2:numDir-1
    for l1=2:numDir-1
      mRadialIdx=[(l1-1)*shift+1  l1*shift+1 ;1 cubeSize ; (l2-1)*shift+1  l2*shift+1  ];

      F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
      A=A+F{c,l}{l2,l1};
    end
   end

  mRadial=ones(shift,cubeSize,shift+1,dataClass);
  mRadial( [1 shift],1:cubeSize,1: shift+1)=.5;
  mRadial(1: shift,1: cubeSize,[1 shift+1])=.5;
  mRadial(shift,1: cubeSize,1)=1/3;

  l2=1;
  l1=numDir;
  mRadialIdx=[(l1-1)*shift+1  l1*shift ;1 cubeSize ; (l2-1)*shift+1  l2*shift+1  ];

  F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
    A=A+F{c,l}{l2,l1};

  mRadial=ones(shift,cubeSize,shift+1,dataClass);
  mRadial( [1 shift],1:cubeSize,1: shift+1)=.5;
  mRadial(1: shift,1: cubeSize,[1 shift+1])=.5;

  for l2=2:numDir-1
     l1=numDir;
      mRadialIdx=[(l1-1)*shift+1  l1*shift ;1 cubeSize ; (l2-1)*shift+1  l2*shift+1  ];

      F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
      A=A+F{c,l}{l2,l1};
  end
  mRadial=ones(shift+1,cubeSize,shift+1,dataClass);
  mRadial( [1 shift+1],1:cubeSize,1: shift+1)=.5;
  mRadial(1: shift+1,1: cubeSize,[1 shift+1])=.5;
  for l2=2:numDir-1
     l1=1;
      mRadialIdx=[(l1-1)*shift+1  l1*shift+1 ;1 cubeSize ; (l2-1)*shift+1  l2*shift+1  ];

      F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
      A=A+F{c,l}{l2,l1};
   end

  mRadial=ones(shift+1,cubeSize,shift,dataClass);
  mRadial( [1 shift+1],1:cubeSize,1: shift)=.5;
  mRadial(1: shift+1,1: cubeSize,[1 shift])=.5;
  mRadial(1,1: cubeSize,shift)=1/3;

  l1=1;
  l2=numDir;
  mRadialIdx=[(l1-1)*shift+1  l1*shift+1 ;1 cubeSize ; (l2-1)*shift+1  l2*shift ];

  F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
  A=A+F{c,l}{l2,l1};

  mRadial=ones(shift+1,cubeSize,shift,dataClass);
  mRadial( [1 shift+1],1:cubeSize,1: shift)=.5;
  mRadial(1: shift+1,1: cubeSize,[1 shift])=.5;

  for l1=2:numDir-1
    l2=numDir;
    mRadialIdx=[(l1-1)*shift+1  l1*shift+1 ;1 cubeSize ; (l2-1)*shift+1  l2*shift ];

    F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
    A=A+F{c,l}{l2,l1};
  end
  mRadial=ones(shift+1,cubeSize,shift+1,dataClass);
  mRadial( [1 shift+1],1:cubeSize,1: shift+1)=.5;
  mRadial(1: shift+1,1: cubeSize,[1 shift+1])=.5;

  for l1=2:numDir-1
  l2=1;
     mRadialIdx=[(l1-1)*shift+1  l1*shift+1 ;1 cubeSize ; (l2-1)*shift+1  l2*shift+1 ];

      F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
      A=A+F{c,l}{l2,l1};
  end

  mRadial=ones(shift,cubeSize,shift,dataClass);
  mRadial( [1 shift],1:cubeSize,1: shift)=.5;
  mRadial(1: shift,1: cubeSize,[1 shift])=.5;
  mRadial( shift,1: cubeSize,shift)=1/3;
  l1=numDir;
  l2=numDir;
  mRadialIdx=[(l1-1)*shift+1  l1*shift ;1 cubeSize ; (l2-1)*shift+1  l2*shift ];

  F{c,l}{l2,l1}=PolarToRec(mRadialIdx,mRadial,cubeSize,P{c},PF{c} );
  A=A+F{c,l}{l2,l1};
end

for c=1:3
   for l2=1:numDir
     for l1=1:numDir
       F{c,l}{l2,l1}=F{c,l}{l2,l1}./A;
        F{c,l}{l2,l1}=squeeze(real(fftshift(ifftn(fftshift(F{c,l}{l2,l1})))));
     end
   end
end
end
Temp=F(1,:);
F(1,:)=F(2,:);
F(2,:)=Temp;

author: thomaskuestner


text:
function Conn = hlp_absvalsq(Conn,methods,force,verb,absonly)
%
% Compute the square of the absolute value of a complex number. If
% 'force'=false then real-valued data is returned unmodified
%
% Inputs:
%
%   Conn        Connectivity structure as computed by est_mvarConnectivity().
%   methods     Cell array of strings denoting which connectivity measures
%               (fields of Conn) to transform.
%   force       force each measure to be squared, irregardless of whether
%               or not it is complex
% Outputs:
%
%   Conn        Transformed connectivity structure.
%
% See Also: pop_est_mvarConnectivity()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift/
%
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    if nargin<3
        force = 0;
    end

    if nargin<4
        verb = 1;
    end

    if nargin<5
        absonly = false;
    end

    if ~iscell(methods)
        error('methods must be a cell matrix of strings');
    end

    for cnd=1:length(Conn)
        for i=1:length(methods)
            if ~isfield(Conn(cnd),methods{i})
                if verb, fprintf('Conn(%d).%s does not exist. Skipping...\n',cnd,methods{i}); end
                continue;
            end
            if isreal(Conn(cnd).(methods{i})) && ~force
               if verb, fprintf('Conn(%d).%s is real, skipping (use ''force'' argument to override)\n',cnd,methods{i}); end
                continue;
            end

            if absonly
                Conn(cnd).(methods{i}) = abs(Conn(cnd).(methods{i}));
            else
                Conn(cnd).(methods{i}) = abs(Conn(cnd).(methods{i})).^2;
            end

%             Conn(cnd).(methods{i}) = Conn(cnd).(methods{i}).*conj(Conn(cnd).(methods{i}));
        end
    end

repo: REST


text:
function f = rhsfn(u)

  mu = -1-2*i;
  f = mu*u;

repo: Aero-matlab


text:
function h=mybar(f,s,xstr,ystr,col,lineWidth,w,we)
%
% mybar(f,[s],[xstr],[ystr],[color],[lineWidth],[barwidth],[errorbarwidth])
%
%  Draws a bar graph with heights determined by f
%  and error-bars determined by s.
%
%  if input matrices are 2d, the ROWS are plotted as
%  separate sub-bars in different colors.
%
%  xlabels set to 'xstr' (columns)
%  legend set to 'ystr' (rows)
%
%  color is a string of colors, one for each subbar
%  (default is 'brgmyc')
%
%  lineWidth (default 2)
%
%  barwidth (default 0.65)
%  errorbarWidth (default 0.3)
%
%  Example:
%    f = rand(5,1)+1;
%    s = rand(5,1)/2;
%    xstr = str2mat('one','two','three','four','five');
%    ystr = '';
%    mybar(f,s,xstr,ystr,'g');
%  7/11/95  gmb Wrote it.
%  4/11/96  gmb Revised it to plot 2d data sets
%  11/19/97 gmb Converted it to Matlab version 5.0
%  12/30/97 djh Fixed bug in xlabels, for some reason it was
%               putting the 1st label under the 2nd bar, and so on.
%  1/24/97  djh Added lineWidth as optional input arg
%  3/15/01  dbr Made legend size allocation a bit more intelligent.
%
%  05/2004  arw Modified so that if you pass in a list of colors for a
%  single data vector, it applies those colors to the individual bars.
%           EXAMPLE:
%          a=rand(4,1);b=rand(4,1)/5; h=mybar(a,b,[],[],'rgby');
% OR:      mybar(a,b,[],[],[1 0 0;0 1 0;0 0 1;1 1 0]);
%
%  06/2005  ras now allows col, xstr, and ystr to be entered as cells

%%%%%% parse inputs %%%%%
if nargin==1, s=zeros(size(f));                                 end
if ~exist('col','var') | isempty(col), col='brgmyck';           end
if ~exist('lineWidth','var') | isempty(lineWidth), lineWidth=2; end
if ~exist('w','var') | isempty(w), w=0.65;                      end
if ~exist('we','var') | isempty(we), we=0.3;                    end

if(size(f,1)==1),
    % row vector instead of column -- be flexible:
    f=f';
    s=conj(s');
end

if (isreal(s)),   s = s+sqrt(-1)*s;                             end

% if col is a cell, convert to N x 3 color table
if iscell(col),   col = cell2ColorOrder(col);                   end
% size(col)
% size(f)

% ensure col spans all columns in f -- repeat colors if needed
if ischar(col)
    while length(col) < max(size(f))
        col = [col,col];
    end
else
    while size(col,1) < max(size(f))
        col = [col; col];
    end
end

sizeShift=.01; % creates space between bars in histogram ADDED LGA 10-11-05

nsubbars=size(f,2);
we=we/nsubbars;

heldstate=ishold;
plot(0,0)
set(gca,'XLim',[1-w,size(f,1)+w]);

yhilim = max(max(f+s))*1.2;
if ~isnan(yhilim)
    %set(gca,'YLim',[0,yhilim]);
end

subx=linspace(-w/2,w/2,nsubbars+1);

newplot(gca)
hold off;
hIndex=1;
for subbar=1:nsubbars
    x=[];
    y=[];
    for i=1:size(f,1);
        tempx= [i+subx(subbar),i+subx(subbar),(i+subx(subbar+1))-sizeShift,(i+subx(subbar+1))-sizeShift];
        tempy= [0,f(i,subbar),f(i,subbar),0];

        % I can find a way to add some spacing here is I want!!! - LGA
        % 101105

        if size(f,2)==1
            % single column -- each bar a diff't color
            if ischar(col)
                thisCol=col(i);
            else
                thisCol=col(i,:);

            end

            h{hIndex}=fill(tempx',tempy',thisCol);
            hold on;

        else
            % multiple columns -- each sub-bar a diff't color
            if ischar(col)
                thisCol=col(subbar);
            else
                thisCol=col(subbar,:);
            end

            h{hIndex}=fill(tempx',tempy',thisCol);
            hold on;

        end
        hIndex=hIndex+1;
        x=[x';tempx]';
        y=[y';tempy]';
    end
    %   if ((length(col)==size(f,1)) & (size(f,2)==1))
    %      % do somethign!
    %   end
    %
    %   h(subbar)={fill(x,y,col(subbar))};
    %
    a=line(get(gca,'XLim'),[0,0]);
    set(a,'Color','k');

    errx=mean([subx(subbar),subx(subbar+1)]);

    for i=1:size(f,1)

        line([i+subx(subbar),i+subx(subbar),i+subx(subbar+1)-sizeShift,i+subx(subbar+1)-sizeShift], ...
            [0,f(i,subbar),f(i,subbar),0],'Color','k','lineWidth',2);

        g=line([i+errx-sizeShift/2,i+errx-sizeShift/2],[f(i,subbar)-imag(s(i,subbar)),f(i,subbar)+real(s(i,subbar))],...
            'Color','k','lineWidth',lineWidth);
        g=line([i+errx-we/2-sizeShift/2,i+errx+we/2-sizeShift/2],[f(i,subbar)-imag(s(i,subbar)),f(i,subbar)-imag(s(i,subbar))],...
            'Color','k','lineWidth',lineWidth);
        g=line([i+errx-we/2-sizeShift/2,i+errx+we/2-sizeShift/2],[f(i,subbar)+real(s(i,subbar)),f(i,subbar)+real(s(i,subbar))],...
            'Color','k','lineWidth',lineWidth);
    end
    %     for i=1:size(f,1)
    %         line([i+subx(subbar),i+subx(subbar),i+subx(subbar+1),i+subx(subbar+1)], ...
    %             [0,f(i,subbar),f(i,subbar),0],'Color','k','lineWidth',lineWidth);
    %         g=line([i+errx,i+errx],[f(i,subbar)-imag(s(i,subbar)),f(i,subbar)+real(s(i,subbar))],...
    %             'Color','k','lineWidth',lineWidth);
    %         g=line([i+errx-we/2,i+errx+we/2],[f(i,subbar)-imag(s(i,subbar)),f(i,subbar)-imag(s(i,subbar))],...
    %             'Color','k','lineWidth',lineWidth);
    %         g=line([i+errx-we/2,i+errx+we/2],[f(i,subbar)+real(s(i,subbar)),f(i,subbar)+real(s(i,subbar))],...
    %             'Color','k','lineWidth',lineWidth);
    %     end
    %hold on
end 					%subbars

% Xlabels
set(gca,'xLimMode','manual');
set(gca,'XTick',[1:size(f,1)]);
if (exist('xstr')) & ~isempty(xstr)
    if ~isempty(xstr)
        set(gca,'XTickLabel',xstr);
    end
end

% Legend
if (exist('ystr')) & ~isempty(ystr)

    %heighten the graph
    ylim=get(gca,'Ylim');
    nLines = length(ystr);
    yMult = 1 + nLines/4; % Each legend line uses about 1/8th of the window height
    set(gca,'YLim',[ylim(1),ylim(2)*yMult]);
    bs=8;

    posgca=get(gca,'Position');
    posgcf=get(gcf,'Position');
    dx=diff(get(gca,'XLim'))/(posgca(3)*posgcf(3));
    dy=diff(get(gca,'YLim'))/(posgca(4)*posgcf(4));

    bw=bs*dx; 				%box width
    bh=bs*dy; 				%box height

    yspace=3.5*bh; 			%vertical spacing

    xlim=get(gca,'Xlim');
    ylim=get(gca,'Ylim');

    xc=xlim(1)+50*dx;
    yc=ylim(2)-30*dy;
    bx=[xc-bw,xc+bw,xc+bw,xc-bw,xc-bw];
    by=[yc+bh,yc+bh,yc-bh,yc-bh,yc+bh];

    textx=xc+bw+10*dx;
    texty=yc;

    %draw box and text
    for i=1:nLines
        if iscell(ystr)
            label = ystr{i};
        else
            label = ystr(i,:);
        end
        if size(col, 1) == 1
            fill(bx,by,col(i));
        else
            fill(bx,by,col(i,:));
        end
        line(bx,by,'lineWidth',lineWidth,'Color','k');
        text(textx,texty,label);
        by=by-yspace;
        texty=texty-yspace;

    end
end
set(gca,'XLim',[1-w,size(f,1)+w]);

if (heldstate==0)
    hold off
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [] = anim8_DIC_images_corr_faces_n_n(IMset,DIC_2Dpair_results,varargin)
%% function for plotting 2D-DIC results imported from Ncorr in step 2
% called inside plotNcorrPairResults
% plotting the images chosen for stereo DIC (2 views) with the
% triangular faces results plotted on top, colored as their correlation
% coefficient.
% on the left side the images from the reference camera (reference image and current images), and on the right side the
% images from the deformed camera
% requirements: GIBBON toolbox
%
% calling options:
% [] = anim8_DIC_images(IMset,DIC_2Dpair_results);
% [] = anim8_DIC_images(IMset,DIC_2Dpair_results,CorCoeffCutOff,CorCoeffDispMax);
%
% INPUT:
% * IMset - a 2nX1 cell array containing 2n grayscale images. The first n
% images are from camera A (the "reference" camera), and the last n images
% are from camera B (the "deformed" camera). The first image in the set is
% considered as the reference image, on which the reference grid of points
% is defined, and all the correlated points and consequent displacements
% and strains, are relative to this image.
% * DIC_2Dpair_results - containig the correlated points, correlation
% coefficients, faces..
% * optional: CorCoeffCutOff - - maximal correlation coefficient to plot
% points
% * optional: CorCoeffDispMax - maximal correlation coefficient in colorbar

%%
nCur=numel(IMset); % number of frames

nSteps=nCur/2; %Number of animation steps
if  rem(nSteps,1)~=0
    error('Number of images in the set should be even');
end

%%
Points=DIC_2Dpair_results.Points;
CorCoeffVec=DIC_2Dpair_results.CorCoeffVec;
nCamRef=DIC_2Dpair_results.nCamRef;
nCamDef=DIC_2Dpair_results.nCamDef;
nImages=DIC_2Dpair_results.nImages;
F=DIC_2Dpair_results.Faces;

nVars = length(varargin);
switch nVars
    case 2
        CorCoeffCutOff=varargin{1};
        if isnan(CorCoeffCutOff)
            CorCoeffCutOff=max(max([CorCoeffVec{:}]));
        end
        CorCoeffDispMax=varargin{2};
        if isnan(CorCoeffDispMax)
            CorCoeffDispMax=max(max([CorCoeffVec{:}]));
        end
    case 1
        CorCoeffCutOff=varargin{1};
        if isnan(CorCoeffCutOff)
            CorCoeffCutOff=max(max([CorCoeffVec{:}]));
        end
        CorCoeffDispMax=max(max([CorCoeffVec{:}]));
    case 0
        CorCoeffCutOff=max(max([CorCoeffVec{:}]));
        CorCoeffDispMax=max(max([CorCoeffVec{:}]));
    otherwise
        error('Wrong number of input arguments');
end

%%
hf=cFigure;
hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

% Ref
ii=1;
subplot(1,2,1)
CFcorr=mean(CorCoeffVec{1}(F),2);
CFcorr(CFcorr>CorCoeffCutOff)=NaN;
hp1=imagesc(repmat(IMset{ii},1,1,3)); hold on;
hp2=gpatch(F,Points{1},CFcorr,'k',0.4);
set(hp2,'EdgeAlpha',0.4);
colormap jet
pbaspect([size(IMset{ii},2) size(IMset{ii},1) 1])
hs1=title(['Ref (Cam ' num2str(nCamRef) ' frame ' num2str(1) ')']);
hc1=colorbar;
caxis([0 CorCoeffDispMax])
title(hc1, 'Corr-Coeff')
axis off

% Cur
ii=nImages+1;
subplot(1,2,2)
CFcorr=mean(CorCoeffVec{1+nCur/2}(F),2);
CFcorr(CFcorr>CorCoeffCutOff)=NaN;
hp3=imagesc(repmat(IMset{ii},1,1,3)); hold on
hp4=gpatch(F,Points{ii},CFcorr,'k',0.4); %colormap jet
colormap jet
pbaspect([size(IMset{ii},2) size(IMset{ii},1) 1])
hs2=title(['Cur ' num2str(ii) ' (Cam ' num2str(nCamDef) ' frame ' num2str(1) ')']);
hc2=colorbar;
caxis([0 CorCoeffDispMax])
title(hc2, 'Corr-Coeff')
axis off

drawnow

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

for ii=1:nImages
    Pnow1=Points{ii};
    Pnow2=Points{ii+nImages};

    cNow1=CorCoeffVec{ii};
    cNow1(cNow1>CorCoeffCutOff)=NaN;
    cNow2=CorCoeffVec{ii+nImages};
    cNow2(cNow2>CorCoeffCutOff)=NaN;

    TitleNow1=['Cur ' num2str(ii) ' (Cam ' num2str(nCamRef) ' frame ' num2str(ii) ')'];
    TitleNow2=['Cur ' num2str(ii) ' (Cam ' num2str(nCamDef) ' frame ' num2str(ii) ')'];

   %Set entries in animation structure
    animStruct.Handles{ii}=[hp1,hp3,hp2,hp2,hp4,hp4,hs1,hs2]; %Handles of objects to animate
    animStruct.Props{ii}={'CData','CData','Vertices','CData','Vertices','CData','String','String'}; %Properties of objects to animate
    animStruct.Set{ii}={repmat(IMset{ii},1,1,3),repmat(IMset{ii+nImages},1,1,3),Pnow1,cNow1,Pnow2,cNow2,TitleNow1,TitleNow2}; %Property values for to set in order to animate

end

anim8(hf,animStruct);

addColorbarLimitsButton(hf);
addColormapButton(hf);
addEdgeColorButton(hf);
addFaceAlphaButton(hf);
addLightButton(hf);
addAmbientStrengthButton(hf);

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% Modified by Rana Odabas 2018
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function test_libltfat_gabframediag
[~,~,enuminfo]=libltfatprotofile;
LTFAT_FIRWIN = enuminfo.LTFAT_FIRWIN;

a = 14;
gl = 34;
M = 64;
g = zeros(gl,1);
d = zeros(a,1);
gPtr = libpointer('doublePtr',g);
dPtr = libpointer('doublePtr',d);

calllib('libltfat','ltfat_firwin_d',LTFAT_FIRWIN.LTFAT_HANN,gl,gPtr);

calllib('libltfat','ltfat_gabframediag_d',gPtr,gl,a,M,a,dPtr);

d =gabframediag(gPtr.Value,a,M,lcm(a,M));
d(1:a)-dPtr.Value

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function tapas_hgf_plotTraj(r)
% Plots the estimated or generated trajectories for the HGF perceptual model
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_hgf_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Optional plotting of standard deviations (true or false)
plotsd = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Time axis
if size(r.u,2) > 1 && ~isempty(find(strcmp(fieldnames(r.c_prc),'irregular_intervals'))) && r.c_prc.irregular_intervals
    t = r.u(:,end)';
else
    t = ones(1,size(r.u,1));
end

ts = cumsum(t);
ts = [0, ts];

% Number of levels
l = length(r.p_prc.p)/5;

% Upper levels
for j = 1:l-1

    % Subplots
    subplot(l,1,j);

    if plotsd == true
        upperprior = r.p_prc.mu_0(l-j+1) +sqrt(r.p_prc.sa_0(l-j+1));
        lowerprior = r.p_prc.mu_0(l-j+1) -sqrt(r.p_prc.sa_0(l-j+1));
        upper = [upperprior; r.traj.mu(:,l-j+1)+sqrt(r.traj.sa(:,l-j+1))];
        lower = [lowerprior; r.traj.mu(:,l-j+1)-sqrt(r.traj.sa(:,l-j+1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mu_0(l-j+1); r.traj.mu(:,l-j+1)], 'b', 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu_0(l-j+1), 'ob', 'LineWidth', 2); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu_', num2str(l-j+1)]);
end

% Input level
subplot(l,1,l);

if plotsd == true
    upperprior = r.p_prc.mu_0(1) +sqrt(r.p_prc.sa_0(1));
    lowerprior = r.p_prc.mu_0(1) -sqrt(r.p_prc.sa_0(1));
    upper = [upperprior; r.traj.mu(:,1)+sqrt(r.traj.sa(:,1))];
    lower = [lowerprior; r.traj.mu(:,1)-sqrt(r.traj.sa(:,1))];

    plot(0, upperprior, 'or', 'LineWidth', 1);
    hold all;
    plot(0, lowerprior, 'or', 'LineWidth', 1);
    fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
         'r', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
end
plot(ts, [r.p_prc.mu_0(1); r.traj.mu(:,1)], 'r', 'LineWidth', 2);
hold all;
plot(0, r.p_prc.mu_0(1), 'or', 'LineWidth', 2); % prior
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0.6 0]); % inputs
if ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    plot(ts(2:end), r.y(:,1), '.', 'Color', [1 0.7 0]); % responses
    title(['Response y (orange), input u (green), and posterior expectation of x_1 ', ...
           '(red) for \rho=', num2str(r.p_prc.rho), ', \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ...
           ', \pi_u=', num2str(r.p_prc.pi_u)], ...
          'FontWeight', 'bold');
    ylabel('y, u, \mu_1');
else
    title(['Input u (green) and posterior expectation of x_1 ', ...
           '(red) for \rho=', num2str(r.p_prc.rho), ', \kappa=', ...
           num2str(r.p_prc.ka), ', \omega=', num2str(r.p_prc.om), ...
           ', \pi_u=', num2str(r.p_prc.pi_u)], ...
          'FontWeight', 'bold');
    ylabel('u, \mu_1');
end
xlim([0 ts(end)]);
xlabel({'Trial number', ' '}); % A hack to get the relative subplot sizes right
hold off;

author: translationalneuromodeling

repo: tapas


text:
function [out] = recharge_7(p1,fin)
%recharge_7

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Constant recharge limited by incoming flux
% Constraints:  -
% @(Inputs):    p1   - maximum recharge rate [mm/d]
%               fin  - incoming flux [mm/d]

out = min(p1,fin);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function cellInfo = Auto_Cmd_Bus(varargin)
% AUTO_CMD_BUS returns a cell array containing bus object information
%
% Optional Input: 'false' will suppress a call to Simulink.Bus.cellToObject
%                 when the MATLAB file is executed.
% The order of bus element attributes is as follows:
%   ElementName, Dimensions, DataType, SampleTime, Complexity, SamplingMode, DimensionsMode, Min, Max, DocUnits, Description

suppressObject = false;
if nargin == 1 && islogical(varargin{1}) && varargin{1} == false
    suppressObject = true;
elseif nargin > 1
    error('Invalid input argument(s) encountered');
end

cellInfo = { ...
  { ...
    'Auto_Cmd_Bus', ...
    '', ...
    '', ...
    'Auto', ...
    '-1', {...
{'timestamp', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'p_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad/s'), sprintf('rate x command in body frame')}; ...
{'q_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad/s'), sprintf('rate y command in body frame')}; ...
{'r_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad/s'), sprintf('rate z command in body frame')}; ...
{'phi_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), sprintf('roll command')}; ...
{'theta_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), sprintf('pitch command')}; ...
{'psi_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'psi_rate_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad/s'), sprintf('yaw rate command')}; ...
{'x_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'y_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'z_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'lat_cmd', 1, 'int32', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'lon_cmd', 1, 'int32', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'alt_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'u_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), sprintf('velocity x command in control frame')}; ...
{'v_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), sprintf('velocity y command in control frame')}; ...
{'w_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), sprintf('velocity z command in control frame')}; ...
{'ax_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'ay_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'az_cmd', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'throttle_cmd', 1, 'uint16', -1, 'real', 'Sample', 'Fixed', [], [], '', sprintf('throttle command')}; ...
{'frame', 1, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', sprintf('Coordinate Frame:\n0:FRAME_GLOBAL_NED\n1:FRAME_LOCAL_FRD\n2:FRAME_BODY_FRD')}; ...
{'reserved', 1, 'uint8', -1, 'real', 'Sample', 'Fixed', [], [], '', ''}; ...
{'cmd_mask', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], '', sprintf('Type mask for auto command:\n  1: p_cmd valid\n  2: q_cmd valid\n  3: r_cmd valid\n  4: phi_cmd valid\n  5: theta_cmd valid\n  6: psi__cmd_valid\n  7: psi_rate_cmd_valid\n  8: x_cmd valid\n  9: y_cmd valid\n10: z_cmd valid\n11: lat_cmd valid\n12: lon_cmd valid\n13: alt_cmd valid\n14: u_cmd valid\n15: v_cmd valid\n16: w_cmd valid\n17: ax_cmd valid\n18: ay_cmd valid\n19: ax_cmd valid\n20: throttle_cmd valid')}; ...
    } ...
  } ...
}';

if ~suppressObject
    % Create bus objects in the MATLAB base workspace
    Simulink.Bus.cellToObject(cellInfo)
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function [kip] = lbf2kip(lbf)
% Convert force from pounds-force to kip.
% Chad A. Greene 2012
kip = lbf*.001;


text:
function [gx] = g_dummy(x,P,u,in)
gx = 1./(1+exp(-x./P));

author: MBB-team

repo: VBA-toolbox


text:
close all
clear

%------------------------------------------------------------------------------------------
nb = 1000;%Total number of initial states
nmpcd=zeros(1,nb);%Store computing time for one time-step for NMPC-DCBF
impcd=zeros(1,nb);%Store computing time for one time-step for iMPC-DCBF
nmpcinf=zeros(1,1);%Store infeasible rate for one time-step for NMPC-DCBF
impcinf=zeros(1,1);%Store infeasible rate for one time-step for iMPC-DCBF
i=24;%i is number of horozon
ts = 1;%Total time steps
load(gamma1_0p6\InitialStateData');%Load the generated 1000 initial states
[a, b, c, d]=compare(i,ts,nb,InitialMat);%Involves NMPC-DCBF and iMPC-DCBF
nmpcd(1,:)=a;
impcd(1,:)=b;
nmpcinf(1,1)=c;
impcinf(1,1)=d;
nmpcplot1=nmpcd(1,:);
nmpcplot1(nmpcplot1==-1)=[];%Eliminate infeasible one time-step trajectories
impcplot1=impcd(1,:);
impcplot1(impcplot1==-1)=[];%Eliminate infeasible one time-step trajectories

save('timecom24','nmpcplot1','impcplot1');
save('feasibility24','nmpcinf','impcinf');

distnmpc = fitdist(nmpcplot1','Normal');
distimpc = fitdist(impcplot1','Normal');
mu1=distnmpc.mu;%Get mean of sample of computing time for NMPC-DCBF
mu2=distimpc.mu;%Get mean of sample of computing time for iMPC-DCBF
sigma1=distnmpc.sigma;%Get variance of sample of computing time for NMPC-DCBF
sigma2=distimpc.sigma;%Get variance of sample of computing time for iMPC-DCBF

%Initialize atmosphere parameters
function [tnmpc, timpc, ratiotnmpc, ratiotimpc]=compare(N11,ttsim,samplen,InitialMat)
tnmpc=[];
timpc=[];
for i=1:samplen
N1=N11;
tsim=ttsim;
xini1=InitialMat(1,i);
yini1=InitialMat(2,i);
thetaini1=InitialMat(3,i);
vini1=InitialMat(4,i);
x01=[xini1;yini1;thetaini1;vini1];
x02=[xini1 yini1 thetaini1 vini1];
t1 = nmpcdcbf(x01, N1, tsim);
t2 = impcdcbf(x02, N1, tsim);
tindex=[N11 i];
disp(tindex);
tnmpc=[tnmpc t1];%Computing time for NMPC-DCBF
timpc=[timpc t2];%Computing time for iMPC-DCBF
end
nnmpc1 = length(tnmpc);
nimpc1 = length(timpc);
tnmpcs = tnmpc;
timpcs = timpc;
tnmpcs(tnmpcs==-1)=[];
timpcs(timpcs==-1)=[];
nnmpc2 = length(tnmpcs);
nimpc2 = length(timpcs);
ratiotnmpc = (nnmpc1-nnmpc2)/nnmpc1;%Infeasible rate for NMPC-DCBF
ratiotimpc = (nimpc1-nimpc2)/nimpc1;%Infeasible rate for iMPC-DCBF
end

function t1 = nmpcdcbf(x00, N1, ttsim)
%% General Flags
run_nmpc_dcbf_one = true;
run_nmpc_dcbf_multiple = true;

%% Setup and Parameters
x0 = x00;
dt = 0.1;
time_total = ttsim *dt;
P = zeros(4,4);%Weight matrix P
P(1,1) = 10; P(2,2) = 10; P(3,3) = 10;P(4,4) = 10;
Q = zeros(4,4);%Weight matrix Q
Q(1,1) = 10; Q(2,2) = 10; Q(3,3) = 10;Q(4,4) = 10;
R = zeros(4,4);%Weight matrix R
R(1,1) = 1; R(2,2) = 1;R(3,3) = 1000;R(4,4) = 1000;
%Variables range as below
xmin = [-10; -10; -10;-10];
xmax = [10; 10; 10;10];
umin = [-7; -5;-inf;-inf];
umax = [7; 5;inf;inf];

%% Discrete-time unicycle model
system.dt = dt;
system.A = [1 0 0 0;
    0 1 0 0;
    0 0 1 0;
    0 0 0 1];
system.B = [x0(4,1)*cos(x0(3,1))*dt;
    x0(4,1)*sin(x0(3,1))*dt;
    0;
    0];
system.C =[0 0 0 0;
    0 0 0 0;
    1*dt 0 0 0;
    0 1*dt 0 0];
system.xl = xmin;
system.xu = xmax;
system.ul = umin;
system.uu = umax;

%% NMPC-DCBF parameters
params_nmpc_dcbf.Q = Q;
params_nmpc_dcbf.R = R;
params_nmpc_dcbf.P = P;
params_nmpc_dcbf.gamma1 = 0.6;
params_nmpc_dcbf.gamma2 = 0.6;

%% Obstacle
obs.pos1 = [0; 0];
obs.r1 = 1;

%% Simulate NMPC-DCBF with other N
params_nmpc_dcbf.N = N1;
if run_nmpc_dcbf_one
    fprintf('Run NMPC-DCBF-24\n');
    controller_nmpc_dcbf_multiple24 = NMPCDCBF1(x0, system, params_nmpc_dcbf);%Highest order mcbf=2
    controller_nmpc_dcbf_multiple24.obs = obs;
    controller_nmpc_dcbf_multiple24.sim(time_total);
end
t1=controller_nmpc_dcbf_multiple24.tt;
end

function t2=impcdcbf(x00, N1, ttsim)
t2=0;
xo = 0;
yo = 0;
r = 1;
N = N1;%Number of horizon
K1 = 1000;%Maximum iteration times, jmax
K = 1000;%Maximum iteration times, jmax
dt = 0.1;
gamma1 = 6;
gamma2 = 6;
nsim = ttsim;

tic;

% Constraints

umin = [-7; -5; -inf; -inf;];
umax = [7; 5; Inf; Inf];
xmin = [-10; -10; -10; -10];
xmax = [ 10;  10;  10;  10];

% Objective function
Q = diag([10 10 10 10]);
QN = Q;
R = 1 * eye(4);
R(3,3) = 1000;
R(4,4) = 1000;
% Initial and reference states
x0 = x00;
xr = [3; 0.01; 0; 0];
ur = [0; 0; 1; 1];

% Dynamic system initialization
BB = [0 0 0 0;0 0 0 0;dt 0 0 0;0 dt 0 0];

% Convex MPC frame
[nx, nu] = size(BB);
% Initialize states set (x00,x01,,,x0N）
x_0 = [];
x0 = x0';
u_0 = zeros(nu, N);
u0 = zeros(nu, 1);% we may need to change this for a better warm-up start
impc = zeros(nx, nsim + 1);
impc(:, 1) = x0;
x0new = x0;
for i=1 : (N+1)
    x_0 = [x_0 x0new];
    x0new = [dt*x0new(4)*cos(x0new(3))+x0new(1);dt*x0new(4)*sin(x0new(3))+x0new(2);dt*u0(1)+x0new(3);dt*u0(2)+x0new(4)];
end

abc = tangentline(x_0, r);
Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);%First order CBF constraints
lcbf = -abc(3, :);
lcbf(1) = [];
lcbf = lcbf';
ucbf = inf * ones(N, 1);
A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);%Second order CBF constraints
lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
ucbf2 = inf * ones(N-1, 1);

% Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
% - quadratic objective
P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
% - linear objective
q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
% - linear dynamics
Ax = getax(x_0, dt, N, nx);
Bu = kron([sparse(1, N); speye(N)], BB);
Cx = getcx(x_0, u_0, dt, N, nx);
Aeq = [Ax, Bu];
leq = [-x0; -Cx];
ueq = leq;
% - input and state constraints
Aineq = speye((N+1)*nx + N*nu);
lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

% - OSQP constraints
A = [Aeq; Aineq; Acbf; A2cbf];
l = [leq; lineq; lcbf; lcbf2];
u = [ueq; uineq; ucbf; ucbf2];

% Create an OSQP object
prob = osqp;
% Setup workspace
prob.setup(P, q, A, l, u, 'warm_start', true);
% Solve
res = prob.solve();

% Check solver status
if ~strcmp(res.info.status, 'solved')
%    error('OSQP did not solve the problem!')
        t2=-1;
        return
end
ctrlx = res.x(1:(N+1)*nx);
ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);

for j = 1 : (K1-1)%Iterations for the first time-step
    storagex = ctrlx;%store the state variables and inputs in order to compare them with the next iteration
    storageu = ctrlu;
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
end % Move for the first step
ctrl = ctrlu(1:nu);
x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
ctrlu = rewrite(ctrlu, nu, N);
t2=t2+toc;
return

x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
u_0 = transu(ctrlu, nu, N);
impc(:, 2) = x0;
storagex = ctrlx;
storageu = ctrlu;

for i = 1 : (nsim-1)%Iterations for the left time steps

   for j = 1 : K
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    x0 = ctrlx(1:nx);
    testx0 = [testx0 x0];
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
    storagex = ctrlx;
    storageu = ctrlu;
   end
   ctrl = ctrlu(1:nu);
   x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
   ctrlu = rewrite(ctrlu, nu, N);
   x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
   u_0 = transu(ctrlu, nu, N);
   impc(:, i+2) = x0;
end
end

% Linerize the CBF constraints (get a, b, c for lines)
function abc = tangentline(xy, r)% x and y from initialize states set, abc are coeeficients for linear equation a*x+b*y+c=0
[xx, ~] = size(xy);%xx=2,yy=N+1
xy(xx,:) = []; % this part should be changed for other case
xy((xx-1),:) = []; % this part should be changed other case
[xx, yy] = size(xy);%xx=2,yy=N+1
xyjiao = zeros(xx, yy);%intersection points
for i = 1 : xx
    for j = 1 : yy
        xyjiao(i, j) = r * xy(i, j) * (1 / (xy(:, j)' * xy(:, j)))^(0.5);%calculate coordinates of intersection points
    end
end
cc = -r^2 * ones(1, yy);
abc = [xyjiao; cc];
end

% Get CBF constraints matrix
function Acbf = cbfmat(abc, nx, nu, dt, gamma, x0)
[~, yy] = size(abc);
Acbfx = zeros((yy-1), yy*nx);
Acbfu = zeros((yy-1), (yy-1)*nu);
for i = 1 : (yy-1)
    Acbfx(i, (i*nx)+1) = abc(1, (i+1));
    Acbfx(i, (i*nx)+2) = abc(2, (i+1));
end
for i = 1 : (yy-1)
    Acbfu(i, ((i-1)*nu+3)) = - (1 - dt * gamma)^(i) * (abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
Acbf = [Acbfx Acbfu];
end

% Transfer vector x into matrix
function res = trans(x0, vector, nxx, nyy)%nxx=nx,nyy=N+1
 res = zeros(nxx, nyy);
 res(:,1) = x0;
 for i = 1 : (nyy -1)
     res(:,i+1)= vector(((i)*nxx+1):(i+1)*nxx);
 end
end

% Transfer vector u into matrix
function resu = transu(vector, nxx, nyy)%nxx=nu,nyy=N
 resu = zeros(nxx, nyy);
 for i = 1 : (nyy)
     resu(:,i)= vector(((i-1)*nxx+1):(i)*nxx);
 end
end

% Rewrite u vector
function reu = rewrite(vector, nu, N)
append = vector((N-1)*nu+1:N*nu);
vector(1:nu) = [];
reu = [vector;append];
end

% Get new x_0
function x_0 = newinit(x0, ctrlu, nx, nu, N, dt)
 x_0 = zeros(nx, N+1);
 x_0(:, 1) = x0;
 for i=1 : N
    u0 = ctrlu((i-1)*nu+1:i*nu);
    x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*u0(1)+x0(3);dt*u0(2)+x0(4)];
    x_0(:, i + 1) = x0;
 end
end

% Get AA matrix
function AA = getaa(x0, dt)
 AA = [1 0 -x0(4)*sin(x0(3))*dt cos(x0(3))*dt;0 1 x0(4)*cos(x0(3))*dt sin(x0(3))*dt;0 0 1 0;0 0 0 1];
end
% Get CC matrix
function CC = getcc(x0, x1, u0, dt)
 CC = [x0(4)*sin(x0(3))*x0(3)*dt-x0(4)*cos(x0(3))*dt+x1(1)-x0(1);-x0(4)*cos(x0(3))*x0(3)*dt-x0(4)*sin(x0(3))*dt+x1(2)-x0(2);-u0(1)*dt+x1(3)-x0(3);-u0(2)*dt+x1(4)-x0(4)];
end
% Get Ax matrix
function Ax = getax(x_0, dt, N, nx)
x0 = x_0(:,1);
AA = getaa(x0, dt);
Ax = kron(speye(N+1), -speye(nx)) + kron(sparse(diag(ones(N, 1), -1)), AA);
for i = 1 : (N-1)
   x0 = x_0(:,i+1);
   AA = getaa(x0, dt);
   Ax(nx*(i+1)+1:nx*(i+1)+nx,nx*i+1:nx*i+nx) = AA;
end
end
% Get Cx matrix
function Cx = getcx(x_0, u_0, dt, N, nx)
Cx = zeros(N*nx, 1);
for i = 1 : N
    u0 = u_0(:,i);
    x0 = x_0(:,i);
    x1 = x_0(:,i+1);
    CC = getcc(x0, x1, u0, dt);
    Cx((i-1)*nx+1:(i-1)*nx+nx) = CC;
end
end
% Get A2cbf
function A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0)
[~, yy] = size(abc);
Acbfx2 = zeros((yy-2), yy*nx);
Acbfx22 = zeros((yy-2), yy*nx);
Acbfu2 = zeros((yy-2), (yy-1)*nu);
for i = 1 : (yy-2)
    Acbfx2(i, (i*nx)+1) = (gamma1-1/dt)*abc(1, (i+1));
    Acbfx2(i, (i*nx)+2) = (gamma1-1/dt)*abc(2, (i+1));
    Acbfx2(i, ((i+1)*nx)+1) = (1/dt)*abc(1, (i+2));
    Acbfx2(i, ((i+1)*nx)+2) = (1/dt)*abc(2, (i+2));
end
for i = 1 : (yy-2)
    Acbfx22(i, (1*nx)+1) = -(1-dt*gamma2)^(i)/dt*abc(1, (1+1));
    Acbfx22(i, (1*nx)+2) = -(1-dt*gamma2)^(i)/dt*abc(2, (1+1));
end
Acbfx2 = Acbfx2 + Acbfx22;
for i = 1 : (yy-2)
    Acbfu2(i, ((i-1)*nu+4)) = - (1 - dt * gamma2)^(i)*(gamma1-1/dt)*(abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
A2cbf = [Acbfx2 Acbfu2];
end
% Get lcbf2
function lcbf2 = getlcbf2(abc, dt, gamma1, gamma2)
[~, yy] = size(abc);
lcbf2 = zeros((yy-2),1);
for i = 1 : (yy-2)
    lcbf2(i, 1) = -abc(3, (i+2))/dt-(gamma1-1/dt)*abc(3, (i+1))+(1 - dt * gamma2)^(i)*abc(3, (1+1))/dt;
end
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
function [dS,f]=mtdspectrumc(data,phi,params)
% Multi-taper frequency derivative of the spectrum - continuous process
%
% Usage:
%
% [dS,f]=mtdspectrumc(data,phi,params)
% Input:
%   Note that all times can be in arbitrary units. But the units have to be
%   consistent. So, if E is in secs, win, t have to be in secs, and Fs has to
%   be Hz. If E is in samples, so are win and t, and Fs=1. In case of spike
%   times, the units have to be consistent with the units of data as well.
%       data        (in form samples x channels/trials or a single vector) -- required
%       phi         (angle for evaluation of derivative) -- required.
%                       e.g. phi=[0,pi/2] gives the time and frequency derivatives
%       params: structure with fields tapers, pad, Fs, fpass, trialave
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           trialave (average over trials/channels when 1, don't average when 0) - optional. Default 0
% Output:
%       dS       (spectral derivative in form phi x frequency x channels/trials if trialave=0 or
%                in form phi x frequency if trialave=1)
%       f        (frequencies)

if nargin < 2; error('Need data and angle'); end;
if nargin < 3; params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear err params
data=change_row_to_column(data);
N=size(data,1);
nfft=max(2^(nextpow2(N)+pad),N);
[f,findx]=getfgrid(Fs,nfft,fpass);
tapers=dpsschk(tapers,N,Fs); % check tapers
K=size(tapers,2);
J=mtfftc(data,tapers,nfft,Fs);
J=J(findx,:,:);
A=sqrt(1:K-1);
A=repmat(A,[size(J,1) 1]);
A=repmat(A,[1 1 size(J,3)]);
S=squeeze(mean(J(:,1:K-1,:).*A.*conj(J(:,2:K,:)),2));
if trialave; S=squeeze(mean(S,2));end;
nphi=length(phi);
for p=1:nphi;
    dS(p,:,:)=real(exp(i*phi(p))*S);
end;
dS=squeeze(dS);
dS=change_row_to_column(dS);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [score] = EMRtest(x,model)
% [score] = EMRtest(x,model): Efficient Manifold Ranking for out-of-sample
%                             retrieval
% Input:
%       - x:  the query point, row vector.
%       - model:  the model learned by EMR
%
% Output:
%       - score: the ranking scores for each point in the database
%
% Usage:
%
%      See: http://www.zjucadcg.cn/dengcai/Data/Examples.html#EMR
%
%Reference:
%
%	 Bin Xu, Jiajun Bu, Chun Chen, Deng Cai, Xiaofei He, Wei Liu, Jiebo
%	 Luo, "Efficient Manifold Ranking for Image Retrieval",in Proceeding of
%	 the 34th International ACM SIGIR Conference on Research and
%	 Development in Information Retrieval (SIGIR), 2011, pp. 525-534.
%
%   version 2.0 --Feb./2012
%   version 1.0 --Sep./2010
%
%   Written by Bin Xu (binxu986 AT gmail.com)
%              Deng Cai (dengcai AT gmail.com)

r = model.r;
a = model.a;
p = size(model.landmarks,1);

% Z construction
D = EuDist2(x,model.landmarks);
[dump,idx] = sort(D);
dump = dump(1:r)/dump(r);
dump = 0.75 * (1 - dump.^2);
z = sparse(idx(1:r),1,dump,p,1);

Z = [model.Z z];
Z = Z';

nSmp =size(Z,1);

y0 = zeros(nSmp,1);
y0(end) = 1;

% Efficient Ranking
feaSum = full(sum(Z,1));
D = Z*feaSum';
D = max(D, 1e-12);
D = D.^(-.5);
H = spdiags(D,0,nSmp,nSmp)*Z;

C = speye(p);
A = H'*H-(1/a)*C;

tmp = H'*y0;
tmp = A\tmp;
score = y0 - H*tmp;

score(end) = [];

author: ZJULearning

repo: MatlabFunc


text:
function bic = ml_bic(ll,num_para,num_data)
%ML_BIC  Bayesian information criterion
%
%   input -----------------------------------------------------------------
%
%       o ll        : (1 x 1), log-likelihood.
%
%       o num_para  : (1 x 1), number of parameters.
%
%       o num_data  : (1 x 1), number of datapoints.
%
%   output ----------------------------------------------------------------
%
%       o bic       : (1 x 1)

 bic =  -2 * ll + num_para * log(num_data);

end

author: epfl-lasa