text:
function x = r8sp_cg ( n, nz_num, row, col, a, b, x )

%*****************************************************************************80
%
%% R8SP_CG uses the conjugate gradient method on an R8SP system.
%
%  Discussion:
%
%    The R8SP storage format stores the row, column and value of each nonzero
%    entry of a sparse matrix.
%
%    It is possible that a pair of indices (I,J) may occur more than
%    once.  Presumably, in this case, the intent is that the actual value
%    of A(I,J) is the sum of all such entries.  This is not a good thing
%    to do, but I seem to have come across this in MATLAB.
%
%    The R8SP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP
%    ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
%
%    The matrix A must be a positive definite symmetric matrix.
%
%    The method is designed to reach the solution after N computational
%    steps.  However, roundoff may introduce unacceptably large errors for
%    some problems.  In such a case, calling the routine again, using
%    the computed solution as the new starting estimate, should improve
%    the results.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Frank Beckman,
%    The Solution of Linear Equations by the Conjugate Gradient Method,
%    in Mathematical Methods for Digital Computers,
%    edited by John Ralston, Herbert Wilf,
%    Wiley, 1967,
%    ISBN: 0471706892,
%    LC: QA76.5.R3.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer NZ_NUM, the number of nonzero elements in
%    the matrix.
%
%    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and
%    column indices of the nonzero elements.
%
%    Input, real A(NZ_NUM), the nonzero elements of the matrix.
%
%    Input, real B(N), the right hand side vector.
%
%    Input/output, real X(N).
%    On input, an estimate for the solution, which may be 0.
%    On output, the approximate solution vector.
%
  b = b(:);
  x = x(:);
%
%  Initialize
%    AP = A * x,
%    R  = b - A * x,
%    P  = b - A * x.
%
  ap = r8sp_mv ( n, n, nz_num, row, col, a, x );

  r(1:n,1) = b(1:n,1) - ap(1:n,1);
  p(1:n,1) = b(1:n,1) - ap(1:n,1);
%
%  Do the N steps of the conjugate gradient method.
%
  for it = 1 : n
%
%  Compute the matrix*vector product AP=A*P.
%
    ap = r8sp_mv ( n, n, nz_num, row, col, a, p );
%
%  Compute the dot products
%    PAP = P*AP,
%    PR  = P*R
%  Set
%    ALPHA = PR / PAP.
%
    pap = p' * ap;
    pr = p' * r;

    if ( pap == 0.0 )
      return
    end

    alpha = pr / pap;
%
%  Set
%    X = X + ALPHA * P
%    R = R - ALPHA * AP.
%
    x(1:n,1) = x(1:n,1) + alpha * p(1:n,1);
    r(1:n,1) = r(1:n,1) - alpha * ap(1:n,1);
%
%  Compute the vector dot product
%    RAP = R*AP
%  Set
%    BETA = - RAP / PAP.
%
    rap = r' * ap;

    beta = - rap / pap;
%
%  Update the perturbation vector
%    P = R + BETA * P.
%
    p(1:n,1) = r(1:n,1) + beta * p(1:n,1);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%POINTS2D_DEMO  One-line description here, please.
%
%   output = points2d_demo(input)
%
%   Example
%   points2d_demo
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2020-01-08,    using Matlab 9.7.0.1247435 (R2019b) Update 2
% Copyright 2020 INRA - Cepia Software Platform.

%% generate data

rng(42);
pts = randn([100 2]) * 15 + 50;

figure; hold on; axis([0 100 0 100]);
drawPoint(pts, 'k+');
print(gcf, 'points2d.png', '-dpng');

centPts = centroid(pts);

%% bounding box

bbox = boundingBox(pts);

figure; hold on; axis([0 100 0 100]);
drawPoint(pts, 'k+');
drawBox(bbox, 'color', 'b', 'linewidth', 2);
print(gcf, 'points2d_bbox.png', '-dpng');

%% equivalent ellipse

elli = equivalentEllipse(pts);

figure; hold on; axis([0 100 0 100]);
drawPoint(pts, 'k+');
drawEllipse(elli, 'color', 'b', 'linewidth', 2);
print(gcf, 'points2d_ellipse.png', '-dpng');

%% equivalent ellipse

hull = convexHull(pts);

figure; hold on; axis([0 100 0 100]);
drawPoint(pts, 'k+');
drawPolygon(hull, 'color', 'b', 'linewidth', 2);
print(gcf, 'points2d_hull.png', '-dpng');

%% everything together

figure; hold on; axis([0 100 0 100]);
hp = drawPoint(pts, 'color', 'k', 'marker', 'o', 'linewidth', 2);
hc = drawPoint(centPts, 'color', 'b', 'marker', '*', 'linewidth', 2, 'MarkerSize', 10);
hb = drawBox(bbox, 'color', [0 0 .7], 'linewidth', 2);
he = drawEllipse(elli, 'color', [.7 0 0], 'linewidth', 2);
hh = drawPolygon(hull, 'color', [0 .7 0], 'linewidth', 2);
legend({'Points', 'Centroid', 'BoundingBox', 'Equiv. Ellipse', 'Conv. Hull'}, 'Location', 'NorthEast');
print(gcf, 'points2d_demo.png', '-dpng');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%DAGNN.UNIFORMSCALINGGRIDGENERATOR  Generate an iso-tropic scaling + translation
%   grid for bilinear resampling.
%   This layer maps 1 x 1 x 3 x N transformation parameters corresponding to
%   the scale and translation y and x respectively, to 2 x Ho x Wo x N
%   sampling grids compatible with dagnn.BlilinearSampler.

% (c) 2016 Ankush Gupta
classdef UniformScalingGridGenerator < dagnn.Layer

 properties
     Ho = 0;
     Wo = 0;
 end

  properties (Transient)
    % the grid --> this is cached
    % has the size: [2 x HoWo]
    xxyy ;
  end

  methods

    function outputs = forward(obj, inputs, ~)
      % input is a 1x1x3xN TENSOR corresponding to:
      % [  s 0 ty ]
      % [  0 s tx ]
      %
      % OUTPUT is a 2xHoxWoxN grid

      % reshape the tfm params into matrices:
      T = inputs{1};
      % check shape:
      sz_T = size(T);
      assert(all(sz_T(1:3) == [1 1 3]), 'transforms have incorrect shape');
      nbatch = size(T,4);
      S = reshape(T(1,1,1,:), 1,1,nbatch); % x,y scaling
      t = reshape(T(1,1,2:3,:), 2,1,nbatch); % translation
      % generate the grid coordinates:
      useGPU = isa(T, 'gpuArray');
      if isempty(obj.xxyy)
        obj.initGrid(useGPU);
      end
      % transform the grid:
      g = bsxfun(@times, obj.xxyy, S); % scale
      g = bsxfun(@plus, g, t); % translate
      g = reshape(g, 2,obj.Ho,obj.Wo,nbatch);
      outputs = {g};
    end

    function [derInputs, derParams] = backward(obj, inputs, ~, derOutputs)
      dY = derOutputs{1};
      useGPU = isa(dY, 'gpuArray');
      nbatch = size(dY,4);

      % create the gradient buffer:
      dA = zeros([1,1,3,nbatch], 'single');
      if useGPU, dA = gpuArray(dA); end

      dY  = reshape(dY, 2,obj.Ho*obj.Wo,nbatch);
      % gradient wrt the linear part:
      dA(1,1,1,:) = reshape(obj.xxyy,1,[]) * reshape(dY, [],nbatch);
      % gradient wrt translation (or bias):
      dA(1,1,2:3,:) = reshape(sum(dY,2),1,1,2,[]);

      derInputs = {dA};
      derParams = {};
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      nBatch = inputSizes{1}(4);
      outputSizes = {[2, obj.Ho, obj.Wo, nBatch]};
    end

    function obj = UniformScalingGridGenerator(varargin)
      obj.load(varargin);
      % get the output sizes:
      obj.Ho = obj.Ho;
      obj.Wo = obj.Wo;
      obj.xxyy = [];
    end

    function obj = reset(obj)
      reset@dagnn.Layer(obj) ;
      obj.xxyy = [] ;
    end

    function initGrid(obj, useGPU)
      % initialize the grid:
      % this is a constant
      xi = linspace(-1, 1, obj.Ho);
      yi = linspace(-1, 1, obj.Wo);
      [yy,xx] = meshgrid(yi,xi);
      xxyy = [xx(:), yy(:)]' ; % 2xM
      if useGPU
        xxyy = gpuArray(xxyy);
      end
      obj.xxyy = xxyy ;
    end

  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [cm3] = in32cm3(in3)
% Convert volume from cubic inches to cubic centimeters.
% Chad Greene 2012
cm3 = in3*16.387064;

author: Sable

repo: mcbench-benchmarks


text:
function roidb = roidb_from_ilsvrc14(imdb,varargin)
% roidb = roidb_from_ilsvrc14(imdb)
%   Builds an regions of interest database from imdb image
%   database. Uses precomputed selective search boxes available
%   in the R-CNN data package.
%
%   Inspired by Andrea Vedaldi's MKL imdb and roidb code.

ip = inputParser;
ip.addRequired('imdb', @isstruct);
ip.addParameter('exclude_difficult_samples',       false,   @islogical);
ip.addParameter('with_selective_search',           false,  @islogical);
ip.addParameter('with_edge_box',                   false,  @islogical);
ip.addParameter('with_self_proposal',              false,  @islogical);
ip.addParameter('rootDir',                         '.',    @ischar);
ip.addParameter('extension',                       '',     @ischar);
ip.addParameter('roidb_name_suffix',               '',     @isstr);
ip.addParameter('regions_file_sp',               '',     @isstr);

ip.parse(imdb, varargin{:});
opts = ip.Results;

cacheDir = fullfile(opts.rootDir, 'imdb' , 'cache' , 'ilsvrc');
if ~exist(cacheDir, 'dir')
    mkdir(cacheDir);
end

try
    flip = imdb.flip;
catch
    flip = false;
end

if isempty(opts.roidb_name_suffix)
   opts.roidb_name_suffix = '';
end
if flip == false
    cache_file = [cacheDir '/roidb_' imdb.name opts.roidb_name_suffix '_unflip.mat'];
else
    cache_file = [cacheDir '/roidb_' imdb.name opts.roidb_name_suffix '_flip.mat'];
end

%%
try
   load(cache_file);
