text:
classdef MMF6 < PROBLEM
% <multi> <real> <multimodal>
% Multi-modal multi-objective test function

%------------------------------- Reference --------------------------------
% C. Yue, B. Qu, and J. Liang, A multi-objective particle swarm optimizer
% using ring topology for solving multimodal multiobjective Problems, IEEE
% Transactions on Evolutionary Computation, 2018, 22(5): 805-817.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        POS;    % Pareto optimal set for IGDX calculation
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            obj.D = 2;
            obj.lower    = [1,-1];
            obj.upper    = [3,2];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            temp = X(:,2)<=0 | X(:,2)<=1 & (X(:,1)<=7/6|8/6<X(:,1)&X(:,1)<=9/6|10/6<X(:,1)&X(:,1)<=11/6|13/6<X(:,1)&X(:,1)<=14/6|15/6<X(:,1)&X(:,1)<=16/6|17/6<X(:,1));
            y    = zeros(size(X,1),1);
            y(temp)  = X(temp,2) - sin(6*pi*abs(X(temp,1)-2)+pi);
            y(~temp) = X(~temp,2) - 1 - sin(6*pi*abs(X(~temp,1)-2)+pi);
            PopObj(:,1) = abs(X(:,1)-2);
            PopObj(:,2) = 1 - sqrt(PopObj(:,1)) + 2*y.^2;
        end
        %% Generate Pareto optimal solutions
        function R = GetOptimum(obj,N)
            % Generate points in Pareto optimal set
            obj.POS(:,1) = linspace(1,3,N/2)';
            obj.POS(:,2) = sin(6*pi*abs(obj.POS(:,1)-2)+pi);
            obj.POS = [obj.POS;obj.POS(:,1),obj.POS(:,2)+1];
            % Generate points on Pareto front
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R(:,1) = linspace(0,1,100)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case 'IGDX'
                    score = feval(metName,Population,obj.POS);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopDec = Population.decs;
            temp1  = PopDec(:,2)<=0 | PopDec(:,2)<=1 & (PopDec(:,1)<=7/6|8/6<PopDec(:,1)&PopDec(:,1)<=9/6|10/6<PopDec(:,1)&PopDec(:,1)<=11/6|13/6<PopDec(:,1)&PopDec(:,1)<=14/6|15/6<PopDec(:,1)&PopDec(:,1)<=16/6|17/6<PopDec(:,1));
            temp2  = PopDec(:,1)<=2;
            Draw(Population(temp1&temp2).objs,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 .5 .5],'Markeredgecolor',[1 .2 .2],{'\it f\rm_1','\it f\rm_2',[]});
            Draw(Population(temp1&~temp2).objs+0.05,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[.5 .5 1],'Markeredgecolor',[.2 .2 1]);
            Draw(Population(~temp1&temp2).objs+0.1,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[.5 1 .5],'Markeredgecolor',[.2 1 .2]);
            Draw(Population(~temp1&~temp2).objs+0.15,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 .5 1],'Markeredgecolor',[1 .2 1]);
            Draw(obj.PF,'-','LineWidth',1,'Color',[1 .2 .2]);
            Draw(obj.PF+0.05,'-','LineWidth',1,'Color',[.2 .2 1]);
            Draw(obj.PF+0.1,'-','LineWidth',1,'Color',[.2 1 .2]);
            Draw(obj.PF+0.15,'-','LineWidth',1,'Color',[1 .2 1]);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function bitsHitV = cumbitor(bitsHitV)
%"cumbitor"
%   Find cumulative bitor of all values in bitsHitV, using a "vectorized"
%   binary tree algorithm to reduce calls to bitor.
%
%JRA 12/15/04
%
%Usage:
%   function bitsHitV = cumbitor(bitsHitV)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

while length(bitsHitV) >= 2
    len = length(bitsHitV);
    if mod(len, 2) == 1
        last = bitsHitV(end);
        bitsHitV = bitor(bitsHitV(1:(len-1)/2), bitsHitV((len-1)/2+1:end-1));
        bitsHitV(1) = bitor(bitsHitV(1), last);
    else
        bitsHitV = bitor(bitsHitV(1:len/2), bitsHitV(len/2+1:end));
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function w = drfiLearnSaliencyFusionWeight( train_dir, gt_dir, num_segmentation, is_resize )
    % Assume that all training images are placed under train_dir.
    % The saliency maps of i-th segmentation are under the folder "i" (e.g.,
    % saliency maps of 3rd segmentation are under the folder "3").
    % Detailed introduction on learning the saliency fusion weight can be
    % found in our supplementary material.

    M = num_segmentation;

    % Resize all training images to a fixed size 200*200
    sub_dir_list = dir(fullfile(train_dir, '*'));

    ind = [];
    for m = 1 : length(sub_dir_list)
        if strcmp(sub_dir_list(m).name, '.') || strcmp(sub_dir_list(m).name, '..')
            ind = [ind, m];
            continue;
        end
    end

    % Remove '.' and '..'
    sub_dir_list(ind) = [];

    normh = 200;
    normw = 200;

    % Resize
    if is_resize
        for m = 1 : length(sub_dir_list)
            image_list = dir(fullfile(train_dir, sub_dir_list(m).name, '*.png'));
            sub_dir_name = sub_dir_list(m).name;

            parfor n = 1 : length(image_list)
                image = imread(fullfile(train_dir, sub_dir_name, image_list(n).name));

                image = imresize(image, [normh, normw]);

                imwrite(image, fullfile(train_dir, sub_dir_name, image_list(n).name));

%                 if mod(jx, 500) == 0
%                     fprintf( 'sub_dir: %s, jx: %d\n', sub_dir_name, jx );
%                 end
            end

            fprintf( '%d / %d\n', m, length(sub_dir_list) );
        end
    end

    image_list = dir(fullfile(train_dir, sub_dir_list(end).name, '*.png'));
    num_image = length(image_list);

    % prepare H and f
    H = zeros(M, M);
    f = zeros(M, 1);

    for ii = 1 : M * M
        [m, n] = ind2sub([M, M], ii);
        sub_dir_name_n = sub_dir_list(n).name;
        sub_dir_name_m = sub_dir_list(m).name;
        if m >= n
            temp = zeros(1, num_image);
            parfor k = 1 : num_image
                image_name = image_list(k).name;
                Akm = im2double(imread(fullfile(train_dir, sub_dir_name_m, image_name)));
                Akn = im2double(imread(fullfile(train_dir, sub_dir_name_n, image_name)));

                if size(Akm, 3) > 1
                    Akm = rgb2gray( Akm );
                end

                if size(Akn, 3) > 1
                    Akn = rgb2gray( Akn );
                end

                Nk = 1;%size(Akm, 1) * size(Akm, 2);
                temp(k) = sum(sum(Akm .* Akn)) / Nk;
                % fprintf( 'ix: %d, jx: %d, n: %d\n', ix, jx, n );
            end
            H(m, n) = 2 * sum( temp );
        else
            H(m, n) = H(n, m);
        end

        fprintf( 'Computing H, m: %d, n: %d\n', m, n );
    end
    H = H / num_image;
    save( 'H.mat', 'H' );
%     load( 'H.mat' );

    for m = 1 : M
        temp = zeros(1, num_image);
        sub_dir_name = sub_dir_list(m).name;
        parfor k = 1 : num_image
            image_name = image_list(k).name;
            Akm = im2double(imread(fullfile(train_dir, sub_dir_name, image_name)));
            if size(Akm, 3) > 1
                Akm = rgb2gray(Akm);
            end

            A = imread(fullfile(gt_dir, image_name));
            A = imresize(A, [normh, normw]);
            A = im2double( A );
            if size(A, 3) > 1
                A = rgb2gray(A);
            end

            A( A > 0.5 ) = 1.0;
            A( A < 0.5 ) = 0;

            % f(ix) = f(ix) - 2 * sum(sum(A .* Ani));
            % temp = temp - 2 * sum(sum(A .* Ani));
            Nk = 1;%size(A, 1) * size(A, 2);
            temp(k) = - 2 * sum(sum(A .* Akm)) / Nk;
        end
        f(m) = sum( temp );
        fprintf( 'comupting f, m: %d\n', m );
    end
    f = f / num_image;
    save( 'f.mat', 'f' );

    % Solve the quadratic programming problem
    w_init = ones(M, 1) / M;

    Aeq = ones(1, M);
    beq = 1;

    lb = zeros(M, 1);
    ub = ones(M, 1);

    opt = optimset( 'Algorithm', 'interior-point-convex' );
    w = quadprog(H, f, [], [], Aeq, beq, lb, ub, w_init, opt );

    w( w < 1e-6 ) = 0;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function ll = fgplvmSequenceLogLikelihood(model, X, Y, varargin)

% FGPLVMSEQUENCELOGLIKELIHOOD Log-likelihood of a sequence for the GP-LVM.
% FORMAT
% DESC returns the log probability of a given latent sequence and an
% associated observed data sequence under the posterior prediction
% induced by the training data for a given GP-LVM model.
% ARG model : the model for which the sequence prediction will be
% made.
% ARG X : the latent sequence for which the posterior distribution
% will be evaluated.
% ARG Y : the observed data sequence for which the posterior
% distribution will be evaluated.
% ARG P1, P2, P3 ... : optional additional arguments to be passed
% to the dynamics' model sequence log likelihood.
% RETURN ll : the sequence log likelihood.
%
% SEEALSO : fgplvmCreate, fgplvmOptimiseSequence, fgplvmSequenceObjective
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
%
% MODIFICATIONS : Carl Henrik Ek, 2007

% FGPLVM
if(nargin<3)
  error('This function requires at least two arguments.');
end

