text:
% THE SPIDER Version 1.71 (24/7/2006)  -- Current Repository
%
% Basic library objects.
%   data         - Storing input data and output results
%   data_global  - Implementation of data object that limits memory overhead
%   algorithm    - Generic algorithm object
%   group        - Groups sets of objects together (algorithms or data)
%   loss         - Evaluates loss functions
%   get_mean     - Takes mean loss over groups of algs
%   chain        - Builds chains of objects: output of one to input of another
%   param        - To train and test different hyperparameters of an object
%   cv           - Cross validation using objects given data
%   kernel       - Evaluates and caches kernel functions
%   distance     - Evaluates and caches distance functions
%
% Statistical Tests objects.
%   wilcoxon     - Wilcoxon test of statistical significance of results
%   corrt_test   - Corrected resampled t-test - for dependent trials
%
% Dataset objects.
%   spiral       - Spiral dataset generator.
%   toy          - Generator of dataset with only a few relevant features
%   toy2d        - Simple 2d Gaussian problem generator
%   toyreg       - Linear Regression with o outputs and n inputs
%
% Pre-Processing objects
%   normalize    - Simple normalization of data
%   map          - General user specified mapping function of data
%
% Density Estimation objects.
%   parzen       - Parzen's windows kernel density estimator
%   indep        - Density estimator which assumes feature independence
%   bayes        - Classifer based on density estimation for each class
%   gauss        - Normal distribution density estimator
%
% Pattern Recognition objects.
%   svm          - Support Vector Machine (svm)
%   c45          - C4.5 for binary or multi-class
%   knn          - k-nearest neighbours
%   platt        - Conditional Probability estimation for margin classifiers
%   mksvm        - Multi-Kernel LP-SVM
%   anorm        - Minimize the a-norm in alpha space using kernels
%   lgcz         - Local and Global Consistent Learner
%   bagging	     - Bagging Classifier
%   adaboost     - ADABoost method
%   hmm          - Hidden Markov Model
%   loom         - Leave One Out Machine
%   l1           - Minimize l1 norm of w for a linear separator
%   kde          - Kernel Dependency Estimation: general input/output machine
%   dualperceptron       - Kernel Perceptron
%   ord_reg_perceptron   - Ordinal Regression Perceptron (Shen et al.)
%   splitting_perceptron - Splitting Perceptron (Shen et al.)
%   budget_perceptron    - Sparse, online Pereceptron (Crammer et al.)
%   randomforest - Random Forest Decision Trees         WEKA-Required
%   j48          - J48 Decision Trees for binary        WEKA-Required
%
% Multi-Class and Multi-label objects.
%   one_vs_rest  - Voting method of one against the rest (also for multi-label)
%   one_vs_one   - Voting method of one against one
%   mc_svm       - Multi-class Support Vector Machine by J.Weston
%   c45          - C4.5 for binary or multi-class
%   knn          - k-nearest neighbours
%
% Feature Selection objects.
%   feat_sel     - Generic object for feature selection + classifier
%   r2w2_sel     - SVM Bound-based feature selection
%   rfe          - Recursive Feature Elimination (also for the non-linear case)
%   l0           - Dual zero-norm minimization (Weston, Elisseeff)
%   fsv          - Primal zero-norm based feature selection (Mangasarian)
%   fisher       - Fisher criterion feature selection
%   mars         - selection algorithm of Friedman (greedy selection)
%   clustub      - Multi-class feature selection using spectral clustering
%   mutinf       - Mutual Information for feature selection.
%
% Regression objects.
%   svr          - Support Vector Regression
%   gproc        - Gaussian Process Regression
%   relvm_r      - Relevance vector machine
%   multi_rr     - (possibly multi-dimensional) ridge regression
%   mrs          - Multivariate Regression via Stiefel Constraints
%   knn          - k-nearest neighbours
%   multi_reg    - meta method for independent multiple output regression
%   kmp          - kernel matching pursuit
%   kpls         - kernel partial least squares
%   lms          - least mean squared regression [now obselete due to multi_rr]
%   rbfnet       - Radial Basis Function Network (with moving centers)
%   reptree      - Reduced Error Pruning Tree       WEKA-Required
%
% Model Selection objects.
%   gridsel      - select parameters from a grid of values
%   r2w2_sel     - Selecting SVM parameters by generalization bound
%   bayessel     - Bayessian parameter selection
%
% Unsupervised objects.
%   one_class_svm - One class SVM
%   kmeans       - K means clustering
%   kvq          - Kernel Vector Quantization
%   kpca         - Kernel Principal Components Analysis
%   ppca         - Probabilistic Principal Component Analysis
%   nmf          - Non-negative Matrix factorization
%   spectral     - Spectral clustering
%   mrank        - Manifold ranking
%   ppca         - Probabilistic PCA
%
% Reduced Set and Pre-Image objects.
%   pmg_mds      - Calculate Pre-Images based on multi-dimensional scaling
%   pmg_rr       - Calculate Pre-Images based on learning and ridge regression
%   rsc_burges   - Bottom Up Reduced Set; calculates reduced set based on gradient descent
%   rsc_fp       - Bottom Up Reduced Set; calculates reduced set for rbf with fixed-point iteration schemes
%   rsc_mds      - Top Down Reduced Set; calculates reduced set with multi-dimensional scaling
%   rsc_learn    - Top Down Reduced Set; calculates reduced set with ridge regression
%   rss_l1       - Reduced Set Selection via L1 penalization
%   rss_l0       - Reduced Set Selection via L0 penalization
%   rss_mp       - Reduced Set Selection via matching pursuit
%
% See also help on: demos, train, test.
% Most of these algorithms are available as extra download on the website:
% http://www.kyb.tuebingen.mpg.de/bs/people/spider/index.html

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = plotcoeffs(f, varargin)
%PLOTCOEFFS   Display coefficients graphically.
%   PLOTCOEFFS(F) plots the absolute values of the coefficients underlying the
%   representation of a CHEBFUN F on a semilogy scale. If F is an array-valued
%   CHEBFUN or has breakpoints, then a curve is plotted for each FUN of each
%   component (column) of F.
%
%   PLOTCOEFFS(F, S) allows further plotting options, such as linestyle,
%   linecolor, etc, as in the standard MATLAB manner. If S contains a string
%   'LOGLOG', the coefficients will be displayed on a log-log scale.
%
%   H = PLOTCOEFFS(F) returns a column vector of handles to lineseries objects.
%
%   As of Chebfun version 5, the coefficients plotted by 'plotcoeffs' are either
%   Chebyshev coefficients for Chebyshev-based chebfuns or Fourier coefficients
%   for periodic chebfuns ('trigfuns').
%
% See also CHEBFUN/PLOT.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Deal with an empty input:
if ( isempty(f) )
    if ( nargout == 1 )
        varargout{1} = plot([]);
    end
    return
end

dom = f(1).domain([1,end]);

% We can only plot the coefficients of one CHEBFUN at a time:
if ( any(cellfun(@(f) isa(f, 'chebfun'), varargin)) )
    error('CHEBFUN:CHEBFUN:plotcoeffs:multipleChebfuns', ...
        'Calls of the form PLOTCOEFFS(F, ''b'',  G, ''r'') are not supported.');
end

% Store the hold state of the current axis:
holdState = ishold;

% Grab some colours:
col = [];
if ( nargin > 1 && ischar(varargin{1}) && numel(varargin{1}) < 4 )
    % (Note. LineStyle definitions have at most 4 characters. Other plotting
    % parameters have more than characters.)
    col = regexp(varargin{1}, '[bgrcmykw]', 'match');
    if ( numel(col) > 1 )
        error('CHEBFUN:plotcoeffs:color', ...
            'Error in color/linetype argument.');
    elseif ( ~isempty(col) )
        col = col{:};
%         varargin(1) = []; % Don't remove as we need marker information.
    end
end
if ( isempty(col) )
    colIdx = find(cellfun(@(arg) all(ischar(arg)) && ...
        any(strfind(arg, 'color')), varargin));
    if ( colIdx )
        col = varargin{colIdx+1};
        varargin(colIdx+(0:1)) = [];
    else
        col = get(0, 'DefaultAxesColorOrder');
        col = repmat(col, ceil(numColumns(f)/7), 1);
    end
end

% Deal with 'LogLog' and input:
doLogLog = any(cellfun(@(s) strcmpi(s, 'loglog'), varargin));

% Convert to a cell array for easy handling:
f = mat2cell(f);

% Initialise the output:
h = cell(1, numel(f));

% Loop over the columns in a quasimatrix:
for k = 1:numel(f)
    % Colours for this columns:
    if ( ischar(col) )
        colk = col;
    else
        colk = col(k,:);
    end
    % Call the column version:
    h{k} = columnPlotCoeffs(f{k}, colk, varargin{:}, 'domain', dom);
end

% Return hold state to what it was before:
if ( ~holdState )
    hold off
end

% Set xScale to logarithmic if requested:
if ( doLogLog )
    set(gca, 'xScale', 'Log');
end

% Give an output if one was requested:
if ( nargout > 0 )
    try
        h = cell2mat(h);
    catch
        % shrug
    end
    varargout = {h};
end

end

function h = columnPlotCoeffs(f, col, varargin)
numFuns = numel(f.funs);

% Initialise handle storage:
h = zeros(numFuns, 1);