catch
    addpath(fullfile(imdb.details.devkit_path, 'evaluation'));

    roidb.name = imdb.name;
    % wsh  regions_file = fullfile('data', 'selective_search_data', [roidb.name '.mat']);
    regions = [];
    if opts.with_selective_search
        fprintf('Loading SS region proposals...');
        regions = load_proposals(regions_file_ss, regions);
        fprintf('done\n');
    end
    if opts.with_edge_box
        regions = load_proposals(regions_file_eb, regions);
    end
    if opts.with_self_proposal
        regions = load_proposals(opts.regions_file_sp, regions);
        regions = renameStructField(regions,'roi', 'boxes');
        regions = renameStructField(regions,'list', 'images');

    end

    if isempty(regions)
        fprintf('Warrning: no ADDITIONAL windows proposal is loaded!\n');
        regions.boxes = cell(length(imdb.image_ids), 1);
        if flip
            regions.images = imdb.image_ids(1:2:end);
        else
            regions.images = imdb.image_ids;
        end
    end

    hash = make_hash(imdb.details.meta_det.synsets);
    WNIDs = {imdb.details.meta_det.synsets.WNID};
    if ~flip

        for i = 1:length(imdb.image_ids)
            if ~mod(i,100)
              tic_toc_print('roidb (%s): %d/%d\n', roidb.name, i, length(imdb.image_ids));
            end
            try
                anno_file = fullfile(imdb.details.bbox_path, [imdb.image_ids{i} '.xml']);
                if strfind(anno_file, '_flip'),
                  flipthat = true;
                  anno_file = strrep(anno_file, '_flip', '');
                else
                  flipthat = false;
                end
                if ~ispc,         anno_file(anno_file == '\') = '/'; end
                if strcmp(imdb.name(1:12),'ilsvrc15_vid')
                  rec = VOCreadrecTrackxml(anno_file, hash);
                else
                  rec = VOCreadrecxml(anno_file, hash);
                end
            catch
                warning('GT(xml) file empty/broken: %s\n', imdb.image_ids{i});

                rec = [];
            end
          if ~isempty(regions)
            if numel(regions.images) ~= numel(imdb.image_ids) %smaller imdb due to subsampling -> find right regions
              j = find(strcmp(regions.images,[imdb.image_ids{i} '.' imdb.extension])) ;
            else
              j = i;
            end
                [~, image_name1] = fileparts(imdb.image_ids{i});
                [~, image_name2] = fileparts(regions.images{j});
                assert(strcmp(image_name1, image_name2));
          end
            roidb.rois(i) = attach_proposals(rec, regions.boxes{j}, WNIDs, opts.exclude_difficult_samples, flipthat);
        end

    else
        % flip case
        for i = 1:length(imdb.image_ids)/2

            tic_toc_print('roidb (%s): %d/%d\n', roidb.name, i, length(imdb.image_ids)/2);
            try
                anno_file = fullfile(imdb.details.bbox_path, [imdb.image_ids{2*i-1} '.xml']);
                if strcmp(imdb.name(1:12),'ilsvrc15_vid')
                  rec = VOCreadrecTrackxml(anno_file, hash);
                else
                  rec = VOCreadrecxml(anno_file, hash);
                end
            catch
                warning('GT(xml) file empty/broken: %s\n', imdb.image_ids{2*i-1});
                rec = [];
            end
            if ~isempty(regions)
                [~, image_name1] = fileparts(imdb.image_ids{i*2-1});
                [~, image_name2] = fileparts(regions.images{i});
                assert(strcmp(image_name1, image_name2));
            end
            roidb.rois(i*2-1) = attach_proposals(rec, regions.boxes{i*2-1}, WNIDs, opts.exclude_difficult_samples, false);
            roidb.rois(i*2) = attach_proposals(rec, regions.boxes{i*2}, WNIDs, opts.exclude_difficult_samples, true);
        end
    end

    rmpath(fullfile(imdb.details.devkit_path, 'evaluation'));

    fprintf('Saving roidb to cache...');
    save(cache_file, 'roidb', '-v7.3');
    fprintf('done\n');
end

% ------------------------------------------------------------------------
function rec = attach_proposals(ilsvrc_rec, boxes, WNIDs, exclude_difficult_samples, flip)

if size(boxes,2) > 4
  boxes = boxes(:,1:4);
end

if isfield(ilsvrc_rec, 'objects') && ~isempty(ilsvrc_rec.objects)
    if exclude_difficult_samples
        valid_objects = ~cat(1, ilsvrc_rec.objects(:).difficult);
    else
        valid_objects = 1:length(ilsvrc_rec.objects(:));
    end
    gt_boxes = cat(1, ilsvrc_rec.objects(valid_objects).bbox);

    %%% ============ NOTE ==============
    % coordinate starts from 0 in ilsvrc
    gt_boxes = gt_boxes + 1;

    if flip
        gt_boxes(:, [1, 3]) = ilsvrc_rec.imgsize(1) + 1 - gt_boxes(:, [3, 1]);
    end
    all_boxes = cat(1, gt_boxes, boxes);
    gt_classes = zeros(length(valid_objects),1);
    for i = 1:length(valid_objects)
      gt_classes(i) = find(strcmp(WNIDs,ilsvrc_rec.objects(valid_objects(i)).class));
    end
    num_gt_boxes = size(gt_boxes, 1);

else
    gt_boxes = [];
    all_boxes = boxes;
    gt_classes = [];
    num_gt_boxes = 0;
end

num_boxes = size(boxes, 1);
rec.boxes = single(all_boxes);
rec.feat = [];
rec.class = uint8(cat(1, gt_classes, zeros(num_boxes, 1)));
if ~isempty(ilsvrc_rec.objects) && isfield(ilsvrc_rec.objects,'trackid')
  rec.trackids = uint8(cat(1, ilsvrc_rec.objects(:).trackid));
  rec.occluded = logical(cat(1, ilsvrc_rec.objects(:).occluded));
  rec.generated = logical(cat(1, ilsvrc_rec.objects(:).generated));
else
  rec.trackids = [];
  rec.occluded = [];
  rec.generated =[];
end
rec.gt = cat(1, true(num_gt_boxes, 1), false(num_boxes, 1));
rec.overlap = zeros(num_gt_boxes+num_boxes, length(WNIDs), 'single');
for i = 1:num_gt_boxes
    rec.overlap(:, gt_classes(i)) = ...
        max(rec.overlap(:, gt_classes(i)), boxoverlap(all_boxes, gt_boxes(i, :)));
end

% ------------------------------------------------------------------------
function regions = load_proposals(proposal_file, regions)
% ------------------------------------------------------------------------
if isempty(regions)
    regions = load(proposal_file);
else
    regions_more = load(proposal_file);
    if ~all(cellfun(@(x, y) strcmp(x, y), regions.images(:), regions_more.images(:), 'UniformOutput', true))
        error('roidb_from_ilsvrc: %s is has different images list with other proposals.\n', proposal_file);
    end
    regions.boxes = cellfun(@(x, y) [double(x); double(y)], regions.boxes(:), regions_more.boxes(:), 'UniformOutput', false);
end

function rec = VOCreadrecTrackxml(path,hash)

x=VOCreadxml(path);
x=x.annotation;

rec.folder=x.folder;
rec.filename=x.filename;
rec.source.database=x.source.database;

rec.size.width=str2double(x.size.width);
rec.size.height=str2double(x.size.height);

rec.imgname=[x.folder x.filename];
rec.imgsize=str2double({x.size.width x.size.height});
rec.database=rec.source.database;

if isfield(x,'object')
    for i=1:length(x.object)
        rec.objects(i)=xmlobjtopas(x.object(i),hash);
    end
else
    rec.objects = [];
end

function p = xmlobjtopas(o,hash)

p.class=o.name;

p.label= get_class2node( hash, p.class );

p.bbox=str2double({o.bndbox.xmin o.bndbox.ymin o.bndbox.xmax o.bndbox.ymax});

p.trackid = str2double(o.trackid);
p.occluded = str2double(o.occluded);
p.generated = str2double(o.generated);

author: feichtenhofer

repo: Detect-Track


text:
function   y = sigmoid(x)
% legacy code
s = warning ('on');
warning ('*** The function `sigmoid` is now deprecated. Please see `VBA_sigmoid` for an alternative.')
warning (s);

% fallback
y = VBA_sigmoid (x);

author: MBB-team

repo: VBA-toolbox


text:
function y = horzcat( varargin )
y = cat( 2, varargin{:} );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [x_nml_wave,UB] = normalize_UB( H_wave,y_wave,x,N,l,u)
%UNTITLED9 Summary of this function goes here
%   Detailed explanation goes here
x_cplx = x(1:N)+j*x(N+1:2*N);
x_temp = x_cplx./abs(x_cplx);
for ii = 1:N
    delta = u(ii)-l(ii);
    if delta <=pi
        x_nml(ii,1) = x_temp(ii);
    else
        if ((angle(x_temp(ii))<=u(ii))&&(angle(x_temp(ii))>=l(ii)))
            x_nml(ii,1) = x_temp(ii);
        else
            d_l = abs(x_cplx(ii)-exp(j*l(ii)));
            d_u = abs(x_cplx(ii)-exp(j*u(ii)));
           if d_l <= d_u
                x_nml(ii,1) = exp(j*l(ii));
           else
                x_nml(ii,1) = exp(j*u(ii));
           end
        end
    end
end
x_nml_wave = [real(x_nml);imag(x_nml)];
UB = objval_func(x_nml_wave,H_wave,y_wave);
end

author: yuanhao-cui

repo: Must-Reading-on-ISAC

path: github-repos/MATLAB/yuanhao-cui-Must-Reading-on-ISAC/Must-Reading-on-ISAC-34cd6615c52ebca121428a979e756c608b195040/Codes/Fan2018TSP/Codes for DFRC Waveform Design/Constant Modulus/normalize_UB.m


text:
function quad = piecewise_linear_product_quad ( a, b, f_num, f_x, f_v, ...
  g_num, g_x, g_v, quad_num )

%*****************************************************************************80
%
%% PIECEWISE_LINEAR_PRODUCT_QUAD: estimate piecewise linear product integral.
%
%  Discussion:
%
%    We are given two piecewise linear functions F(X) and G(X) and we wish
%    to estimate the value of the integral
%
%      INTEGRAL = Integral ( A <= X <= B ) F(X) * G(X) dx
%
%    The functions F(X) and G(X) are defined as tables of coordinates X and
%    values V.  A piecewise linear function is evaluated at a point X by
%    evaluating the interpolant to the data at the endpoints of the interval
%    containing X.
%
%    It must be the case that A <= B.
%
%    It must be the case that the node coordinates F_X(*) and G_X(*) are
%    given in ascending order.
%
%    It must be the case that:
%
%      F_X(1) <= A and B <= F_X(F_NUM)
%      G_X(1) <= A and B <= G_X(G_NUM)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the limits of integration.
%
%    Input, integer F_NUM, the number of nodes for F.
%
%    Input, real F_X(F_NUM), the node coordinates for F.
%
%    Input, real F_V(F_NUM), the nodal values for F.
%
%    Input, integer G_NUM, the number of nodes for G.
%
%    Input, real G_X(G_NUM), the node coordinates for G.
%
%    Input, real G_V(G_NUM), the nodal values for G.
%
%    Input, integer QUAD_NUM, the number of quadrature points.
%
%    Output, real QUAD, an estimate for the integral of F(X) * G(X)
%    from A to B.
%
  quad = 0.0;

  f_left = 1;
  g_left = 1;

  a2 = a;
  a2 = max ( a2, f_x(1) );
  a2 = max ( a2, g_x(1) );

  b2 = b;
  b2 = min ( b2, f_x(f_num) );
  b2 = min ( b2, g_x(g_num) );

  for i = 1 : quad_num

    xq =  ( (                2 * i - 1 ) * b2 ...
          + ( 2 * quad_num - 2 * i + 1 ) * a2 )  ...
          / ( 2 * quad_num             );

    f_left = r8vec_bracket3 ( f_num, f_x, xq, f_left );

    fq = f_v(f_left) + ( xq - f_x(f_left) ) * ( f_v(f_left+1) - f_v(f_left) ) ...
      / ( f_x(f_left+1) - f_x(f_left) );

    g_left = r8vec_bracket3 ( g_num, g_x, xq, g_left );

    gq = g_v(g_left) + ( xq - g_x(g_left) ) * ( g_v(g_left+1) - g_v(g_left) ) ...
      / ( g_x(g_left+1) - g_x(g_left) );

    quad = quad + fq * gq;

  end

  quad = quad * ( b - a ) / quad_num;

  return
end


text:
function value = p11_i4 ( action, name, value )

%*****************************************************************************80
%
%% P11_I4 sets or gets I4 parameters for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P11_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P11_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P11_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P11_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P11_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P11_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P11_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P11_I4 - Fatal error!' );

  end

  return
end


text:
filename='Cantilever_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function f = minus(f, g)
%-   CHEBFUN minus.
%   F - G subtracts CHEBFUNs F and G, or a scalar from a CHEBFUN if either F or
%   G is a scalar.
%
%   H = MINUS(F, G) is called for the syntax 'F - G'.

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

f = plus(f, uminus(g));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% minimum distance from nose of fly to ellipse of any other fly within 30
% degrees
function [data,units] = compute_dnose2ell_anglerange(trx,n,anglerange_str,varargin)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);

% name of closestfly stat
fn = sprintf('closestfly_nose2ell_angle_%sto%s',anglerange_str{:});

% access closestfly to ensure that dnose2ell_anglerange is computed
for i = 1:numel(flies),
  trx(flies(i)).(fn);
end
%compute_closestfly_nose2ell_anglerange(trx,n,anglerange,varargin{:});

% anglerange_deg = round(anglerange*180/pi);
% fn = sprintf('dnose2ell_angle_%dto%d',anglerange_deg(1),anglerange_deg(2));
% fn = strrep(fn,'-','min');

% name of dnose2ell stat
fn = sprintf('dnose2ell_angle_%sto%s',anglerange_str{:});

for i1 = 1:nflies,
  fly1 = flies(i1);
  data{i1} = trx(fly1).(fn);
end
units = parseunits('mm');

author: kristinbranson

repo: JAABA


text:
% OtStop Stop using the Omnitrans interface.
%
% Usage:
%    OtStop
%
% Releases the Omnitrans resources allocated with a call to OtStart().

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * 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
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% 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 OWNER 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.

% 2006-02-28 Created (Johan Meijdam, VORtech)

function OtStop()
   Ot('stop');
end


text:
function [wordS, start, stop] = word(stringS,n)
%function [wordS start stop] = word(stringS,n)
%Returns the n'th blank delimited word in string of words. n
%must be positive. If there are fewer than n words in the string,
%the null string is returned.  The beginning and ending indices are
%given by 'start' and 'stop'.
%Inspired by the REXX function 'word'.
%
%ex:
%   c=word('a and b',3)
%
%LM: 7 March 02, JOD.
%
% 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/>.