logTwoPi = log(2*pi);
if model.isMissingData
  [mu, covarSigma]  = gpPosteriorMeanCovar(model, X);
  for i = 1:model.d
    missing = true;
    if ~any(isnan(Y(:, i)))
      U = jitChol(covarSigma{i});
      missing = false;
    end
    Ydiff = (Y(:, i)-mu(:, i));
    ll = 0;
    if missing
      ind = find(~isnan(Ydiff));
      if length(ind) ~= 0
        U = jitChol(covarSigma{i}(ind, ind));
      end
    else
      ind = [1:size(Ydiff, 1)]';
    end
    if length(ind) ~= 0
      UinvYdiff = U'\Ydiff(ind);
      logDet = logdet([], U);
      ll = ll + logDet + (UinvYdiff'*UinvYdiff);
      ll = ll + logTwoPi*length(ind);
    end
  end
else
  [mu, covarSigma, factors] = gpPosteriorMeanCovar(model, X);
  missing = true;
  if ~any(isnan(Y))
    U = jitChol(covarSigma);
    missing = false;
  end
  Ydiff = (Y-mu);
  ll =0;
  for i = 1:model.d
    if missing
      ind = find(~isnan(Ydiff(:, i)));
      if length(ind) ~= 0
        U = jitChol(covarSigma(ind, ind));
      end
    else
      ind = [1:size(Ydiff, 1)]';
    end
    if length(ind) ~= 0
      UinvYdiff = U'\Ydiff(ind, i);
      logDet = logdet([], U);
      ll = ll + logDet + log(factors(i)) + (UinvYdiff'*UinvYdiff)/factors(i);
      ll = ll + logTwoPi*length(ind);
    end
  end
end
ll = -0.5*ll;

if isfield(model, 'dynamics') & ~isempty(model.dynamics)
  % A dynamics model is being used.
  feval = str2func([model.dynamics.type 'SequenceLogLikelihood']);
  if isfield(model, 'dynamicsBalancing') & ~isempty(model.dynamicsBalancing)
    ll = ll + model.dynamicsBalancing*feval(model.dynamics, X, varargin{:});
  else
    ll = ll + feval(model.dynamics, X, varargin{:});
  end
elseif isfield(model, 'prior') &  ~isempty(model.prior)
  for i = 1:size(X, 1)
    ll = ll + priorLogProb(model.prior, X(i, :));
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function str = makeValidFieldName(str)
% From MATLAB doc: field names must begin with a letter, which may be
% followed by any combination of letters, digits, and underscores.
% Invalid characters will be converted to underscores, and the prefix
% "x0x[Hex code]_" will be added if the first character is not a letter.
    isoct=exist('OCTAVE_VERSION','builtin');
    pos=regexp(str,'^[^A-Za-z]','once');
    if(~isempty(pos))
        if(~isoct)
            str=regexprep(str,'^([^A-Za-z])','x0x${sprintf(''%X'',unicode2native($1))}_','once');
        else
            str=sprintf('x0x%X_%s',toascii(str(1)),str(2:end));
        end
    end
    if(isempty(regexp(str,'[^0-9A-Za-z_]', 'once' ))) return;  end
    if(~isoct)
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        pos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(pos)) return; end
        str0=str;
        pos0=[0 pos(:)' length(str)];
        str='';
        for i=1:length(pos)
            str=[str str0(pos0(i)+1:pos(i)-1) sprintf('_0x%X_',str0(pos(i)))];
        end
        if(pos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end

author: AvaisP


text:
function points = randomPointInBox3d(box, N, varargin)
%RANDOMPOINTINBOX3D Generate random point(s) within a 3D box.
%
%   PTS = randomPointInBox3d(BOX)
%   Generate a random point within the 3D box BOX. The result is a 1-by-3
%   row vector.
%
%   PTS = randomPointInBox3d(BOX, N)
%   Generates N points within the box. The result is a N-by-3 array.
%
%   BOX has the format:
%   BOX = [XMIN XMAX YMIN YMAX ZMIN ZMAX].
%
%   Example
%     % draw points within a box
%     box = [10 40 20 60 30 50];
%     pts =  randomPointInBox3d(box, 500);
%     figure(1); hold on;
%     drawBox3d(box);
%     drawPoint3d(pts, '.');
%     axis('equal');
%     axis([0 100 0 100 0 100]);
%     view(3);
%
%   See also
%   points3d, boxes3d

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-06-27, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

if nargin < 2
    N = 1;
end

% extract box bounds
xmin = box(1);
ymin = box(3);
zmin = box(5);

% compute size of box
dx = box(2) - xmin;
dy = box(4) - ymin;
dz = box(6) - zmin;

% compute point coordinates
points = [rand(N, 1)*dx+xmin , rand(N, 1)*dy+ymin , rand(N, 1)*dz+zmin];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function varargout = plot3(f, g, h, varargin)
%PLOT3   Plot for FUN objects in 3-D space.
%   PLOT3() is a three-dimensional analogue of PLOT().
%
%   PLOT3(X, Y, Z), where Z, Y, and Z are three CHEBTECH objects, plots a line
%   in 3-space. X, Y, and Z may be array-valued, but must have the same number
%   of columns.
%
%   Various line types, plot symbols, and colors may be obtained with PLOT3(X,
%   Y, Z, S) where s is a 1, 2 or 3 character string made from the characters
%   listed under the PLOT command.
%
%   H1 = PLOT3(X, Y, Z) returns a column vector of handles to lineseries
%   objects, one handle per plotted line (in the case of vector-valued CHEBTECH
%   objects). [H1, H2] returns a second vector of column handles, this time for
%   each of the marker plots.
%
% See also PLOT, PLOTDATA.

% 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

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

%%
% We can only plot real data in this way:
f = real(f);
g = real(g);
h = real(h);

% Get the data for plotting from PLOTDATA():
data = plotData(f, g, h);

%%
% Plot the curve:
h1 = plot3(data.xLine, data.yLine, data.zLine);
set(h1, 'Marker', 'none')
hold on

%%
% Plot the points:
h2 = plot3(data.xPoints, data.yPoints, data.zPoints, varargin{:});

% Change the style accordingly:
set(h2,'LineStyle', 'none')
if ( all(strcmp(get(h2, 'Marker'), 'none')) && (length(f) < 257) )
    set(h2, 'Marker', 'o')
end

%%

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

% Give an output if one was requested:
if ( nargout > 0 )
    varargout{1} = h1;
    varargout{2} = h2;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function x = intercept_model(nvols_per_run, varargin)
% Build design matrix X for intercepts
% given vector of session lengths [s1 s2 s3] in images
%
% :Usage:
% ::
%
%     x = intercept_model(nvols_per_run, [indx of dummy scans in each session])
%
% :Examples:
% ::
%
%    nvols_per_run = [166 166 144 137];
%    x = intercept_model(nvols_per_run);
%
%    x = intercept_model(repmat(166, 1, 5));
%
%    Xi = intercept_model(EXPT.FIR.nruns, 1:2);
%
% ..
%    tor modified april 07: separate column for each run
% ..

    if ~isvector(nvols_per_run), error('nvols_per_run must be a vector of the number of volumes for each run'); end
    nvols_per_run = nvols_per_run(:)'; %ensure row vector

    npoints = sum(nvols_per_run); % number of time points total
    nruns = length(nvols_per_run); % number of runs

    x = zeros(npoints, nruns);

    st = cumsum([1 nvols_per_run]);
    en = st(2:end) - 1; % ending values
    st = st(1:end-1); % starting values

    for i = 1:nruns
        x(st(i):en(i), i) = 1;
    end

    if length(varargin) > 0 && ~isempty(varargin{1})
        % Model dummy regressors
        x = [x model_dummy(npoints, st', varargin{1})];
    end
end

function x = model_dummy(npoints, st, wh)
    % dummy regressors for first scan or two (at least, that's the intended
    % use)
    % separate column for each run!

    k = length(wh) .* length(st);
    x = zeros(npoints, k);

    dosamecol = 0;

    if dosamecol
        for i = 1:k
            ind = st(r) + wh(i) - 1; % which to filter
            ind(ind < 1) = []; % if negative numbers (end of runs), this makes it ok
            x(ind, i) = 1;
        end

    else

        colindx = 1;

        for r = 1:length(st) % for each run
            for i = 1:length(wh)

                ind = st(r) + wh(i) - 1; % which to filter
                ind(ind < 1) = []; % if negative numbers (end of runs), this makes it ok

                x(ind, colindx) = 1;
                colindx = colindx + 1;
            end
        end

    end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% test for data interpolation using geodesic distances

rep  = 'data/';
rep  = '';

name = 'mm';
name = 'cavern';
name = 'disk';
name = 'toto';
name = 'boat';
name = 'cartoon';

clear options;
options.null = 0;

%% compute a binary shape
n = 300;
if strcmp(name, 'toto')
    M0 = load_image([rep name],n);
    M = sum(M0,3);
    M = 1-(M==M(1));
    % compute the constraint map
    CM = zeros(n) + Inf; CM(M==0) = -Inf;
    options.constraint_map = CM; % constraint the propagation inside the shape
    W = ones(n);
    mode = 'interpolation';
else
    n0 = [];
    if strcmp(name, 'disk')
        n0 = n;
    end
    if strcmp(name, 'cartoon')
        n0 = 400;
    end
    M = load_image([rep name],n0);
    M = rescale(crop(M,n)); M0 = M;
    M = rescale(sum(M,3));
    W = compute_edge_energy(M, 1.2);
    W = 1./rescale(W,0.001,1);
    mode = 'colorization';
end

% RGB2COL matrix
CM = rand(3); CM(:,1) = 1;
[CM,R] = qr(CM); CM(:,1) = 1/3; CM = CM';

%% pick some points
b = 1;
points = [];
while b==1
    clf;
    hold on;
    imagesc(rescale(M)); axis image; axis off;
    if not(isempty(points))
        plot_scattered(points(end:-1:1,:));
    end
    colormap gray(256); axis ij;
    hold off;
    [y,x,b] = ginput(1);
    if b==1
        points(:,end+1) = [x;y];
    end
end
npoints = size(points,2);
% random colors (sqrt to make it more colorfull
f = rand(3,npoints);
f = ( f ./ repmat(sum(f,1),[3, 1]) ).^(1);
% use fixed palette
fc = [1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1; 0.2 0.7 1; 1 0.2 0.7; 0.7 1 0.2; 0.2 1 0.7]';
sel =  floor(rand(npoints,1)*6)+1;
f = fc(:,sel);

if strcmp(mode, 'colorization')
    % random CbCr
    f = CM*f; f = f(2:end,:);
end

%% perform interpolation
if strcmp(mode, 'colorization')
    sigma_list = 0.001/n;
    alpha_list = 3;
else
    % the smaller, the more interpolation is performed
    % the higher, the more diffusion it is
    sigma_list = [25 25 25 0.01 0.01 0.01]/n;
    % the higher, the more "voronoi"-like it is
    alpha_list = [1 3 50 1 3 50];
end

nsigma = length(sigma_list);
nrows = 1;
if nsigma>3
    nrows = 2;
end

clf;
for i=1:nsigma
    options.sigma = sigma_list(i);
    options.alpha = alpha_list(i);
    [A,G] = perform_geodesic_interpolation(W,points,f,options);

    if strcmp(mode, 'colorization')
        % cat the luminosity component
        A = cat(3, M,A);
        A = reshape(A, [n^2 3])';
        A = ( (CM^(-1)) * A )';
        A = reshape( A, [n n 3] );
    end
    A = clamp(A);
    subplot(nrows,ceil(nsigma/nrows), i);
    hold on;
    imagesc(rescale(A)); axis image; axis off; axis ij;
    plot_scattered(points(end:-1:1,:));
    hold off;
    title(['\alpha=' num2str(options.alpha) ', \sigma=' num2str(options.sigma)]);
end

author: alecjacobson

repo: gptoolbox


text:
function valstr = m2json(val)
    if isstruct(val)
        valstr = struct2json(val);
    elseif iscell(val)
        valstr = cell2json(val);
    elseif isa(val, 'numeric')
        sz = size(val);
        if length(find(sz>1))>1 % 2D or higher array
            valstr = '';
            for i = 1:sz(1)
                valsubstr = [sprintf('%.15g, ', val(i,:))];
                valsubstr = valsubstr(1:(end-2));
                valstr = [valstr ', [' valsubstr ']'];
            end
            valstr = valstr(3:end); % trail leading commas
        else
            valstr = [sprintf('%.15g, ', val)];
            valstr = valstr(1:(end-2));
        end
        if length(val)>1
            valstr = ['[' valstr ']'];
        elseif length(val) == 0
            valstr = '[]';
        end
        valstr = strrep(valstr, 'Inf', 'null');
        valstr = strrep(valstr, 'NaN', 'null');
    elseif ischar(val)
        if size(val, 1) > 1
            valstr = cell2json(cellstr(val));
        else
            val = checkescape(val); %add escape characters
            valstr = ['"' val '"'];
        end
    elseif islogical(val)
        if val
            valstr = 'true';
        else
            valstr = 'false';
        end
    else
        valstr = ''; % wtf is it?
    end

author: goGPS-Project


text:
% Mex file interface to the lp_solve 5.5 toolkit. Please see
% reference guide for more information.
%
% mxlpsolve is a low-level interface to the `lp_solve toolkit`. It may be called
% directly, or may be used to build higher level functions for solving
% various kinds of linear and mixed-integer linear programs. It uses an
% integer handle to point to a linear program.
%
%     return = mxlpsolve('add_column', lp, [column])
%
%     return = mxlpsolve('add_columnex', lp, [column])
%
%     return = mxlpsolve('add_constraint', lp, [row], constr_type, rh)
%
%     return = mxlpsolve('add_constraintex', lp, [row], constr_type, rh)
%
%     return = mxlpsolve('add_SOS', lp, name, sostype, priority, [sosvars], [weights])
%
%     return = mxlpsolve('column_in_lp', lp, [column])
%
%     mxlpsolve('default_basis', lp)
%
%     return = mxlpsolve('del_column', lp, column)
%
%     return = mxlpsolve('del_constraint', lp, del_row)
%
%     lp_handle = mxlpsolve('copy_lp', lp)
%
%     mxlpsolve('delete_lp', lp)
%
%     mxlpsolve('dualize_lp', lp)
%
%     mxlpsolve('free_lp', lp)
%
%     return = mxlpsolve('get_anti_degen', lp)
%
%     [bascolumn] = mxlpsolve('get_basis', lp {, nonbasic})
%
%     return = mxlpsolve('get_basiscrash', lp)
%
%     return = mxlpsolve('get_bb_depthlimit', lp)
%
%     return = mxlpsolve('get_bb_floorfirst', lp)
%
%     return = mxlpsolve('get_bb_rule', lp)
%
%     return = mxlpsolve('get_bounds_tighter', lp)
%
%     return = mxlpsolve('get_break_at_value', lp)
%
%     name = mxlpsolve('get_col_name', lp, column)
%
%     [names] = mxlpsolve('get_col_name', lp)
%
%     [column, return] = mxlpsolve('get_column', lp, col_nr)
%
%     [column, return] = mxlpsolve('get_columnex', lp, col_nr)
%
%     return = mxlpsolve('get_constr_type', lp, row)
%
%     return = mxlpsolve('get_constr_value', lp, row {, primsolution})
%
%     [constr_type] = mxlpsolve('get_constr_type', lp)
%
%     [constr, return] = mxlpsolve('get_constraints', lp)
%
%     [duals, return] = mxlpsolve('get_dual_solution', lp)
%
%     return = mxlpsolve('get_epsb', lp)
%
%     return = mxlpsolve('get_epsd', lp)
%
%     return = mxlpsolve('get_epsel', lp)
%
%     return = mxlpsolve('get_epsint', lp)
%
%     return = mxlpsolve('get_epsperturb', lp)
%
%     return = mxlpsolve('get_epspivot', lp)
%
%     lp_handle = mxlpsolve('get_handle', lp_name)
%
%     return = mxlpsolve('get_improve', lp)
%
%     return = mxlpsolve('get_infinite', lp)
%
%     return = mxlpsolve('get_lowbo', lp, column)
%
%     [return] = mxlpsolve('get_lowbo', lp)
%
%     return = mxlpsolve('get_lp_index', lp, orig_index)
%
%     name = mxlpsolve('get_lp_name', lp)
%
%     value = mxlpsolve('get_mat', lp, row, col)
%
%     [matrix, return] = mxlpsolve('get_mat', lp {, sparse})
%
%     return = mxlpsolve('get_max_level', lp)
%
%     return = mxlpsolve('get_maxpivot', lp)
%
%     return = mxlpsolve('get_mip_gap', lp, absolute)
%
%     return = mxlpsolve('get_nameindex', lp, name, isrow)
%
%     return = mxlpsolve('get_Ncolumns', lp)
%
%     return = mxlpsolve('get_negrange', lp)
%
%     return = mxlpsolve('get_nonzeros', lp)
%
%     return = mxlpsolve('get_Norig_columns', lp)
%
%     return = mxlpsolve('get_Norig_rows', lp)
%
%     return = mxlpsolve('get_Nrows', lp)
%
%     return = mxlpsolve('get_obj_bound', lp)
%
%     [row_vec, return] = mxlpsolve('get_obj_fn', lp)
%
%     [row_vec, return] = mxlpsolve('get_obj_fun', lp)
%
%     return = mxlpsolve('get_objective', lp)
%
%     name = mxlpsolve('get_objective_name', lp)
%
%     return = mxlpsolve('get_orig_index', lp, lp_index)
%
%     name = mxlpsolve('get_origcol_name', lp, column)
%
%     [names] = mxlpsolve('get_origcol_name', lp)
%
%     name = mxlpsolve('get_origrow_name', lp, row)
%
%     [names] = mxlpsolve('get_origrow_name', lp)
%
%     return = mxlpsolve('get_pivoting', lp)
%
%     return = mxlpsolve('get_presolve', lp)
%
%     return = mxlpsolve('get_presolveloops', lp)
%
%     [pv, return] = mxlpsolve('get_primal_solution', lp)
%
%     return = mxlpsolve('get_print_sol', lp)
%
%     return = mxlpsolve('get_rh', lp, row)
%
%     [rh] = mxlpsolve('get_rh', lp)
%
%     return = mxlpsolve('get_rh_range', lp, row)
%
%     [rh_ranges] = mxlpsolve('get_rh_range', lp)
%
%     [row, return] = mxlpsolve('get_row', lp, row_nr)
%
%     [row, return] = mxlpsolve('get_rowex', lp, row_nr)
%
%     name = mxlpsolve('get_row_name', lp, row)
%
%     [names] = mxlpsolve('get_row_name', lp)
%
%     return = mxlpsolve('get_scalelimit', lp)
%
%     return = mxlpsolve('get_scaling', lp)
%
%     [objfrom, objtill, objfromvalue, objtillvalue, return] = mxlpsolve('get_sensitivity_obj', lp)
%
%     [objfrom, objtill, objfromvalue, objtillvalue, return] = mxlpsolve('get_sensitivity_objex', lp)
%
%     [duals, dualsfrom, dualstill, return] = mxlpsolve('get_sensitivity_rhs', lp)
%
%     [duals, dualsfrom, dualstill, return] = mxlpsolve('get_sensitivity_rhsex', lp)
%
%     return = mxlpsolve('get_simplextype', lp)
%
%     [obj, x, duals, return] = mxlpsolve('get_solution', lp)
%
%     return = mxlpsolve('get_solutioncount', lp)
%
%     return = mxlpsolve('get_solutionlimit', lp)
%
%     return = mxlpsolve('get_status', lp)
%
%     return = mxlpsolve('get_statustext', lp, statuscode)
%
%     return = mxlpsolve('get_timeout', lp)
%
%     return = mxlpsolve('get_total_iter', lp)
%
%     return = mxlpsolve('get_total_nodes', lp)
%
%     return = mxlpsolve('get_upbo', lp, column)
%
%     [upbo] = mxlpsolve('get_upbo', lp)
%
%     return = mxlpsolve('get_var_branch', lp, column)
%
%     [var_branch] = mxlpsolve('get_var_branch', lp)
%
%     return = mxlpsolve('get_var_dualresult', lp, index)
%
%     return = mxlpsolve('get_var_primalresult', lp, index)
%
%     return = mxlpsolve('get_var_priority', lp, column)
%
%     [var_priority] = mxlpsolve('get_var_priority', lp)
%
%     [var, return] = mxlpsolve('get_variables', lp)
%
%     return = mxlpsolve('get_verbose', lp)
%
%     return = mxlpsolve('get_working_objective', lp)
%
%     return = mxlpsolve('has_BFP', lp)
%
%     return = mxlpsolve('has_XLI', lp)
%
%     return = mxlpsolve('is_add_rowmode', lp)
%
%     return = mxlpsolve('is_anti_degen', lp, testmask)
%
%     return = mxlpsolve('is_binary', lp, column)
%
%     [binary] = mxlpsolve('is_binary', lp)
%
%     return = mxlpsolve('is_break_at_first', lp)
%
%     return = mxlpsolve('is_constr_type', lp, row, mask)
%
%     return = mxlpsolve('is_debug', lp)
%
%     return = mxlpsolve('is_feasible', lp, [values] {, threshold})
%
%     return = mxlpsolve('is_free', lp, column)
%
%     return = mxlpsolve('is_unbounded', lp, column)
%
%     [free] = mxlpsolve('is_free', lp)
%
%     [free] = mxlpsolve('is_unbounded', lp)
%
%     return = mxlpsolve('is_infinite', lp, value)
%
%     return = mxlpsolve('is_int', lp, column)
%
%     [int] = mxlpsolve('is_int', lp)
%
%     return = mxlpsolve('is_integerscaling', lp)
%
%     return = mxlpsolve('is_maxim', lp)
%
%     return = mxlpsolve('is_nativeBFP', lp)
%
%     return = mxlpsolve('is_nativeXLI', lp)
%
%     return = mxlpsolve('is_negative', lp, column)
%
%     [negative] = mxlpsolve('is_negative', lp)
%
%     return = mxlpsolve('is_piv_mode', lp, testmask)
%
%     return = mxlpsolve('is_piv_rule', lp, rule)
%
%     return = mxlpsolve('is_presolve', lp, testmask)
%
%     return = mxlpsolve('is_scalemode', lp, testmask)
%
%     return = mxlpsolve('is_scaletype', lp, scaletype)
%
%     return = mxlpsolve('is_semicont', lp, column)
%
%     [semicont] = mxlpsolve('is_semicont', lp)
%
%     return = mxlpsolve('is_SOS_var', lp, column)
%
%     [SOS_var] = mxlpsolve('is_SOS_var', lp)
%
%     return = mxlpsolve('is_trace', lp)
%
%     return = mxlpsolve('is_use_names', lp, isrow)
%
%     versionstring = mxlpsolve('lp_solve_version')
%
%     lp_handle = mxlpsolve('make_lp', rows, columns)
%
%     mxlpsolve('print_constraints', lp {, columns})
%
%     return = mxlpsolve('print_debugdump', lp, filename)
%
%     mxlpsolve('print_duals', lp)
%
%     mxlpsolve('print_lp', lp)
%
%     mxlpsolve('print_objective', lp)
%
%     mxlpsolve('print_scales', lp)
%
%     mxlpsolve('print_solution', lp {, columns})
%
%     mxlpsolve('print_str', lp, str)
%
%     mxlpsolve('print_tableau', lp)
%
%     [handle_vec] = mxlpsolve('print_handle')
%
%     lp_handle = mxlpsolve('read_freeMPS', filename {, verbose})
%
%     lp_handle = mxlpsolve('read_lp_file', filename {, verbose {, lp_name}})
%
%     lp_handle = mxlpsolve('read_lp', filename {, verbose {, lp_name}})
%
%     lp_handle = mxlpsolve('read_LP', filename {, verbose {, lp_name}})
%
%     lp_handle = mxlpsolve('read_mps', filename {, verbose})
%
%     lp_handle = mxlpsolve('read_MPS', filename {, verbose})
%
%     return = mxlpsolve('read_params', lp, filename {, options})
%
%     return = mxlpsolve('reset_params', lp)
%
%     lp_handle = mxlpsolve('read_XLI', xliname, modelname {, dataname {, options {, verbose}}})
%
%     return = mxlpsolve('set_add_rowmode', lp, turnon)
%
%     mxlpsolve('set_anti_degen', lp, anti_degen)
%
%     return = mxlpsolve('set_basis', lp, [bascolumn], nonbasic)
%
%     mxlpsolve('set_basiscrash', lp, mode)
%
%     mxlpsolve('set_basisvar', lp, basisPos, enteringCol)
%
%     mxlpsolve('set_bb_depthlimit', lp, bb_maxlevel)
%
%     mxlpsolve('set_bb_floorfirst', lp, bb_floorfirst)
%
%     mxlpsolve('set_bb_rule', lp, bb_rule)
%
%     return = mxlpsolve('set_BFP', lp, filename)
%
%     return = mxlpsolve('set_binary', lp, column, must_be_bin)
%
%     return = mxlpsolve('set_binary', lp, [must_be_bin])
%
%     return = mxlpsolve('set_bounds', lp, column, lower, upper)
%
%     return = mxlpsolve('set_bounds', lp, [lower], [upper])
%
%     mxlpsolve('set_bounds_tighter', lp, tighten)
%
%     mxlpsolve('set_break_at_first', lp, break_at_first)
%
%     mxlpsolve('set_break_at_value', lp, break_at_value)
%
%     return = mxlpsolve('set_col_name', lp, column, name)
%
%     return = mxlpsolve('set_col_name', lp, [names])
%
%     return = mxlpsolve('set_column', lp, col_no, [column])
%
%     return = mxlpsolve('set_columnex', lp, col_no, [column])
%
%     return = mxlpsolve('set_constr_type', lp, row, con_type)
%
%     return = mxlpsolve('set_constr_type', lp, [con_type])
%
%     mxlpsolve('set_debug', lp, debug)
%
%     mxlpsolve('set_epsb', lp, epsb)
%
%     mxlpsolve('set_epsd', lp, epsd)
%
%     mxlpsolve('set_epsel', lp, epsel)
%
%     mxlpsolve('set_epsint', lp, epsint)
%
%     mxlpsolve('set_epslevel', lp, epslevel)
%
%     mxlpsolve('set_epsperturb', lp, epsperturb)
%
%     mxlpsolve('set_epspivot', lp, epspivot)
%
%     return = mxlpsolve('set_free', lp, column)
%
%     return = mxlpsolve('set_unbounded', lp, column)
%
%     mxlpsolve('set_improve', lp, improve)
%
%     mxlpsolve('set_infinite', lp, infinite)
%
%     return = mxlpsolve('set_int', lp, column, must_be_int)
%
%     return = mxlpsolve('set_int', lp, [must_be_int])
%
%     return = mxlpsolve('set_lowbo', lp, column, value)
%
%     return = mxlpsolve('set_lowbo', lp, [values])
%
%     return = mxlpsolve('set_lp_name', lp, name)
%
%     return = mxlpsolve('set_mat', lp, [matrix])
%
%     return = mxlpsolve('set_mat', lp, row, column, value)
%
%     mxlpsolve('set_maxim', lp)
%
%     mxlpsolve('set_maxpivot', max_num_inv)
%
%     mxlpsolve('set_minim', lp)
%
%     mxlpsolve('set_mip_gap', lp, absolute, mip_gap)
%
%     mxlpsolve('set_negrange', negrange)
%
%     return = mxlpsolve('set_obj', lp, column, value)
%
%     return = mxlpsolve('set_obj', lp, [values])
%
%     mxlpsolve('set_obj_bound', lp, obj_bound)
%
%     return = mxlpsolve('set_obj_fn', lp, [row])
%
%     return = mxlpsolve('set_obj_fnex', lp, [row])
%
%     return = mxlpsolve('set_outputfile', lp, filename)
%
%     mxlpsolve('set_pivoting', lp, pivoting)
%
%     mxlpsolve('set_preferdual', lp, dodual)
%
%     mxlpsolve('set_presolve', lp, do_presolve {, maxloops})
%
%     mxlpsolve('set_print_sol', lp, print_sol)
%
%     return = mxlpsolve('set_rh', lp, row, value)
%
%     return = mxlpsolve('set_rh', lp, [values])
%
%     return = mxlpsolve('set_rh_range', lp, row, deltavalue)
%
%     return = mxlpsolve('set_rh_range', lp, [deltavalues])
%
%     mxlpsolve('set_rh_vec', lp, [rh])
%
%     return = mxlpsolve('set_row', lp, row_no, [row])
%
%     return = mxlpsolve('set_rowex', lp, row_no, [row])
%
%     return = mxlpsolve('set_row_name', lp, row, name)
%
%     return = mxlpsolve('set_row_name', lp, [names])
%
%     mxlpsolve('set_scalelimit', lp, scalelimit)
%
%     mxlpsolve('set_scaling', lp, scalemode)
%
%     return = mxlpsolve('set_semicont', lp, column, must_be_sc)
%
%     return = mxlpsolve('set_semicont', lp, [must_be_sc])
%
%     mxlpsolve('set_sense', lp, maximize)
%
%     mxlpsolve('set_simplextype', lp, simplextype)
%
%     mxlpsolve('set_solutionlimit', lp, simplextype)
%
%     mxlpsolve('set_timeout', lp, sectimeout)
%
%     mxlpsolve('set_trace', lp, trace)
%
%     return = mxlpsolve('set_upbo', lp, column, value)
%
%     return = mxlpsolve('set_upbo', lp, [values])
%
%     mxlpsolve('set_use_names', lp, isrow, use_names)
%
%     return = mxlpsolve('set_var_branch', lp, column, branch_mode)
%
%     return = mxlpsolve('set_var_branch', lp, [branch_mode])
%
%     return = mxlpsolve('set_var_weights', lp, [weights])
%
%     mxlpsolve('set_verbose', lp, verbose)
%
%     return = mxlpsolve('set_XLI', lp, filename)
%
%     result = mxlpsolve('solve', lp)
%
%     return = mxlpsolve('time_elapsed', lp)
%
%     mxlpsolve('unscale', lp)
%
%     return = mxlpsolve('write_freemps', lp, filename)
%
%     return = mxlpsolve('write_freeMPS', lp, filename)
%
%     return = mxlpsolve('write_lp', lp, filename)
%
%     return = mxlpsolve('write_LP', lp, filename)
%
%     return = mxlpsolve('write_mps', lp, filename)
%
%     return = mxlpsolve('write_MPS', lp, filename)
%
%     return = mxlpsolve('write_params', lp, filename {, options})
%
%     return = mxlpsolve('write_XLI', lp, filename {, options {, results}})
%

disp('mxlpsolve driver not found !!!');
disp('Check if mxlpsolve.dll is on your system and in a directory known to MATLAB.');
disp('Press enter to see the paths where MATLAB looks for the driver.');
pause
path
disp('A path can be added via the menu: File, Set Path');
error('mxlpsolve.dll not found');

author: opencobra

repo: cobratoolbox


text:
function pdf = gompertz_pdf ( x, a, b )

%*****************************************************************************80
%
%% GOMPERTZ_PDF evaluates the Gompertz PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) = B * A**X / exp ( B * ( A**X - 1 ) / log ( A ) )
%
%    for
%
%      0.0 <= X
%      1.0 <  A
%      0.0 <  B
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Johnson, Kotz, and Balakrishnan,
%    Continuous Univariate Distributions, Volume 2, second edition,
%    Wiley, 1994, pages 25-26.
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameters of the PDF.
%    1 < A, 0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < 0.0 )

    pdf = 0.0;

  elseif ( 1.0 < a )

    pdf = exp ( log ( b ) + x * log ( a ) - ( b / log ( a ) ) * ( a^x - 1.0 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = duplicator_fh(duplication_indices,xdim,w)
%
% This factory creates a function handle to an MV2DF, which represents the
% function:
%
%    y = x(duplication_indices)
%

if nargin==0
    test_this();
    return;
end

map = @(x) x(duplication_indices);

%xdim = max(duplication_indices);
ydim = length(duplication_indices);
c = zeros(1,ydim);
r = zeros(1,ydim);
at = 1;
for i=1:xdim
    ci = find(duplication_indices==i);
    n = length(ci);
    r(at:at+n-1) = i;
    c(at:at+n-1) = ci;
    at = at + n;
end
reverse = sparse(r,c,1,xdim,ydim);

transmap = @(y) reverse*y;

f = @(w) linTrans(w,map,transmap);

if exist('w','var') && ~isempty(w)
    f = f(w);
end

end

function test_this()

dup = [ 1 3 1 3];
x = [ 1 2 3 4];

f = duplicator_fh(dup,length(x));

y = f(x),

test_MV2DF(f,x);

end

author: nesl


text:
% READEDF - read eeg data in EDF format.
%
% Usage:
%    >> [data,header] = readedf(filename);
%
% Input:
%    filename - file name of the eeg data
%
% Output:
%    data   - eeg data in (channel, timepoint)
%    header - structured information about the read eeg data
%      header.length - length of header to jump to the first entry of eeg data
%      header.records - how many frames in the eeg data file
%      header.duration - duration (measured in second) of one frame
%      header.channels - channel number in eeg data file
%      header.channelname - channel name
%      header.transducer - type of eeg electrods used to acquire
%      header.physdime - details
%      header.physmin - details
%      header.physmax - details
%      header.digimin - details
%      header.digimax - details
%      header.prefilt - pre-filterization spec
%      header.samplerate - sampling rate
%
% Author: Jeng-Ren Duann, CNL/Salk Inst., 2001-12-21

% Copyright (C) Jeng-Ren Duann, CNL/Salk Inst., 2001-12-21
%
% 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.

% 03-21-02 editing header, add help -ad

function [data,header] = readedf(filename);

if nargin < 1
    help readedf;
    return;
end

fp = fopen(filename,'r','ieee-le');
if fp == -1,
  error('File not found ...!');
  return;
end

hdr = setstr(fread(fp,256,'uchar')');
header.length = str2num(hdr(185:192));
header.records = str2num(hdr(237:244));
header.duration = str2num(hdr(245:252));
header.channels = str2num(hdr(253:256));
header.channelname = setstr(fread(fp,[16,header.channels],'char')');
header.transducer = setstr(fread(fp,[80,header.channels],'char')');
header.physdime = setstr(fread(fp,[8,header.channels],'char')');
header.physmin = str2num(setstr(fread(fp,[8,header.channels],'char')'));
header.physmax = str2num(setstr(fread(fp,[8,header.channels],'char')'));
header.digimin = str2num(setstr(fread(fp,[8,header.channels],'char')'));
header.digimax = str2num(setstr(fread(fp,[8,header.channels],'char')'));
header.prefilt = setstr(fread(fp,[80,header.channels],'char')');
header.samplerate = str2num(setstr(fread(fp,[8,header.channels],'char')'))./header.duration;

fseek(fp,header.length,-1);
data = fread(fp,'int16');
fclose(fp);

data = reshape(data,header.duration*header.samplerate(1),header.channels,header.records);
temp = [];
for i=1:header.records,
  temp = [temp data(:,:,i)'];
end
data = temp;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function a = log10(a)
%LOG10        Gradient logarithm  log10(a)
%

% written  12/06/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 05/22/09     S.M. Rump  improved log10
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  wng = warning;
  warning off

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  if isa(a.x,'intval')
    INTLAB_STDFCTS_LOG10_ = getappdata(0,'INTLAB_STDFCTS_LOG10_');
    clog10 = infsup(INTLAB_STDFCTS_LOG10_.INF,INTLAB_STDFCTS_LOG10_.SUP);
    ax = clog10 ./ full(a.x(:));
  else
    ax = 1./( log(10) * full(a.x(:)) );
  end
  a.x = log10(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  warning(wng)

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function cond = conex4_condition ( )

%*****************************************************************************80
%
%% CONEX4_CONDITION returns the L1 condition of the CONEX4 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real COND, the L1 condition number.
%
  a_norm = 33.0;
  b_norm = 136.0;
  cond = a_norm * b_norm;

  return
end


text:
%
% n=erlangbinv(p,rho)
%
% This function finds the smallest n such the Erlang B probability that
% a system with n servers, no waiting line, Poisson arrival rate lambda,
% service rate (per server) mu, and intensity rho=lambda/mu will have
% a probability <=p of having all servers busy.
%
% The Erlang B probability is given by
%
%  B=(rho^m/m!)/(sum(rho^k/k!),k=0..m)
%
% We use a recurrence relation which is more accurate than direct evaluation
% of the formula.  This recurrence relation is a "folk theorem".  The author
% would appreciate a reference to its first publication.  The recurrence is
%
%  B(0,rho)=1
%
%  B(n,rho)=(rho*B(n-1,rho)/n)/(1+rho*B(n-1,rho)/n)
%
% This routine simply loops through the recursion until B is <= p, and
% then returns n.
%
function n=erlangbinv(p,rho)
%
% Start the recursion with B=1.
%
B=1;
%
% Loop, iterating the recursion until the probability is <= p.
%
n=1;
while (1 == 1),
  B=((rho*B)/n)/(1+rho*B/n);
  if (B <= p),
    return;
  end;
  n=n+1;
end;

author: Sable

repo: mcbench-benchmarks


text:
function v = barywts(n)
%BARYWTS   Barycentric weights for Chebyshev points of 2nd kind.
%   BARYWTS(N) returns the N barycentric weights for polynomial interpolation on
%   a Chebyshev grid of the 2nd kind. The weights are normalised so that they
%   have infinity norm equal to 1 and the final entry is positive.
%
% See also BARY, CHEBPTS.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% See Thm. 5.2 of Trefethen, Approximation Theory and Approximation Practice,
% SIAM, 2013 for more information.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( n == 0 )                      % Special case (no points)
    v = [];
elseif ( n == 1 )                  % Special case (single point)
    v = 1;
else                               % General case
    v = [ones(n-1,1) ; .5];        % Note v(end) is positive.
    v(end-1:-2:1) = -1;
    v(1) = .5*v(1);
end

end


text:
function [w, g_sqr] = adagrad(w, g_sqr, grad, opts, lr)
%ADAGRAD
%   Example AdaGrad solver, for use with CNN_TRAIN and CNN_TRAIN_DAG.
%
%   Set the initial learning rate for AdaGrad in the options for
%   CNN_TRAIN and CNN_TRAIN_DAG. Note that a learning rate that works for
%   SGD may be inappropriate for AdaGrad; the default is 0.001.
%
%   If called without any input argument, returns the default options
%   structure.
%
%   Solver options: (opts.train.solverOpts)
%
%   `epsilon`:: 1e-10
%      Small additive constant to regularize variance estimate.
%
%   `rho`:: 1
%      Moving average window for variance update, between 0 and 1 (larger
%      values result in slower/more stable updating). This is similar to
%      RHO in AdaDelta and RMSProp. Standard AdaGrad is obtained with a RHO
%      value of 1 (use total average instead of a moving average).
%
%   A possibly undesirable effect of standard AdaGrad is that the update
%   will monotonically decrease to 0, until training eventually stops. This
%   is because the AdaGrad update is inversely proportional to the total
%   variance of the gradients seen so far.
%   With RHO smaller than 1, a moving average is used instead. This
%   prevents the final update from monotonically decreasing to 0.

% Copyright (C) 2016 Joao F. Henriques.
% 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).

if nargin == 0  % Return the default solver options
  w = struct('epsilon', 1e-10, 'rho', 1) ;
  return ;
end

g_sqr = g_sqr * opts.rho + grad.^2 ;

w = w - lr * grad ./ (sqrt(g_sqr) + opts.epsilon) ;

repo: LapSRN


text:
function prediction = predict(model,Xt,n, simfct, args)
% Iterative prediction of a trained LS-SVM NARX model (in recurrent mode)
%
% >> Yp = predict({Xw,Yw,type,gam,sig2}, Xt, nb)
% >> Yp = predict(model,                 Xt, nb)
%
% The model needs to be trained using Xw, Yw which is the result of
% windowize or windowizeNARX. The number of time lags for the model
% is determined by the dimension of the input, or if not
% appropriate, by the number of given starting values.
%
% By default, the model is evaluated on the past points using
% simlssvm. However, if one wants to use this procedure for other
% models, this default can be overwritten by your favorite training
% function. This function (denoted by simfct) has to follow the following syntax:
%
% >> simfct(model,inputs,arguments)
%
% thus:
%
% >> Yp = predict(model, Xt, nb, simfct)
% >> Yp = predict(model, Xt, nb, simfct, arguments)
%
%
% Full syntax
%
%     1. Using the functional interface for the LS-SVMs:
%
% >> Yp = predict({Xw,Yw,type,gam,sig2,kernel,preprocess}, Xt, nb)
%
%       Outputs
%         Yp            : nb x m matrix with the predictions
%       Inputs
%         Xw            : N x d matrix with the inputs of the training data
%         Yw            : N x w matrix with the outputs of the training data
%         type          : 'function estimation' ('f') or 'classifier' ('c')
%         gam           : Regularization parameter
%         sig2          : Kernel parameter (bandwidth in the case of the 'RBF_kernel')
%         kernel(*)     : Kernel type (by default 'RBF_kernel')
%         preprocess(*) : 'preprocess' or 'original' (by default)
%         Xt            : nb*d matrix of the starting points for the prediction
%         nb(*)         : Number of outputs to predict
%
%
%     2. Using the object oriented interface with LS-SVMs:
%
% >> Yp = predict(model, Xt, nb)
%
%       Outputs
%         Yp    : nb x m matrix with the predictions
%       Inputs
%         Xt    : nb x d matrix of the starting points for the prediction
%         nb(*) : Number of outputs to predict
%
%
%     3. Using another model:
%
% >> Yp = predict(model, Xt, nb, simfct, arguments)
%
%       Outputs
%         Yp           : nb x m matrix with the predictions
%       Inputs
%         Xt           : nb x d matrix of the starting points for the prediction
%         nb           : Number of outputs to predict
%         simfct       : Function used to evaluate a test point
%         arguments(*) : Cell with the extra arguments passed to simfct
%
% See also:
%   windowize, trainlssvm, simlssvm.

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

if size(Xt,2)~=1 & size(Xt,1)~=1,
  error('prediction only implemented for the one-dimensional autonomous case');
end
eval('model = initlssvm(model{:});',' ');
eval('xdim = model.x_dim;','xdim = max(size(Xt)); ');
eval('n; alle=0;','n=length(Xt)-xdim; alle=1;');
Xt = Xt(1:xdim);
Xt = reshape(Xt,length(Xt),1);
eval('simfct;','simfct=''simlssvm'';');
eval('model = trainlssvm(model);');

xdelays = length(Xt);
prediction = zeros(xdelays+n,1);
prediction(1:xdelays) = Xt;

% closed loop
eval('for i=1:n,  prediction(xdelays+i) = feval(simfct,model, prediction(i-1+(1:xdelays))'',args); end',...
     'for i=1:n,  prediction(xdelays+i) = feval(simfct,model,prediction(i-1+(1:xdelays))''); end');

if ~alle,
  prediction = prediction(xdelays+1:end);
end

author: peterhcharlton

repo: RRest


text:
function cdf = cosine_cdf ( x, a, b )

%*****************************************************************************80
%
%% COSINE_CDF evaluates the Cosine CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, the parameter of the PDF.
%    0.0 < B.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= a - pi * b )

    cdf = 0.0;

  elseif ( x <= a + pi * b )

    y = ( x - a ) / b;

    cdf = 0.5 + ( y + sin ( y ) ) / ( 2.0 * pi );

  elseif ( a + pi * b < x )

    cdf = 1.0;

  end

  return
end


text:
function test_suite = test_vectorAngle3d
%TESTVECTORANGLE3D  One-line description here, please.
%
%   output = testVectorAngle3d(input)
%
%   Example
%   testVectorAngle3d
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-11-16,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testOrthogonalUnitVectors(testCase) %#ok<*DEFNU>

v1 = [1 0 0];
v2 = [0 1 0];
v3 = [0 0 1];
exp = pi/2;

angle1 = vectorAngle3d(v1, v2);
testCase.assertEqual(exp, angle1);

angle2 = vectorAngle3d(v1, v3);
testCase.assertEqual(exp, angle2);

angle3 = vectorAngle3d(v2, v3);
testCase.assertEqual(exp, angle3);

angle1 = vectorAngle3d(v1, -v2);
testCase.assertEqual(exp, angle1);

angle2 = vectorAngle3d(v1, -v3);
testCase.assertEqual(exp, angle2);

angle3 = vectorAngle3d(v2, -v3);
testCase.assertEqual(exp, angle3);

function testOrthogonalVectors(testCase)

v1 = [3 0 0];
v2 = [0 4 0];
v3 = [0 0 5];
exp = pi/2;

angle1 = vectorAngle3d(v1, v2);
testCase.assertEqual(exp, angle1);

angle2 = vectorAngle3d(v1, v3);
testCase.assertEqual(exp, angle2);

angle3 = vectorAngle3d(v2, v3);
testCase.assertEqual(exp, angle3);

function testParallelVectors(testCase)

v1 = [3 0 0];
v2 = [5 0 0];
exp = 0;

angle1 = vectorAngle3d(v1, v2);
testCase.assertEqual(exp, angle1);

v1 = [3 4 5]*7;
v2 = [3 4 5]*11;
exp = 0;

angle1 = vectorAngle3d(v1, v2);
testCase.assertEqual(exp, angle1);

function testSingleByArray(testCase)

v0 = [7 0 0];
v1 = [3 0 0];
v2 = [0 4 0];
v3 = [0 0 5];
vecs = cat(1, v1, v2, v3);

exp = [0; pi/2; pi/2];
angles = vectorAngle3d(v0, vecs);
testCase.assertEqual(exp, angles);

function testArrayByArray(testCase)

v0 = [7 0 0];
v1 = [3 0 0];
v2 = [0 4 0];
v3 = [0 0 5];
vecs1 = cat(1, v0, v0, v0);
vecs2 = cat(1, v1, v2, v3);

exp = [0; pi/2; pi/2];
angles = vectorAngle3d(vecs1, vecs2);
testCase.assertEqual(exp, angles);


text:
classdef ECGtask_ECG_delineation < ECGtask

% ECGtask for ECGwrapper (for Matlab)
% ---------------------------------
%
% Description:
%
% Abstract class for defining ECGtask interface
%
% Adding user-defined QRS detectors:
% A QRS detector that has the following interface can be added to the task:
%
%     [positions_single_lead, position_multilead] = your_ECG_delineator( ECG_matrix, ECG_header, progress_handle, payload_in);
%
% where the arguments are:
%    + ECG_matrix, is a matrix size [ECG_header.nsamp ECG_header.nsig]
%    + ECG_header, is a struct with info about the ECG signal, such as:
%         .freq, the sampling frequency
%         .desc, description about the signals.
%    + progress_handle, is a handle to a waitbar object, that can be used
%          to track the progress within your function. See the
%          documentation about this class in this kit.
%    + payload_in, is a user data variable allowed to be sent each call to
%          your function. It is sent, via the payload property of this
%          class, for example:
%
%         this_ECG_wrappers.ECGtaskHandle.payload = your_variable;
%         this_ECG_wrappers.ECGtaskHandle.payload = {your_var1 your_var2};
%         this_ECG_wrappers.ECGtaskHandle.payload = load(cached_filenames);
%
%          In the context of delineation, it is thought to be a user
%          corrected, or "gold quality" QRS location, in order to
%          improve the wave delineation quality. If "payload_in" is a
%          struct, this function will automatically filter and time-shift
%          all QRS detection fields started with the string "corrected_".
%          For this purpose, QRScorrector task, automatically appends this
%          string to eachm anually reviewed QRS location series.
%
% the output of your function must be:
%    + positions_single_lead, a cell array size ECG_header.nsig with the
%          QRS sample locations found in each lead.
%    + position_multilead, a numeric vector with the QRS locations
%          calculated using multilead rules.
%
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 18/2/2013
% Last update: 18/2/2013

    properties(GetAccess = public, Constant)
        name = 'ECG_delineation';
        target_units = 'ADCu';
        doPayload = true;
        cAnnotationFields = { 'Pon' 'P' 'Poff' 'QRSon' 'qrs' 'Q' 'R' 'S' 'QRSoff' 'Ton' 'T' 'Toff' };

    end

    properties( GetAccess = public, SetAccess = private)
        % if user = memory;
        % memory_constant is the fraction respect to user.MaxPossibleArrayBytes
        % which determines the maximum input data size.
        memory_constant = 0.3;

        started = false;

    end

    properties( Access = private, Constant)

        cQRSdelineators = {'all-delineators' 'wavedet' };

    end

    properties( Access = private )

        delineators2do
        bWFDBdelineators
        tmp_path_local
        command_sep
        WFDB_bin_path
        lead_idx = [];
        lead_names

    end

    properties
        progress_handle
        only_ECG_leads = false;
        delineators = 'all-delineators';
        wavedet_config
        payload
        tmp_path
    end

    methods

        function obj = ECGtask_ECG_delineation(obj)

            % wavedet configuration structure
            obj.wavedet_config.setup.wavedet.QRS_detection_only = false;

        end

        function Start(obj, ECG_header, ECG_annotations)

            if( obj.only_ECG_leads )
                obj.lead_idx = get_ECG_idx_from_header(ECG_header);
            else
%                 'all-leads'
                obj.lead_idx = 1:ECG_header.nsig;
            end

            % lead names desambiguation
            str_aux = regexprep(cellstr(ECG_header.desc), '\W*(\w+)\W*', '$1');
            obj.lead_names = regexprep(str_aux, '\W', '_');

            [str_aux2, ~ , aux_idx] = unique( obj.lead_names );
            aux_val = length(str_aux2);

            if( aux_val ~= ECG_header.nsig )
                for ii = 1:aux_val
                    bAux = aux_idx==ii;
                    aux_matches = sum(bAux);
                    if( sum(bAux) > 1 )
                        obj.lead_names(bAux) = strcat( obj.lead_names(bAux), repmat({'v'}, aux_matches,1), cellstr(num2str((1:aux_matches)')) );
                    end
                end
            end

            obj.lead_names = regexprep(obj.lead_names, '\W*(\w+)\W*', '$1');
            obj.lead_names = regexprep(obj.lead_names, '\W', '_');

            if( strcmpi('all-delineators', obj.delineators) )
                obj.delineators2do = obj.cQRSdelineators(2:end);
            else
                if( ischar(obj.delineators) )
                    obj.delineators2do = {obj.delineators};
                else
                    obj.delineators2do = obj.delineators;
                end
            end

            % local path required to avoid network bottlenecks in distributed filesystems
            if( isunix() && exist('/scratch/', 'dir') )
                str_username = getenv('USER');
                obj.tmp_path_local = ['/scratch/' str_username filesep];
                if( ~exist(obj.tmp_path_local, 'dir') )
                    if(~mkdir(obj.tmp_path_local))
                        obj.tmp_path_local = '/scratch/';
                    end
                end
                obj.tmp_path = tempdir;
            else
                if( isempty(obj.tmp_path) )
                    obj.tmp_path = tempdir;
                    obj.tmp_path_local = tempdir;
                end
            end

            obj.started = true;

        end

        function payload_out = Process(obj, ECG, ECG_start_offset, ECG_sample_start_end_idx, ECG_header, ECG_annotations, ECG_annotations_start_end_idx )

            payload_out = [];

            if( ~obj.started )
                obj.Start(ECG_header);
                if( ~obj.started )
                    cprintf('*[1,0.5,0]', 'Task %s unable to be started for %s.\n', obj.name, ECG_header.recname);
                    return
                end
            end

            % payload property is used in this task to input an external QRS
            % detector, or manually corrected detections.
            aux_struct = get_QRS_from_payload(obj.payload, ECG_start_offset, ECG_sample_start_end_idx);

            cant_QRSdelineators = length(obj.delineators2do);

            for ii = 1:cant_QRSdelineators

                this_delineator = obj.delineators2do{ii};

                [this_delineator, this_delineator_name] = strtok(this_delineator, ':');
                if( isempty(this_delineator_name) )
                    this_delineator_name = this_delineator;
                else
                    this_delineator_name = this_delineator_name(2:end);
                end

                cprintf( 'Blue', [ 'Processing ECG delineator ' this_delineator_name '\n' ] );

                %% perform ECG delineation

                switch( this_delineator )

                    case 'wavedet'
                    %% Wavedet delineation

                        obj.progress_handle.checkpoint(['Processing ' this_delineator])

                        try

                            [position_multilead, positions_single_lead] = wavedet_interface(ECG, ECG_header, aux_struct, obj.lead_idx, obj.wavedet_config, ECG_sample_start_end_idx, ECG_start_offset, obj.progress_handle);

                            for jj = 1:length(obj.lead_idx)
                                payload_out.wavedet.(obj.lead_names{obj.lead_idx(jj)}) = positions_single_lead(jj);
                            end

                            payload_out.wavedet.multilead = position_multilead;

                        catch aux_ME

                            strAux = sprintf('Wavedet failed in recording %s\n', ECG_header.recname);
                            strAux = sprintf('%s\n', strAux);

                            report = getReport(aux_ME);
                            fprintf(2, '%s\nError report:\n%s', strAux, report);

                        end

                    case 'user'
                        %% user-defined delineator

                        obj.progress_handle.checkpoint(['Processing ' this_delineator_name])

                        try

                            if( exist(this_delineator_name) == 2 )

%                                 ud_func_pointer = eval(['@' this_delineator_name]);
                                ud_func_pointer = str2func(this_delineator_name);

                                obj.progress_handle.checkpoint([ 'User defined function: ' this_delineator_name])

                                ECG_header_aux = trim_ECG_header(ECG_header, obj.lead_idx);

                                [positions_single_lead, position_multilead] = ud_func_pointer( double(ECG(:,obj.lead_idx)), ECG_header_aux, obj.progress_handle, aux_struct);

                                % filter and offset delineation
                                for jj = 1:length(obj.lead_idx)

                                    % filter heartbeats within range
                                    bAux = positions_single_lead(jj).qrs >= ECG_sample_start_end_idx(1) & positions_single_lead(jj).qrs <= ECG_sample_start_end_idx(2);

                                    aux_struct = [];
                                    for fn = rowvec(fieldnames(positions_single_lead(jj)))
                                        aux_val = positions_single_lead(jj).(fn{1});
                                        if( any(strcmpi(obj.cAnnotationFields, fn{1})) )
                                            aux_struct.(fn{1}) = aux_val(bAux) + ECG_start_offset - 1;
                                        else
                                            aux_struct.(fn{1}) = aux_val;
                                        end
                                    end

                                    payload_out.(this_delineator_name).(obj.lead_names{obj.lead_idx(jj)}) = aux_struct;

                                end

                                if( ~isempty(position_multilead) )

                                    % filter heartbeats within range
                                    bAux = position_multilead.qrs >= ECG_sample_start_end_idx(1) & position_multilead.qrs <= ECG_sample_start_end_idx(2);

                                    aux_struct = [];
                                    for fn = rowvec(fieldnames(position_multilead))
                                        aux_val = position_multilead.(fn{1});
                                        if( any(strcmpi(obj.cAnnotationFields, fn{1})) )
                                            aux_struct.(fn{1}) = aux_val(bAux) + ECG_start_offset - 1;
                                        else
                                            aux_struct.(fn{1}) = aux_val;
                                        end
                                    end

                                    position_multilead = aux_struct;

                                    payload_out.(this_delineator_name).multilead = position_multilead;

                                end

                            else
                                disp_string_framed(2, sprintf('Function "%s" is not reachable in path.', this_delineator_name));
                                fprintf(1, 'Make sure that exist(%s) == 2\n',this_delineator_name);
                            end

                        catch aux_ME

                            disp_string_framed(2, sprintf('Delineator "%s" failed in recording %s lead %s', this_delineator_name, ECG_header.recname, ECG_header.desc(jj,:) ) );

                            report = getReport(aux_ME);
                            fprintf(2, 'Error report:\n%s', report);

                        end

                end

            end

        end

        function payload = Finish(obj, payload, ECG_header)

        end

        function payload = Concatenate(obj, plA, plB)

            if( isempty(plA) )

                payload = plB;

            else
                for this_ECG_delineator = rowvec(fieldnames(plA))

                    this_ECG_delineator = this_ECG_delineator{1};

                    if( isfield(plB, this_ECG_delineator)  )

                        for this_lead = rowvec(fieldnames(plA.(this_ECG_delineator)))

                            this_lead = this_lead{1};

                            if( isfield(plB.(this_ECG_delineator), this_lead)  )

                                for fn = rowvec(fieldnames(plA.(this_ECG_delineator).(this_lead) ))

                                    if( isfield(plB.(this_ECG_delineator).(this_lead), fn{1}) )
                                        payload.( this_ECG_delineator).(this_lead).(fn{1}) = [ colvec(plA.( this_ECG_delineator).(this_lead).(fn{1}) ); colvec( plB.(this_ECG_delineator).(this_lead).(fn{1}) ) ];
                                    else
                                        error('ECGtask_ECG_delineation:BadTMPfiles', ['Results from ' this_ECG_delineator '.' this_lead  '.' fn{1} ' not found. Skipping concatenation.'])
                                    end

                                end
                            else
                                error('ECGtask_ECG_delineation:BadTMPfiles', ['Results from ' this_ECG_delineator '.' this_lead ' not found. Skipping concatenation.'] )
                            end
                        end

                    else
                        error('ECGtask_ECG_delineation:BadTMPfiles', ['Results from ' this_ECG_delineator ' not found. Skipping concatenation.'] )
                    end
                end
            end

        end

        %% property restriction functions

        function set.delineators(obj,x)
            if( (ischar(x) || iscellstr(x)) )
                x = cellstr(x);
                aux_val = colvec(intersect(obj.cQRSdelineators, x));
                aux_idx = find(cellfun(@(a)(~isempty(strfind(a, 'user:'))), x));
                obj.delineators = [aux_val; colvec(x(aux_idx)) ];
            else
                warning('ECGtask_ECG_delineation:BadArg', 'Invalid delineators.');
            end
        end

        function set.tmp_path_local(obj,x)
            if( ischar(x) )
                if(exist(x, 'dir'))
                    obj.tmp_path_local = x;
                else
                    if(mkdir(x))
                        obj.tmp_path_local = x;
                    else
                        warning('ECGtask_ECG_delineation:BadArg', ['Could not create ' x ]);
                    end
                end

            else
                warning('ECGtask_ECG_delineation:BadArg', 'tmp_path_local must be a string.');
            end
        end

        function set.tmp_path(obj,x)
            if( ischar(x) )
                if(exist(x, 'dir'))
                    obj.tmp_path = x;
                else
                    if(mkdir(x))
                        obj.tmp_path = x;
                    else
                        warning('ECGtask_ECG_delineation:BadArg', ['Could not create ' x ]);
                    end
                end

            else
                warning('ECGtask_ECG_delineation:BadArg', 'tmp_path_local must be a string.');
            end
        end

    end

    methods ( Access = private )

    end

end

author: marianux

repo: ecg-kit

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


text:
% Chapter 11: Interior-point methods
%
%  log_utility_flow.m - Section 11.8.4: Network rate optimization
help Contents

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% this program is designed to induce yearly variability for the data after
% induced monthly variability by FFT using a simple linear relationship

function [yearly_corrected_precip]=yearly_precip_correction(monthly_corrected_precip)
% load the yearly precipitation after FFT
load('Pnew_precip');
P_FFT=Pnew_precip';

load(monthly_corrected_precip);

n=length(P_FFT); % years
m=length(monthly_corrected_precip);  % days

% calculate yearly precip
% the fourth column of Y is precip
j=1;
Z=zeros(n,1);
for i=1:365:m
    Z(j,1)=sum(monthly_corrected_precip(i:i+365-1,4));
    %Z(j,1)=sum(gP(i:i+365-1,1));
    j=j+1;
end

% calculate the ratio of yearly precip between those after FFT and
% generated by weather generator
for i=1:n
    P_ratio(i,1)=P_FFT(i,1)/Z(i,1);
end

% extend the yearly precip ratio to daily scale,the data in each
% year are the same
P_extent=zeros(m,1);
j=1;
for i=1:365:m
    P_extent(i:i+365-1,1)=P_ratio(j,1);
    j=j+1;
end

% adjust the daily precip generated by WG using above ratios
P_adjust=zeros(size(monthly_corrected_precip));
P_adjust(:,1:3)=monthly_corrected_precip(:,1:3);
for i=1:m
    P_adjust(i,4)=monthly_corrected_precip(i,4)*P_extent(i,1);
end
yearly_corrected_precip=P_adjust;

mm=find(yearly_corrected_precip(:,4)<0);
yearly_corrected_precip(mm,4)=0;

save('yearly_corrected_precip','yearly_corrected_precip')


text:
function excelrecord2fcn(T, n, endo, stringdates1, Y, Ytilde, pref, EPStilde)
% script excelrecord2
% records the information contained in the worksheets 'actual fitted' and 'resids' of the excel spreadsheet 'results.xls'

% first compute the cell for actual/ fitted

% create the cell that will be saved on excel
afcell={};

% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);

% loop over variables (horizontal dimension)
for ii=1:n
% create cell of actual/fitted for variable ii
temp=['actual and fitted: ' endo{ii,1}];
af_i=[temp {''} {''};{''} {''} {''};{''} {'sample'} {'fitted'};stringdates1 num2cell(Y(:,ii)) num2cell(Ytilde(:,ii))];
afcell=[afcell af_i vertspace];
end

% trim
afcell=afcell(:,1:end-1);

% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),afcell,'actual fitted','B2');
end

% then compute the cell for the residuals

% create the cell that will be saved on excel
horzspace=repmat({''},1,n);
rescell=[{'residuals'} horzspace;{''} horzspace;{''} endo';stringdates1 num2cell(EPStilde)];

% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),rescell,'resids','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
% GetSegyTraceHeaderInfo : Returns a array of a SEGY Y TraceHeader value
%
% Call :
% [value]=GetSegyHeaderInfo(SegyTraceHeaders,header)
%
% header is a header value like 'cdp','dt','TraceNumber'
%

%
% (C) 2001-2004, Thomas Mejer Hansen, tmh@gfy.ku.dk/thomas@cultpenguin.com
%
%
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
function [value]=GetSegyTraceHeaderInfo(SegyTraceHeaders,header)

  try
    value=[SegyTraceHeaders.(header)];
  catch
    SegymatVerbose(['Could not get  ',header,' header information from supplied SegyTraceHeaders']);
  end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function g = wangPriorGradient(prior, x)

% WANGPRIORGRADIENT Gradient wrt x of the Wang prior.

% PRIOR

% Compute gradient of prior
g = -prior.M./x;

author: SheffieldML


text:
% absolute change in area
function [data,units] = compute_absdarea(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  data{i} = abs(trx(fly).darea);
end
units = parseunits('mm^2/s');

author: kristinbranson

repo: JAABA


text:
% DEMO_VALIDATION_CSP - The demo exemplifies how to use the crossvalidation
%  function to validate CSP-based classification. The important issue here
%  is that the CSP analysis has to be performed WITHIN the cross-validation
%  on each training set. The matrix of spatial filters that is obtained from
%  the training set needs to be transfered to the test set. If you would like
%  to know more about such valiation issues, see [Lemm et al, Neuroimage
%  2011].
%  In the function crossvalidation, there is the possibility to specify
%  processing chains separately for training and test set. Variables obtained
%  from training data may be transfered to the test data. Each step in the
%  processing chain is an application of a 'proc_*' function that transforms
%  the features.

file= fullfile(BTB.DataDir, 'demoMat', 'VPkg_08_08_07', ...
               'calibration_motorimageryVPkg');
[cnt, mrk]= file_loadMatlab(file);
mrk= mrk_selectClasses(mrk, [1 2]);

[filt_b,filt_a]= butter(5, [9 13]/cnt.fs*2);
cnt= proc_filt(cnt, filt_b, filt_a);
fv= proc_segmentation(cnt, mrk, [750 3750]);

proc.train= {{'CSPW', @proc_cspAuto, 3}
             @proc_variance
             @proc_logarithm
            };
proc.apply= {{@proc_linearDerivation, '$CSPW'}
             @proc_variance
             @proc_logarithm
            };

crossvalidation(fv, {@train_RLDAshrink, 'Gamma',0}, ...
                'SampleFcn', {@sample_chronKFold, 8}, ...
                'Proc', proc)

author: bbci


text:
function variance = log_series_variance ( a )

%*****************************************************************************80
%
%% LOG_SERIES_VARIANCE returns the variance of the Logarithmic Series PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the PDF.
%    0.0 < A < 1.0.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  alpha = - 1.0 / log ( 1.0 - a );

  variance = a * alpha * ( 1.0 - alpha * a ) / ( 1.0 - a )^2;

  return
end


text:
function SQL_Store(writeWhat,dbname)
% SQL_Store 	Upload data to the mySQL database.
%
% Uploads data in the HCTSA.mat file in the current directory back into the
% mySQL database. Should be done to store the result new computations done by
% TS_Compute.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check Inputs
% ------------------------------------------------------------------------------
if nargin < 1
	writeWhat = 'null'; % 'nullerror'
    % find all nulls in the database and write over them if there are values in local files
end
if ~ismember(writeWhat,{'null','error','nullerror'})
    error('Unknown specifier ''%s''',writeWhat)
end
if nargin < 3
	dbname = '';
end

% ------------------------------------------------------------------------------
%% Open a mySQL database connection
% ------------------------------------------------------------------------------
% Much faster to use java connector than the Matlab database toolbox
% (exec commands are slow; update commands are even slower)
[dbc, dbname] = SQL_OpenDatabase(dbname,0,0);

% ------------------------------------------------------------------------------
%% Load local files
% ------------------------------------------------------------------------------
%% Read in information from local files
fid = 1; % haha no more logging option...
loadTimer = tic;
fprintf(fid,'Loading data from HCTSA.mat...');
load('HCTSA.mat')
fprintf(fid,' Done in %s.\n',BF_TheTime(toc(loadTimer)));
clear loadTimer

% ------------------------------------------------------------------------------
% Check that this file was actually retrieved from the database (such that these
% ts_ids and op_ids represent items stored in the linked mySQL database)
% ------------------------------------------------------------------------------
if ~fromDatabase
    theLoc = which('HCTSA.mat');
    error(['It looks like:\n%s\nwas generated using TS_init and does not match up ' ...
            'with ID entries a mySQL database...'],theLoc);
end

% ------------------------------------------------------------------------------
% Check that TS_CalcTime exists
% ------------------------------------------------------------------------------
% It will not exist if you've (by default) not retrieved calculation time data
% from the database, and not generated it from the computation
if ~exist('TS_CalcTime','var')
    error('No calculation time data found in HCTSA.mat');
end

% ------------------------------------------------------------------------------
%% Preliminary definitions
% ------------------------------------------------------------------------------
numTS = height(TimeSeries); % Number of time series
numOps = height(Operations); % Number of operations
ts_id_loc = TimeSeries.ID; % tsids in local file
op_id_loc = Operations.ID; % opids in local file
ts_ids_string = BF_cat(ts_id_loc,',');
op_ids_string = BF_cat(op_id_loc,',');

% ------------------------------------------------------------------------------
%% Check that nothing has been deleted since the data was loaded...
% ------------------------------------------------------------------------------
% Time series
SelectString = sprintf('SELECT COUNT(ts_id) FROM TimeSeries WHERE ts_id IN (%s)',ts_ids_string);
numTS_db = mysql_dbquery(dbc,SelectString);
numTS_db = numTS_db{1};

% Operations
SelectString = sprintf('SELECT COUNT(op_id) FROM Operations WHERE op_id IN (%s)',op_ids_string);
numOps_db = mysql_dbquery(dbc,SelectString);
numOps_db = numOps_db{1};

if (numOps == numOps_db) && (numTS == numTS_db)
    fprintf(1,'All local time series and operation ids still exist in %s. This is good.\n',dbname);
else
    if numTS_db < numTS
        fprintf(1,'There are %u time series that no longer match the database',(numTS-numTS_db));
    end

    if numOps_db < numOps
    	fprintf(1,'There are %u operations that no longer match the database',(numOps-numOps_db));
    end
    error(['It could be dangerous to write back to a changed database. ' ...
                'You should start a SQL_Retrieve from scratch.'])
end

% ------------------------------------------------------------------------------
%% Find the elements that are empty in the database
%       (and hopefully full in the local file)
% ------------------------------------------------------------------------------
% Parts of calculated subsection that are empty in storage
fprintf(1,'Retrieving %s elements from the Results table in %s...',writeWhat,dbname);

switch writeWhat
case 'null'
    % collect nulls in the database
    SelectString = sprintf(['SELECT ts_id, op_id FROM Results WHERE ts_id IN (%s)' ...
    					' AND op_id IN (%s) AND QualityCode IS NULL'],ts_ids_string,op_ids_string);
case 'nullerror'
    % Collect all NULLS and previous errors
    SelectString = sprintf(['SELECT ts_id, op_id, QualityCode FROM Results WHERE ts_id IN (%s)' ...
    					' AND op_id IN (%s) AND (QualityCode IS NULL OR QualityCode = 1)'], ...
        					ts_ids_string,op_ids_string);
case 'error'
    % Collect all previous errors (assume done a SQL_Retrieve using 'error' input)
    SelectString = sprintf(['SELECT ts_id, op_id FROM Results WHERE ts_id IN (%s)' ...
    					' AND op_id IN (%s) AND QualityCode = 1'], ...
        					ts_ids_string,op_ids_string);
end

retrievalTimer = tic; % Time the retrieval (should be fast)
[qrc,emsg] = mysql_dbquery(dbc,SelectString);
if ~isempty(emsg)
    fprintf(1,'\n'); error('Error selecting %s elements from %s',writeWhat,dbname);
elseif isempty(qrc)
    fprintf(1,'\nNo %s elements in this range in the database anymore! Nothing to write.\n',writeWhat);
    SQL_CloseDatabase(dbc); return
else
	fprintf(1,' Retrieved %u entries in %s\n',length(qrc),BF_TheTime(toc(retrievalTimer)));
end
clear retrievalTimer % Stop timing

ts_id_db = vertcat(qrc{:,1}); % ts_ids (in op_id pairs) of empty database elements in this ts_id/op_id range
op_id_db = vertcat(qrc{:,2}); % op_ids (in ts_id pairs) of empty database elements in this ts_id/op_id range
numWrite = length(ts_id_db);  % Number of database elements to attempt to write back to

% ------------------------------------------------------------------------------
% Give user feedback about what database writing will occur
% ------------------------------------------------------------------------------
switch writeWhat
case 'null'
    fprintf(1,['There are %u NULL entries in Results.\nWriting calculated ' ...
                    'elements of TS_DataMat to %s...\n'],numWrite,dbname);
case 'error'
    fprintf(1,['There are %u entries in Results (all previous errors) ' ...
                    'that are being written to %s...\n'],numWrite,dbname);
    fprintf(1,['Previous results stored as errors in the database WILL NOT ' ...
                                    'be overwritten with newer errors\n']);
case 'nullerror'
    q_db = qrc(:,3); % empties (NULL) and fatal error (1)
    q_db(cellfun(@isempty,q_db)) = {0}; % turn NULLs to 0s
    q_db = vertcat(q_db{:}); % turn cells to a numeric vector
    % so now nulls in database are labeled 0, and previous errors are labeled 1
    fprintf(1,['There are %u entries in Results (either NULL or previous errors) ' ...
                    'that are being written to %s...\n'],numWrite,dbname);
    fprintf(1,['Note that previous results stored as errors in the database WILL NOT ' ...
                                'be overwritten with newer errors\n']);
    fprintf(1,'However, NULLS will be written over with any result from the local files\n');
end

localIndex = zeros(numWrite,2);
localIndex(:,1) = arrayfun(@(x)find(ts_id_loc == x,1),ts_id_db); % Indices of rows in local file for each entry in the database
localIndex(:,2) = arrayfun(@(x)find(op_id_loc == x,1),op_id_db); % Indicies of columns in local file for each entry in the database
updateMe = zeros(numWrite,1); % Label iterations that should be written to the database

% ------------------------------------------------------------------------------
% Begin writing each element to the database, one at a time, using UPDATE commands
% ------------------------------------------------------------------------------

writeBackTimer = tic; % Time how long this takes to give user feedback
numReports = 3; % the number of time remaining updates to provide the user

for i = 1:numWrite

    % Retrieve the elements
    TS_DataMat_ij = TS_DataMat(localIndex(i,1),localIndex(i,2));
    TS_Quality_ij = TS_Quality(localIndex(i,1),localIndex(i,2));
    TS_CalcTime_ij = TS_CalcTime(localIndex(i,1),localIndex(i,2));

    switch writeWhat
    case 'null'
        if isfinite(TS_DataMat_ij)
            updateMe(i) = 1; % There is a value in TS_DataMat -- write it back to the NULL entry in the database
        end
    case 'error'
        if isfinite(TS_DataMat_ij) && TS_Quality_ij~=1
            updateMe(i) = 1; % There is a now a non-error value in TS_DataMat previously returned an error (in the database)
        end
    case 'nullerror'
        if isfinite(TS_DataMat_ij) && (q_db(i)==0 || TS_Quality_ij~=1)
            updateMe(i) = 1; % there is a value in TS_DataMat -- write it to the entry in the database
        end
		% (i) Has been calculated and a value stored in TS_DataMat (isfinite()), and
		% (ii) Either the database entry is NULL or we didn't get an error (prevents writing errors over errors)
    end

    if updateMe(i)

        if isnan(TS_CalcTime_ij) % happens when there is an error in the code
            TS_CalcTime_string = 'NULL';
        else
            TS_CalcTime_string = sprintf('%f',TS_CalcTime_ij);
        end

        % I can't see any simple way around running lots of single UPDATE commands (for each entry)
    	updateString = sprintf(['UPDATE Results SET Output = %19.17g, QualityCode = %u, CalculationTime = %s ' ...
						'WHERE ts_id = %u AND op_id = %u'],TS_DataMat_ij,TS_Quality_ij, ...
						TS_CalcTime_string,ts_id_db(i),op_id_db(i));
        [~,emsg] = mysql_dbexecute(dbc, updateString);
        if ~isempty(emsg)
            SQL_CloseDatabase(dbc) % close the database connection before calling the error...
        	error('Error storing (ts_id,op_id) = (%u,%u) to %s??!!\n%s\n', ...
                			TimeSeries.ID(localIndex(i,1)),Operations.ID(localIndex(i,2)),dbname,emsg);
        end
    end

    % Give user feedback on how long is remaining:
    if (i==50) && (i < numWrite/numReports) % give an initial estimate:
        fprintf(1,['Based on the first 50 retrievals, this is taking ' ...
                'approximately %s per entry to write to the database.\n'],BF_TheTime(toc(writeBackTimer)));
		fprintf(1,'Approximately %s remaining...\n',BF_TheTime(toc(writeBackTimer)/i*(numWrite-i)));
    elseif mod(i,floor(numWrite/numReports))==0 % Give numReports more time updates...
		fprintf(1,['Approximately %s remaining -- %u (/ %u possible) entries have been'  ...
			' written to %s...\n'],BF_TheTime(toc(writeBackTimer)/i*(numWrite-i)),sum(updateMe),i,dbname);
	end
end

% ------------------------------------------------------------------------------
% Finished writing back to the database!
% ------------------------------------------------------------------------------

fprintf(1,['Successfully wrote %u new calculation results ' ...
            '(/ %u) to the Results table of %s in %s.\n'],...
            sum(updateMe),numWrite,dbname,BF_TheTime(toc(writeBackTimer)));

if any(~updateMe) % Some entries were not written to the database
    fprintf(1,['%u entries were not written (previously-calculated errors) and remain ' ...
                            'awaiting calculation in the database.\n'],sum(~updateMe));
end

SQL_CloseDatabase(dbc) % Close the database connection

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function planC = deriveProstateCTVfromPTV(ptvStructNum, planC)
% function deriveProstateCTVfromPTV(ptvStructNum, planC)
%
% APA, 08/24/2012

if ~exist('planC','var')
    global planC
end

global stateS

indexS = planC{end};

scanNum = getStructureAssociatedScan(ptvStructNum,planC);

gtvStructNum = length(planC{indexS.structures}) + 1;

gtvStructS = newCERRStructure(scanNum, planC, gtvStructNum);

rasterSegments = [];

for slcNum = 1:length(planC{indexS.scan}(scanNum).scanInfo)

    % Calculate centroid for this slice
    [rasterSegs, planC] = getRasterSegments(ptvStructNum,planC,slcNum);

    if ~isempty(rasterSegs)

        meanY = mean(rasterSegs(:,2));

        %rasterSegmentsPosterior = rasterSegs(rasterSegs(:,2) <= meanY,:);
        %rasterSegmentsAnterior  = rasterSegs(rasterSegs(:,2) > meanY,:);

        rasterSegmentsShrinked_10_6 = structMargin(rasterSegs, [10 6]/10, scanNum, planC);
        rasterSegmentsShrinked_10_6 = structDiff(rasterSegs,rasterSegmentsShrinked_10_6, scanNum, planC);

        rasterSegmentsShrinked_10 = structMargin(rasterSegs, 10/10, scanNum, planC);
        rasterSegmentsShrinked_10 = structDiff(rasterSegs,rasterSegmentsShrinked_10, scanNum, planC);

        patientPosition = planC{indexS.scan}.scanInfo(1).DICOMHeaders.PatientPosition;
        if isequal(patientPosition,'HFP') || isequal(patientPosition,'FFP')
            rasterSegmentsAnterior = rasterSegmentsShrinked_10(rasterSegmentsShrinked_10(:,2) <= meanY, :);
            rasterSegmentsPosterior = rasterSegmentsShrinked_10_6(rasterSegmentsShrinked_10_6(:,2) > meanY, :);
        else
            rasterSegmentsPosterior = rasterSegmentsShrinked_10_6(rasterSegmentsShrinked_10_6(:,2) <= meanY, :);
            rasterSegmentsAnterior = rasterSegmentsShrinked_10(rasterSegmentsShrinked_10(:,2) > meanY, :);
        end

        rasterSegments = [rasterSegments; rasterSegmentsAnterior; rasterSegmentsPosterior];

    end

end

gtvStructS.rasterSegments = rasterSegments;
gtvStructS.rasterized = 1;

contourS = rasterToPoly(rasterSegments, scanNum, planC);

gtvStructS.contour = contourS;

gtvStructS.structureName = 'PROSTATE_CTV';

planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, gtvStructS, gtvStructNum);

planC = updateStructureMatrices(planC, gtvStructNum);

% Refresh View
if ~isempty(stateS) && isfield(stateS,'handle') && isfield(stateS.handle,'CERRSliceViewer') && isnumeric(stateS.handle.CERRSliceViewer)
    stateS.structsChanged = 1;
    CERRRefresh
end

author: cerr

repo: CERR


text:
function t=convertLabelsToTable(L,n)
%convertLabelsToTable Converts labeled stack of images to the table of cluster-sizes.
%   t=convertLabelsToTable(L,n)
%   L is a labelled image, n number of clustrers - output of the function
%   [L,n] = bwlabelStack(bw,ncon);
%   t is the table (#frames X #clusters) containing sizes of the individual
%   clusters in each frame.
sizeZ=size(L,3);
maxClust = max(n);
t=zeros(sizeZ,maxClust);
for frame=1:sizeZ
    for cluster = 1:maxClust
        clustTmp=L(:,:,frame)==cluster;
        t(frame,cluster)=sum(clustTmp(:));
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
v=ver('matlab'); Release = v.Release;
switch Release
    case '(R2017a)'
        vivado = '2016.2';
    case '(R2017b)'
        vivado = '2017.4';
    case '(R2018b)'
        vivado = '2017.4';
end

if ispc
    hdlsetuptoolpath('ToolName', 'Xilinx Vivado', ...
        'ToolPath', ['C:\Xilinx\Vivado\',vivado,'\bin\vivado.bat']);
elseif isunix
    hdlsetuptoolpath('ToolName', 'Xilinx Vivado', ...
        'ToolPath', ['/opt/Xilinx/Vivado/',vivado,'/bin/vivado']);
end

author: analogdevicesinc


text:
function s = lp_e(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15,i16,...
             i17,i18,i19,i20)
%
%  Auxiliary routine for the execution of strings with MATLAB
%  expressions, which depend on numerical values.
%
%  Input:
%
%    i1,...,i20 strings or numerical values.
%
%  Output:
%
%    s          string.
%
%
%  LYAPACK 1.0 (Thilo Penzl, November 1999)

ni = nargin;

if ni>=1, if ischar(i1), s = i1; else, s = num2str(i1); end, end
if ni>=2, if ischar(i2), s = [s,i2]; else, s = [s,num2str(i2)]; end, end
if ni>=3, if ischar(i3), s = [s,i3]; else, s = [s,num2str(i3)]; end, end
if ni>=4, if ischar(i4), s = [s,i4]; else, s = [s,num2str(i4)]; end, end
if ni>=5, if ischar(i5), s = [s,i5]; else, s = [s,num2str(i5)]; end, end
if ni>=6, if ischar(i6), s = [s,i6]; else, s = [s,num2str(i6)]; end, end
if ni>=7, if ischar(i7), s = [s,i7]; else, s = [s,num2str(i7)]; end, end
if ni>=8, if ischar(i8), s = [s,i8]; else, s = [s,num2str(i8)]; end, end
if ni>=9, if ischar(i9), s = [s,i9]; else, s = [s,num2str(i9)]; end, end
if ni>=10, if ischar(i10), s = [s,i10]; else, s = [s,num2str(i10)]; end, end
if ni>=11, if ischar(i11), s = [s,i11]; else, s = [s,num2str(i11)]; end, end
if ni>=12, if ischar(i12), s = [s,i12]; else, s = [s,num2str(i12)]; end, end
if ni>=13, if ischar(i13), s = [s,i13]; else, s = [s,num2str(i13)]; end, end
if ni>=14, if ischar(i14), s = [s,i14]; else, s = [s,num2str(i14)]; end, end
if ni>=15, if ischar(i15), s = [s,i15]; else, s = [s,num2str(i15)]; end, end
if ni>=16, if ischar(i16), s = [s,i16]; else, s = [s,num2str(i16)]; end, end
if ni>=17, if ischar(i17), s = [s,i17]; else, s = [s,num2str(i17)]; end, end
if ni>=18, if ischar(i18), s = [s,i18]; else, s = [s,num2str(i18)]; end, end
if ni>=19, if ischar(i19), s = [s,i19]; else, s = [s,num2str(i19)]; end, end
if ni>=20, if ischar(i20), s = [s,i20]; else, s = [s,num2str(i20)]; end, end
if ni>=21, error('Too many input arguments.'); end