% Call plotcoeffs at the tech level:
for j = 1:numFuns
    h(j) = plotcoeffs(f.funs{j}, varargin{:}, 'color', col); hold on
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function K = squared_levenstein(kern,dat1,dat2,ind1,ind2,kerParam)
%-------- convert matrix to cell array ----
D= [];
Xf = get_x(dat1);
Xs = get_x(dat2);
X1 = {};
X2 = {};

for i = 1:size(Xf,1)
  tmp = Xf(i,:);
  X1{i} = tmp(tmp > 0);
end

for i = 1:size(Xs,1)
  tmp = Xs(i,:);
  X2{i} = tmp(tmp > 0);
end

%----- get distance matrix -----------
for i = 1:length(ind1)
    disp(['D row ' num2str(i)])
    for j = 1:length(ind2)
        D(i,j) = d(X1{i},X2{j});
    end
end
%---get kernel matrix from distance matrix--
K = [];
m = size(D);

for i = 1:length(ind1)
    disp(['D->K row ' num2str(i)])
    for j = 1:length(ind2)
        si = length(X1{i});
        sj = length(X2{j});
        p = (si + sj + D(i,j))/2;
        K(i,j) = cos(2*atan(sqrt(((p-si)*(p-sj))/(p*(p-D(i,j))))));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = d(a,b)
    %----------initializing-------------
    C = zeros(length(a)+1,length(b)+1);
    C(:,1) = [0:length(a)]';
    C(1,:) = [0:length(b)];
    %------compute distance-------------
    for i = 2:size(C,1)
        for j = 2:size(C,2)
            delta = 1-abs(sign(a(i-1)-b(j-1)));
            C(i,j) = min([C(i-1,j)+1,C(i,j-1)+1,C(i-1,j-1) + 1 - delta]);
        end
    end

    s = C(end,end);


text:
function value = i4_rise ( x, n )

%*****************************************************************************80
%
%% I4_RISE computes the rising factorial function [X]^N.
%
%  Discussion:
%
%    [X]^N = X * ( X + 1 ) * ( X + 2 ) * ... * ( X + N - 1 ).
%
%    Note that the number of ways of arranging N objects in M ordered
%    boxes is [M]^N.  (Here, the ordering of the objects in each box matters).
%    Thus, 2 objects in 2 boxes have the following 6 possible arrangements:
%
%      -|12, 1|2, 12|-, -|21, 2|1, 21|-.
%
%    Moreover, the number of non-decreasing maps from a set of
%    N to a set of M ordered elements is [M]^N / N!.  Thus the set of
%    nondecreasing maps from (1,2,3) to (a,b,c,d) is the 20 elements:
%
%      aaa, abb, acc, add, aab, abc, acd, aac, abd, aad
%      bbb, bcc, bdd, bbc, bcd, bbd, ccc, cdd, ccd, ddd.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the argument of the rising factorial function.
%
%    Input, integer N, the order of the rising factorial function.
%    If N = 0, RISE = 1, if N = 1, RISE = X.  Note that if N is
%    negative, a "falling" factorial will be computed.
%
%    Output, integer VALUE, the value of the rising factorial function.
%
  value = 1;

  arg = x;

  if ( 0 < n )

    for i = 1 : n
      value = value * arg;
      arg = arg + 1;
    end

  elseif ( n < 0 )

    for i = -1 : -1 : n
      value = value * arg;
      arg = arg - 1;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [f] = spm_cost_SHC_fx(x,v,P)
% equations of motion for foraging problem using SHCs
% problem
% FORMAT [f] = spm_cost_SHC_fx(x,v,P)
%
% x   - hidden states (x.x, x.v x.q and x.a)
% v   - exogenous inputs
% P   - parameters
%
% The parameters associate increases in some physiological states x.q with
% positions in physical space, encoded by radial basis functions x.a
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_cost_SHC_fx.m 3757 2010-03-08 11:41:53Z guillaume $

% location and radius of attractors A
%--------------------------------------------------------------------------
global A; X   = A.x;

% gradient of Hamiltonian (G) is determined by the attractor state x.a
%--------------------------------------------------------------------------
[m,i] = max(x.a);
G     = x.x - X(:,i);

% motion of physical states
%--------------------------------------------------------------------------
f   = x;
f.x = x.v;
f.v = -G*8 - x.v*4;

% motion of physiological states (using basis functions of position)
%--------------------------------------------------------------------------
for i = 1:size(X,2)
    b(i,1) = norm(x.x - X(:,i)) < A.d;
end

f.q = P'*b - x.q/2;
f.a = P*(x.q < A.u) - b*4 - sum(x.a);

% flow
%--------------------------------------------------------------------------
dt  = 1/8;
f.x = f.x*dt;
f.v = f.v*dt;
f.q = f.q*dt;
f.a = f.a*dt;

author: spm


text:
classdef SPX_Modulator < handle

    methods(Static)

        function [ outputSequence ] = modulate_bits_with_signals( inputSequence, s1, s0)
            %MODULATEBINARYSEQUENCE Modulates a sequence of bits to create output
            % The size of signal for each bit
            [N, L] = size(s1);
            % If s0 is not defined, we will assume s0 to be all zeros
            if ~exist('s0','var')
              s0=zeros(size(s1));
            end
            if ~iscolumn(s1)
                error('s1 must be a column vector.');
            end
            if ~iscolumn(s0)
                error('s0 must be a column vector.');
            end
            if ~isvector(inputSequence)
                error('inputSequence must be a vector.');
            end
            % Number of bits
            E = length(inputSequence);
            % Combine the two signals into a matrix
            s = [s0 s1];
            % The row number for each output
            r  = 1:N;
            outputSequence  = s(r, inputSequence+1);
            outputSequence = reshape(outputSequence, N*E,1);
        end

        function [ transmittedSequence ] = modulate_bits_with_gaussian_noise(...
            transmittedBits, N, sigma1, sigma0)
            %MODULATE_BITS_WITH_GAUSSIAN_NOISE Generates Gaussian noise samples with different
            %variances for the transmitted bits

            %   Detailed explanation goes here

            % Number of transmitted bits
            B = length(transmittedBits);
            % We create space for transmitted sequence
            transmittedSequence = zeros(B * N , 1);
            for i=1:B
                pos = (i-1)*N;
                bit = transmittedBits(i);
                if bit
                    transmittedSequence(pos + (1:N))  = sigma1 * randn(N,1);
                else
                    transmittedSequence(pos + (1:N))  = sigma0 * randn(N,1);
                end
            end

        end

        function [ outputSequence ] = modulate_mary_symbols_with_signals( inputSequence, signals )
            %MODULATE_MARY_SYMBOLS_WITH_SIGNALS Modulates a sequence of M-ary symbols to create
            %an output signal
            %
            % Input:
            %   inputSequence: Sequence of M-ary symbols.
            %       We assume that the symbols are in [1,M]
            %   signals: signals to be used for each symbol
            %       Every signal has length N[as a column vector]
            %       and there are M such signals.

            % The size of signal and number of signals
            [N, M] = size(signals);
            % Number of symbols
            E = length(inputSequence);
            % The modulation process
            % First generate signal for each symbol
            % column by column
            outputSequence  = signals(1:N, inputSequence);
            % Then reshape the output to a column vector
            outputSequence = reshape(outputSequence, N*E,1);
        end

    end

end

author: indigits

repo: sparse-plex


text:
function n = order( val, base )
%Order of magnitude of number for specified base. Default base is 10.
%order(0.002) will return -3., order(1.3e6) will return 6.
%Author Ivar Smith

if nargin < 2
    base = 10;
end
n = floor(log(abs(val))./log(base));

author: Sable

repo: mcbench-benchmarks


text:
function rswVector = rotateVectorFromEciToRsw(eciVector, rVect, vVect)
    R = normVector(rVect);
    W = normVector(crossARH(rVect, vVect));
    S = normVector(crossARH(W,R));

    ECI2RSWRotMat = [R,S,W];

    rswVector = ECI2RSWRotMat \ eciVector;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Defining A Sensor Mask By Opposing Corners Example
%
% This example demonstrates how to define a sensor mask using the grid
% coordinates of two opposing corners of a rectangle, and then use this for
% the detection of the pressure field generated by an initial pressure
% distribution within a two-dimensional homogeneous propagation medium. It
% builds on the Homogeneous Propagation Medium and Using A Binary Sensor
% Mask examples.
%
% author: Bradley Treeby
% date: 21st August 2014
% last update: 21st August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

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

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = 128;           % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500; 	% [m/s]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [Pa]
disc_x_pos = 50;    % [grid points]
disc_y_pos = 50;    % [grid points]
disc_radius = 8;    % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 3; % [Pa]
disc_x_pos = 80;    % [grid points]
disc_y_pos = 60;    % [grid points]
disc_radius = 5;    % [grid points]
disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

source.p0 = disc_1 + disc_2;

% define the first rectangular sensor region by specifying the location of
% opposing corners
rect1_x_start = 25;
rect1_y_start = 31;
rect1_x_end   = 30;
rect1_y_end   = 50;

% define the second rectangular sensor region by specifying the location of
% opposing corners
rect2_x_start = 71;
rect2_y_start = 81;
rect2_x_end   = 80;
rect2_y_end   = 90;

% assign the list of opposing corners to the sensor mask
sensor.mask = [rect1_x_start, rect1_y_start, rect1_x_end, rect1_y_end;...
               rect2_x_start, rect2_y_start, rect2_x_end, rect2_y_end].';

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% query the size of the output data - note that when using a sensor mask
% defined using opposing corners, the output data is returned as a
% structure
size(sensor_data(1).p)
size(sensor_data(2).p)

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