%Get the number of blank-delimited words:
num_words=words(stringS);

len = length(stringS);

if num_words<n
   wordS='';
   start = '';
   stop  = '';
else
   index=1;
   indices_begin=[];
   indices_end=[];
   beginning=0;
   ending=1;
   wordNum = 0;
   while index<= len
      char=stringS(index);
      if ~isspace(char) && beginning==0
        indices_begin=[indices_begin,index];
        ending=0;
        beginning=1;
      end
      if isspace(char) && ending==0
        indices_end=[indices_end,index-1];
        beginning=0;
        ending=1;
        wordNum = wordNum + 1;
        if wordNum == n
          index = len + 1; %finish
        end
      end
      index=index+1;
   end
   %special case: if the last character was not a space,
   %then we need to count it as ending a word:
   pos = len;
   if ~isspace(stringS(pos))
     indices_end=[indices_end,pos];
   end
   wordS=stringS(indices_begin(n):indices_end(n));
   start = indices_begin(n);
   stop  = indices_end(n);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [] = plx_ad_info(filename)
% plx_ad_info(filename): prints .plx file a/d channel info
%
% plx_ad_info(filename)
%
% INPUT:
%   filename - .plx file name. Will use file open dialog if filename is empty string
%
% OUTPUT:
%   (none)

if nargin ~= 1
    error 'expected 1 input argument';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    pl2 = PL2GetFileIndex(filename);
    PL2Print(pl2.AnalogChannels);
    return
end

[n, names] = plx_adchan_names(filename);
[n, raw] = plx_ad_chanmap(filename);
[n, samplecounts] = plx_adchan_samplecounts(filename);
[n, freqs] = plx_adchan_freqs(filename);
fprintf(' %18.18s', 'RawChannelNum');
fprintf(' %18.18s', 'Name');
fprintf(' %18.18s', 'SampleCount');
fprintf(' %18.18s', 'SamplesPerSecond');
fprintf('\n');

for i=1:n
    fprintf(' %18d', raw(i));
    fprintf(' %18.18s', names(i,:));
    fprintf(' %18d', samplecounts(i));
    fprintf(' %18.10g', freqs(i));
    fprintf('\n');
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function data_obj = import_onsets(data_obj, design_file)
% Import onsets from a csv or Excel file and attach them to the fmri_glm_design_matrix object

            % DesginmatrixTable INPUT
            % TRY readtable() -- may be more flexible than xlsread.
            xlsread(DesginmatrixTable)

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [U,data,SS,R] = arap(V,F,b,bc,varargin)
  % ARAP Solve for the as-rigid-as-possible deformation according to various
  % manifestations including:
  %   (1) "As-rigid-as-possible Surface Modeling" by [Sorkine and Alexa 2007]
  %   (2) "A local-global approach to mesh parameterization" by [Liu et al.
  %     2010] or "A simple geometric model for elastic deformation" by [Chao et
  %     al.  2010]
  %   (3) Adapted version of "As-rigid-as-possible Surface Modeling" by
  %     [Sorkine and Alexa 2007] presented in section 4.2 of or "A simple
  %     geometric model for elastic deformation" by [Chao et al.  2010]
  %
  % U = arap(V,F,b,bc) given a rest mesh (V,F) and list of constraint vertex
  % indices (b) and their new postions (bc) solve for pose mesh positions (U),
  % using default choice for 'Energy'
  %
  % U = arap(V,F,b,bc,'ParameterName','ParameterValue',...)
  %
  % Inputs:
  %   V  #V by dim list of rest domain positions
  %   F  #F by {3|4} list of {triangle|tetrahedra} indices into V
  %   b  #b list of indices of constraint (boundary) vertices
  %   bc  #b by dim list of constraint positions for b
  %   Optional:
  %     'Energy'
  %       followed by a string specifying which arap energy definition to use.
  %       One of the following:
  %         'spokes'  "As-rigid-as-possible Surface Modeling" by [Sorkine and
  %           Alexa 2007], rotations defined at vertices affecting incident
  %           edges
  %         'elements'  "A local-global approach to mesh parameterization" by
  %           [Liu et al.  2010] or "A simple geometric model for elastic
  %           deformation" by [Chao et al.  2010], rotations defined at
  %           elements (triangles or tets)
  %         'spokes-and-rims'  Adapted version of "As-rigid-as-possible Surface
  %           Modeling" by [Sorkine and Alexa 2007] presented in section 4.2 of
  %           or "A simple geometric model for elastic deformation" by [Chao et
  %           al.  2010], rotations defined at vertices affecting incident
  %           edges and opposite edges
  %     'V0' #V by dim list of initial guess positions
  %       dim by dim by #C list of linear transformations initial guesses,
  %       optional (default is to use identity transformations)
  %     'Data' see output
  %     'Groups'
  %       followed by #V list of group indices (1 to k) for each vertex, such
  %       that vertex i is assigned to group G(i)
  %     'Tol'
  %       stopping critera parameter. If variables (linear transformation matrix
  %       entries) change by less than 'tol' the optimization terminates,
  %       default is 0.75 (weak tolerance)
  %     'MaxIter'
  %       max number of local-global iterations, default is 100
  %     'Dynamic'
  %        #V by dim list of external forces
  %     'TimeStep'
  %        scalar time step value
  %     'Vm1'
  %        #V by dim positions at time t-1
  %     'Tikhonov' followed by constant Tikhonov regularization parameter
  %       alpha:
  %       http://en.wikipedia.org/wiki/Tikhonov_regularization#Relation_to_probabilistic_formulation
  %     'Youngs'  followed by young's modulus (only for dynamics)
  %     'Flat' followed by whether to add the constraint that Z=0 {false}
  %     'RemoveRigid' followed by whether to add a constraint that places an
  %       arbitrary point at the origin and another along the x-axis {false}
  %     'Aeq'/'Beq'  followed by #Aeq by dim*#V linear equality constraints
  %        matrix and righthand sides respectively {[]}/{[]}
  % Outputs:
  %   U  #V by dim list of new positions
  %   data  struct of reusable data
  %     .CSM dim*n by dim*n sparse matrix containing special laplacians along the
  %       diagonal so that when multiplied by repmat(U,dim,1) gives covariance
  %       matrix elements, can be used to speed up next time this function is
  %       called, see function definitions
  %
  % Known issues: 'Flat',true + 'Energy','elements' should only need a 2D
  % rotation fit, but this does 3D to stay general (e.g. if one were to use
  % 'spokes' then the edge-set cannot be pre-mapped to a common plane)
  %
  % See also: takeo_arap
  %

  % parse input
  G = [];

  % number of vertices
  n = size(V,1);
  assert(isempty(b) || max(b) <= n,'b should be in bounds [1,size(V,1)]');
  assert(isempty(b) || min(b) >= 1,'b should be in bounds [1,size(V,1)]');

  indices = 1:n;
  max_iterations = 100;
  tol = 0.001;
  interior = indices(~ismember(indices,b));
  U = [];
  Vm1 = [];
  youngs = 1.2e9;
  % default is Sorkine and Alexa style local rigidity energy
  energy = [];
  % default is no external forces
  fext = [];
  % defaults is unit time step
  h = 1;
  % Tikhonov regularization alpha
  alpha_tik = 0;
  % flatten/parameterization
  flat = false;
  % remove rigid transformation invariance
  remove_rigid = false;
  G = [];
  debug = false;
  data = [];
  Aeq = [];
  Beq = [];
  ref_data = [];
  weight = [];

  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Energy','V0','Data','Tikhonov','Groups','Tol', ...
     'MaxIter','Dynamic','TimeStep','Vm1','RemoveRigid','Debug', ...
     'Aeq','Beq','Flat','Ref','Weight','Youngs'}, ...
    {'energy','U','data','alpha_tik','G','tol','max_iterations','fext', ...
    'h','Vm1','remove_rigid','debug','Aeq','Beq','flat','ref_data','weight','youngs'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if isempty(energy)
    switch size(F,2)
    case 4
      energy = 'elements';
    case 3
      energy = 'spokes-and-rims';
    end
  end
  if isempty(fext)
    dynamic = false;
    fext = zeros(size(V));
  else
    dynamic = true;
  end

  if isempty(data)
    if strcmp(energy,'elements') && numel(G) ~= size(F,1) && numel(G) == n
      % groups are defined per vertex, convert to per face using mode
      data.G = mode(G(F),2);
    else
      data.G = G;
    end
  end

  if isempty(data.G)
    k = n;
  else
    k = max(data.G);
  end

  ref_map = [];
  if isempty(ref_data)
    if flat
      [ref_V,ref_F,ref_map] = plane_project(V,F);
      assert(strcmp(energy,'elements'),'flat only makes sense with elements');
    else
      ref_V = V;
      ref_F = F;
    end
  else
    ref_V = ref_data{1};
    ref_F = ref_data{2};
    ref_map = ref_data{3};
  end
  dim = size(ref_V,2);

  if isempty(bc)
    bc = sparse(0,dim);
  end

  assert(dim == size(bc,2));
  assert(size(Aeq,1) == size(Beq,1));
  assert((size(Aeq,2) == 0) || (size(Aeq,2) == dim*size(V,1)));

  if isempty(U)
    if(dim == 2)
      U = laplacian_mesh_editing(V,F,b,bc);
    else
      U = V;
    end
    U = U(:,1:dim);
  end
  assert(n == size(U,1));
  assert(dim == size(U,2));

  if dynamic
    V0 = U(:,1:dim);
    if isempty(Vm1)
      Vm1 = V0;
    end
    %% This parameter seem good for a 0.1m tall hotdog with h=0.0333
    % Larger --> more rigid
    dyn_alpha = h^2*youngs;
    % Smaller --> more damped
    mom = 1e10;
    M = massmatrix(ref_V,ref_F);
    if ~isempty(ref_map)
      M = ref_map*M*ref_map';
    end
    DQ = 0.5*mom*M;
    Dl =     mom*M*(-2*V0 + Vm1) - h^2*fext;
  else
    dyn_alpha = 1;
    DQ = sparse(size(V,1),size(V,1));
    Dl = sparse(size(V,1),dim);
  end

  if ~isfield(data,'L') || isempty(data.L)
    if isempty(weight)
      data.L = cotmatrix(ref_V,ref_F);
    else
      wA = doublearea(ref_V,ref_F);
      wG = grad(ref_V,ref_F);
      data.L = -0.5*wG'*( repmat(wA.*weight,dim,1) .* wG);
    end

    if ~isempty(ref_map)
      data.L = ref_map*data.L*ref_map';
    end
  end

  rr.b = cell(dim,1);
  rr.bc = cell(dim,1);
  if ~isfield(data,'rr')
    data.rr = [];
    if ~isfield(data.rr,'preF')
      data.rr.preF = cell(dim,1);
    end
  end
  if ~isfield(data,'preF')
    data.preF = [];
  end
  if remove_rigid
    if ~isempty(b)
      warning('RemoveRigid`s constraints are not typically wanted if |b|>0');
    end
    % the only danger is picking two points which end up mapped very close to
    % each other
    [~,f] = farthest_points(V,dim);
    for c = 1:dim
      rr.b{c} = f([1:dim-(c-1)])';
      rr.bc{c} = zeros(numel(rr.b{c}),1);
    end
    % Immediately remove rigid transformation from initial guess
    U = bsxfun(@minus,U,U(f(1),:));
    % We know that f(1) is at origin
    switch dim
    case 3
      % rotate about y-axis so that f(2) has (x=0)
      theta = atan2(U(f(2),1),U(f(2),3));
      R = axisangle2matrix([0 1 0],theta);
      U = U*R;
      % rotate about x-axis so that f(2) has (y=0)
      theta = atan2(U(f(2),3),U(f(2),2))+pi/2;
      R = axisangle2matrix([1 0 0],theta);
      U = U*R;
      % rotate about z-axis so that f(3) has (x=0)
      theta = atan2(U(f(3),2),U(f(3),1))+pi/2;
      R = axisangle2matrix([0 0 1],theta);
      U = U*R;
    case 2
      % rotate so that f(2) is on y-axis (x=0)
      theta = atan2(U(f(2),2),U(f(2),1))+pi/2;
      R = [cos(theta) -sin(theta);sin(theta) cos(theta)];
      U = U*R;
    otherwise
      error('Unsupported dimension');
    end
  end

  all = [interior(:);b(:)]';

  % cholesky factorization
  %cholL = chol(-L(interior,interior),'lower');
  % Why lu and not cholesky?
  %[luL,luU,luP,luQ,luR] = lu(-L(interior,interior));

  %R = repmat(eye(dim,dim),[1 1 n]);

  if ~isfield(data,'ae') || isempty(data.ae)
    data.ae = avgedge(V,F);
  end

  % build covariance scatter matrix used to build covariance matrices we'll
  % later fit rotations to
  if ~isfield(data,'CSM') || isempty(data.CSM)
    %assert(size(ref_V,2) == dim);
    data.CSM = ...
      covariance_scatter_matrix(ref_V,ref_F,'Energy',energy);
    if ~isempty(ref_map)
      data.CSM = data.CSM * repdiag(ref_map',dim);
    end

    % if there are groups then condense scatter matrix to only build
    % covariance matrices for each group
    if ~isempty(G)
      G_sum = group_sum_matrix(data.G,k);
      %CSM = [G_sum sparse(k,n); sparse(k,n) G_sum] * CSM;
      data.CSM = repdiag(G_sum,dim) * data.CSM;
    end
  end

  % precompute rhs premultiplier
  if ~isfield(data,'K') || isempty(data.K)
    [~,data.K] = arap_rhs(ref_V,ref_F,[],'Energy',energy,'Weight',weight);
    if ~isempty(ref_map)
      data.K = repdiag(ref_map,dim) * data.K;
    end
  end

  % initialize rotations with identies (not necessary)
  R = repmat(eye(dim,dim),[1 1 size(data.CSM,1)/dim]);

  iteration = 1;
  U_prev = V;
  data.energy = inf;
  while true

    if iteration > max_iterations
      if debug
        fprintf('arap: Iter (%d) > max_iterations (%d)\n',iteration,max_iterations);
      end
      break;
    end

    if iteration > 1
      change = max(abs(U(:)-U_prev(:)));
      if debug
        fprintf('arap: iter: %d, change: %g, energy: %g\n', ...
          iteration,change/data.ae,data.energy);
      end
      if change/data.ae <tol
        if debug
          fprintf('arap: change (%g) < tol*ae (%g * %g)\n',change,tol,data.ae);
        end
        break;
      end
    end

    U_prev = U;

    % energy after last global step
    U(b,:) = bc;
    %E = arap_energy(V,F,U,R);
    %[1 E]

    % compute covariance matrix elements
    S = zeros(size(data.CSM,1),dim);
    S(:,1:dim) = data.CSM*repmat(U,dim,1);
    % dim by dim by n list of covariance matrices
    SS = permute(reshape(S,[size(data.CSM,1)/dim dim dim]),[2 3 1]);
    % fit rotations to each deformed vertex
    R = fit_rotations(SS,'SinglePrecision',false);

    % energy after last local step
    U(b,:) = bc;
    %E = arap_energy(V,F,U,R);
    %[2 E]

    % This is still the SLOW way of building the right hand side, the entire
    % step of building the right hand side should collapse into a
    % #handles*dim*dim+1 by #groups matrix times the rotations at for group

    % distribute group rotations to vertices in each group
    if ~isempty(data.G)
      R = R(:,:,data.G);
    end

    U(b,:) = bc;

    %B = arap_rhs(V,F,R);
    Rcol = reshape(permute(R,[3 1 2]),size(data.K,2),1);
    Bcol = data.K * Rcol;
    B = reshape(Bcol,[size(Bcol,1)/dim dim]);

    zQ = -0.5*dyn_alpha*data.L+DQ+alpha_tik*speye(size(data.L));
    zL = -dyn_alpha*B+Dl;
    % RemoveRigid requires to solve each indepently
    if remove_rigid
      assert(isempty(Aeq),'Linear equality constraints not supported')
      for c = 1:dim
        eff_b = [b rr.b{c}];
        eff_bc = [bc(:,c);rr.bc{c}];
        [U(:,c),data.rr.preF{c}] = min_quad_with_fixed( ...
          zQ, ...
          zL(:,c),eff_b,eff_bc,[],[],data.rr.preF{c});
      end
    else
      if isempty(Aeq)
        cellfun(@(v) assignin('base',v,evalin('caller',v)),{'zQ','zL'});
        [U,data.preF] = min_quad_with_fixed( ...
          zQ, ...
          zL,b,bc,Aeq,Beq,data.preF);
      else
        assert(dim == size(bc,2));
        assert(dim == size(B,2));
        % repeat diagonal to solve for all coordinates together
        mQ = repdiag(zQ,dim);
        mL = reshape(zL,[],1);
        mb = reshape(bsxfun(@plus,b(:),size(V,1)*(0:dim-1)),[],1)';
        mbc = bc(:);
        [U,data.preF] = min_quad_with_fixed(mQ,mL,mb,mbc,Aeq,Beq,data.preF);
        U = reshape(U,[],dim);
      end
    end
    energy_prev = data.energy;
    data.energy = trace(U'*(zQ)*U+U'*(zL))+trace(V'*(0.5*dyn_alpha*data.L*V));
    if data.energy > energy_prev
      if debug
        fprintf('arap: energy (%g) increasing (over %g) (iter %d)\n', ...
          data.energy,energy_prev,iteration);
      end
      break;
    end

    %U(interior,:) = -L(interior,interior) \ (B(interior,:) + L(interior,b)*bc);
    %U(interior,:) = -L(interior,all)*L(all,interior) \ (L(interior,all)*B(all,:) + L(interior,all)*L(all,b)*bc);
    %U(interior,:) = luQ*(luU\(luL\(luP*(luR\(B(interior,:)+L(interior,b)*bc)))));
    %U(interior,:)=cholL\((B(interior,:)+L(interior,b)*bc)'/cholL)';
    iteration = iteration + 1;
  end
  U = U(:,1:dim);
end

author: alecjacobson

repo: gptoolbox


text:
function data = create_synthetic_dataset(data)
% create_synthetic_dataset creates test data for running nldr algorithms.
%
% inputs:
%    data      a struct describing the test data
%              .dataset the number of the example, see code for more infos
%              .n       the number of data points (default=400)
%              .state   the initial state for the random numbers (default=0)
%              .noise   the variance of Gaussian noise to add (default=0)
%              other options for some of the data sets (see code)
%              alternatively, data = 1 chooses the dataset directly,
%              the number of points defaults to 1000
%
% outputs:
%    data      a struct containing .x the generated data, each column is
%              a data point, and other stuff:
%              .z     the "correct" embedding
%              .e     some random noise of same dimensionality
%              .x_noisefree  the noisefree version of .x, i.e.
%                     .x = .xnoise_free + sqrt(.noise) * .e
%
% Adapted from create.m, originally written by
% (c) Stefan Harmeling, 2006
% using the examples of the original LLE and ISOMAP code.

if ~isfield(data, 'dataset'),
  number = data;
  clear data
  data.dataset = number;
end
if ~isfield(data, 'n'), data.n = 400; end
if ~isfield(data, 'noise'), data.noise = 0.0; end
if ~isfield(data, 'state'), data.state = 0; end

% set the randomness
rand('state', data.state);
randn('state', data.state);

data.typ = 'data';
switch data.dataset
 case 0 % "swiss roll with hole"
  data.name = 'swiss roll with hole';
  n = data.n;
  a = 1;   % swiss roll goes from a*pi to b*pi
  b = 4;
  y = rand(2,n);
  % punch a rectangular hole at the center
  l1 = 0.05; l2 = 0.15;
  y = y - 0.5;
  ok = find((abs(y(1,:))>l1) | (abs(y(2,:))>l2));
  i = length(ok);
  y(:, 1:i) = y(:, ok);
  while (i<n)
    p = rand(2,1) - 0.5;
    if (abs(p(1))>l1) || (abs(p(2))>l2)
      i = i + 1;
      y(:,i) = p;
    end
  end
  y = y + 0.5;
  tt = (b-a)*y(1,:) + a;
  tt = pi*tt;
  height = 21*y(2,:);
  data.col = tt;
  data.x = [tt.*cos(tt); height; tt.*sin(tt)];
  data.z = [tt; height]; % the ground truth
  data.az = -4;
  data.el = 13;

 case -1 % "swiss roll" dataset extracted from LLE's swissroll.m
  data.name = 'uniform swiss roll';
  n = data.n;
  a = 1;   % swiss roll goes from a*pi to b*pi
  b = 4;
  y = rand(2,n);
  data.z = y;  % the ground truth
  switch 1
   case 1
    % uniform distribution along the manifold (in data space)
    tt = sqrt((b*b-a*a)*y(1,:)+a*a);
   case 2
%    error('do not use this case')
    % nonuniform distribution along the manifold (in data space)
    tt = (b-a)*y(1,:) + a;
  end
  tt = pi*tt;
  % now tt should go from a*pi to b*pi
  height = 21*y(2,:);
  data.col = tt;
  data.x = [tt.*cos(tt); height; tt.*sin(tt)];
  data.az = -4;
  data.el = 13;

 case 1 % "swiss roll (uniform in embedding space)"
  % dataset extracted from LLE's swissroll.m
  data.name = 'classic swiss roll';
  n = data.n;
  a = 1;   % swiss roll goes from a*pi to b*pi
  b = 4;
  y = rand(2,n);
  tt = (b-a)*y(1,:) + a;
  tt = pi*tt;
  height = 21*y(2,:);
  data.col = tt;
  data.x = [tt.*cos(tt); height; tt.*sin(tt)];
  data.z = [tt; height]; % the ground truth
  data.az = -4;
  data.el = 13;

 case 11 % "undersampled swiss roll"
  % dataset extracted from LLE's swissroll.m
  data.name = 'undersampled swiss roll';
  data.n = 100;
  n = data.n;
  a = 1;   % swiss roll goes from a*pi to b*pi
  b = 4;
  y = rand(2,n);
  tt = (b-a)*y(1,:) + a;
  tt = pi*tt;
  height = 21*y(2,:);
  data.col = tt;
  data.x = [tt.*cos(tt); height; tt.*sin(tt)];
  data.z = [tt; height]; % the ground truth
  data.az = -4;
  data.el = 13;

 case 12 % "swiss roll"
  % dataset extracted from LLE's swissroll.m
  data.name = 'classic swiss roll';
  data.n = 400;
  n = data.n;
  a = 1;   % swiss roll goes from a*pi to b*pi
  b = 4;
  y = rand(2,n);
  tt = (b-a)*y(1,:) + a;
  tt = pi*tt;
  height = 21*y(2,:);
  data.col = tt;
  data.x = [tt.*cos(tt); height; tt.*sin(tt)];
  data.z = [tt; height]; % the ground truth
  data.az = -4;
  data.el = 13;

 case 2 % "scurve" dataset extracted from LLE's scurve.m
  data.name = 'scurve';
  n = data.n;
  % I added 'ceil' and 'floor' to account for the case that n is odd
  angle = pi*(1.5*rand(1,ceil(n/2))-1); height = 5*rand(1,n);
  data.x = [[cos(angle), -cos(angle(1:floor(n/2)))]; height;[ sin(angle), 2-sin(angle)]];
  data.col = [angle, 1.5*pi + angle];
  data.z = [angle, 1.5*pi+angle; height]; % the ground truth

 case 3 % "square" dataset, a uniformly sampled 2D square randomly
         % rotated into higher dimensions
  data.name = 'square';
  n = data.n;
  d = 2;     % intrinsic dimension
  % optional parameter for dataset==3
  % data.D      dimension of the data
  if ~isfield(data, 'D'), data.D = 3; end
  % generate random rotation matrix
  D = data.D;
  A = randn(D, D);
  options.disp = 0;
  [R, dummy] = eigs(A*A', d, 'LM', options);
  tt = rand(d, n);
  data.col = tt(1,:);
  data.x = R*tt;
  data.z = tt;   % the ground truth
  data.az = 7;
  data.el = 40;

 case 4 % spiral: two dimensional "swiss roll"
  data.name = 'spiral';
  n = data.n;
  tt = (3*pi/2)*(1+2*rand(1, n));
  data.col = tt;
  data.x = [tt.*cos(tt); tt.*sin(tt)];
  data.z = tt; % the ground truth

 case -4 % spiral: two dimensional "swiss roll"
  data.name = 'noisy spiral';
  n = data.n;
  tt = (3*pi/2)*(1+2*rand(1, n));
  data.col = tt;
  data.x = [tt.*cos(tt); tt.*sin(tt)];
  data.x = data.x + randn(size(data.x));
  data.z = tt; % the ground truth

 case 5 % hole: a dataset with a hole
  data.name = 'hole';
  n = data.n;
  data.x = rand(2,n) - 0.5;
  % punch a rectangular hole at the center
  l1 = 0.2; l2 = 0.2;
  ok = find((abs(data.x(1,:))>l1) | (abs(data.x(2,:))>l2));
  i = length(ok);
  data.x(:, 1:i) = data.x(:, ok);
  while (i<n)
    p = rand(2,1) - 0.5;
    if (abs(p(1))>l1) || (abs(p(2))>l2)
      i = i + 1;
      data.x(:,i) = p;
    end
  end
  data.col = data.x(2,:);
  data.z = data.x;

 case 6 % P : taken from Saul's slides
  % note that for k=20, isomap and lle work fine which is very different
  % from the plots that Saul showed in his slides.
  data.name = 'P';
  load x
  x(2,:) = 500-x(2,:);
  data.x = x;
  data.z = x;
  data.col = data.z(2,:);
  data.n = size(x, 2);

 case 7 % fishbowl: uniform in data space
  gamma = 0.8;
  data.name = 'fishbowl (uniform in data space)';
  n = data.n;
  data.x = rand(3,n)-0.5;
  %project all data onto the surface of the unit sphere
  data.x = data.x ./ repmat(sqrt(sum(data.x.*data.x, 1)), [3 1]);
  ok = find(data.x(3,:) < gamma);
  i = length(ok);
  data.x(:, 1:i) = data.x(:, ok);
  while (i < n)
    p = rand(3,1)-0.5;
    p = p / sqrt(p'*p);
    if (p(3) < gamma)
      i = i+1;
      data.x(:, i) = p;
    end
  end
  % the projection on the plane works as follows:
  % start a beam from (0,0,1) through each surface point on the sphere
  % and look where it hits the xy plane.
  data.z = data.x(1:2,:) ./ repmat(1-data.x(3,:), [2 1]);
  data.col = data.x(3,:);
  data.az = -18;
  data.el = 16;
 case 8 % fishbowl: uniform in embedding space
  data.name = 'fishbowl (uniform in embedding space)';
  n = data.n;
  data.z = rand(2, n) - 0.5;
  % keep the disc
  ok = find(sum(data.z .* data.z) <= 0.25);
  i = length(ok);
  data.z(:, 1:i) = data.z(:, ok);
  while (i < n)
    p = rand(2,1) - 0.5;
    if (p'*p <= 0.25)
      i = i + 1;
      data.z(:, i) = p;
    end
  end
  gamma = 0.8;  % same role/parameter as in case 7
  data.z = 2*sqrt((1+gamma)/(1-gamma))*data.z;
  % project the disc onto the sphere
  alpha = 2 ./ (1 + sum(data.z .* data.z, 1));
  data.x = [repmat(alpha, [2 1]).*data.z; zeros(1, n)];
  data.x(3,:) = 1-alpha;
  data.col = data.x(3,:);
  data.az = -18;
  data.el = 16;

 case 9  % a gaussian blob
  data.name = 'gaussian blob';
  n = data.n;
  data.x = randn(3,n);
  data.z = data.x(2:3,:);
  data.col = data.x(3,:);

end

data.D = size(data.x, 1);  % dimensionality of the data
% finally generate noise
data.e = randn(size(data.x));
data.x_noisefree = data.x;  % the noise free data
data.x = data.x_noisefree + sqrt(data.noise)*data.e;

% precalculate the distanzmatrix
data.distances = distanz(data.x);

repo: gspbox


text:
function F = fuzzysysfcn(inmf,outmf,vrange,op)
%FUZZYSYSFCN Fuzzy system function.
%   F = FUZZYSYSFCN(INMF,OUTMF,VRANGE,OP) creates a fuzzy system
%   function, F, corresponding to a set of rules and output membership
%   functions. INMF is an M-by-N matrix of input membership function
%   handles. M is the number of rules, and N is the number of fuzzy
%   system inputs. OUTMF is a cell array containing output membership
%   functions. numel(OUTMF) can be either M or M + 1. If it is M + 1,
%   then the "extra" output membership function is used for an
%   automatically computed "else rule." VRANGE is a two-element vector
%   specifying the valid range of input values for the output membership
%   functions. OP is a function handle specifying how to combine the
%   antecedents for each rule. OP can be either @min or @max. If OP is
%   omitted, then @min is used.
%
%  The output, F, is a function handle that computes the fuzzy system's
%  output, given a set of inputs, using the syntax:
%  out = F(Z1,Z2,Z3,...,ZN);
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

if nargin < 4
   op = @min;
end

% The lambda functions are independent of the inputs Z1,Z2,...,ZN, so
% they can be computed in advance.
L = lambdafcns(inmf,op);

F = @fuzzyOutput;

   %-------------------------------------------------------------------%
   function out = fuzzyOutput(varargin)
      Z = varargin;
      % The implication functions and aggregation functions have to be
      % computed separately for each input value.  Therefore we have to
      % loop over each input value to determine the corresponding output
      % value. Zk is a cell array that will be used to pass scalar
      % values for each input (Z1,Z2,...,ZN) to IMPLFCNS.
      Zk = cell(1,numel(Z));
      % Initialize the array of output values to be the same size as the
      % first input,Z{1}.
      out = zeros(size(Z{1}));
      for k = 1:numel(Z{1})
         for p = 1:numel(Zk)
            Zk{p} = Z{p}(k);
         end
         Q = implfcns(L,outmf,Zk{:});
         Qa = aggfcn(Q);
         out(k) = defuzzify(Qa,vrange);
      end
   end

end

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
function varargout = slogfrac(varargin)
%logfrac
%
% y = SLOGFRAC(x)
%
% Computes/declares log(1 + x(1)/x(2))

switch class(varargin{1})
    case 'double'
        x = varargin{1};
        % Safe version with defined negative values (helps fmincon when
        % outside feasible region)
        if all(x==0)
            varargout{1} = log(2);% ?definition...
        elseif x(1)==0
            varargout{1} = log(1);
        else
            aux = 1 + (x(1)./(x(2)+sqrt(eps)));
            if aux <= 0
                varargout{1} = -15;
            else
                varargout{1} = log(1 + (x(1)./(x(2)+sqrt(eps))));
            end

        end

    case 'sdpvar'
        if min(size(varargin{1}))>1
            error('SLOGFRAC only defined for vector arguments');
        else
            varargout{1} = yalmip('define',mfilename,varargin{1});
        end

    case 'char'

        operator = CreateBasicOperator('callback');
        operator.range = [-inf inf];
        operator.domain = [-inf inf];
        operator.bounds = @bounds;
        operator.derivative = @(x) ([1./(x(1)+x(2)+sqrt(eps));1./(x(1)+x(2)+sqrt(eps))-(x(2)+sqrt(eps)).^-1]);

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

function [L,U] = bounds(xL,xU)

% Derive bounds on x1/x2
z = [xL(1)./xL(2) xU(1)./xL(2) xL(1)./xU(2) xU(1)./xU(2)];
fL = min(z);
fU = max(z);

if fL <= -1
    L = -inf;
else
    L = log(1 + fL);
end
if fU <= -1
    U = -inf;
else
    U = log(1 + fU);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% DEMCMU35SEQUENCEOPTIMISE

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'cmu35gplvm';

% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);
skel = acclaimReadSkel('35.asf');
[tmpchan, skel] = acclaimLoadChannels('35_01.amc', skel);

% REMOVE LEG TEST
% Indices associated with right leg.
legInd = [8:14];
% Indices associated with upper body.
bodyInd = [21:50];

% Load saved model.
capName = dataSetName;
capName(1) = upper(capName(1));

for experimentNo = 1:3;
  load(['dem' capName num2str(experimentNo) '.mat']);
  for missing = 1:2
    if missing==1
      type = 'Leg';
      missingInd = legInd;
    else
      type = 'Body';
      missingInd = bodyInd;
    end
    if exist(['dem' capName 'Yvals' type num2str(experimentNo) '.mat'], 'file')==2
      load(['dem' capName 'Yvals' type num2str(experimentNo) '.mat'])
    else
      disp('Reconstructing ...');
      startInd = 63;

      Ytest(startInd:end, missingInd) = NaN;
      model = gpComputeAlpha(model);
      ll = [];
      for j = 1:size(model.X_u, 1)
        for i = 1:size(Ytest, 1);
          ll(i, j) = fgplvmPointLogLikelihood(model, model.X_u(j, :), ...
                                              Ytest(i, :));
        end
      end
      [void, ind] = max(ll, [], 2);
      Xinit = model.X_u(ind, :);
      Xpred = fgplvmOptimiseSequence(model, Xinit, Ytest, 1, 1000);
      Xpred = Xpred(startInd:end, :);
      Ytest = Ytest(startInd:end, :);
      Ypred = gpOut(model, Xpred);

      save(['dem' capName 'Yvals' type num2str(experimentNo) '.mat'], ...
           'Xpred', 'Ypred');
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_suite = test_grFindGeodesicPath
%TEST_GRFINDGEODESICPATH  One-line description here, please.
%
%   output = test_grFindGeodesicPath(input)
%
%   Example
%   testGrFindGeodesicPath
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-22,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

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

[nodes, edges] = createTestGraph01;
path = grFindGeodesicPath(nodes, edges, 3, 8);
exp = [4 1 3 10];
testCase.assertEqual(exp, path);

function testGraph03(testCase)

[nodes, edges] = createTestGraph03;
path = grFindGeodesicPath(nodes, edges, 1, 12);
exp = [1 2 5 7 11 13];
testCase.assertEqual(exp, path);

author: mattools

repo: matGeom


text:
%WAVE_BASES  1D Wavelet functions Morlet, Paul, or DOG
%
%  [DAUGHTER,FOURIER_FACTOR,COI,DOFMIN] = ...
%      wave_bases(MOTHER,K,SCALE,PARAM);
%
%   Computes the wavelet function as a function of Fourier frequency,
%   used for the wavelet transform in Fourier space.
%   (This program is called automatically by WAVELET)
%
% INPUTS:
%
%    MOTHER = a string, equal to 'MORLET' or 'PAUL' or 'DOG'
%    K = a vector, the Fourier frequencies at which to calculate the wavelet
%    SCALE = a number, the wavelet scale
%    PARAM = the nondimensional parameter for the wavelet function
%
% OUTPUTS:
%
%    DAUGHTER = a vector, the wavelet function
%    FOURIER_FACTOR = the ratio of Fourier period to scale
%    COI = a number, the cone-of-influence size at the scale
%    DOFMIN = a number, degrees of freedom for each point in the wavelet power
%             (either 2 for Morlet and Paul, or 1 for the DOG)
%
%----------------------------------------------------------------------------
%   Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo
%   University of Colorado, Program in Atmospheric and Oceanic Sciences.
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made.  This
%   routine is provided as is without any express or implied warranties
%   whatsoever.
%----------------------------------------------------------------------------
function [daughter,fourier_factor,coi,dofmin] = ...
	wave_bases(mother,k,scale,param);

mother = upper(mother);
n = length(k);

if (strcmp(mother,'MORLET'))  %-----------------------------------  Morlet
	if (param == -1), param = 6.;, end
	k0 = param;
	expnt = -(scale.*k - k0).^2/2.*(k > 0.);
	norm = sqrt(scale*k(2))*(pi^(-0.25))*sqrt(n);    % total energy=N   [Eqn(7)]
	daughter = norm*exp(expnt);
	daughter = daughter.*(k > 0.);     % Heaviside step function
	fourier_factor = (4*pi)/(k0 + sqrt(2 + k0^2)); % Scale-->Fourier [Sec.3h]
	coi = fourier_factor/sqrt(2);                  % Cone-of-influence [Sec.3g]
	dofmin = 2;                                    % Degrees of freedom
elseif (strcmp(mother,'PAUL'))  %--------------------------------  Paul
	if (param == -1), param = 4.;, end
	m = param;
	expnt = -(scale.*k).*(k > 0.);
	norm = sqrt(scale*k(2))*(2^m/sqrt(m*prod(2:(2*m-1))))*sqrt(n);
	daughter = norm*((scale.*k).^m).*exp(expnt);
	daughter = daughter.*(k > 0.);     % Heaviside step function
	fourier_factor = 4*pi/(2*m+1);
	coi = fourier_factor*sqrt(2);
	dofmin = 2;
elseif (strcmp(mother,'DOG'))  %--------------------------------  DOG
	if (param == -1), param = 2.;, end
	m = param;
	expnt = -(scale.*k).^2 ./ 2.0;
	norm = sqrt(scale*k(2)/gamma(m+0.5))*sqrt(n);
	daughter = -norm*(i^m)*((scale.*k).^m).*exp(expnt);
	fourier_factor = 2*pi*sqrt(2./(2*m+1));
	coi = fourier_factor/sqrt(2);
	dofmin = 1;
else
	error('Mother must be one of MORLET,PAUL,DOG')
end

return

author: grinsted

repo: wavelet-coherence


text:
function [events, times] = get_events_of_multi_phase_appliance(input_params)

    dataset = input_params.dataset;
    household = input_params.household;
    evaluation_days = input_params.evaluation_days;
    granularity = input_params.granularity;
    filteringMethod = input_params.filtering_method;
    filtLength = input_params.filtLength;
    edgeThreshold = input_params.filtLength;
    plevelMinLength = input_params.plevelMinLength;
    eventPowerStepThreshold = input_params.eventPowerStepThreshold;
    maxEventDuration = input_params.maxEventDuration;

    events = {};
    times = {};
    for phase = [1:3]
        % get real, apparent and reactive (distoritve and translative
        % component) power
        power = getPower(dataset, household, evaluation_days, granularity, phase);

        % apply filter to normalized apparent power and get edges
        function_handle = str2func(filteringMethod);
        normalized_apparent_power_filtered = function_handle(power.normalized_apparent, filtLength);
        [rows, cols] = find(abs(diff(normalized_apparent_power_filtered)) > edgeThreshold);
        edges = sparse(rows, cols, ones(length(rows),1), 1, size(normalized_apparent_power_filtered,2)-1);

        % get power levels (period between two edges with similar power values)
        [plevel] = getPowerLevelsStartAndEndTimes(edges, plevelMinLength);
        if isempty(plevel.startidx)
            continue;
        end

        % get characteristics of power levels
        plevel = getPowerLevelProperties(plevel, power, plevelMinLength);

        % generate event vectors by taking the diffference between two consecutive power levels
        event_vecs = zeros(length(plevel.startidx)-1, 4);
        eventIsValid = zeros(length(plevel.startidx), 1);
        numOfEvents = 0;
        for i = 1:length(plevel.startidx)-1
               if abs(plevel.mean.end(i,1) - plevel.mean.start(i+1,1)) > eventPowerStepThreshold && plevel.startidx(i+1) - plevel.endidx(i) < maxEventDuration
                    eventIsValid(i) = 1;
                    numOfEvents = numOfEvents + 1;
                    event_vecs(numOfEvents, 1:3) = plevel.mean.start(i+1, :) - plevel.mean.end(i, :);
                    max_std_true_power = max(plevel.std(i,1), plevel.std(i+1,1));
                    max_std_reactive_power = max(plevel.std(i,2), plevel.std(i+1,2));
                    oscillationTerm = norm([max_std_true_power, max_std_reactive_power]);
                    event_vecs(numOfEvents, 4) = oscillationTerm;
               end
        end
        event_vecs = event_vecs(1:numOfEvents, :);
        timeOfEvents = plevel.endidx(eventIsValid==1)';

        events{end+1} = event_vecs;
        times{end+1} = timeOfEvents;
    end
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function s=shuffle(x,c)
%Syntax: s=shuffle(x,c)
%______________________
%
% Makes c shuffled surrogates of a time series x.
%
% s is the shuffled time series.
% x is the original time series.
% c is the number of surrogates.
%
%
% References:
%
% Theiler J, Galdrikian B, Longtin A, Eubank S, Farmer D J (1992): Using
% Surrogate Data to Detect Nonlinearity in Time Series. In Nonlinear Modeling
% and Forecasting, eds. Casdagli M & Eubank S. 163-188. Addison-Wesley
%
% Theiler J, Eubank S,Galdrikian B, Longtin A,  Farmer D J (1992): Testing
% for nonlinearity in time series: the method of surrogate data. Physica D
% 58: 77-94
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 12 Apr 2001.

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % N is the time series length
   N=length(x);
end

if nargin<2 | isempty(c)==1
   c=1;
else
   % c must be scalar
   if sum(size(c))>2
      error('c must be scalar.');
   end
   % c must be greater or equal than 1
   if c<1
      error('c must be greater or equal than 1.');
   end
end

for i=1:c

    % Make a random vector with indices up to n
    s1=randperm(N);

    % Shuffle the original x with respect to s1
    s(:,i)=x(s1);
end

author: Sable

repo: mcbench-benchmarks


text:
function test94()

% Test for fixedrankNewGHquotientfactory geometry (low rank completion)
%
% Test CG scheme

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:

clc; close all;
%     reset(RandStream.getDefaultStream);
%     randnfoo = randn(1, 1); %#ok<NASGU>

m = 500;
n = 500;
r = 5;
A = randn(m, r);
B = randn(n, r);
C = A*B';

problem.M = fixedrankNewGHquotientfactory(m, n, r);

df = problem.M.dim();
p = 4*df/(m*n);
mask = rand(m, n) <= p;

symm = @(M) .5*(M+M');

problem.cost = @cost;
    function f = cost(X)
        f = 0.5*(norm(mask.*(X.G*X.H' - C), 'fro')^2);
    end

problem.grad = @grad;
    function g = grad(X)
        G = X.G;
        H = X.H;
        GtG = X.G'*X.G;
        HtH = X.H'*X.H;
        invGtG = eye(size(GtG)) / GtG;
        invHtH = eye(size(HtH)) / HtH;

        R = mask.^2 .* (G*H' - C);
        g = struct('G', (R*H)*invHtH, ...
            'H', (R'*G)*invGtG );
    end

problem.hess = @hess;
    function Hess = hess(X, eta)

        S = mask.*( X.G*X.H' - C);
        S_star  = mask.*(eta.G*X.H' + X.G*eta.H');

        GtG = X.G'*X.G;
        HtH = X.H'*X.H;
        invGtG = eye(size(GtG)) / GtG;
        invHtH = eye(size(HtH)) / HtH;

        Hess.G = S_star*X.H*invHtH;
        Hess.H = S_star'*X.G*invGtG;

        ShH = S*eta.H;
        ShG = S'*eta.G;
        SH = S*X.H;
        SG = S'*X.G;
        gradG = SH*invHtH;
        gradH = SG*invGtG;

        Hess.G = Hess.G + ShH*invHtH - 2*SH*(invHtH * symm(eta.H'*X.H) * invHtH);
        Hess.H = Hess.H + ShG*invGtG - 2*SG*(invGtG * symm(eta.G'*X.G) * invGtG);

        % I still need a correction factor for the non-constant metric
        Hess.G = Hess.G + gradG*symm(eta.H'*X.H)*invHtH + eta.G*symm(gradH'*X.H)*invHtH - X.G*symm(eta.H'*gradH)*invHtH;
        Hess.H = Hess.H + gradH*symm(eta.G'*X.G)*invGtG + eta.H*symm(gradG'*X.G)*invGtG - X.H*symm(eta.G'*gradG)*invGtG;

        Hess = problem.M.proj(X, Hess);

    end

%
%     checkgradient(problem);
%     drawnow;
%     pause;
%     checkhessian(problem);
%     drawnow;
%     pause;
%
%     problem = rmfield(problem, 'hess');

% [U S V] = svds(mask.*C, r);
% G0 = U*(S.^0.5);
% H0 = V*(S.^0.5);
G0 = randn(m, r);
H0 = randn(n, r);

X0 = struct('G', G0, 'H', H0 );

options.statsfun = @statsfun;

    function stats = statsfun(problem, x, stats)
        stats.RMSE = nan;%norm(C - x.M*x.N', 'fro')/sqrt(m*n);
    end

% options.linesearch = @linesearch_adaptive;

options.maxiter = inf;
options.maxinner = 30;
options.maxtime = 120;
%     options.mininner = problem.M.dim();
options.tolgradnorm = 1e-2;
options.Delta_bar = min(m, n) * r;
options.Delta0 = options.Delta_bar / 8;

%     options.useRand = true;
% fprintf('----------- TR -----------\n');
% [Xopt costopt info] = trustregions(problem, X0, options);

fprintf('----------- Steepest descent with NW -----------\n');
options.linesearch = @linesearch;
[Xopt costopt infos_SD] = steepestdescent(problem, X0, options);

fprintf('----------- Steepest descent with adaptive -----------\n');
options.linesearch = @linesearch_adaptive;
[Xopt costopt infos_SD_adaptive] = steepestdescent(problem, X0, options);

options.beta_type = 'P-R'; % Other options are 'steep' and 'FR'
options.orth_value = 0.1; Inf;

fprintf('----------- CG algorithm with default -----------\n');
options.linesearch = @linesearch;
[Xopt costopt infos_CG] = conjugategradient(problem, X0, options);

fprintf('----------- CG algorithm with adaptive -----------\n');
options.linesearch = @linesearch_adaptive;
[Xopt costopt infos_CG_adaptive] = conjugategradient(problem, X0, options);

%     keyboard;

% subplot(3, 1, 1);
% semilogy([info.iter], [info.cost], '.-');
% subplot(3, 1, 2);
% semilogy([info.iter], [info.gradnorm], '.-');
% subplot(3, 1, 3);
% semilogy([info.iter], [info.RMSE], '.-');

% Cost versus iterations
fs = 20;
figure;
semilogy(0:length([infos_SD.cost])-1,[infos_SD.cost],'Color','b','LineWidth',2);
hold on;
semilogy(0:length([infos_SD_adaptive.cost])-1,[infos_SD_adaptive.cost],'Color','r','LineWidth',2);
semilogy(0:length([infos_CG.cost])-1,[infos_CG.cost],'--','Color','b','LineWidth',2);
semilogy(0:length([infos_CG_adaptive.cost])-1,[infos_CG_adaptive.cost],'--','Color','r','LineWidth',2);

hold off;
ax1 = gca;
set(ax1,'FontSize',fs);
xlabel(ax1,'Number of iterations','FontSize',fs);
ylabel(ax1,'Cost ','FontSize',fs);
mincost = min([infos_SD.cost]);
maxcost = max([infos_SD.cost]);
axis([get(gca,'XLim') mincost maxcost])
legend('SD NW','SD adaptive', 'CG NW','CG adaptive');
legend 'boxoff';
box off;

end

author: NicolasBoumal

repo: manopt


text:
function y=atanSlope3(x)
%%ATANSLOPE3 Evaluate the function -3*(atan(x)-x)/x^3 avoiding problems at
%            the point x=0, where the limit equals 1. This function arises
%            in range enclosure algorithms for interval arithmetic.
%
%INPUTS: x A vector of matrix of real values at which the function should
%          be evaluated.
%
%OUTPUTS: y The value(s) of the function -3*(atan(x)-x)/x^3 at the point(s)
%           in x.
%
%For values of x having magnitude over 1e-2 or for complex values, the
%function is evaluated as written. For smaller values, a Taylor series
%expansion is used to avoid the singularity near the origin.
%
%This is one of the slope function in Table 10.6 of Section 10.6.2 of [1].
%This file implements the function for non-intervals. See the Interval
%class for an implementation for Interval arithmetic.
%
%REFERENCES:
%[1] IEEE Standard for Interval Arithmetic," in IEEE Std 1788-2015,
%    pp.1-97, June 30 2015.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Allocate space
y=zeros(size(x));

%For large enough values, use the standard function.
sel=abs(x)>1e-2|~isreal(x);
x1=x(sel);
y(sel)=-3*(atan(x1)-x1)./x1.^3;

%For smaller magnitude values, use an eigth-order Taylor series expansion.
%The relative error at 10^-2 is on the order of 1e-23, which is less than
%the precision of a double floating point number at that value. The series
%is given in Horner form.
x1=x(~sel);
y(~sel)=1+x1.^2.*(-(3/5)+x1.^2.*(3/7+x1.^2.*(-(1/3)+(3*x1.^2)/11)));

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
global I hi

I=imread('smile.png'); % load image

hi=imshow(I,'InitialMagnification',600);
hold on;
tol=2; % tolerance
r=[20;30]; % start point of flood
zz_flood_fill_movie(r,tol); % make flood fill


text:
%% Demo 2: With missing data

% Add path
addpath(genpath('../nricp'));
addpath(genpath('../../toolbox_graph/toolbox_graph/'));
addpath(genpath('../../icp'));

% load data
load ('../data/faceSource.mat');
load ('../data/faceTargetMissing.mat');

% Specify that surface normals are available and can be used.
Options.useNormals = 1;

% Specify that the source deformations should be plotted.
Options.plot = 1;

% Perform non-rigid ICP
[pointsTransformed, X] = nricp_landmarks(Source, TargetMissing, Options);

author: RhythmJnh

repo: Non-rigid-ICP


text:
function y = geomean( varargin )

%GEOMEAN   Geometric mean.
%   GEOMEAN(X) = GEO_MEAN(X) = PROD(X).^(1/LENGTH(X)). We have replaced this
%   function with GEO_MEAN to better match our function naming conventions.
%   Please start using it instead.

warning( 'CVX:Renamed', [ ...
    'The function "geomean" has been renamed "geo_mean". Please start\n', ...
    'using the new name. The old name will be removed in a future release.' ], 1 );

y = geo_mean( varargin{:} );

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

author: cvxr

repo: CVX


text:
function imOut=compareImagesBySide(imIn,showlines,shownumbers,frac)
% Compares images by putting parts of the imags side by side.
% imOut=compareImagesBySide(imIn,showlines,shownumbers,frac)
% imIn - stack of images to be compared
% showlines - if true (default) then line is drawn to separate individual sections.
% shownumbers - if true (default) then number to individual section is
% displayed.
% frac - vector specifying fractions of the individual sections (default: frac=size(Im(1))/numberOfImages

if ~exist('showlines','var')
    showlines = 1;
end

if ~exist('shownumbers','var')
    shownumbers = 1;
end

sizeIm=size(imIn);
numberOfImages=sizeIm(3);
if ~exist('frac','var')
    frac=repmat(floor(sizeIm(1)/numberOfImages),1,size(imIn,3));
end

imOut = imIn(0:frac(1)-1,:,0);
for ii=1:numberOfImages-1
    startInd=frac(ii);
    endInd=sum(frac(1:ii+1))-1;
    imOut = cat(1,imOut,imIn(startInd:endInd,:,ii));
end

colorvalue = 'w';

if showlines
dipshow(double(imOut))
for ii=1:numberOfImages-1
    startInd=(ii)*frac(ii)-1;
    line([startInd,startInd],[0,sizeIm(2)],'color',colorvalue,'linewidth',1)
end
end

if shownumbers
for ii=1:numberOfImages
    startInd=(ii)*frac(ii)-1-0.5*frac(ii);
    text(startInd,sizeIm(2)/10,num2str(ii),'color',colorvalue,'fontsize',20)
end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function T = getTransViaRotateGivenCenter(theta_vector,center,rotatorLength)
%T = getTransViaRotateGivenCenter(theta_vector,win_size)
%   T: m*1 t_concord
%   theta_vector: m*1 in degree! Make sure abs(theta)<180
%   center: 1*2 indicates the x-y coordinates of rotation center
%   rotatorLength: scalar, length of the ratator. Set this variable to
%       avoid numerical problems.

warning('This function might get non-affine tform');
assert(size(theta_vector,1)==1 || size(theta_vector,2)==1);
assert(size(center,1)==1 || size(theta_vector,2)==1);
assert(length(center)==2);
assert(length(rotatorLength)==1);
theta_vector = theta_vector(:);
assert(all(abs(theta_vector) < 180));

pose_src = zeros(length(theta_vector),4);
pose_src(:,[1 3]) = repmat(center(:)',length(theta_vector),1);
pose_src(:,2) = center(1) + rotatorLength * cosd(theta_vector);
pose_src(:,4) = center(2) + rotatorLength * sind(theta_vector);
pose_dst = [center(1) center(1)+rotatorLength center(2) center(2)];

T = getTransToSpecific(pose_src,pose_dst);

end

author: zhusz


text:
%%
% Test for power diagrams

%%
% Test for power diagrams

addpath('power_bounded/');

rep = 'results/power-diagrams/';
[~,~] = mkdir(rep);

% number of sites
N = 50;
% x and y coordinate of the sites
sigma = .1;
xy = sigma*randn(N, 2)+1/2;
xy(1,:) = 1/2; % center points

% the bounding box in clockwise order
bb = [0,0; 0,1; 1,1; 1,0];

lambda = [1 5 10 100];
lambda = linspace(1,1.1,5);

for i=1:length(lambda)% weights
    w = ones(N,1);w(1) = lambda(i);
    % get the power diagram
    [V,C] = power_bounded(xy(:,1),xy(:,2), w, bb);
    % draw the resulted power diagram
    plot_power(xy,V,C,bb);
    saveas(gcf, [rep 'power-' num2str(i), '.eps'], 'epsc');
    drawnow;
end

author: optimaltransport


text:
function EnableDisableGuiPlotBttns(handles)
global maingui
if isempty(handles)
    return
end

if ~isempty(maingui.dataTree.currElem.GetRawData())
    set(handles.radiobuttonPlotRaw, 'enable','on')
else
    set(handles.radiobuttonPlotRaw, 'enable','off')
end
if ~isempty(maingui.dataTree.currElem.GetDod())
    set(handles.radiobuttonPlotOD, 'enable','on')
else
    set(handles.radiobuttonPlotOD, 'enable','off')
end
if ~isempty(maingui.dataTree.currElem.GetDc())
    set(handles.radiobuttonPlotConc, 'enable','on')
else
    set(handles.radiobuttonPlotConc, 'enable','off')
end

raw_checked  = get(handles.radiobuttonPlotRaw, 'value');
if strcmp(get(handles.radiobuttonPlotRaw, 'enable'), 'on')
    raw_enable  = true;
else
    raw_enable  = false;
end
if strcmp(get(handles.radiobuttonPlotOD, 'enable'), 'on')
    OD_enable  = true;
else
    OD_enable  = false;
end
if strcmp(get(handles.radiobuttonPlotConc, 'enable'), 'on')
    Conc_enable  = true;
else
    Conc_enable  = false;
end

if ~isempty(maingui.dataTree.groups(1).CondNames) > 0  % If there is stim
    iCondGrp = get(handles.popupmenuConditions, 'value');
    CondName = maingui.dataTree.groups(1).CondNames{iCondGrp};
    if ~isempty(maingui.dataTree.currElem.GetDodAvg(CondName)) || ~isempty(maingui.dataTree.currElem.GetDcAvg(CondName))
        set(handles.checkboxPlotHRF, 'enable','on');
        if ~isa(maingui.dataTree.currElem, 'RunClass')
            set(handles.checkboxPlotHRF, 'value',1);
            if ~isempty(maingui.dataTree.currElem.GetDcAvg())
                set(handles.radiobuttonPlotConc, 'value',1);
            elseif ~isempty(maingui.dataTree.currElem.GetDodAvg())
                set(handles.radiobuttonPlotOD, 'value',1);
            end
        end
    else
        set(handles.checkboxPlotHRF, 'enable','off');
        set(handles.checkboxPlotHRF, 'value',0);
    end
elseif raw_enable && raw_checked
    set(handles.checkboxPlotHRF, 'enable','off');
    set(handles.checkboxPlotHRF, 'value',0);
elseif ~OD_enable && ~Conc_enable
    set(handles.checkboxPlotHRF, 'enable','off');
    set(handles.checkboxPlotHRF, 'value',0);
else
    set(handles.checkboxPlotHRF, 'enable','off');
    set(handles.checkboxPlotHRF, 'value',0);
end

if isa(maingui.dataTree.currElem, 'RunClass')
    if ~OD_enable && ~Conc_enable
        set(handles.radiobuttonPlotRaw, 'value',1)
    end
end

author: BUNPC


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [phase] = getPhase(household, applianceID, dataset)

    % returns the phase of the specified (house, appliance)-pair

    phase_matrix = getPhaseMatrix(dataset);
    phase = phase_matrix(applianceID, household);

end


text:
%ByteCodes Bytecode enumeration for the brick
%
% Notes::
% - bytecodes can be found in the EV3 documentation and source code
% (bytecodes.h)

classdef ByteCodes < uint8
    enumeration
        %VM
        Error (0) % = 0x00,
        Nop (1) % = 0x01,
        ProgramStop (2) % = 0x02,
        ProgramStart (3) % = 0x03,
        ObjectStop (4) % = 0x04,
        ObjectStart (5) % = 0x05,
        ObjectTrig (6) % = 0x06,
        ObjectWait (7) % = 0x07,
        Return (8) % = 0x08,
        Call (9) % = 0x09,
        ObjectEnd (10) % = 0x0A,
        Sleep (11) % = 0x0B,
        ProgramInfo (12) % = 0x0C,
        Label(13) % = 0x0D,
        Probe(14) % = 0x0E,
        Do(15) % = 0x0F,

        %Move
        InitBytes (47) %= 0x2F,

        %Branch
        Jr (64) %= 0x40,

        %VM
        Info (124) %= 0x7C,
        String (125) %= 0x7D,
%        MemoryWrite = 0x7E,
%        MemoryRead = 0x7F,

        %UI
        UIFlush (128) %= 0x80,
        UIRead (129) %= 0x81,
        UIWrite (130) %= 0x82,
        UIButton (131) %= 0x83,
        UIDraw (132) %= 0x84,

        %Timer
        TimerWait (133) %= 0x85,
        TimerReady (134) %= 0x86,
        TimerRead (135) %= 0x87,

        %Sound
        Sound (148) %= 0x94,
        SoundTest (149) %= 095,
        SoundReady (150) % = 0x96,

%       %Input
%       InputSample = 0x97,
        InputDeviceList (152) % = 0x98,
        InputDevice (153) %= 0x99,
        InputRead (154) % = 0x9a,
%       InputTest = 0x9b,
        InputReady (156) % = 0x9c,
        InputReadSI (157) % = 0x9d,
%       InputReadExt = 0x9e,
%       InputWrite = 0x9f,
%
%       %Output
%       OutputGetType = 0xa0,
        OutputSetType (161) % = 0xa1,
        OutputReset (162) % = 0xa2,
        OutputStop (163) % = 0xA3,
        OutputPower (164) % = 0xA4,
        OutputSpeed (165) % = 0xA5,
        OutputStart	(166) % = 0xA6,
        OutputPolarity (167) % = 0xA7,
        OutputRead (168) % = 0xA8,
        OutputTest (169) % = 0xA9,
        OutputReady (170) % = 0xAA,
%       OutputPosition = 0xAB, NOT IMPLEMENTED IN THE CODE
        OutputStepPower (172) % = 0xAC,
        OutputTimePower (173) % = 0xAD,
        OutputStepSpeed (174) % = 0xAE,
        OutputTimeSpeed (175) % = 0xAF,
        OutputStepSync (176) % = 0xB0,
        OutputTimeSync (177) % = 0xB1,
        OutputClrCount (178) % = 0xB2,
        OutputGetCount (179) % = 0xB3,
%
%       %Memory
%        File = 0xC0,
%        Array = 0xc1,
%        ArrayWrite = 0xc2,
%        ArrayRead = 0xc3,
%        ArrayAppend = 0xc4,
%        MemoryUsage = 0xc5,
%        FileName = 0xc6,
%

         %Communications
         COMGet (211) % = 0xD3,
         COMSet (212) % = 0xD4,

%       %Mailbox
         MailboxOpen (216) % = 0xD8,
         MailboxWrite (217) % = 0xD9,
         MailboxRead (218) % = 0xDA,
         MailboxTest (219) % = 0xDB,
         MailboxReady (220) % = 0xDC,
         MailboxClose (221) % = 0xDD,
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function varargout = angleSort3d(pts, varargin)
%ANGLESORT3D Sort 3D coplanar points according to their angles in plane
%
%   PTS2 = angleSort3d(PTS);
%   Considers all points are located on the same plane, and sort them
%   according to the angle on plane. PTS is a [Nx2] array. Note that the
%   result depend on plane orientation: points can be in reverse order
%   compared to expected. The reference plane is computed besed on the 3
%   first points.
%
%   PTS2 = angleSort3d(PTS, PTS0);
%   Computes angles between each point of PTS and PT0. By default, uses
%   centroid of points.
%
%   PTS2 = angleSort3d(PTS, PTS0, PTS1);
%   Specifies the point which will be used as a start.
%
%   [PTS2, I] = angleSort3d(...);
%   Also return in I the indices of PTS, such that PTS2 = PTS(I, :);
%
%   See also:
%   points3d, angles3d, angleSort
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2005-11-24
% Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY :
%   04/01/2007: remove unused variables

% default values
pt0     = mean(pts, 1);
pt1     = pts(1,:);

if length(varargin)==1
    pt0 = varargin{1};
elseif length(varargin)==2
    pt0 = varargin{1};
    pt1 = varargin{2};
end

% create support plane
plane   = createPlane(pts(1:3, :));

% project points onto the plane
pts2d   = planePosition(pts, plane);
pt0     = planePosition(pt0, plane);
pt1     = planePosition(pt1, plane);

% compute origin angle
theta0  = atan2(pt1(2)-pt0(2), pt1(1)-pt0(1));
theta0  = mod(theta0 + 2*pi, 2*pi);

% translate to reference point
n       = size(pts, 1);
pts2d   = pts2d - repmat(pt0, [n 1]);

% compute angles
angle   = atan2(pts2d(:,2), pts2d(:,1));
angle   = mod(angle - theta0 + 4*pi, 2*pi);

% sort points according to angles
[angle, I] = sort(angle); %#ok<ASGLU>

% format output
if nargout<2
    varargout{1} = pts(I, :);
elseif nargout==2
    varargout{1} = pts(I, :);
    varargout{2} = I;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function filtered_image = gaussianbpf(I,d0,d1)
	% Butterworth Bandpass Filter
	% This simple  function was written for my Digital Image Processing course
	% at Eastern Mediterranean University taught by
	% Assoc. Prof. Dr. Hasan Demirel
	% for the 2010-2011 Spring Semester
	% for the complete report:
	% http://www.scribd.com/doc/51981950/HW4-Frequency-Domain-Bandpass-Filtering
	%
	% Written By:
	% Leonardo O. Iheme (leonardo.iheme@cc.emu.edu.tr)
	% 24th of March 2011
	%
	%   I = The input grey scale image
	%   d0 = Lower cut off frequency
	%   d1 = Higher cut off frequency
	%
	% The function makes use of the simple principle that a bandpass filter
	% can be obtained by multiplying a lowpass filter with a highpass filter
	% where the lowpass filter has a higher cut off frquency than the high pass filter.
	%
	% Usage GAUSSIANBPF(I,DO,D1)
	% Example
	% ima = imread('grass.jpg');
	% ima = rgb2gray(ima);
	% filtered_image = gaussianbpf(ima,30,120);
	% Gaussian Bandpass Filter
	%
	% biafra ahanonu
	% updated: 2013.11.09 [13:20:43]
	% changelog
		% 2013.11.09 [13:21:06] updated function so it no longer converts to uint8 as this causes problems. Also, changed filtered_image = fftI + filter3.*fftI; to filtered_image = filter3.*fftI;
		% 2014.06.03 - updated to do binary fft

	f = double(I);
	[nx ny] = size(f);
	% f = uint8(f);
	fftI = fft2(f,2*nx-1,2*ny-1);
	fftI = fftshift(fftI);

	% subplot(2,2,1)
	% imshow(f,[]);
	% title('Original Image')

	% subplot(2,2,2)
	% fftshow(fftI,'log')
	% title('Fourier Spectrum of Image')
	% Initialize filter.
	% filter1 = ones(2*nx-1,2*ny-1);
	% filter2 = ones(2*nx-1,2*ny-1);
	% filter3 = ones(2*nx-1,2*ny-1);
	% % ======
	% % Gaussian
	% for i = 1:2*nx-1
	%     for j =1:2*ny-1
	%         dist = ((i-(nx+1))^2 + (j-(ny+1))^2)^.5;
	%         % Use Gaussian filter.
	%         filter1(i,j) = exp(-dist^2/(2*d1^2));
	%         filter2(i,j) = exp(-dist^2/(2*d0^2));
	%         filter3(i,j) = 1.0 - filter2(i,j);
	%         filter3(i,j) = filter1(i,j).*filter3(i,j);
	%     end
	% end
	% ======
	% binary
	imageSize = size(fftI);
	ci = [round(imageSize(1)/2), round(imageSize(2)/2), 10];     % center and radius of circle ([c_row, c_col, r])
	[xx,yy] = ndgrid((1:imageSize(1))-ci(1),(1:imageSize(2))-ci(2));
	filter3 = logical((xx.^2 + yy.^2)<(ci(3)^2));
	% imagesc(mask)
	% ======
	% Update image with passed frequencies
	fftI_filtered = filter3.*fftI;
	% subplot(2,2,3)
	% fftshow(filter3,'log')
	% title('Frequency Domain Filter Function Image')
	filtered_image = ifftshift(fftI_filtered);
	filtered_image = ifft2(filtered_image,2*nx-1,2*ny-1);
	filtered_image = real(filtered_image(1:nx,1:ny));
	filtered_image = double(filtered_image);
	% filtered_image = uint8(filtered_image);
	% ======
	% imagesc(filter3)
	subplot(2,2,1)
	fftshow(fftI_filtered,'log')
	subplot(2,2,2)
	fftshow(fftI,'log')
	subplot(2,2,3)
	imagesc(I)
	subplot(2,2,4)
	imshow(filtered_image,[])
	title('Bandpass Filtered Image')

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [S, V, D] = MySVDs(A, R)
[m, n] = size(A);
if 10*m < n
    AAT = A*A';
    [S, V, D] = svd(AAT);
    V = diag(V) .^ 0.5;
    tol = max(size(A)) * eps(max(V));
    R = min([R, sum(V > tol)]);
    V = V(1:R);
    S = S(:,1:R);
    D = A'*S*diag(1./V);
    V = diag(V);
    return;
end
if m > 10*n
    [S, V, D] = MySVD(A');
    mid = D;
    D = S;
    S = mid;
    return;
end
[S,V,D] = svd(A);
S = S(:, 1:R);
V = V(1:R, 1:R);
D = D(:, 1:R);

author: andrewssobral


text:
function pass = test_real( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e2*pref.techPrefs.chebfuneps;

% Example 1
f = ballfun(@(x,y,z)z);
V = real(ballfunv(f,f,f));
g = ballfun(@(x,y,z)z);
exact = ballfunv(g,g,g);
pass(1) = norm(V-exact)<tol;

% Example 2
f = ballfun(@(x,y,z)y+1i*z);
V = real(ballfunv(f,f,f));
g = ballfun(@(x,y,z)y);
exact = ballfunv(g,g,g);
pass(2) = norm(V-exact)<tol;

% Example 3
f1 = ballfun(@(x,y,z)x);
f2 = ballfun(@(x,y,z)1i*z);
f3 = ballfun(@(x,y,z)cos(y)+1i*sin(x));
V = real(ballfunv(f1,f2,f3));
g1 = ballfun(@(x,y,z)x);
g2 = ballfun(@(x,y,z)0);
g3 = ballfun(@(x,y,z)cos(y));
exact = ballfunv(g1,g2,g3);
pass(3) = norm(V-exact)<tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun


text:
function Y = spdiag(V,K)
%SPDIAG Sparse diagonal matrices.
%   SPDIAG(V,K) when V is a vector with N components is a sparse square
%   matrix of order N+ABS(K) with the elements of V on the K-th diagonal.
%   K = 0 is the main diagonal, K > 0 is above the main diagonal and K < 0
%   is below the main diagonal.
%
%   SPDIAG(V) is the same as SPDIAG(V,0) and puts V on the main diagonal.
%
%   See also DIAG, SPDIAGS.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  June 2008

if (nargin<2)
  K = 0;
end

n = length(V) + abs(K);

if (K>0)
  i = 1:length(V);
  j = K+1:n;
elseif (K<0)
  i = -K+1:n;
  j = 1:length(V);
else
  i = 1:n;
  j = 1:n;
end

Y = sparse(i,j,V(:),n,n);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [V, converged, i] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt)
%FDPF  Solves the power flow using a fast decoupled method.
%   [V, CONVERGED, I] = FDPF(YBUS, SBUS, V0, BP, BPP, REF, PV, PQ, MPOPT)
%   solves for bus voltages given the full system admittance matrix (for
%   all buses), the complex bus power injection vector (for all buses),
%   the initial vector of complex bus voltages, the FDPF matrices B prime
%   and B double prime, and column vectors with the lists of bus indices
%   for the swing bus, PV buses, and PQ buses, respectively. The bus voltage
%   vector contains the set point for generator (including ref bus)
%   buses, and the reference angle of the swing bus, as well as an initial
%   guess for remaining magnitudes and angles. MPOPT is a MATPOWER options
%   vector which can be used to set the termination tolerance, maximum
%   number of iterations, and output options (see MPOPTION for details).
%   Uses default options if this parameter is not given. Returns the
%   final complex voltages, a flag which indicates whether it converged
%   or not, and the number of iterations performed.
%
%   See also RUNPF.

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

%% default arguments
if nargin < 7
    mpopt = mpoption;
end

%% options
tol     = mpopt.pf.tol;
max_it  = mpopt.pf.fd.max_it;
lu_vec  = have_feature('lu_vec');

%% initialize
converged = 0;
i = 0;
V = V0;
Va = angle(V);
Vm = abs(V);

%% set up indexing for updating V
npv = length(pv);
npq = length(pq);

%% evaluate initial mismatch
mis = (V .* conj(Ybus * V) - Sbus(Vm)) ./ Vm;
P = real(mis([pv; pq]));
Q = imag(mis(pq));

%% check tolerance
normP = norm(P, inf);
normQ = norm(Q, inf);
if mpopt.verbose > 1
    fprintf('\niteration     max mismatch (p.u.)  ');
    fprintf('\ntype   #        P            Q     ');
    fprintf('\n---- ----  -----------  -----------');
    fprintf('\n  -  %3d   %10.3e   %10.3e', i, normP, normQ);
end
if normP < tol && normQ < tol
    converged = 1;
    if mpopt.verbose > 1
        fprintf('\nConverged!\n');
    end
end

%% reduce B matrices
Bp = Bp([pv; pq], [pv; pq]);
Bpp = Bpp(pq, pq);

%% factor B matrices
if lu_vec
    [Lp,  Up,  pp,  qp ] = lu(Bp,  'vector');
    [Lpp, Upp, ppp, qpp] = lu(Bpp, 'vector');
    [junk, iqp ] = sort(qp);
    [junk, iqpp] = sort(qpp);
    % [~, iqp ] = sort(qp);
    % [~, iqpp] = sort(qpp);
else
    [Lp, Up, Pp] = lu(Bp);
    [Lpp, Upp, Ppp] = lu(Bpp);
end

%% do P and Q iterations
while (~converged && i < max_it)
    %% update iteration counter
    i = i + 1;

    %%-----  do P iteration, update Va  -----
    if lu_vec
        dVa = -( Up \  (Lp \ P(pp)) );
        dVa = dVa(iqp);
    else
        dVa = -( Up \  (Lp \ (Pp * P)));
    end

    %% update voltage
    Va([pv; pq]) = Va([pv; pq]) + dVa;
    V = Vm .* exp(1j * Va);

    %% evalute mismatch
    mis = (V .* conj(Ybus * V) - Sbus(Vm)) ./ Vm;
    P = real(mis([pv; pq]));
    Q = imag(mis(pq));

    %% check tolerance
    normP = norm(P, inf);
    normQ = norm(Q, inf);
    if mpopt.verbose > 1
        fprintf('\n  P  %3d   %10.3e   %10.3e', i, normP, normQ);
    end
    if normP < tol && normQ < tol
        converged = 1;
        if mpopt.verbose
            fprintf('\nFast-decoupled power flow converged in %d P-iterations and %d Q-iterations.\n', i, i-1);
        end
        break;
    end

    %%-----  do Q iteration, update Vm  -----
    if lu_vec
        dVm = -( Upp \ (Lpp \ Q(ppp)) );
        dVm = dVm(iqpp);
    else
        dVm = -( Upp \ (Lpp \ (Ppp * Q)) );
    end

    %% update voltage
    Vm(pq) = Vm(pq) + dVm;
    V = Vm .* exp(1j * Va);

    %% evalute mismatch
    mis = (V .* conj(Ybus * V) - Sbus(Vm)) ./ Vm;
    P = real(mis([pv; pq]));
    Q = imag(mis(pq));

    %% check tolerance
    normP = norm(P, inf);
    normQ = norm(Q, inf);
    if mpopt.verbose > 1
        fprintf('\n  Q  %3d   %10.3e   %10.3e', i, normP, normQ);
    end
    if normP < tol && normQ < tol
        converged = 1;
        if mpopt.verbose
            fprintf('\nFast-decoupled power flow converged in %d P-iterations and %d Q-iterations.\n', i, i);
        end
        break;
    end
end

if mpopt.verbose
    if ~converged
        fprintf('\nFast-decoupled power flow did not converge in %d iterations.\n', i);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% Copyright (C) 2014-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym acos (@var{x})
%% Symbolic acos function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = acos (x)
%%   @result{} y = (sym) acos(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = acos(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('acos', x);
end

%!error acos (sym(1), 2)
%!assert (isequaln (acos (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = acos(x);
%! f2 = acos(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = acos(A);
%! f2 = acos(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = acos (d);
%! f = acos (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function wysiwyg
    %WYSIWYG -- this function is called with no args and merely
    %       changes the size of the figure on the screen to equal
    %       the size of the figure that would be printed,
    %       according to the papersize attribute.  Use this function
    %       to give a more accurate picture of what will be
    %       printed.
    %       Dan(K) Braithwaite, Dept. of Hydrology U.of.A  11/93

    unis = get(gcf,'units');
    ppos = get(gcf,'paperposition');
    set(gcf,'units',get(gcf,'paperunits'));
    pos = get(gcf,'position');
    pos(3:4) = ppos(3:4);
    set(gcf,'position',pos);
    set(gcf,'units',unis);

author: CelsoReyes