% plot a time trace from the centre of each rectangle
figure;
plot(squeeze(sensor_data(1).p(end/2, end/2, :)), 'r-');
hold on;
plot(squeeze(sensor_data(2).p(end/2, end/2, :)), 'b-');
legend('Trace from centre of first rectangle', 'Trace from centre of second rectangle');
xlabel('Time Index');
ylabel('Pressure');
axis tight;

repo: Aero-matlab


text:
clear;
clc
close all;

%% EXAMPLE 2.3 Calculation of Gravity at various latitudes and heights
lat = deg2rad(45);
lon = 0;
hgt = 10000;

[~, ~, ~, ~, gravity]= ch_earth(lat, lon, hgt);
gravity

author: yandld


text:
clc;
close all
%% which one to show?
idx = 24;
coordIndices = [1,2,3];

imgFig1 = figure(1);
set(imgFig1, 'Position', [100 100 1400 900]) % [1 1 width height]

subplot(2,2,1);
imagesc(imgMat(:,:,:,idx)); axis off image;
subplot(2,2,2);
imagesc(instanceMaskMat(:,:,:,idx)); axis off image;
subplot(2,2,3);
A = (predInstanceMaskMat0(:,:,coordIndices,idx) + 1) / 2;
imagesc(A);  axis off image;
%% 3D surface
subplot(2,2,4);

r = 1;
[x,y,z] = sphere(50);
x0 = 0; y0 = 0; z0 = 0;
x = x*r + x0;
y = y*r + y0;
z = z*r + z0;

% figure
lightGrey = 0.7*[1 1 1]; % It looks better if the lines are lighter
surface(x,y,z, 'FaceColor', 'none', 'EdgeColor',lightGrey)
hold on
%% points
hold on;
points = reshape(predInstanceMaskMat0(:,:,coordIndices,idx), [], 3);
points = points';
pointsColor = points - (-1);%min(points(:)); % -1
pointsColor = pointsColor ./ 2;%max(pointsColor(:));
for i = 1:size(points,2)
    plot3( points(1,i), points(2,i), points(3,i), 's', 'MarkerSize',3, 'MarkerFaceColor', pointsColor(:,i)', 'MarkerEdgeColor', pointsColor(:,i)');
end
hold off;

axis off square
view([1 1 0.75]) % adjust the viewing angle
zoom(1.4)

%% save result
if flagSaveFig
    export_fig( sprintf('%s/%04d_visualization_single.jpg', saveFolder, i) );
end

author: aimerykong


text:
function plotData(x, y)
%PLOTDATA Plots the data points x and y into a new figure

plot(x, y, 'rx', 'MarkerSize', 10); % Plot the data
% Hint: You can use the 'rx' option with plot to have the markers
%       appear as red crosses. Furthermore, you can make the
%       markers larger by using plot(..., 'rx', 'MarkerSize', 10);

ylabel('Profit in $10,000s'); % Set the y ? axis label
xlabel('Population of City in 10,000s'); % Set the x ? axis label

figure; % open a new figure window

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function y = bwblkslv(L,b) %#ok
% BWBLKSLV Solves block sparse upper-triangular system.
%    y = bwblkslv(L,b) yields the same result as
%              y(L.perm,:) = L.L'\b
%    However, BWBLKSLV is faster than the built-in operator "\",
%    because it uses dense linear algebra and loop-unrolling on
%    supernodes.
%
%    Typical use, with X sparse m x m positive definite and b is m x n:
%            L = sparchol(symbchol(X),X);
%            y = bwblkslv(L,fwblkslv(L,b));
%    Then y solves X*y=b.
%
% See also symbchol, fwblkslv, mldivide, mrdivide

% This file is part of CholTool 1.00
% Copyright (C) 1998 Jos F. Sturm
% CRL, McMaster University, Canada.
% Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
i_fish=8;
% also load fish 8 in GUI
% change alpha_max to 1 in 'DrawCellsOnAnatProj.mat'

CR = getappdata(hfig,'CellResp');

%% plot basic stats of time series
TS = CR(:,1:3000);
TS_max = max(TS,[],2);
TS_avr = mean(TS,2);
figure;
hist(TS_max,0:0.1:5)
figure
hist(TS_avr,0:0.01:1)

%% color-code based on average activity
cIX_plot = 1:size(CR,1);
TS_avr = mean(TS,2);

TSmax = 0.3;
TSmin = min(TS_avr);
TS_avr(TS_avr>TSmax) = TSmin;

score = (TS_avr-TSmin)/(TSmax-TSmin);
gIX_plot = ceil(score*64);
gIX_plot(gIX_plot==0)=1;
wIX = score;%.^2; % set transparency weight, and enhance contrast

numK = max(gIX_plot);
clrmap = parula(numK);

setappdata(hfig,'isRefAnat',0);
I = LoadCurrentFishForAnatPlot(hfig,cIX_plot,gIX_plot,clrmap,wIX);
DrawCellsOnAnat(I);

DrawCustomColorbar(clrmap,[TSmin,TSmax],2,gcf);
%%
savedir = GetOutputDataDir;
filename = fullfile(savedir,'Zstack.tif');
WriteZstack(hfig,filename,cIX_plot,gIX_plot,clrmap);

author: xiuyechen

repo: FishExplorer


text:
function results = vl_test_slic(varargin)
% VL_TEST_SLIC
vl_test_init ;

function s = setup()
s.im = im2single(vl_impattern('roofs1')) ;

function test_slic(s)
segmentation = vl_slic(s.im, 10, 0.1) ;

author: jianxiongxiao

repo: ProfXkit


text:
% STD_READLIMOFILE- Read limo file into eeglab structure
%
% Usage:
%   >> [MeasureData, parameters, MeasureRange1, MeasureRange2] = std_readfilelimo(filestruct,'key', val)
%
% Inputs:
% filestrcuct  - Structure with the following fields:
%                datasets : Index of he datasets in datasetinfo (real)
%                trials   : index of the trials to pull out     (cell array)
%                value    : Cell array with the values of the independent variables of the set.
%                case     : name of the subject /case
%                filebase : Full Path and filename of the
%
% Optional inputs:
%   'measure'     - ['itcbeta1' | 'itcbeta2' |'itcr2r' |'itcr2f' |'itcr2p' |
%                    'erpbeta1' | 'erpbeta2' |'erpr2r' |'erpr2f' |'erpr2p' |
%                    'erspbeta1'| 'erspbeta2'|'erspr2r'|'erspr2f'|'erspr2p'|
%                    'specbeta1'| 'specbeta2'|'specr2r'|'specr2f'|'specr2p ]
%   'timelimits'  - [min max] ERSP time (latency in ms) range of interest
%   'freqlimits'  - [min max] ERSP frequency range (in Hz) of interest
%   'channels'    - [cell or integer] channel labels - for instance
%                   { 'cz' 'pz' }
%
% Outputs:
%   MeasureData    - the multi-channel or multi-component data. The size of this
%                    output depends on the number of dimension (for ERSP or ERP
%                    for instance). The last dimension is channels or components.
%   parameters     - structure containing parameters saved with the data file.
%   MeasureRange1  - time points (ERP, ERSP) or frequency points (spectrum)
%   MeasureRange2  - frequency points (ERSP, ITCs)
%
%
%
% See also:
%
%
% Authors: Arnaud Delorme, Ramon Martinez-Cancino
%
% Copyright (C) 2014  Ramon Martinez-Cancino,INC, SCCN
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
function [MeasureData, parameters, MeasureRange1, MeasureRange2] = std_readfilelimo(filestruct, varargin)

parameters = [];
try
    options = varargin;
    if ~isempty( varargin ),
        for i = 1:2:numel(options)
            g.(options{i}) = options{i+1};
        end
    else g = []; end
catch
    error('std_readfilelimo() error: calling convention {''key'', value, ... } error');
end

try, g.channels;                catch, g.channels          = [];        end; %
try, g.components;              catch, g.components        = [];        end; %
try, g.measure;                 catch, g.measure           = [];        end; %
try, g.timelimits;              catch, g.timelimits        = [];        end; %
try, g.triallimits;             catch, g.triallimits       = [];        end; %
try, g.freqlimits;              catch, g.freqlimits        = [];        end; %
try, g.getparamonly;            catch, g.getparamonly      = 'off';     end; %  work on this

%  Checking entries
%..........................................................................

% Indetify if components or channels
%..........................................................................
if isempty(g.components)
    chanoric = 'Channels';
else
    chanoric = 'Components';
end

% Identify type of analysis for limo stuff path
%..........................................................................
if strcmp(g.measure(1:3),'erp') && strcmp(chanoric,'Channels')
    type = 'Channels_Time';
elseif strcmp(g.measure(1:3),'erp') && strcmp(chanoric,'Components')
    type = 'Components_Time';
elseif strcmp(g.measure(1:3),'spe') && strcmp(chanoric,'Channels')
    type = 'Channels_Frequency';
elseif strcmp(g.measure(1:3),'spe') && strcmp(chanoric,'Components')
    type = 'Components_Frequency';
elseif strcmp(g.measure(1:3),'ers') && strcmp(chanoric,'Channels')
    type = 'Channels_Time-Frequency';
elseif strcmp(g.measure(1:3),'ers') && strcmp(chanoric,'Components')
    type = 'Components_Time-Frequency';
end

% Study pathname
%..........................................................................
filebasetemp = filestruct.filebase;
[path, name, tmp] = fileparts(filebasetemp);
[path2study,tmp,tmp]= fileparts(path);

% Getting design index in STUDY
%..........................................................................
designindx = name(7:9);
if isempty(str2num(designindx))
    designindx = designindx(1:2);
    if isempty(str2num(designindx))
        designindx = designindx(1);
        if isempty(str2num(designindx))
            error('std_readfilelimo(): Fail to detect index of the design in STUDY');
        end
    end
end

% Creating path to LIMO.mat and loading it
%..........................................................................
filelist        = dir(path2study);
fileindxtmp     = strmatch('LIMO_', {filelist.name});
limofoldername  = filelist(fileindxtmp(1)).name;
path2limofiles  = [path2study filesep limofoldername filesep filestruct.case filesep 'GLM' num2str(designindx) '_' type];
structmp        = load ([path2limofiles filesep 'LIMO.mat']);
limostruct      = structmp.LIMO;

% Identifying condition
%..........................................................................
dims = []; cond_indxs = [];
condtemp = {limostruct.data.studyinfo.value};

for i =  1 : length(condtemp)
    dims(i)       = length(condtemp{i});
    tmp = condtemp{i};
    for j = 1: dims(i)
        tmp{j} = num2str(tmp{j});
    end
    condtemp{i} = tmp;
    cond_indxs(i) = find(strcmp(condtemp{i},num2str(filestruct.value{i})));
end
% Sub2ind stuff
command = 'sub2ind(dims';
for i = 1: length(dims)
    command = [command ',cond_indxs(' num2str(i) ')'];
end
command = [command ')'];
CondIndx = eval(command);

% Checking time and set TimeIndxVec if type = TIME
%..........................................................................
if all([~isempty(g.timelimits),isfield(limostruct.data,'sampling_rate')])
    if any([g.timelimits(1) < limostruct.data.start, g.timelimits(2) > limostruct.data.end])
        error('std_readfilelimo(): Time indices provided are out of bonds');
    else
        TimeVec = limostruct.data.start: 1000/limostruct.data.sampling_rate:limostruct.data.end
        [tmp,TimeIndx(1)] = min(abs(TimeVec - g.timelimits(1)));
        [tmp,TimeIndx(2)] = min(abs(TimeVec - g.timelimits(2)));
        TimeIndxVec = TimeIndx(1):TimeIndx(2);
    end
elseif all([isempty(g.timelimits),isfield(limostruct.data,'sampling_rate')]) % will use all the timepoints
    TimeIndxVec = 1:length(limostruct.data.start: 1000/limostruct.data.sampling_rate:limostruct.data.end);
end

% Checking frequencies and set FreqIndxVec if type = Frequency
%..........................................................................
 if all([~isempty(g.freqlimits),isfield(limostruct.data,'freqlist')])
     if any([g.freqlimits(1) < limostruct.data.lowf, g.freqlimits(2) > limostruct.data.highf])
         error('std_readfilelimo(): Frequency indices provided are out of bonds');
     else
         [tmp,FreqIndx(1)] = min(abs(limostruct.data.freqlist - g.freqlimits(1)));
         [tmp,FreqIndx(2)] = min(abs(limostruct.data.freqlist - g.freqlimits(2)));
         FreqIndxVec = FreqIndx(1):FreqIndx(2);
     end
 elseif all([isempty(g.freqlimits),isfield(limostruct.data,'freqlist')]) % will use all frequencies
     FreqIndxVec = 1 : length(limostruct.data.freqlist);
 end

 % Checking frequencies and time and set FreqIndxVec and TimeIndxVec if type = Time-Frequency
 %..........................................................................
 if all([isfield(limostruct.data,'tf_times'), isfield(limostruct.data,'tf_freq')])

     %1- Time
     if ~isempty(g.timelimits)
         if any([g.timelimits(1) < limostruct.data.start, g.timelimits(2) > limostruct.data.end])
             error('std_readfilelimo(): Time indices provided are out of bonds');
         else
             [tmp,TimeIndx(1)] = min(abs(limostruct.data.tf_times - g.timelimits(1)));
             [tmp,TimeIndx(2)] = min(abs(limostruct.data.tf_times - g.timelimits(2)));
             TimeIndxVec = TimeIndx(1):TimeIndx(2);
         end
     else % will use all timepoints
         TimeIndxVec = 1:length(limostruct.data.tf_times);
     end

     %2- Frequency
     if ~isempty(g.freqlimits)
         if any([g.freqlimits(1) < limostruct.data.lowf, g.freqlimits(2) > limostruct.data.highf])
             error('std_readfilelimo(): Frequency indices provided are out of bonds');
         else
             [tmp,FreqIndx(1)] = min(abs(limostruct.data.tf_freqs - g.freqlimits(1)));
             [tmp,FreqIndx(2)] = min(abs(limostruct.data.tf_freqs - g.freqlimits(2)));
             FreqIndxVec = FreqIndx(1):FreqIndx(2);
         end
     else % will use all frequencies
         FreqIndxVec = 1:length(limostruct.data.tf_freqs);
     end

 end

% Checking channels or ic to load
%..........................................................................
% 1- Channels
if isempty(g.channels) && strcmp(chanoric, 'Channels')
    ChanIcIndx = 1:length({limostruct.data.chanlocs.labels});
elseif ~isempty(g.channels) && strcmp(chanoric, 'Channels')
     for i = 1 : length(g.channels)
         if isempty(find(strcmp(lower(g.channels(i)),lower({limostruct.data.chanlocs.labels}))))
             error('std_readfilelimo(): Invalid channel label');
         else
             ChanIcIndx = find(strcmp(lower(g.channels(i)),lower({limostruct.data.chanlocs.labels})));
         end
     end
end

% 2- Component
if isempty(g.components) && strcmp(chanoric, 'Components')
    ChanIcIndx = 1:length({limostruct.data.chanlocs.labels});
elseif ~isempty(g.channels) && strcmp(chanoric, 'Components')
     for i = 1 : length(g.channels)
         if isempty(find(strcmp(g.channels(i),{limostruct.data.chanlocs.labels})))
             error('std_readfilelimo(): Invalid channel label');
         else
             ChanIcIndx = find(strcmp(g.channels(i),{limostruct.data.chanlocs.labels}));
         end
     end
end

% Retrieving data
%..........................................................................
if isfield(limostruct.data,'sampling_rate')
    TimeVec = limostruct.data.start: 1000/limostruct.data.sampling_rate:limostruct.data.end
end

switch g.measure
    % beta1
    %..................................................................
    case 'erpbeta1'
        datatmp = load([path2limofiles filesep 'Betas.mat']); % PUT THIS OUT OF THE LOOP
        data = datatmp.Betas;
        MeasureData   = data(ChanIcIndx,TimeIndxVec,CondIndx);   % Chan/IC x Time x Condition

    case 'specbeta1'
        datatmp = load([path2limofiles filesep 'Betas.mat']);
        data = datatmp.Betas;
        MeasureData   = data(ChanIcIndx,FreqIndxVec,CondIndx);   % Chan/IC x Freq x Condition

    case {'erspbeta1','itcbeta1'}
        datatmp = load([path2limofiles filesep 'Betas.mat']);
        data = datatmp.Betas;
        MeasureData   = data(FreqIndxVec, TimeIndxVec,CondIndx); % Freq x Time x Condition

    % beta2
    %..................................................................
    case {'erpbeta2','erspbeta2','specbeta2','itcbeta2'}
        error('std_readfilelimo(): Invalid option in this realease');

    % erpr2
    %..................................................................
    case 'erpr2r'
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(ChanIcIndx,TimeIndxVec,1);

     case 'erpr2f'
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(ChanIcIndx,TimeIndxVec,2);

    case 'erpr2p'
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(ChanIcIndx,TimeIndxVec,3);

    % erspr2 & itcr2
    %..................................................................
    case {'erspr2r','itcr2r'}
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(FreqIndxVec,TimeIndxVec,1);

     case {'erspr2f','itcr2f'}
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(FreqIndxVec,TimeIndxVec,2);

    case {'erspr2p','itcr2p'}
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(FreqIndxVec,TimeIndxVec,3);

    %spec
    %..................................................................
    case {'specr2r'}
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(ChanIcIndx,FreqIndxVec,1);

     case {'specr2f'}
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(ChanIcIndx,FreqIndxVec,2);

    case {'specr2p'}
        datatmp = load([path2limofiles filesep 'R2.mat']);
        data = datatmp.R2;
        MeasureData   = data(ChanIcIndx,FreqIndxVec,3);
end

% Measures ranges
%..........................................................................
switch g.measure(1:3)
    case 'erp'
        MeasureRange1 = TimeVec(TimeIndxVec);
        MeasureRange2 = [];
    case {'ers','itc'}
        MeasureRange1 = limostruct.data.tf_times(TimeIndxVec);
        MeasureRange2 = limostruct.data.tf_freqs(FreqIndxVec);
    case 'spe'
        MeasureRange1 = limostruct.data.freqlist(FreqIndxVec);
        MeasureRange2 = [];
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [sol, infos] = gsp_ml_rls_tv(G, xl, y, k, tau,lambda, A, At, param)
%GSP_ML_RLS_TV Manifold Learning regularized least square with TV regularization
%   Usage: sol = gsp_ml_rls_tv(G, xl, y, k, tau,lambda);
%          sol = gsp_ml_rls_tv(G, xl, y, k, tau, lambda, A, At);
%          sol = gsp_ml_rls_tv(G, xl, y, k, tau, lambda, A, At, param);
%          [sol, infos] = gsp_ml_rls_tv(...)
%
%   Input parameters:
%       G       : Graph
%       xl      : labeled points
%       y       : labels
%       k       : kernel
%       tau     : regularization parameters
%       lambda  : regularization parameters
%       A       : Operator
%       At      : Adoint operator
%       param   : Optional parameters
%   Output parameters:
%       sol     : solution of the problem (kernel coefficients)
%       infos   : convergence info
%
%   *param* is a structure of optional argument given to the solver
%   gradient_descent. Please see the function gradient descent for more
%   information.
%
%   In *param*, you also have to set an upperbound for the operator A as
%   param.nu!
%
%   This function solves the following problem:
%
%   ..  argmin_alpha  || A (K alpha) - y ||_2^2
%                       + tau *alpha^T K alpha
%                       + lambda || L K alpha ||_TVG
%
%   If tau is set to zero, then the following problem is solved
%
%   ..  argmin_alpha alpha^T K alpha
%                      + lambda || L K alpha ||_TVG
%                      s. t.  A (K alpha) = y
%

% Author: Nathanael Perraudin
% Date  : 8 decembre 2014

if nargin<7
    A = @(x) x;
end

if nargin<8
    At = A;
end

if nargin<9
    param = struct;
end

if ~isfield(param, 'tol'), param.tol = 1e-6; end
if ~isfield(param, 'nu'), param.nu = 1; end
if ~isfield(param, 'verbose'), param.verbose = 1; end

N = size(xl,2);

% Evaluate the kernel on the data points
K = gsp_rkhs_evaluate(k,xl);
nu = norm(K);

alpha_in = zeros(N,size(y,2));

if ~isfield(G,'D');
    G = gsp_adj2vec(G);
end
if ~isfield(G,'lmax');
    G = gsp_estimate_lmax(G);
end

paramtv.verbose = param.verbose - 1;
ftv.eval = @(x) sum(lambda*gsp_norm_tv(G,K*x));
ftv.prox = @(x,T) gsp_prox_tv(x,lambda*T,G,paramtv);

if tau >0

    fp.eval = @(x) tau * sum(gsp_norm_tik(K,x));
    fp.grad = @(x) 2*tau*K*x;

    ffid.eval = @(x) norm(A(K*x)-y,'fro')^2;
    ffid.grad = @(x) 2*K'*At(A(K*x)-y);

    ftot.grad = @(x) ffid.grad(x) + fp.grad(x);
    ftot.eval = @(x) ffid.eval(x) + fp.eval(x);

    param.gamma = 0.5/(tau*nu+nu^2*param.nu^2+lambda*nu*G.lmax);

    [sol,infos] = forward_backward(alpha_in, ftv, ftot, param);
else

    fp.eval = @(x)  x'*K*x;
    fp.grad = @(x) 2*K*x;

    paramproj.A = @(x) A(K*x);
    paramproj.At = @(x) K'*At(x);
    paramproj.nu = nu^2*param.nu^2;
    paramproj.tight = 0;
    paramproj.verbose = param.verbose-1;
    paramproj.y = y;
    paramproj.maxit = 50;
    ffid.eval = @(x) eps;
    ffid.prox = @(x,T) proj_b2(x,T,paramproj);

    param.gamma = 0.5/(nu+lambda*nu*G.lmax);

    [sol,infos] = generalized_forward_backward(alpha_in, {ffid,ftv}, fp, param);

end

end

repo: gspbox


text:
function target = myConvolve(kConv, data, kernel, stride, task)
% Matlab interfaces of alex fast convolution kernels.
% There are are total 3 kinds of operations : given the two of (images,
% filters, hidacts) to compute the remaining one. And for each operation,
% there are two kernels : one for the second layer which is optimized
% for one channel input and one for the other layers.

% forward:  (data,kernel) == (images, filters)
% weight:   (data,kernel) == (images, hidacts)
% backward: (data,kernel) == (hidacts, filters)

% interfaces for cuda code kFunctions.cu.

if strcmp(task,'forward')
    numColors = size(data,5);
    if numColors == 1
        numFilters = size(kernel,1); filter_batch_size = 16;
        filter_batch_num = numFilters / filter_batch_size;
        assert(filter_batch_num == floor(filter_batch_num));

        numImages = size(data,1); imgSizeX = size(data,2); imgSizeY = size(data,3); imgSizeZ = size(data,4);
        filterSize =  size(kernel,2);

        paddingStart = 0; moduleStride = stride; imgStride = numImages;
        numModulesX = (imgSizeX - filterSize) / stride + 1; numModulesY = (imgSizeY - filterSize) / stride + 1; numModulesZ = (imgSizeZ - filterSize) / stride + 1;
        target = zeros(numImages, numModulesX, numModulesY, numModulesZ, numFilters, 'single');
        for b = 1 : filter_batch_num
            filterPerThread = 4; imagePerThread = 1;
            kConv.ThreadBlockSize = [32, 4];
            kConv.GridSize = [ceil(numImages/(32 * imagePerThread)), numModulesX * numModulesY * numModulesZ * filter_batch_size / (filterPerThread * 4) ]; % filterPerThread, BlockSize.Y
            target_gpu = feval(kConv,...
               target(:,:,:,:,(b-1)*filter_batch_size+1 : b*filter_batch_size), data, kernel((b-1)*filter_batch_size+1 : b*filter_batch_size,:,:,:,:),...
               numImages, filter_batch_size, imgSizeZ, imgSizeY, imgSizeX, filterSize, ...
               paddingStart, moduleStride, numModulesZ, numModulesY, numModulesX, imgStride);
           target(:,:,:,:,(b-1)*filter_batch_size+1 : b*filter_batch_size) = gather(target_gpu);
        end
    else
        numImages = size(data,1); imgSizeX = size(data,2); imgSizeY = size(data,3); imgSizeZ = size(data,4);
        numFilters = size(kernel,1); filterSize =  size(kernel,2);

        paddingStart = 0; moduleStride = stride; imgStride = numImages; numGroups = 1;
        numModulesX = (imgSizeX - filterSize) / stride + 1; numModulesY = (imgSizeY - filterSize) / stride + 1; numModulesZ = (imgSizeZ - filterSize) / stride + 1;

        filterPerThread = 8; imagePerThread = 1;

        kConv.ThreadBlockSize = [32, 4];
        kConv.GridSize = [ceil(numImages/(32 * imagePerThread)), numModulesX * numModulesY * numModulesZ * numFilters / (filterPerThread * 4) ]; % filterPerThread, BlockSize.Y

        target = zeros(numImages, numModulesX, numModulesY, numModulesZ, numFilters, 'single');

        target_gpu = feval(kConv,...
           target, data, kernel,...
           numImages, numFilters, imgSizeZ, imgSizeY, imgSizeX, filterSize, ...
           paddingStart, moduleStride, numModulesZ, numModulesY, numModulesX, imgStride, numColors, numGroups);

        target = gather(target_gpu);
    end

elseif strcmp(task,'weight')
    numColors = size(data,5);
    if numColors == 1
        numImages = size(data,1); imgSizeX = size(data,2); imgSizeY = size(data,3); imgSizeZ = size(data,4);
        numModulesX = size(kernel,2); numModulesY = size(kernel,3); numModulesZ = size(kernel,4); numFilters = size(kernel,5);
        paddingStart = 0; moduleStride = stride; imgStride = numImages; partialSum = numModulesX * numModulesY * numModulesZ;
        filterSize = imgSizeX - stride * (numModulesX - 1);

        pixelsPerThread = 5; preLoadCases = 32; scaleOutput = 1 ./ (numImages * partialSum);

        kConv.ThreadBlockSize = [16, 8];
        kConv.GridSize = [numFilters*numModulesX*numModulesY*numModulesZ/partialSum/16, ceil(filterSize^3 /(8*pixelsPerThread))];

        target = zeros(numFilters, filterSize, filterSize, filterSize, numColors, 'single');
        target_gpu = feval(kConv, ....
           target, data, kernel,...
           numImages, numFilters, numModulesZ, numModulesY, numModulesX, imgSizeZ, imgSizeY, imgSizeX, ...
           filterSize, paddingStart, moduleStride, imgStride, partialSum, scaleOutput);

        target = gather(target_gpu);
    else
        numImages = size(data,1); imgSizeX = size(data,2); imgSizeY = size(data,3); imgSizeZ = size(data,4);
        numFilters = size(kernel,5); numModulesX = size(kernel,2); numModulesY = size(kernel,3); numModulesZ = size(kernel,4);
        paddingStart = 0; moduleStride = stride; imgStride = numImages; partialSum = numModulesX * numModulesY * numModulesZ;
        filterSize = imgSizeX - stride * (numModulesX - 1);

        preLoadCases = 32; filtersPerThread = 2; colorsPerThread = 8;
        scaleOutput = 1 ./ (numImages * partialSum); numGroups = 1;

        kConv.ThreadBlockSize = [16, 8];
        kConv.GridSize = [numFilters*numModulesX*numModulesY*numModulesZ/partialSum/16/filtersPerThread, ceil(filterSize^3 / 8) * (numColors / colorsPerThread)];

        target = zeros(numFilters, filterSize, filterSize, filterSize, numColors, 'single');
        target_gpu = feval(kConv, ....
           target, data, kernel,...
           numImages, numFilters, numModulesZ, numModulesY, numModulesX, imgSizeZ, imgSizeY, imgSizeX, ...
           filterSize, paddingStart, moduleStride, imgStride, numColors, numGroups, partialSum, scaleOutput);

        target = gather(target_gpu);
    end
elseif strcmp(task,'backward')
    numColors = size(kernel,5);
    if numColors == 1
        numImages = size(data,1); numModulesX = size(data,2); numModulesY = size(data,3); numModulesZ = size(data,4); moduleStride = stride;
        numFilters = size(kernel,1); filterSize = size(kernel,2);

        imgSizeX = stride * (numModulesX - 1) + filterSize; imgSizeY = stride * (numModulesY - 1) + filterSize; imgSizeZ = stride * (numModulesZ - 1) + filterSize;
        paddingStart = 0; imgsPerThread = 2;

        kConv.ThreadBlockSize = [16, 16];
        kConv.GridSize = [ceil(numImages/(imgsPerThread *16)), imgSizeZ * ceil(imgSizeY/4) * ceil(imgSizeX/4)];

        target = zeros(numImages, imgSizeX, imgSizeY, imgSizeZ, numColors, 'single');
        target_gpu = feval(kConv,....
            target, data, kernel,...
            numModulesZ, numModulesY, numModulesX, numImages, numFilters, filterSize, ...
            imgSizeZ, imgSizeY, imgSizeX, paddingStart, moduleStride);

        target = gather(target_gpu);
    else
        numImages = size(data,1); numModulesX = size(data,2); numModulesY = size(data,3); numModulesZ = size(data,4); moduleStride = stride;
        numFilters = size(kernel,1); filterSize = size(kernel,2);

        imgSizeX = stride * (numModulesX - 1) + filterSize; imgSizeY = stride * (numModulesY - 1) + filterSize; imgSizeZ = stride * (numModulesZ - 1) + filterSize;
        paddingStart = 0; numGroups = 1;

        colorsPerThread = 4; imgsPerThread = 1;

        kConv.ThreadBlockSize = [32, 4];
        kConv.GridSize = [ceil(numImages/(imgsPerThread * 32)) * (numColors / (4 * colorsPerThread)), imgSizeZ * imgSizeY * imgSizeX];

        target = zeros(numImages, imgSizeX, imgSizeY, imgSizeZ, numColors, 'single');
        target_gpu = feval(kConv,....
            target, data, kernel,...
            numModulesZ, numModulesY, numModulesX, numImages, numFilters, filterSize, ...
            imgSizeZ, imgSizeY, imgSizeX, paddingStart, moduleStride, numColors, numGroups);

        target = gather(target_gpu);
    end
end

author: zhirongw


text:
%% ARDRONE TEST

% Author: James A. Douthwaite

classdef ARdrone_LQR < ARdrone_prev
%% INITIALISE THE AGENT SPECIFIC PARAMETERS
    % LQR SPECIFIC PARAMETERS
    properties

    end
    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function this = ARdrone_LQR(varargin)

            % Call the super class
            this@ARdrone_prev(varargin);                                         % Create the super class 'agent'

            % DYNAMIC CONSTRAINTS
            this.v_nominal = 2.0;
            this.v_max = 3;

            % Append control parameter
            %obj.DYNAMICS.Q = diag([1 1 1 1 1 1 2 1 1 1 1 1])*1.5E1;
            this.DYNAMICS.Q = diag([1 1 1 1 1 1 1 1 1 1 1 1])*1.5E1;
            this.DYNAMICS.R = diag(ones(4,1))*1E-4;
            this.DYNAMICS.N = zeros(size(this.DYNAMICS.SS.B));        % Terminal input penalisation

            % Calculate the linear feedback from the LQR
            this.DYNAMICS.K_lqr = this.CreateController_LQR(...
                this.DYNAMICS.SS,...
                this.DYNAMICS.Q,...
                this.DYNAMICS.R,...
                this.DYNAMICS.N);

            % //////////////// Check for user overrides ///////////////////
            % - It is assumed that overrides to the properties are provided
            %   via the varargin structure.
            [this] = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Setup
        % - The same as any 2D/3D object
        % Main
        function this = main(this,ENV,varargin)
            % INPUTS:
            % varargin - Cell array of inputs
            % >dt      - The timestep
            % >objects - The detectable objects cell array of structures
            % OUTPUTS:
            % obj      - The updated project

            % //////////// CHECK FOR NEW INFORMATION UPDATE ///////////////
            % UPDATE THE AGENT WITH THE NEW ENVIRONMENTAL INFORMATION
            this = this.GetAgentUpdate(ENV,varargin{1});

            % /////////////////// WAYPOINT TRACKING ///////////////////////
            % Get desired heading
            desiredHeadingVector = this.GetTargetHeading();
            desiredVelocity = desiredHeadingVector*this.v_nominal;

            % Express velocity vector in NED coordinates
            if ENV.currentTime <= 1
                desiredVelocity = zeros(3,1);
            end

            p_wp = this.targetWaypoint.position(:,this.targetWaypoint.sampleNum);

            % Convert from 'xyz' to 'ned'
%             desiredVelocity = enu2ned(desiredVelocity);
%             desiredNEDstate = [p_wp;zeros(3,1);zeros(6,1)];

            % LQR controller update
            this = this.Controller_position(ENV,enu2ned(p_wp));

%             this = this.Controller_velocity(ENV,enu2ned(desiredVelocity));

            % \\\\\\\\\\\\\\\ RECORD THE AGENT-SIDE DATA \\\\\\\\\\\\\\\\\\
            this.DATA.inputNames = {'$\dot{x}$ (m/s)','$\dot{y}$ (m/s)','$\dot{z}$ (m/s)',...
                                   '$\dot{\phi}$ (rad/s)','$\dot{\theta}$ (rad/s)','$\dot{\psi}$ (rad/s)'};
            this.DATA.inputs(1:numel(this.DATA.inputNames),ENV.currentStep) = this.localState(7:end);          % Record the control inputs
        end
    end

    %% /////////////////////// AUXILLARY METHODS //////////////////////////
    % CONTROLLER METHODS
    methods
        % ARdrone controller (local position)
        function [this] = Controller_position(this,ENV,NEDposition)

            % CALCULATE THE NEW STATE REFERENCE
            if ~this.IsIdle()
                X_desired = [NEDposition;zeros(2,1);-this.localState(6);zeros(6,1)];
            else
                % Idle reference
                X_desired = [zeros(5,1);this.localState(6);zeros(6,1)];
            end

            % /////////////////// AIRCRAFT DYNAMICS ///////////////////////
            useLinearModel = false;
            if ~useLinearModel
                [this.localState] = this.UpdateNonLinearPlant(ENV,this.localState,X_desired);
            else
                [this.localState] = this.UpdateLinearPlant(ENV,this.localState,X_desired);
            end

            % \\\\\\\\\\\\\\\\\\\\\ GLOBAL UPDATE \\\\\\\\\\\\\\\\\\\\\\\\\
            this = this.GlobalUpdate_ARdrone(ENV.dt,this.localState);
        end
        % ARdrone controller (local velocity)
        function [this] = Controller_velocity(this,ENV,NEDVelocity)

            % CALCULATE THE NEW STATE REFERENCE
            if ~this.IsIdle()
                X_desired = [zeros(5,1);this.localState(6);NEDVelocity;zeros(3,1)];
            else
                % Idle reference
                X_desired = [zeros(5,1);this.localState(6);zeros(6,1)];
            end

            % /////////////////// AIRCRAFT DYNAMICS ///////////////////////
            useLinearModel = false;
            if ~useLinearModel
                [this.localState] = this.UpdateNonLinearPlant(ENV,this.localState,X_desired);
            else
                [this.localState] = this.UpdateLinearPlant(ENV,this.localState,X_desired);
            end

            % \\\\\\\\\\\\\\\\\\\\\ GLOBAL UPDATE \\\\\\\\\\\\\\\\\\\\\\\\\
            this = this.GlobalUpdate_ARdrone(ENV.dt,this.localState);
        end
    end
    methods (Static)
        % Create an LQR controller (feedback gain)
        function [K_lqr] = CreateController_LQR(SS,Q,R,N)
            % This function designs the LQR controller used to provide
            % error feedback.

            if nargin < 6
                N = zeros(size(SS.B));
            end

            % Get the LQR feedback
            [K_lqr,~,~] = lqr(SS,Q,R,N);

            % IGNORE POSITION FEEDBACK
            K_lqr(:,1:3) = 0;   % Omit position feedback
            % IGNORE XY POSITION FEEDBACK
%             K_lqr(:,1:2) = 0; % Omit XY position feedback
            % IGNORE HEADING FEEDBACK
%             K_lqr(:,6) = 0;
        end
    end
    % //////////////////////// MODELLING/DYNAMICS /////////////////////////
    methods
        % ODE45 - UPDATE NONLINEAR CLOSED LOOP DYNAMICS
        function [X] = UpdateNonLinearPlant(this,TIME,X0,X_desired)
            % This function computes the state update for the current agent
            % using the ode45 function.

            % Input sanity check
            assert(isstruct(TIME),'Expecting a time structure.');

            X = X0; % Default to no change

            % DETERMINE THE INTEGRATION PERIOD
            if TIME.currentTime ~= TIME.timeVector(end)
                % INTEGRATE THE DYNAMICS OVER THE TIME STEP
                [~,Xset] = ode45(@(t,X) this.ARdrone_nonLinear_closedLoop(X,X_desired),...
                    [0 TIME.dt],X0,odeset('RelTol',1e-2,'AbsTol',TIME.dt*1E-2));
                X = Xset(end,:)';
            end
        end
        % ODE45 - UPDATE LINEAR CLOSED LOOP DYNAMICS
        function [X] = updateLinearPlant(this,TIME,X0,X_desired)
            % This function computes the state update for the current agent
            % using the ode45 function.

            % Input sanity check
            assert(isstruct(TIME),'Expecting a time structure.');

            X = X0; % Default to no change

            % Check for the last time-step
            if TIME.currentTime ~= TIME.timeVector(end)
                % INTEGRATE THE DYNAMICS OVER THE TIME STEP
                [~,Xset] = ode45(@(t,X) this.ARdrone_linear_closedLoop(X,X_desired),...
                    [0 TIME.dt],X0,odeset('RelTol',1e-2,'AbsTol',TIME.dt*1E-2));
                X = Xset(end,:)';
            end
        end
    end
    methods
        % CLOSED-LOOP NONLINEAR DYNAMICS
        function [dX] = ARdrone_nonLinear_closedLoop(this,X,X_desired)
            % THE STATE ERROR
            Xerror = X - X_desired;
            % THE NOMINAL INPUT
            [Uss] = this.ARdrone_nominalInput(X);
            % GET THE NOMINAL INPUT + LQR FEEDBACK
            U = Uss - this.DYNAMICS.K_lqr*Xerror;
            % CALL THE OPENLOOP DYNAMICS
            [dX] = this.ARdrone_nonLinear_openLoop(X,U);
        end
        % CLOSED-LOOP LINEAR DYNAMICS
        function [dX] = ARdrone_linear_closedLoop(this,X,X_desired)
            % THE STATE ERROR
            Xerror = X - X_desired;
            % GET THE NOMINAL INPUT + LQR FEEDBACK
            dU = - this.DYNAMICS.K_lqr*Xerror;
            % CALL THE OPENLOOP DYNAMICS
            [dX] = this.ARdrone_linear_openLoop(X,dU);
        end
    end
end
% AGENT STATE VECTOR [x;y;z;psi;the;phi;v;u;w;p;q;r]

repo: OpenMAS


text:
function [vid_ycbcr] = convert_rgb2ycbcr(vid_rgb)

N = length(vid_rgb);

for i = 1: N
    vid_ycbcr{i} = rgb2ycbcr(vid_rgb{i});
end

author: thomas-koehler

repo: SupER


text:
function calcAllFeatureSets_STS(pathWORK,pathMINE,fSetNameType,outcome,setSize,nonTextStruct,textCells,textCellsName,paramAll,freedomMat,baseline,alpha,delta,nBoot)
% -------------------------------------------------------------------------
% function calcAllFeatureSets_STS(pathWORK,pathMINE,fSetNameType,outcome,setSize,nonTextStruct,textCells,textCellsName,paramAll,freedomMat,baseline,alpha,delta,nBoot)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set reduction for a given feature set type
% and for all experiments with different degrees of freedom. See ref. [1]
% for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the STS WORKSPACE directory.
% - pathMINE: Full path to the MINE.jar executable. The executable can be
%             downloaded at: <http://www.exploredata.net/>.
% - fSetNameType: String specifying the name of the type of feature set
%                 (e.g., 'PET', 'SEPARATE', 'FUSED', etc.)
% - outcome: Column vector of size [nInst X 1] specifying the outcome status
%            (1 or 0) for all instances.
% - setSize: Size of the output feature set (typically set to 25 in ref. [1]).
% - nonTextStruct: Structure data for non-texture features. This structure
%                  is of the same format as the one saved as output to
%                  computeAllNonTextureFeatures_STS.m, for example.
% - textCells: Cell vector of organized texture cells for all the
%              different scans. Format: textCells = {cell1, cell2, etc.},
%              where cell1, cell2, etc. are the files saved as output to
%              organizeSeparateTextures_STS.m or organizeFusedTextures_STS.m,
%              for example.
% - textCellsName: Cell vector corresponding to the name of the
%                  corresponding cells in textCells.
%                  (e.g., textCellsName = {'PET','T1','T2FS'} for SEPARATE
%                  scans,  textCellsName = {'PET_T1','PET_T2FS'} for FUSED
%                  scans, as defined in ref. [1])
% - paramAll: Cell vector incorporating all texture extraction parameters
%             tested in textCells. See EXAMPLE below for more details.
% - freedomMat:  Matrix of row vectors of 1's and 0's to specify the degree
%                of freedom on texture extraction parameters for all
%                experiments. For example, for an nth experiment where
%                extraction parameters 1, 2 and 4 in paramAll are allowed
%                to vary, use freedomMat(n,:) = [1,1,0,1].
% - baseline: Vector of numerical values specifying the baseline texture
%             extraction parameters for each entry in paramAll. See EXAMPLE
%             below for more details.
% - alpha: Numerical values specifying the coefficient of the first part of
%          the Gain equation, as defined in ref. [1].
% - delta: Numerical values specifying the coefficient of the second part
%          of the Gain equation, as defined in ref. [1] (third part is set
%          to 0 in this function).
% - nBoot: Number of bootstrap samples to use.
%
% See <https://github.com/mvallieres/radiomics/tree/master/STUDIES/STSstudy/Functions>
% to find computeAllNonTextureFeatures_STS.m, organizeSeparateTextures_STS.m
% and organizeFusedTextures_STS.m. See masterScript_STS.m for a complete
% example of how to utilize the current function.
% -------------------------------------------------------------------------
% OUTPUTS: Feature sets are saved in a folder named 'FSET' in the STS
%          WORKSPACE.
% -------------------------------------------------------------------------
% EXAMPLE:
% MRIinv_cell = {'NoInv','Inv'};
% MRIweight_mat = [1/4,1/3,1/2,2/3,3/4];
% R_mat = [1/2,2/3,1,3/2,2];
% scale_cell = {'pixelW',1,2,3,4,5};
% algo_cell = {'Equal','Lloyd'};
% Ng_mat = [8,16,32,64];
%
% FOR FUSED SCANS
% paramAll = {MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat};
% freedomMat = [1 1 1 1 ; 0 0 0 0 ; 0 1 0 1]; (example)
% baseline = [1 3 3 1 2 3];
%
% FOR SEPARATE SCANS
% paramAll = {R_mat,scale_cell,algo_cell,Ng_mat};
% freedomMat = [1 1 1 1 1 1 ; 0 0 0 0 0 0 ;  0 1 0 1 0 1]; (example)
% baseline = [3 1 2 3];
%
% NOTE: paramAll must always be of the same format and size for SEPARATE
% and FUSED SCANS, with the same ordering of different extraction
% parameters.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd([pathWORK,'/FSET'])

nParamType = size(freedomMat,2);
nFreedom = size(freedomMat,1);
tStart = tic;
for i = 1:nFreedom
    nameSave=['FSET_',fSetNameType,'_'];
    for j = 1:nParamType
        nameSave = [nameSave,num2str(freedomMat(i,j))];
    end
    fprintf(['COMPUTING ',nameSave,' ... '])
    tic
    [fSet] = featureSetReduction(pathMINE,fSetNameType,outcome,setSize,nonTextStruct,textCells,textCellsName,paramAll,freedomMat(i,:),baseline,alpha,delta,nBoot);
    toc
    save(nameSave,'fSet')
end
time = toc(tStart);
fprintf('TOTAL TIME: %.2f seconds\n',time)

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [ccf] = spm_mar2ccf(mar,n)
% Get the cross covariance function from MAR coefficients or structure
% FORMAT [ccf] = spm_mar2ccf(mar,n)
%
% mar   - MAR coefficients or structure (see spm_mar.m)
% n     - number of time bins
%
% ccf   - (2*n + 1,i,j) cross covariance functions between I and J
%
% The mar coefficients are either specified in a cell array (as per
% spm_mar) or as a vector of (positive) coefficients as per spm_Q. The
% former are the negative values of the latter. If mar is a matrix of size
% d*p x d - it is assumed that the (positive) coefficients  run fast over
% lag = p, as per the DCM routines.
%
% see also:
%  spm_ccf2csd.m, spm_ccf2mar, spm_csd2ccf.m, spm_csd2mar.m, spm_mar2csd.m,
%  spm_csd2coh.m, spm_Q.m, spm_mar.m and spm_mar_spectral.m
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mar2ccf.m 6481 2015-06-16 17:01:47Z karl $

% Nyquist
%--------------------------------------------------------------------------
if nargin < 2, n = 128; end

% format coefficients into an array of negative coeficients (cf lag.a)
%--------------------------------------------------------------------------
if isvector(mar) && isnumeric(mar)
    mar = mar(:);
end
if isnumeric(mar)
    d  = size(mar,2);
    p  = size(mar,1)/d;
    for i = 1:d
        for j = 1:d
            for k = 1:p
                lag(k).a(i,j) = -mar((i - 1)*p + k,j);
            end
        end
    end
    mar = lag;
else
    d  = length(mar.lag(1).a);
    p  = length(mar.lag);
end

% covariance of innovations
%--------------------------------------------------------------------------
try
    c = mar.noise_cov;
catch
    c = eye(d,d);
end

% create AR representation and associated convolution kernels
%--------------------------------------------------------------------------
A     = cell(d,d);
B     = cell(d,d);
C     = cell(d,d);
for i = 1:d
    for j = 1:d
        a      = [0; mar.a((1:p) + (i - 1)*p,j)];
        A{i,j} = spdiags(ones(n,1)*a',-(0:p),n,n);
        C{i,j} = speye(n,n)*c(i,j);
    end
    B{i,i}     = speye(n,n);
end
A     = spm_cat(A);
B     = spm_cat(B);
C     = spm_cat(C);
K     = inv(B - A);

% compute cross-covariance matrices and reduces to an array of vectors
%--------------------------------------------------------------------------
CCF   = K*C*K';
ccf   = zeros(n,d,d);
for i = 1:d
    for j = 1:d
        ccf(:,i,j) = full(CCF((1:n) + (i - 1)*n,ceil(n/2) + (j - 1)*n));
    end
end


text:
function dimAdjust = getDimAdjust(L)
%GETDIMADJUST   Adjust dimension of discretization.
%   GETDIMADJUST(L) returns zero unless L is a LINOP.
%
%   A = GETDIMADJUST(L), where L is a LINOP,  returns a matrix of values, A, of
%   the same dimension as L.BLOCKS, which informs discretization methods how
%   the dimension L.DIM of that block must be adjusted (usually increased) so
%   that the highest order derivatives of each of the variables appearing in
%   the LINOP system are discretized at the same dimension.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( ~isa(L, 'linop') )

    % No adjustment
    dimAdjust = 0;

else

    % The input adjustment size of the (j,k) entry is max(diffOrder(:,k))
    dimAdjust = max(getDiffOrder(L), [], 1);
    dimAdjust = max(dimAdjust, 0);
    dimAdjust = repmat(dimAdjust, size(L, 1), 1);

end

end


text:
function [ craft_vel craft_theta a e EoverM GMsun] = predictor(craft_v,angle1, source_planet,target_planet )
%PREDICTOR - ellipse determination function
% the equations in this procedure is based on the
%  This will determine the critical ellipse parametres - mainly a and e,
%  based on the trajectory characteristics of the craft.
% the procedure for determining a and e from craft trajectory is based on:
%'Satellite Orbits and Gravitational Assist for Planets' (2008)-Larry Bogan
%http://www.bogan.ca/astrpages.html

%target planets are indexed , as below.
%1 = sun
%2 = earth
%3 = jupiter
%4 = saturn
%5 = uranus
%6 =neptune

mass = [1.9891e+30,5.9736e+24,1.8983e+27,5.68462313752e+26,8.6810e25,1.0243e26, 750];
m_sun = mass(1);
% planet masses
m_craft = mass(7);
G = 6.67428e-11;
tp = target_planet;
sp = source_planet;
AU=149597870691;
radius = [ (6.955e8/AU) 1 5.2 9.582 20.083 30.1036];

%####################################################################
% Step 1  Determine GMsun
% This determines GM with respect to the sun. This is necesscary as all of
% the tjrajectory ellipses are with respect to the sun as a focus.

GMsun_metres = m_sun * G;

GMsun = GMsun_metres / AU;  % this GMsun in m/s)^2 - being converted to AU

GMsun = GMsun / 1e6; %to go from square metres to square kilometres
% Should be 887 AU/(km/s)^2

%####################################################################
%Step 2 Energy over mass relationship

KE = (m_craft * craft_v^2)/2;

GPE = GMsun * m_craft / radius(sp);
E = (KE - GPE);
EoverM = (KE - GPE)/ m_craft;

% if KE - GPE is negative, it means the orbit is still elliptical around
% the sun - if not it is a hyperbolic orbit and it will leave the solar
% system.

%####################################################################
%Step 3 determine semimajor axis

a = -0.5 * GMsun / EoverM;

%####################################################################
%Step 4 - work out circular velocity, based on using the semi major axis as
%a radius. The radius is converted to metres from AU, and period is
%converted to seconds from years.

P = a^(3/2);

Vc = 2 * pi * a * AU   / P;

Vc3 = Vc / (365 * 86164);
Vcirc = Vc3 / 1000;

%####################################################################
% Step 5 from semimajor axis determine the period
% using Kepeler's third law - already been done above

P = a^(3/2);

%####################################################################
%Step 6 work out the areal velocity -

V_wrt_planet = craft_v * cosd(angle1);

V_wrt_planet_AU = (V_wrt_planet / (AU / 1000)) * 86164 * 365;

A =  radius(sp) * V_wrt_planet_AU /2;

%A = area of ellipse / priod - areal velocity

%####################################################################
%Step 7
% determines eccentricity from Areal Velocity and period of Ellipse

KK = A * P / (pi * a^2);

e = sqrt(1 - KK^2);

%####################################################################
%Step 8 - from semimajor axis and eccentricity, we get apophelion
%and periphelion

r_a = a*(1 + e);
r_p = a*(1 - e);

%####################################################################
%Step 9 determine true anomoly - the angle between periphelion and the
%spacecraft

%1 = sun
%2 = earth
%3 = jupiter
%4 = saturn
%5 = uranus
%6 =neptune
%7 = our spacecraft

V_per = Vcirc * sqrt((1+e) /(1-e));
V_aps = Vcirc * sqrt((1-e) /(1+e));

GSI  = radius(tp) * (( mass(tp)/mass(1))^(2/5));
%
craft_vel = sqrt( (2* EoverM) + (2 * GMsun ./ (radius(tp)- GSI)));

%its reduced by tp.
craft_theta = acosd((a*(1 - e^2)/(radius(tp)- GSI) - 1)/e);
end


text:
function [ b, info ] = cgtsl ( n, c, d, e, b )

%*****************************************************************************80
%
%% CGTSL solves a complex general tridiagonal system.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, complex C(N); the subdiagonal of the matrix in entries C(2:N).
%
%    Input, complex D(N); the diagonal of the tridiagonal matrix.
%
%    Input, complex E(N), the superdiagonal of the matrix in entries E(1:N-1).
%
%    Input, complex B(N), the right hand side.
%
%    Output, complex B(N).  the solution.
%
%    Output, integer INFO.
%    0, normal value.
%    K, if the K-th element of the diagonal becomes exactly zero.  The
%    subroutine returns when this is detected.
%
  info = 0;
  c(1) = d(1);

  if ( 1 <= n-1 )

    d(1) = e(1);
    e(1) = 0.0;
    e(n) = 0.0;

    for k = 1 : n-1

      if ( cabs1 ( c(k) ) <= cabs1 ( c(k+1) ) )

        t      = b(k);
        b(k)   = b(k+1);
        b(k+1) = t;

        t      = c(k);
        c(k)   = c(k+1);
        c(k+1) = t;

        t      = d(k);
        d(k)   = d(k+1);
        d(k+1) = t;

        t      = e(k);
        e(k)   = e(k+1);
        e(k+1) = t;

      end

      if ( cabs1 ( c(k) ) == 0.0 )
        info = k;
        return
      end

      t = -c(k+1) / c(k);
      c(k+1) = d(k+1) + t * d(k);
      d(k+1) = e(k+1) + t * e(k);
      e(k+1) = 0.0;
      b(k+1) = b(k+1) + t * b(k);

    end

  end

  if ( cabs1 ( c(n) ) == 0.0 )
    info = n;
    return
  end
%
%  Back solve.
%
  b(n) = b(n) / c(n);

  if ( 1 < n )

    b(n-1) = ( b(n-1) - d(n-1) * b(n) ) / c(n-1);

    for k = n-2 : -1 : 1
      b(k) = ( b(k) - d(k) * b(k+1) - e(k) * b(k+2) ) / c(k);
    end

  end

  return
end


text:
function nam_mat = tif2mat(nam)
%% convert tiff files into mat files
% inputs:
%   nam: file names
% output:
%   nam_mat: name of the *.mat file
% Author: Pengcheng Zhou, Carnegie Mellon University

%% tiff file information
if isempty(nam)
    nam_mat = [];
    fprintf('empty tiff file');
    return;
else
    [tmp_dir, tmp_file, ~] = fileparts(nam);
    nam_mat = sprintf('%s%s%s.mat', tmp_dir, filesep, tmp_file);
end

info = imfinfo(nam);

if  isfield(info(1), 'ImageDescription') && ~isempty(strfind(info(1).ImageDescription,'ImageJ'))
    junk1=regexp(info(1).ImageDescription,'images=\d*','match');
    if isempty(junk1)
        junk1=regexp(info(1).ImageDescription,'frames=\d*','match');
    end
    junk2=strjoin(junk1);
    T=strread(junk2,'%*s %d','delimiter','=');
else
    T = length(info);   % number of frames
end
d1 = info.Height;   % height of the image
d2 = info.Width;    % width of the image
Ysiz = [d1, d2, T]';

fprintf('CNMF_E is converting TIFF file to *.mat file');
% create a mat file
Tchunk = min(T, round(2^29/d1/d2)); %each chunk uses at most 4GB
Y = smod_bigread2(nam, 1, Tchunk);  %#ok<*NASGU>
save(nam_mat, 'Y', 'Ysiz', '-v7.3');
if Tchunk==T
    return;
else
    data = matfile(nam_mat, 'Writable', true);
    t0 = Tchunk+1;
    while t0<=T
        num2read = min(t0+Tchunk-1, T) - t0 + 1;
        tmpY = smod_bigread2(nam, t0, num2read);
        data.Y(:, :, (1:num2read)+t0-1) = tmpY;
        t0 = t0 + num2read;
    end
end

author: zhoupc