text:
% JM: 2006/07/19
% viewer for data Ic on domain Omega

function varargout = viewImage(Ic,Omega,m,varargin);
% this is just output and hence not commented; fell free to improve (-;

global VIEWER VIEWPARA

% default parameter:
fig       = 1;
name      = 'image';
plots     = 1;
sub       = [1,1,1];
figname   = '';

%% overwrite default parameter
for k=1:1:length(varargin)/2,
  %disp([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
  eval([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
end;

if ~plots, return;  end;          % do nothing, just return;

if nargin == 0,
  fprintf('VIEWER = %s\n',VIEWER);
  VIEWPARA
  return;
end;

if isstr(Ic) & strcmp(Ic,'set'),
%   fprintf('set VIEWPARA\n');
  VIEWPARA = {Omega,m,varargin{:}};
  for j = 1:length(VIEWPARA),
    if strcmp(VIEWPARA{j},'viewer'),
      VIEWER = VIEWPARA{j+1};
      break;
    end;
  end;
%   viewImage;
  return;
end;

% setup figure
fig = figure(fig);
subplot(sub(1),sub(2),sub(3)); cla;

if ~isempty(figname),
  set(fig,'numbertitle','off','name',sprintf('[JM-%d]: %s',fig,figname));
end;

% call the viewer
fh = feval(VIEWER,Ic,Omega,m,VIEWPARA{:});
str = sprintf('%s',name);
title(str);

if nargout>0,  varargout{1} = fh;  end;
return;
%==============================================================================

%==============================================================================
function ih = image2D(B,Omega,m,varargin)

colmap = [];

%% overwrite default parameter
for k=1:1:length(varargin)/2,
  %disp([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
  eval([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
end;

% compute pixel size and generate a regular grid
h  = Omega./m;
z1 = h(1)/2:h(1):Omega(1);
z2 = h(2)/2:h(2):Omega(2);

% we are now ready to use MATLAB's image
ih = image(z1,z2,reshape(B,m)');
axis image; axis xy                     % apply some nice axis
if ~isempty(colmap), colormap(colmap);  end;
return;

%==============================================================================
function ih = imagesc2D(B,Omega,m,varargin)

colmap = [];

%% overwrite default parameter
for k=1:1:length(varargin)/2,
  %disp([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
  eval([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
end;

% compute pixel size and generate a regular grid
h  = Omega./m;
z1 = h(1)/2:h(1):Omega(1);
z2 = h(2)/2:h(2):Omega(2);

% we are now ready to use MATLAB's image
ih = imagesc(z1,z2,reshape(B,m)');
axis image; axis xy                     % apply some nice axis
if ~isempty(colmap), colormap(colmap);  end;
return;

%==============================================================================
function p = volview(I,Omega,m,varargin)

I = reshape(I,m);
% I = flipdim(I,3);
isovalue    = 0;
viewpoint   = [-37.5,30];
facecolor   = .75*[1,1,1];
facealpha   = .8;

%% overwrite default parameter
for k=1:1:length(varargin)/2,
  %disp([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
  eval([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
end;

I(I<isovalue) = 0;

p(1) = patch(isosurface(I,isovalue));

switch computer,
case 'PCWIN',
  set(p(1),...
    'FaceColor',facecolor,...
    'EdgeColor','none');
otherwise,
  set(p(1),...
    'FaceColor',facecolor,...
    'EdgeColor','none',...
    'FaceAlpha',facealpha);
end;

isonormals(I,p(1))

p(2) = patch(isocaps(I,isovalue));

set(p(2),...
  'FaceColor','interp',...
  'EdgeColor','none');

view(viewpoint);
axis equal
axis([1,size(I,2),1,size(I,1),1,size(I,3)]);
camlight; lighting phong
material dull
drawnow
%==============================================================================

%==============================================================================
function ih = montage(I,Omega,m,varargin)

I = reshape(I,m);

threshold = 0;
framesx   = [];
framesy   = [];
direction = 'x';
number    = 'off';
colmap    = [];
%% overwrite default parameter
for k=1:1:length(varargin)/2,
  %disp([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
  eval([varargin{2*k-1},'=varargin{',int2str(2*k),'};']);
end;

if threshold ~= 0,
  I(I<threshold) = 0;
end;

I = squeeze(I);
l = length(size(I));

fig = gcf;
switch l,
  case 1,
    hh = plot(1:size(I,1),I);
  case 2,
    hh = subimage(I,colmap);
    if ~isempty(colmap), colormap(colmap); end;
  case 3,
    [m1,m2,m3] = size(I);

    switch direction,
      case 'x',
        aa = m1; b1 = m2; b2 = m3; I1 = 1:m2; I2 = 1:m3;
      case 'y',
        aa = m2; b1 = m1; b2 = m3; I1 = 1:m1; I2 = 1:m3;
      case 'z',
        aa = m3; b1 = m1; b2 = m2; I1 = 1:m1; I2 = 1:m2;

      otherwise, jmerror(direction);
    end;

    if isempty(framesx) & isempty(framesy),
      framesx = ceil(sqrt(aa));
      framesy = ceil(aa/framesx);
    elseif isempty(framesx),
      framesx = ceil(aa/framesy);
    else
      framesy = ceil(aa/framesx);
    end;

    %   fprintf('image is %dx%dx%d, frames [%d,%d]\n',...
    %     m1,m2,m3,framesx,framesy);

    B = zeros(b1*framesx,b2*framesy);

    a1 = 0; a2 = 0;
    for j=1:aa,
      switch direction,
        case 'x',  P = squeeze(I(j,:,:));
        case 'y',  P = squeeze(I(:,j,:));
        case 'z',  P = squeeze(I(:,:,j));
      end;
      B(a1+I1,a2+I2) = P;
      a2 = a2 + b2;
      if ~rem(j,framesy), a1 = a1 + b1; a2 = 0; end;
    end;

    ih = image(B);
    if ~isempty(colmap), colormap(colmap); end;
    hold on

    if strcmp(number,'on'),
      a1 = 0; a2 = 0;
      for j=1:aa,
        text(a2+10,a1+20,['#',int2str(j)],'fontsize',5,'color','r');
        a2 = a2 + b2;
        if ~rem(j,framesy), a1 = a1 + b1; a2 = 0; end;
      end;
    end;

    for j=0:framesx,
      plot(0.5+[0,framesy*b2],0.5+b1*j*[1,1],'b-');
    end;
    for j=0:framesy,
      plot(0.5+b2*j*[1,1],0.5+[0,framesx*b1],'b-');
    end;

    axis equal
    axis([0 framesy*b2+1, 0, framesx*b1+1])
    axis off
    hold off;
  otherwise,
    jmerror('dimension=?');
end;
% ==============================================================================

%==============================================================================


text:
function plot_z_coord(V,F)
%PLOT_Z_COORD  plot the z-coordinate on the input mesh
%
% plot_z_coord(V,F);
%
% Input:
%  V,F  mesh that is to be plotted

end

author: odedstein

repo: sgi-introduction-course


text:
function [LX, MX, PX] = mppca(X, no_dims, no_analyzers, tol, maxiter, minstd)
%MPPCA Runs EM algorithm and computes local factor analyzers
%
%   [LX, MX, PX] = mppca(X, no_dims, no_analyzers, tol, maxiter, minstd)
%
% Runs EM algorithm to determine coordinates of factor analyzers. The data
% is given in the DxN matrix X. no_dims indicates the number of dimensions that
% the local factor analyzers compute their embedding in. The number of
% factor analyzers that is used is given by no_analyzers. The variable tol indicates
% the tolreance in considering the EM as converged, whereas maxiter
% indicates the maximum number of iterations for the EM algorithm. The
% variable minstd sets the minimum standard deviation of the Gaussians that
% the EM algorithm fits.
% The function returns in LX the lowdimensional representations of X of all
% factor analyzers, in MX the means of the factors analyzers, and in PX the
% noise covariances.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    % Initialize some variables
    [D N] = size(X);                % size of data
    epsilon = 1e-9;                 % regularization parameter

    % Estimate minimum variance allowed
    minvar = minstd ^ 2;

    % Randomly initialize factor analyzers
    mm = mean(X, 2);
    ss = cov(X');
    try                                                 % for small problems
        cc = chol(ss);
        MX = bsxfun(@plus, cc' * randn(D, no_analyzers), mm);
        LX = minstd * randn(D, no_dims, no_analyzers);
        PX = 2 * mean(diag(cc)) * ones(D, 1);
    catch                                               % for large problems (or nearly singular covariance matrices)
        cc = std(X, [], 2);
        MX = bsxfun(@plus, bsxfun(@times, cc, randn(D, no_analyzers)), mm);
        LX = minstd * randn(D, no_dims, no_analyzers);
        PX = 2 * mean(cc) * ones(D, 1);
    end
    clear mm ss cc

    % Compute squared data
    X2 = X .^ 2;

    % Initialize some variables
    const = -D / 2 * log(2 * pi);
    lik = -Inf;
    R  = zeros(no_analyzers, N);
    czz = zeros(no_dims, no_dims, no_analyzers);
    zz  = zeros(no_dims, N, no_analyzers);

    % Run for maxiter iterations at max
    for i=1:maxiter

        % Progress bar
        if rem(i, 10) == 0
            fprintf('.');
        end

        % E step
        pii = 1 ./ PX;
        for k=1:no_analyzers
            l        = LX(:,:,k);                                                                       % select k-th local representation
            ltpi     = bsxfun(@times, pii, l)';
            ltpil    = ltpi * l;
            iltpil   = eye(no_dims) + ltpil;
            cc       = chol(iltpil);
            cci      = inv(cc);
            covz     = cci * cci';                                                                      % compute covariance
            delta    = X - MX(:,k * ones(1, N));
            meanz    = ((eye(no_dims) - ltpil * covz) * ltpi) * delta;
            czz(:,:,k) = covz;
            zz(:,:,k)  = meanz;                                                                         % update local representation
            R(k,:)    = -.5 * (pii' * (delta .* delta) - sum(meanz .* (iltpil * meanz), 1)) - ...       % update reponsibilie
                          sum(log(diag(cc)));
        end

        % Compute responsibilities of datapoints to the clusters
        R = R + const + .5 * sum(log(pi));
        R = exp(bsxfun(@minus, R, max(R, [], 1)));
        R = bsxfun(@rdivide, R, sum(R, 1));

        % Update likelihood of estimation
        oldlik = lik;
        lik = sum(max(R, [], 1) + log(sum(R, 1)));

        % Stop EM after convergence
        if abs(oldlik - lik) < tol
            break;
        end

        % M step
        PX = 0;
        for k=1:no_analyzers                         % Update all factor analyzers
            r    = R(k,:);
            z    = zz(:,:,k);
            rz   = bsxfun(@times, z, r);
            sr   = sum(r);
            srz  = sum(rz, 2);
            srxz = X * rz';
            srx  = X * r';
            m1   = [srxz srx];
            m2   = [sr * czz(:,:,k) + z * rz' srz; srz' sr];
            m1   = m1 / (m2 + (rand(size(m2)) * epsilon));     % random regularization to make sure that INV(m2) does not contain Infs
            LX(:,:,k) = m1(:,1:no_dims);
            MX(:,k) = m1(:,no_dims + 1);
            PX = PX + X2 * r' - sum(LX(:,:,k) .* srxz, 2) - MX(:,k) .* srx;
        end
        PX = max(minvar, PX / N);
        PX(:) = mean(PX);
    end

    % Done
    disp(' ');

repo: cluster


text:
function f = subsasgn(f, index, val)
%SUBSASGN   Chebfun SUBSASGN.
% ( )
%   F(X) = VAL assigns the values VAL at locations specified by X to the
%   CHEBFUN F. SIZE(X, 1) should be equal to LENGTH(VAL) and SIZE(X, 2) should
%   be the number of columns in F. SUBSASGN introduces new breakpoints
%   in F at points in X that were not originally in F.DOMAIN. See DEFINEPOINT
%   for further details.
%
% .
%   CHEBFUN properties are restricted, so F.PROP = VAL has no effect.
%
% {}
%   F{A, B} = G redefines the CHEBFUN F in the interval [A, B] using G. See
%   CHEBFUN/DEFINEINTERVAL for further details.
%
% See also SUBSREF, DEFINEPOINT, DEFINEINTERVAL.

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

% TODO: Error checking. Particularly integer indicies and quasimatrix indexing.

% TODO: Document for array-valued CHEBFUN objects and quasimatrices.

idx = index(1).subs;

switch index(1).type
    case '.'

        % [TODO]: Restrict access to this.
        f = builtin('subsasgn', f, index, val);

    case '()'

        if ( f(1).isTransposed && numel(idx) >= 2 )
            idx(1:2) = idx([2 1]);
        end
        if ( ischar(idx{1}) && strcmp(idx{1}, ':') )
            % Assign a column:
            f = assignColumns(f, idx{2}, val);
        else
            % Define a point value:
            f = definePoint(f, idx{1}, val);
        end

    case '{}'

        % Define an interval:
        f = defineInterval(f, [idx{:}], val);

    otherwise
        error('CHEBFUN:CHEBFUN:subsasgn:UnexpectedType',...
            ['??? Unexpected index.type of ' index(1).type]);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Recursive Least-Squares Algorithm with exponential weighting
%
% From S. Haykin, "Adaptive Filtering Theory (3rd Ed.)", Prentice Hall,
% Chapter 13.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef rls < linear_filter

    properties (GetAccess = 'public', SetAccess = 'private')
        lambda = .99; % forgetting factor
        c = 1E-4; % regularization
    end

    properties (GetAccess = 'public', SetAccess = 'private')
        P = []; % inverse autocorrelation matrix
        w = []; % filter coefficients
    end

    methods

        function obj = rls(parameters) % constructor
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(parameters)'
                    if ismember(fn,fieldnames(obj))
                        obj.(fn{1}) = parameters.(fn{1});
                    end
                end
            end
        end

        function y_est = evaluate(obj,x) % evaluate the algorithm
            if numel(obj.w)>0
                y_est = x*obj.w;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(obj,x,y) % train the algorithm
            if numel(obj.w)==0 % initialize
                m = length(x);
                obj.w = zeros(m,1);
                obj.P = obj.c\eye(m);
            end

            g = obj.P*x'/(obj.lambda+x*obj.P*x'); % gain vector
            err = y - x*obj.w; % instantaneous error
            obj.w = obj.w + g*err; % update filter coefficients
            obj.P = obj.lambda\(obj.P - g*x*obj.P); % update inv. autocorr.
        end

    end
end

repo: kafbox


text:
function varargout = quiver3(varargin)
%QUIVER3  3-D quiver plot of a SPHEREFUNV at data mapped by a SPHEREFUN.
%
%   This is not allowed and returns an error.  This function exists so that
%   the error message is meaningful to a SPHEREFUN user.
%
% See also SPHEREFUNV/QUIVER3.

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

error('CHEBFUN:SPHEREFUN:QUIVER3:notSupported',...
        'quiver3 is not yet supported for SPHEREFUN.');

end


text:
function [N,E] = rentian_scaling(A,XYZ,n)
%RENTIAN_SCALING    Physical Rentian scaling
%
% [N E] = rentian_scaling(A,XYZ,n)
%
% Physical Rentian scaling (or more simply Rentian scaling) is a property
% of systems that are cost-efficiently embedded into physical space. It is
% what is called a "topo-physical" property because it combines information
% regarding the topological organization of the graph with information
% about the physical placement of connections. Rentian scaling is present
% in very large scale integrated circuits, the C. elegans neuronal network,
% and morphometric and diffusion-based graphs of human anatomical networks.
% Rentian scaling is determined by partitioning the system into cubes,
% counting the number of nodes inside of each cube (N), and the number of
% edges traversing the boundary of each cube (E). If the system displays
% Rentian scaling, these two variables N and E will scale with one another
% in loglog space. The Rent's exponent is given by the slope of log10(E)
% vs. log10(N), and can be reported alone or can be compared to the
% theoretical minimum Rent's exponent to determine how cost efficiently the
% network has been embedded into physical space. Note: if a system displays
% Rentian scaling, it does not automatically mean that the system is
% cost-efficiently embedded (although it does suggest that). Validation
% occurs when comparing to the theoretical minimum Rent's exponent for that
% system.
%
% Inputs:
% 	A       MxM adjacency matrix
%           must be unweighted, binary, and symmetric.
% 	XYZ     Vector of node placement coordinates
%           must be Mx3 matrix, where M is the number of nodes.
% 	n       Number of partitions to compute. Each partition is a data
%           point. You want a large enough number to adequately estimate
%           the Rent's exponent.
%
% Outputs:
%	N       nx1 vector of the number of nodes in each of the n partitions.
%	E       nx1 vector of the number of edges crossing the boundary of each
%           partition.
%
% Subsequent Analysis:
%   Rentian scaling plots are then created by: figure; loglog(E,N,'*');
%
%	To determine the Rent's exponent, p, it is important not to use
%	partitions which may be affected by boundary conditions. In Bassett et
%	al. 2010 PLoS CB, only partitions with N<M/2 were used in the
%	estimation of the Rent's exponent. Thus, we can define N_prime =
%	N(find(N<M/2)) and E_prime = E(find(N<M/2)). Next we need to determine
%	the slope of Eprime vs. Nprime in loglog space, which is the Rent's
%	exponent. There are many ways of doing this with more or less
%	statistical rigor. Robustfit in MATLAB is one such option:
%       [b,stats] = robustfit(log10(N_prime),log10(E_prime))
%
%   Then the Rent's exponent is b(1,2) and the standard error of the
%   estimation is given by stats.se(1,2).
%
% Note: n=5000 was used in Bassett et al. 2010 in PLoS CB.
%
%
% Reference:
%   Danielle S. Bassett, Daniel L. Greenfield, Andreas Meyer-Lindenberg,
%   Daniel R. Weinberger, Simon W. Moore, Edward T. Bullmore. Efficient
%   physical embedding of topologically complex information processing
%   networks in brains and computer circuits. PLoS Comput Biol, 2010,
%   6(4):e1000748.
%
%
% Danielle Bassett, UCSB, 2010

% determine the number of nodes in the system
M = numel(XYZ(:,1));
% rescale coordinates so that they are all greater than unity
XYZn = XYZ-repmat(min(XYZ)-1,M,1);
% find the absolute minimum and maximum over all directions
nmax = max(max(XYZn));
nmin = min(min(XYZn));

% initialize variables
count = 0;
N = zeros(n,1);
E = zeros(n,1);

% create partitions, and count the number of nodes inside the partition (N) and the number of % edges traversing the boundary of the partition  (E)
while count<(n+1);
    % define cube end points
    randx = sort((1+nmax-nmin).*rand(2,1),'ascend');
    % find nodes in cube
    L = find(XYZn(:,1)>randx(1) & XYZn(:,1)<randx(2) & XYZn(:,2)>randx(1) & XYZn(:,2)<randx(2) & XYZn(:,3)>randx(1) & XYZn(:,3)<randx(2));
    if ~isempty(L)
        count = count+1;
        % count edges crossing the boundary of the cube
        E(count,1) = sum(sum(A(L,setdiff(1:M,L))));
        % count nodes inside of the cube
        N(count,1) = numel(L);
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Rectangle
% Concrete subclass of <ZeroVolShape.html |Shape|> representing a rectangular
% patch.

%%% Description
% |Rectangle| represents the shape of a rectangle on a plane.  It does not have
% a volume, but it is used to force a user-defined primary grid plane and four
% primary grid lines, which corresponds to the four sides of the rectangle, in
% <maxwell_run.html |maxwell_run|>.

%%% Construction
%  shape = Rectangle(normal_axis, intercept, rect)
%  shape = Rectangle(normal_axis, intercept, rect, dl_max)

% *Input Arguments*
%
% * |normal_axis|: axis normal to the plane of the rectangle.  It should be one
% of |Axis.x|, |Axis.y|, |Axis.z|.
% * |intercept|: location of the rectangle in the |normal_axis| direction.
% * |rect|: bounds of the rectangle in the plane.  If |normal_axis == Axis.y|,
% it should be in the format of |[zmin zmax; xmin xmax]|.
% * |dl_max|: maximum grid size allowed in the rectangle.  It can be either |[dx
% dy dz]| or a single real number |dl| for |dx = dy = dz|.  If unassigned,
% |dl_max = Inf| is used.  In the |normal_axis| direction, |dl_max| is
% meaningless.

%%% Example
%   % Create an instance of Rectangle.
%   shape = Rectangle(Axis.z, 100, [0 100; 0 100]);
%
%   % Use the constructed shape in maxwell_run().
%   [E, H] = maxwell_run({INITIAL ARGUMENTS}, 'OBJ', {'vacuum', 'none', 1.0}, shape, {REMAINING ARGUMENTS});

%%% See Also
% <Plane.html |Plane|>, <Line.html |Line|>, <Point.html |Point|>,
% <ZeroVolShape.html |ZeroVolShape|>, <maxwell_run.html |maxwell_run|>

classdef Rectangle < ZeroVolShape
	% Plane is a Shape for a plane.

	properties
		normal_axis
		intercept
		rect
	end

	methods
        function this = Rectangle(normal_axis, intercept, rect, dl_max)
			chkarg(istypesizeof(normal_axis, 'Axis'), '"normal_axis" should be instance of Axis.');
			chkarg(istypesizeof(intercept, 'real'), '"intercept" should be real.');
			chkarg(istypesizeof(rect, 'real', [Dir.count Sign.count]), ...
				'"rect" should be %d-by-%d array with real elements.', Dir.count, Sign.count);

			bound = NaN(Axis.count, Sign.count);
			bound(normal_axis, :) = [intercept intercept];
			[h, v] = cycle(normal_axis);
			bound([h v], :) = rect;

			box = Box(bound);
			plane = Plane(normal_axis, intercept);
			function level = lsf(x, y, z)
				chkarg(istypeof(x, 'real'), '"x" should be array with real elements.');
				chkarg(istypeof(y, 'real'), '"y" should be array with real elements.');
				chkarg(istypeof(z, 'real'), '"z" should be array with real elements.');
				chkarg(isequal(size(x), size(y), size(z)), '"x", "y", "z" should have same size.');

				% box.lsf(loc{:}) generates -Inf even for points slightly away
				% from the rectangle plane, because of the division by
				% s(normal_dir) == 0.  (Note: 1/0 = Inf.)  This is fine as a
				% behavior of lsf(), but then ZeroVolShape.lsf(loc{:}, true)
				% fails.  To solve this difficulty, modify lsf() not to divide
				% by zero.
				loc = {x, y, z};
				loc{normal_axis} = intercept(ones(size(x)));

				level = min(box.lsf(loc{:}), plane.lsf(x, y, z));  % intersection between box and plane
			end

			lprim = cell(1, Axis.count);
			for w = Axis.elems
				lprim{w} = box.bound(w,:);
			end

			if nargin < 4  % no dl_max
				super_args = {lprim, @lsf};
			else
				dl_max = expand2row(dl_max, Axis.count);
				dl_max(normal_axis) = Inf;  % dl_max is meaningless in normal direction
				super_args = {lprim, @lsf, dl_max};
			end

			this = this@ZeroVolShape(super_args{:});
			this.normal_axis = normal_axis;
			this.intercept = intercept;
			this.rect = rect;
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
classdef MeanSquareErrorNode < GraphNodeCost
    properties

    end

    methods
        function obj = MeanSquareErrorNode(weight)
            obj = obj@GraphNodeCost('MeanSquareError', weight);
        end

        function obj = forward(obj,prev_layers)
            obj = obj.preprocessingForward(prev_layers);

            output = prev_layers{1}.a;
            target = prev_layers{2}.a;
            diff = output - target;
            [D,T,N] = size(output);
            obj.a = real(0.5/T/N * sum( diff(:) .* conj(diff(:)) ));   % support both real and complex numbers
            obj = forward@GraphNodeCost(obj, prev_layers);
        end

        function obj = backward(obj,prev_layers, future_layers)
            input1 = prev_layers{1}.a;
            input2 = prev_layers{2}.a;
            diff = input1 - input2;
            [D,T,N] = size(input1);
            obj.grad{1} = diff / (T*N);
            obj.grad{2} = -obj.grad{1};
            obj = backward@GraphNodeCost(obj, prev_layers, future_layers);
        end

    end

end

author: singaxiong

repo: SignalGraph


text:
function [F,G] = snopt_callback(x,mmodel)

persistent model
if nargin > 1
    model = mmodel;
    return
end

if model.SimpleLinearObjective | model.SimpleQuadraticObjective
    [f,df] = fmincon_fun_liftlayer(x,model);
    [g,geq,dg,dgeq] = fmincon_con_liftlayer(x,model);
else
    [f,df,xevaled] = fmincon_fun_liftlayer(x,model);
    [g,geq,dg,dgeq] = fmincon_con_liftlayer(x,model,xevaled);
end

F = [f;g;geq];
n = length(x);
m = length(g);
p = length(geq);
G = [reshape(df,1,n);dg';dgeq'];
if ~isempty(model.A)
    F = [F;model.A*x - model.b];
    G = [G;model.A];
end
if ~isempty(model.Aeq)
    F = [F;model.Aeq*x - model.beq];
    G = [G;model.Aeq];
end
G = full(G(model.sparsityElements));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [out,im] = pyramid (infile,level,outfile)
% PYRAMID (infile,singvals,outfile)
% Image compression based on Gaussian Pyramid.
%  infile is input file name present in the current directory
%  level is a positive integer (the scaling factor)
%  outfile is output file name which will be created.
%
% It is possble to change the frequency response of the low pass
% gaussian filter changing the parameter "varianza" (see below).
%
% Example of use:
% [out,im]=pyramid('input_image.bmp',3,'output_image.bmp');
% im will be a double cell array with the 3 resampled layers
% and out will be the reconstructed image.
%
% The input image A can be RGB or GRAYSCALE.
%
% RGB case:
% If A is of class double, all values must be in the range [0,1],
% and A must be m-by-n-by-3.
% If A is of class uint16 or uint8, A must be m-by-n-by-3.
%
% GRAYSCALE case:
% If A is of class double, all values must be in the range [0,1],
% and the number of dimensions of A must be 2. If A is of class
% uint16 or uint8, the number of dimensions of A must be 2.
% uint16 or double.
%
%
% References:
%
% For more details concernings the algorithm implemented please visit the following links:
%
%  http://www.cs.ucf.edu/courses/cap6411/lect826h.PDF
%
%  http://www.wisdom.weizmann.ac.il/~maksimf/ex2/g3.html
%
%
% Please contribute if you find this software useful.
% Report bugs to luigi.rosa@tiscali.it
%
%
%*****************************************************************
% Luigi Rosa
% Via Centrale 27
% 67042 Civita di Bagno
% L'Aquila --- ITALY
% email  luigi.rosa@tiscali.it
% mobile +39 340 3463208
% http://utenti.lycos.it/matlab
%*****************************************************************
%
%
if (exist(infile)==2)
   a = imread(infile);
   figure('Name','Input image');
   imshow(a);
else
   warndlg('The file does not exist.',' Warning ');
   im=[];
   out=[];
   return
end

%------------------------------------
% The low-pass filter ---------------
%------------------------------------
% Gaussian Low Pass Filter ----------
%------------------------------------
x=[-16:1:16];
y=[-16:1:16];
dimx=size(x,2);
dimy=size(y,2);
varianza=2.3;
filtro=zeros(dimx,dimy);
for ii=1:dimx
    for jj=1:dimy
        esponente=exp(-(x(ii)^2+y(jj)^2)/(varianza));
        filtro(ii,jj)=esponente;
    end
end
% normalization
filtro=filtro/sum(sum(filtro));
%------------------------------------
%------------------------------------

if isgray(a)

    dvalue = double(a);
    dx     = size(dvalue,1);
    dy     = size(dvalue,2);
    dmin   = min(dx,dy);

    red_p  = floor(log2(dmin));
    red    = min(red_p,level);

    im=cell(red+1,1);

    % Image compression
    im{1}=dvalue;
    for ii=1:red
         % low-pass filter the image
         filtered        = conv2fft(dvalue,filtro,'same');
         % downsampling the image
         sottocampionato = dyaddown(filtered,1,'m');
         % save image
         im{ii+1}        = sottocampionato;
         % next step
         dvalue          = sottocampionato;
    end

    % Image reconstruction
    for ii=1:red
         % upsampling the image
         sovracampionato = dyadup(dvalue,0,'m');
         % low-pass filter the image (there is a scaling factor)
         filtered        = 4*conv2fft(sovracampionato,filtro,'same');
         % next step
         dvalue          = filtered;
    end

    if isa(a,'uint8')
        for ii=2:red+1
             testo=strcat('Image at level-',num2str(ii-1));
             figure('Name',testo);
             imshow(uint8(im{ii}));
             pause(0.6);
        end
        figure('Name','Reconstrunction');
        imshow(uint8(dvalue));
        pause(0.6);
        out=uint8(dvalue);

    end

    if isa(a,'uint16')
        for ii=2:red+1
             testo=strcat('Image at level-',num2str(ii-1));
             figure('Name',testo);
             imshow(uint16(im{ii}));
             pause(0.6);
        end
        figure('Name','Reconstrunction');
        imshow(uint16(dvalue));
        pause(0.6);
        out=uint16(dvalue);

    end

    if isa(a,'double')
        for ii=2:red+1
             testo=strcat('Image at level-',num2str(ii-1));
             figure('Name',testo);
             imshow(uint32(im{ii}));
             pause(0.6);
        end
        figure('Name','Reconstrunction');
        imshow((dvalue));
        pause(0.6);
        out=uint32(dvalue);

   end

    imwrite(out, outfile);
    return;
end

%------------------------------------------------------
if isrgb(a)
    dvalue = double(a);
    dx     = size(dvalue,1);
    dy     = size(dvalue,2);
    dmin   = min(dx,dy);

    red_p  = floor(log2(dmin));
    red    = min(red_p,level);

    im=cell(red+1,1);

    % Image compression
    im{1}=dvalue;
    for ii=1:red
         % low-pass filter the image
         filtered_r        = conv2fft(dvalue(:,:,1),filtro,'same');
         filtered_g        = conv2fft(dvalue(:,:,2),filtro,'same');
         filtered_b        = conv2fft(dvalue(:,:,3),filtro,'same');
         % downsampling the image
         sottocampionato_r = dyaddown(filtered_r,1,'m');
         sottocampionato_g = dyaddown(filtered_g,1,'m');
         sottocampionato_b = dyaddown(filtered_b,1,'m');
         % save image
         clear sottocampionato;
         sottocampionato(:,:,1)=sottocampionato_r;
         sottocampionato(:,:,2)=sottocampionato_g;
         sottocampionato(:,:,3)=sottocampionato_b;
         im{ii+1}        = sottocampionato;
         % next step
         dvalue          = sottocampionato;
    end

    % Image reconstruction
    for ii=1:red
         % upsampling the image
         sovracampionato_r = dyadup(dvalue(:,:,1),0,'m');
         sovracampionato_g = dyadup(dvalue(:,:,2),0,'m');
         sovracampionato_b = dyadup(dvalue(:,:,3),0,'m');
         % low-pass filter the image (there is a scaling factor)
         filtered_r        = 4*conv2fft(sovracampionato_r,filtro,'same');
         filtered_g        = 4*conv2fft(sovracampionato_g,filtro,'same');
         filtered_b        = 4*conv2fft(sovracampionato_b,filtro,'same');
         % next step
         clear dvalue;
         dvalue(:,:,1)          = filtered_r;
         dvalue(:,:,2)          = filtered_g;
         dvalue(:,:,3)          = filtered_b;
    end

    if isa(a,'uint8')
        for ii=2:red+1
             testo=strcat('Image at level-',num2str(ii-1));
             figure('Name',testo);
             imshow(uint8(im{ii}));
             pause(0.6);
        end
        figure('Name','Reconstrunction');
        imshow(uint8(dvalue));
        pause(0.6);
        out=uint8(dvalue);

    end

    if isa(a,'uint16')
        for ii=2:red+1
             testo=strcat('Image at level-',num2str(ii-1));
             figure('Name',testo);
             imshow(uint16(im{ii}));
             pause(0.6);
        end
        figure('Name','Reconstrunction');
        imshow(uint16(dvalue));
        pause(0.6);
        out=uint16(dvalue);

    end

    if isa(a,'double')
        for ii=2:red+1
             testo=strcat('Image at level-',num2str(ii-1));
             figure('Name',testo);
             imshow((im{ii}));
             pause(0.6);
        end
        figure('Name','Reconstrunction');
        imshow((dvalue));
        pause(0.6);
        out=(dvalue);

   end

    imwrite(out, outfile);
    return;

end

author: Sable

repo: mcbench-benchmarks


text:
function box_display_test01 ( )

%*****************************************************************************80
%
%% BOX_DISPLAY_TEST01 plots total degree index sets.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_DISPLAY_TEST01:\n' );
  fprintf ( 1, '  Plot total degree index sets.\n' );
  fprintf ( 1, '\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 1';
  box_display ( m, n, @(x,y)x+y<=0, @(x,y)x+y<=1, title_string );
  print ( '-dpng', 'td1.png' );
  fprintf ( 1, '  Created "td1.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 2';
  box_display ( m, n, @(x,y)x+y<=1, @(x,y)x+y<=2, title_string );
  print ( '-dpng', 'td2.png' );
  fprintf ( 1, '  Created "td2.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 3';
  box_display ( m, n, @(x,y)x+y<=2, @(x,y)x+y<=3, title_string );
  print ( '-dpng', 'td3.png' );
  fprintf ( 1, '  Created "td3.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 4';
  box_display ( m, n, @(x,y)x+y<=3, @(x,y)x+y<=4, title_string );
  print ( '-dpng', 'td4.png' );
  fprintf ( 1, '  Created "td4.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 5';
  box_display ( m, n, @(x,y)x+y<=4, @(x,y)x+y<=5, title_string );
  print ( '-dpng', 'td5.png' );
  fprintf ( 1, '  Created "td5.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 6';
  box_display ( m, n, @(x,y)x+y<=5, @(x,y)x+y<=6, title_string );
  print ( '-dpng', 'td6.png' );
  fprintf ( 1, '  Created "td6.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 7';
  box_display ( m, n, @(x,y)x+y<=6, @(x,y)x+y<=7, title_string );
  print ( '-dpng', 'td7.png' );
  fprintf ( 1, '  Created "td7.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 8';
  box_display ( m, n, @(x,y)x+y<=7, @(x,y)x+y<=8, title_string );
  print ( '-dpng', 'td8.png' );
  fprintf ( 1, '  Created "td8.png".\n' );

  m = 12;
  n = 12;
  title_string = 'Total Degree <= 9';
  box_display ( m, n, @(x,y)x+y<=8, @(x,y)x+y<=9, title_string );
  print ( '-dpng', 'td9.png' );
  fprintf ( 1, '  Created "td9.png".\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = r8_cos ( x )

%*****************************************************************************80
%
%% R8_COS evaluates the cosine of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the cosine of X.
%
  persistent ntsn
  persistent pi2
  persistent pi2rec
  persistent pihi
  persistent pilo
  persistent pirec
  persistent sincs
  persistent xmax
  persistent xsml
  persistent xwarn
%
%  pihi + pilo = pi.  pihi is exactly representable on all machines
%  with at least 8 bits of precision.  whether it is exactly
%  represented depends on the compiler.  this routine is more
%  accurate if it is exactly represented.
%
  pi2 = 1.57079632679489661923132169163975;
  pi2rec = 0.63661977236758134307553505349006;
  pihi = 3.140625;
  pilo = 9.6765358979323846264338327950288E-04;
  pirec = 0.31830988618379067153776752674503;

  if ( isempty ( ntsn ) )

    sincs = [ ...
      -0.374991154955873175839919279977323464, ...
      -0.181603155237250201863830316158004754, ...
      +0.005804709274598633559427341722857921, ...
      -0.000086954311779340757113212316353178, ...
      +0.000000754370148088851481006839927030, ...
      -0.000000004267129665055961107126829906, ...
      +0.000000000016980422945488168181824792, ...
      -0.000000000000050120578889961870929524, ...
      +0.000000000000000114101026680010675628, ...
      -0.000000000000000000206437504424783134, ...
      +0.000000000000000000000303969595918706, ...
      -0.000000000000000000000000371357734157, ...
      +0.000000000000000000000000000382486123, ...
      -0.000000000000000000000000000000336623, ...
      +0.000000000000000000000000000000000256 ]';

    ntsn = r8_inits ( sincs, 15, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( 2.0 * r8_mach ( 3 ) );
    xmax = 1.0 / r8_mach ( 4 );
    xwarn = sqrt ( xmax );

  end

  absx = abs ( x );
  y = absx + pi2;

  if ( xmax < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_COS - Warning!\n' );
    fprintf ( 1, '  No precision because |X| is big.\n' );
    value = 0.0
    return
  end

  if ( xwarn < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_COS - Warning!\n' );
    fprintf ( 1, '  Answer < half precision because |X| is big.\n' );
  end

  value = 1.0;

  if ( absx < xsml )
    return
  end

  xn = r8_aint ( y * pirec + 0.5 );
  n2 = r8_aint ( mod ( xn, 2.0 ) + 0.5 );
  xn = xn - 0.5;
  f = ( absx - xn * pihi ) - xn * pilo;

  xn = 2.0 * ( f * pi2rec ) * ( f * pi2rec ) - 1.0;

  value = f + f * r8_csevl ( xn, sincs, ntsn );

  if ( n2 ~= 0 )
    value = - value;
  end

  if ( value < -1.0 )
    value = -1.0;
  elseif ( 1.0 < value )
    value = 1.0;
  end

  return
end


text:
%Develop by Renoald
%Email:Renoald@live.com
%11/2/2013
%This function can used to view 3D point clouds
%Format of 3D point clouds
%Xi Yi Zi Ri Gi Bi
function ViewPointCloud(file)
%read File
p=dlmread(file);
[i,j]=size(p);
dist=ones(i,1);
if j==6
    %scatter3(p(:,1),p(:,2),p(:,3),'Marker','o')

else
    for i=1 : i
    dist(i,1)=Dist3D(0,0,0,p(i,1),p(i,2),p(i,3));
end
    figure(1);
    %Plot point clouds
    fname=strcat('Number of points:',num2str(i));
    scatter3(p(:,1),p(:,2),p(:,3),1,dist(:,1))
    title(fname);
    xlabel('x axis');
    ylabel('y axis');
    zlabel('z axis');
end


text:
function[]=ylin
% YLIN   Sets y-axis scale to linear.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 J.M. Lilly --- type 'help jlab_license' for details

set(gca,'yscale','linear')

%Update for compatibility with ZOOL
if isappdata(gcf,'ylimitslock')
   if getappdata(gcf,'ylimitslock')
      h=gca;
      hall=axeshandles;
      for i=1:length(hall)
         axes(hall(i))
         set(gca,'yscale','linear')
      end
      axes(h)
   end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%% Example
% An example script, it is used to show how to use the different
% functions of the KUKA iiwa matlab toolbox

% First run the following script in Matlab
% Then start the client on the KUKA iiwa controller

% Mohammad SAFEEA, 24th of April 2018

close all,clear all,clc
warning('off');
ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
else

      %% Get position roientation of end effector

      setBlueOn(t_Kuka); % turn on blue light

          %% move to initial position
jPos={0,pi*20/180,0,-pi*70/180,0,pi*60/180,0}; % initial confuguration
relVel=0.15; % relative velocity
movePTPJointSpace( t_Kuka , jPos, relVel); % point to point motion in joint space

%% Linear relative motion of end effector, relative to base frame
deltaX=0;deltaY=0;deltaZ=100.;
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;

vel=250; % linear velocity of end effector, mm/sec

movePTPLineEefRelBase( t_Kuka , Pos, vel);
Pos{3}=-deltaZ;
movePTPLineEefRelBase( t_Kuka , Pos, vel);

%% Linear relative motion of end effector, relative to EEF frame
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;

vel=150; % linear velocity of end effector, mm/sec

movePTPLineEefRelEef( t_Kuka , Pos, vel);
Pos{3}=-deltaZ;
movePTPLineEefRelEef( t_Kuka , Pos, vel);

setBlueOff(t_Kuka); % turn off blue light

      %% turn off the server
       net_turnOffServer( t_Kuka );

       fclose(t_Kuka);
end
warning('on');

repo: KST-Kuka-Sunrise-Toolbox


text:
function qr=rotqc2qr(qc)
%ROTQC2QR converts a matrix of complex quaternion row vectors into real form
%
% Inputs:
%
%     QC(2m,n)   mxn matrix of complex-valued quaternions
%
% Outputs:
%
%     QR(4m,n)   mxn matrix of real-valued quaternions
%
% The complex-valued quaternion [r+j*b  a+j*c] becomes [r a b c]

%
%      Copyright (C) Mike Brookes 2000-2006
%      Version: $Id: rotqc2qr.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[m,n]=size(qc);
i=(1:2:2*m)-mod(0:m-1,2);
qr=zeros(2*m,n);
qr(i,:)=real(qc);
qr(i+2,:)=imag(qc);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%function chol_comparison
% Test what is the best way to evaluate Cholesky(-like) decomposition.

n = 1000;

% Large and ill-conditioned covariance matrix
C = gpcov(log(100), 1:n, 1:n);

% Regular Cholesky with diagonal addition
tic
Lchol = chol(C + 100*eps*eye(n), 'lower');
Lchol(1:10, 1:10)
toc

% LDL with negative values to zero
tic
[Lldl, Dldl] = ldl(C);
%Dldl(Dldl<0) = 0;
%Lldl = Lldl * sqrt(Dldl);
d = diag(Dldl);
d(d<0) = 0;
Lldl(1:10,1:10)
Lldl = bsxfun(@times, Lldl, sqrt(d(:))');
toc

% My LDL with negative values to zero
tic
[Lmyldl, Dmyldl] = ldl(C);
%Dldl(Dldl<0) = 0;
%Lldl = Lldl * sqrt(Dldl);
d = diag(Dmyldl);
d(d<0) = 0;
Lmyldl(1:10,1:10)
Lmyldl = bsxfun(@times, Lmyldl, sqrt(d(:))');
toc
%return

% Matlab's cholcov based on eigenvalue decomposition
% Note that this solution does not, in general, give triangular matrix!!!
% Also, this is very slow!
tic
Lcholcov = cholcov(C)';
Lcholcov = [Lcholcov, zeros(n, n-cols(Lcholcov))];
toc

% Errors:
error_chol = norm(C - Lchol*Lchol')
error_ldl = norm(C - Lldl*Lldl')
error_myldl = norm(C - Lmyldl*Lmyldl')
%error_ldl = norm(C - Lldl*Dldl*Lldl')
error_cholcov = norm(C - Lcholcov*Lcholcov')

% What other measures could be used?
% How accurate these methods are in solving linear equations?

x = randn(n,1);
y = C * x;

opts.LT = true;
opts.TRANSA = false;
x_chol = linsolve(Lchol, y, opts);
opts.TRANSA = true;
x_chol = linsolve(Lchol, x_chol, opts);

opts.LT = true;
opts.TRANSA = false;
x_ldl = linsolve(Lldl, y, opts);
opts.TRANSA = true;
x_ldl = linsolve(Lldl, x_ldl, opts);

opts.LT = true;
opts.TRANSA = false;
x_myldl = linsolve(Lmyldl, y, opts);
opts.TRANSA = true;
x_myldl = linsolve(Lmyldl, x_myldl, opts);

opts.LT = false;
opts.TRANSA = false;
x_cholcov = linsolve(Lcholcov, y, opts);
opts.TRANSA = true;
x_cholcov = linsolve(Lcholcov, x_cholcov, opts);

error_chol = norm(x - x_chol)
error_ldl = norm(x - x_ldl)
error_myldl = norm(x - x_myldl)
error_cholcov = norm(x - x_cholcov)

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [res] = fiff_pick_types_evoked(orig,meg,eeg,stim,include,exclude)
%
% [res] = fiff_pick_types_evoked(orig,meg,eeg,stim,include,exclude)
%
% Pick desired channels types from evoked-response data
%
% orig      - The original data
% meg       - Include MEG channels
% eeg       - Include EEG channels
% stim      - Include stimulus channels
% include   - Additional channels to include (if empty, do not add any)
% exclude   - Channels to exclude (if empty, do not exclude any)
%
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.10  2008/06/12 19:59:19  msh
%   Fixed the help text
%
%   Revision 1.9  2007/12/23 17:01:24  msh
%   Fixed error handling multiple datasets
%
%   Revision 1.8  2006/04/26 00:14:24  msh
%   Return empty structure from fiff_pick_channels evoked and fiff_pick_types_evoked if there is no match.
%   Accuracy checking was incorrect in mne_add_coil_defs
%
%   Revision 1.7  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.6  2006/04/21 11:33:17  msh
%   Added raw segment reading routines
%
%   Revision 1.5  2006/04/20 23:33:20  msh
%   Separated fiff_pick_channels and fiff_pick_types from the _evoked versions
%
%   Revision 1.4  2006/04/17 11:52:15  msh
%   Added coil definition stuff
%
%   Revision 1.3  2006/04/15 12:21:00  msh
%   Several small improvements
%
%   Revision 1.2  2006/04/14 15:49:49  msh
%   Improved the channel selection code and added ch_names to measurement info.
%
%   Revision 1.1  2006/04/14 00:45:42  msh
%   Added channel picking and fiff_invert_transform
%
%

me='MNE:fiff_pick_types_evoked';

if nargin == 5
    exclude = [];
elseif nargin == 4
    include = [];
    exclude = [];
elseif nargin == 3
    include = [];
    exclude = [];
    stim    = false;
elseif nargin == 2
    include = [];
    exclude = [];
    stim    = false;
    eeg     = false;
elseif nargin ~= 6
    error(me,'Incorrect number of arguments');
end

sel = fiff_pick_types(orig.info,meg,eeg,stim,include,exclude);

if isempty(sel)
    res = [];
    fprintf(1,'Warning : No channels match the selection.\n');
    return;
end

res = orig;
res.info.chs      = res.info.chs(sel);
res.info.ch_names = res.info.ch_names(sel);
%
%   Create the reduced data set
%
for k = 1:length(res.evoked)
    res.evoked(k).epochs = res.evoked(k).epochs(sel,:);
end
res.info.nchan    = length(sel);

return;

end


text:
%%*****************************************************************
%% NTdirfun: compute (dX,dZ), given dy, for the NT direction.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

    function [dX,dy,dZ] = NTdirfun(blk,At,par,Rd,EinvRc,xx,m);

    global solve_ok

    dX = cell(size(blk,1),1); dZ = cell(size(blk,1),1); dy = [];
    if (any(isnan(xx)) | any(isinf(xx)))
       solve_ok = 0;
       fprintf('\n  linsysolve: solution contains NaN or inf.');
       return;
    end
%%
    dy = xx(1:m);
    count = m;
%%
    for p=1:size(blk,1)
       pblk = blk(p,:);
       if strcmp(pblk{1},'l')
          %%dZ{p} = Rd{p} - At{p}*dy;
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),[],[],dy));
          tmp   = par.dd{p}.*dZ{p};
          dX{p} = EinvRc{p} - tmp;
       elseif strcmp(pblk{1},'q')
          %%dZ{p} = Rd{p} - At{p}*dy;
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),[],[],dy));
          tmp = par.dd{p}.*dZ{p} + qops(pblk,qops(pblk,dZ{p},par.ee{p},1),par.ee{p},3);
          dX{p} = EinvRc{p} - tmp;
       elseif strcmp(pblk{1},'s')
          %%dZ{p} = Rd{p} - smat(pblk,At{p}*dy(par.permA(p,:)),par.isspAy(p));
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),par.permA(p,:),par.isspAy(p),dy));
          tmp   = Prod3(pblk,par.W{p},dZ{p},par.W{p},1);
          dX{p} = EinvRc{p}-tmp;
       elseif strcmp(pblk{1},'u');
          n = sum(pblk{2});
          dZ{p} = zeros(n,1);
          dX{p} = xx(count+[1:n]);
          count = count + n;
       end
    end
%%*******************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
nmpcdata=zeros(3,6);%Store data in matrix for nmpc and impc, first row: mean of computing time; second row: standard deviation of computing time; third row: infeasible rate
impcdata=zeros(3,6);
for ii = 1:6
i=ii*4;
filenm1 = ['timecom' num2str(i) '.mat'];
filenm2 = ['feasibility' num2str(i) '.mat'];
load(filenm1);
load(filenm2);
distnmpc = fitdist(nmpcplot1','Normal');
distimpc = fitdist(impcplot1','Normal');
mu1=distnmpc.mu;%Get mean of sample of computing time for NMPC-DCBF
mu2=distimpc.mu;%Get mean of sample of computing time for iMPC-DCBF
sigma1=distnmpc.sigma;%Get variance of sample of computing time for NMPC-DCBF
sigma2=distimpc.sigma;%Get variance of sample of computing time for iMPC-DCBF
nmpcdata(1,ii)=mu1;
impcdata(1,ii)=mu2;
nmpcdata(2,ii)=sigma1;
impcdata(2,ii)=sigma2;
nmpcdata(3,ii)=nmpcinf;
impcdata(3,ii)=impcinf;
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF


text:
classdef PoissonD
%%POISSOND Functions to handle the Poisson distribution.
%Implemented methods are: mean, var, PMF, CDF, momentGenFun, cumGenFun,
%                         rand
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(lambda)
%%MEAN Obtain the mean of the Poisson distribution for the specified mean
%      parameter.
%
%INPUTS: lambda The mean (and variance) of the Poisson distribution.
%
%OUTPUTS: val The mean of the Poisson distribution under consideration.
%
%The Poisson distribution is discussed in Chapter 2.8 of [1].
%
%EXAMPLE:
%We verify the computed mean by comparing it to the sample mean.
% lambda=10;
% meanVal=PoissonD.mean(lambda)
% numSamp=1e5;
% sampMeanVal=mean(PoissonD.rand([1,numSamp],lambda))
%One will find both values are about 10.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=lambda;
end

function val=var(lambda)
%%VAR Obtain the variance of the Poisson distribution for the specified
%     mean parameter
%
%INPUTS: lambda The mean (and variance) of the Poisson distribution,
%
%OUTPUTS: val The variance of the Poisson distribution.
%
%The Poisson distribution is discussed in Chapter 2.8 of [1].
%
%EXAMPLE:
%We verify that the computed variance by comparing it to the sample
%variance.
% lambda=10;
% varVal=PoissonD.var(lambda)
% numSamp=1e5;
% sampVarVal=var(PoissonD.rand([1,numSamp],lambda))
%One will find both values are about 10.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=lambda;
end

function val=PMF(x,lambda)
%%PMF Evaluate the Poisson probability mass function (PMF) at given points
%     with a specified mean parameter.
%
%INPUTS: x The point(s) at which the Poisson PMF is to be evaluated. x is
%          an integer. For nonzero PMF values, x>=0.
%   lambda The mean (and variance) of the Poisson distribution under
%          consideration.
%
%OUTPUTS: val The value(s) of the Poisson PMF with mean lambda.
%
%The Poisson distribution is discussed in Chapter 2.8 of [1].
%
%EXAMPLE:
%Here, we validate the PMF by generating random samples and comparing the
%PMF plot with a histogram of the random samples.
% lambda=6;
% numSamples=1e5;
% figure(1)
% clf
% histogram(PoissonD.rand([numSamples,1],lambda),'BinWidth',1,'Normalization','pdf')
% hold on
% x=0:16;
% vals=PoissonD.PMF(x,lambda);
% stem(x,vals,'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('PDF(x)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the histogram matches well with the plot.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(x));
    sel=(x>=0);

    val(sel)=(lambda.^x(sel)./factorial(x(sel)))*exp(-lambda);
end

function [FVal,pVal]=CDF(x,lambda)
%%CDF Evaluate the cumulative distribution function of the Poisson
%     distribution at desired points.
%
%INPUTS: x The integer point(s) at which the Poisson CDF is evaluated.
%   lambda The mean (and variance) of the Poisson distribution.
%
%OUTPUTS: FVal The CDF of the Poisson distribution evaluated at the desired
%              point(s).
%         pVal The value of the PMF of the Poisson distribution at the
%              desired point(s).
%
%The recursion is from Chapter 4.2 of [1].
%
%EXAMPLE:
%We validate the CDF value by comparing it to a value computed from random
%samples.
% lambda=6;
% x=3;
% numSamples=1e5;
% prob=PoissonD.CDF(x,lambda)
% probSamp=mean(PoissonD.rand([numSamples,1],lambda)<=x)
%One will find the values ot both be about 0.151.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    FVal=zeros(size(x));
    pVal=zeros(size(x));

    x=floor(x);%In case decimals were passed.
    numEls=numel(x);

    for curEl=1:numEls
        p=exp(-lambda);
        F=p;
        for I=1:fix(x(curEl))
            p=p*lambda/I;
            F=F+p;
        end
        FVal(curEl)=F;
        pVal(curEl)=p;
    end
    FVal(x<0)=0;
    pVal(x<0)=0;
end

function momentVal=momentGenFun(lambda,numDerivs,t)
%%MOMENTGENFUN Evaluate the moment generating function (or one of its
%              derivatives) of the Poisson distribution. Taking the kth
%              derivative of the moment generating function and evaluating
%              it at t=0 provides the kth noncentral moment of the
%              distribution.
%
%INPUTS: lambda The mean (and variance) of the Poisson distribution.
%     numDerivs The number of derivatives to take with respect to the
%               argument of the moment generating function. numDerivs>=0.
%             t The numPointsX1 or 1XnumPoints vector of points where the
%               moment generating function should be evaluated. If this
%               parameter is omitted or an empty matrix is passed, the
%               default value of 0 is used.
%
%OUTPUTS: momentVal A numPointsX1 vector of the values of the derivatives
%                   of the moment generating function given at the points
%                   in t or at t=0 if t is omitted.
%
%The moment generating function of a random variable is defined to be
%E(exp(t*x)) where E is the expected value operator, x is the random
%variable and t is a real parameter. It can be shown that the moment
%generating function of the Poisson distribution is
%f^0(t)=exp(lambda*(exp(t)-1))
%By noting the pattern in differentiation (using the chain rule), it can be
%seen that the derivatives of the moment generating function can be
%recurively expressed. Specifically, if f^n(t) is the nth derivative, the
%value for n>0 is
%f^n(t)=lambda*exp(t)*(sum_{k=0}^(n-1) binomial(n-1,k)*f^k(t))
%This function evaluates the moments using the above recursion. The
%binomial terms are computed recurively using the recursion inherent to
%Pascal's triangle.
%
%Note that because the function recursively accumulates values for the
%derivatives, a lot of memory will be used if numDerivs is very large.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

if(nargin<3||isempty(t))
    t=0;
end

numT=length(t);
momentVal=zeros(numT,1);%Allocate Space
for curT=1:numT
    tVal=t(curT);

    fVals=zeros(numDerivs+1,1);
    fVals(1)=exp(lambda*(exp(tVal)-1));%The moment generating function value.

    pascalTriangleRowOld=zeros(numDerivs+1,1);
    pascalTriangleRowOld(1)=1;
    expCoeff=exp(tVal)*lambda;
    for curDeriv=1:numDerivs
        %Get the next row of Pascal's triangle
        pascalTriangleRow=zeros(numDerivs+1,1);
        pascalTriangleRow(1)=1;
        for i=2:(curDeriv-1)
            pascalTriangleRow(i)=pascalTriangleRowOld(i-1)+pascalTriangleRowOld(i);
        end
        pascalTriangleRow(curDeriv)=1;

        fVals(curDeriv+1)=expCoeff*sum(pascalTriangleRow.*fVals);
        pascalTriangleRowOld=pascalTriangleRow;
    end

    momentVal(curT)=fVals(end);
end

end

function cumVal=cumGenFun(lambda,numDerivs,t)
%%CUMGENFUN Evaluate the cumulant generating function (or one of its
%           derivatives) of the Poisson distribution. The cumulant
%           generating function is the natural logarithm of the moment
%           generating function.
%
%INPUTS: lambda The mean (and variance) of the Poisson distribution.
%     numDerivs The number of derivatives to take with respect to the
%               argument of the cumulant generating function, numDerivs>=0.
%             t The numPointsX1 or 1XnumPoints vector of points where the
%               moment generating function should be evaluated. If this
%               parameter is omitted or an empty matrix is passed, the
%               default value of 0 is used.
%
%OUTPUTS: cumVal A numPointsX1 or 1XnumPoints vector of the values of the
%                derivatives of the cumulant generating function given at
%                the points in t or at t=0 if t is omitted.
%
%The moment generating function of a random variable is defined to be
%E(exp(t*x)) where E is the expected value operator, x is the random
%variable and t is a real parameter. The cumulant generating function is
%defined as the natural logarithm of the moment generating function. It can
%be shown that the moment generating function of the Poisson distribution
%is
%E(exp(t*x))=exp(lambda*(exp(t)-1))
%Thus, the cumulant generating function is just
%lambda*(exp(t)-1)
%All derivatives of the cumulant generating function are consequently
%exp(t)*lambda.
%
%September 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

if(nargin<3||isempty(t))
    t=0;
end

if(numDerivs==0)
    cumVal=lambda*(exp(t)-1);
else
    cumVal=exp(t)*lambda;
end

end

function vals=rand(N,lambda)
%%RAND Generate Poisson random variables with a given mean.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%   lambda The mean (and variance) of the Poisson distribution.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated Poisson random variables.
%
%The Poisson random variables are generated according to the simple
%method given in Chapter 4.2 of [1].
%
%The algorithm can handle small and large values of lambda. However, there
%is always a possibility that the value desired might be on the tail
%end of the distribution and be very slow. Thus, the number of iterations
%(the offset from the mean) is set to the maximum of 10 and 10* the
%variance.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N, N];
    else
        dims=N;
    end

    maxIter=max(10,10*lambda);

    vals=zeros(dims);
    for curRow=1:dims(1)
        for curCol=1:dims(2)
            I=fix(lambda);
            [F,p]=PoissonD.CDF(I,lambda);
            U=rand(1);

            if(U<=F)
            %This means that the random number being generated is less than
            %or equal to I and we must search downward.
                while(U<F&&I>-1)
                    F=F-p;
                    p=p*(I/lambda);
                    I=I-1;
                end
                I=I+1;
            else
            %This means that the random number being generated is greater
            %than I and we must search upward starting from I+1.
                curIter=0;
                while(U>F&&curIter<maxIter)
                    p=p*lambda/(I+1);
                    F=F+p;
                    I=I+1;
                    curIter=curIter+1;
                end
            end

            vals(curRow ,curCol)=I;
        end
    end
end
end
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:
% The contact lib by Mathworks, may change to another lib written by
% Professor Hartmut Geyer
addpath(genpath(fullfile(pwd,'contat_lib')));
%parameters for ground plane geometry and contract
ground.stiff = 5e4;
ground.damp = 1000;
ground.height = 0.4;

%robot body size
body.x_length = 0.6;
body.y_length = 0.3;
body.z_length = 0.15;
body.shoulder_size = 0.07;
body.upper_length = 0.16;
body.lower_length = 0.25;
body.foot_radius = 0.035;
body.shoulder_distance = 0.2;
body.max_stretch = body.upper_length + body.lower_length;
body.knee_damping = 0.1;

% parameters for leg control
ctrl.pos_kp = 30;
ctrl.pos_ki = 0;
ctrl.pos_kd = 0;
ctrl.vel_kp = 4.4;
ctrl.vel_ki = 0;
ctrl.vel_kd = 0.0;

% parameters for high level plan
planner.touch_down_height = body.foot_radius; % from foot center to ground

planner.stand_s = 0;
planner.stand_u = 30*pi/180;
planner.stand_k = -60*pi/180;
stance_pos = forward_kinematics(planner.stand_s, planner.stand_u, planner.stand_k, body);
planner.stand_height = stance_pos(3);

planner.flight_height = 1.1*planner.stand_height;% when leg is flying in the air

% gait control parameters
planner.time_circle = 3;
planner.swing_ang = 12/180*pi;
planner.init_shake_ang = 3/180*pi;

%some gait control goals
planner.tgt_body_ang = 0;
planner.tgt_body_vx = 0.3;  % tested 0.15-0.45
planner.Ts = 0.1; % for x_direction leg place
planner.Kv = 0.3;  % for x_dreiction leg place

planner.y_Ts = 0.21; % for y_direction leg place
planner.y_Kv = -0.34;  % for y_direction leg place

% state transition thredsholds
planner.state0_vel_thres = 0.05;
% the sample time of the ctrl block is set to be 0.0025 (400Hz);
planner.state0_trans_thres = 300; %(0.75 seconds)
planner.state0_swing_ang  = 10*pi/180;
planner.state0_swing_T = 1200;
planner.state12_trans_speed = 0.2;
planner.leg_swing_time = 70; %(use 0.25s to generate a motion)

%robot weight
body_weight = 600*body.x_length*body.y_length*body.z_length;
leg_density = 660;
should_weight = leg_density*0.07^3;
upperleg_weight = leg_density*0.04*0.04*body.upper_length;
lowerleg_weight = leg_density*0.04*0.04*body.lower_length;
foot_weight = 1000*4/3*pi*body.foot_radius^3;

total_weight = body_weight + 4*(should_weight+upperleg_weight+lowerleg_weight+foot_weight);

body_inertia = diag([0.151875;0.516375;0.6075]);

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function x = sparsify( x, mode )

global cvx___
narginchk(2,2);
persistent remap

%
% Check mode argument
%

if ~ischar( mode ) || size( mode, 1 ) ~= 1,
    error( 'Second arugment must be a string.' );
end
isobj = strcmp( mode, 'objective' );

pr = cvx___.problems( end );
touch( pr.self, x );
bz = x.basis_ ~= 0;
bs = sum( bz, 1 );
bc = bz( 1, : );
tt = bs > bc + 1;
at = any( tt );
if at,

    %
    % Replace posynomials with log-convex monomials --- that is, unless
    % we are taking the exponential of one. (Not that I know why we
    % would do that!) In that case, we should just use a standard
    % linear replacement and leave it at that.
    %

    if ~isequal( mode, 'exponential' ),
        if isempty( remap ),
            remap = cvx_remap( 'posynomial' );
        end
        t2 = remap( cvx_classify( x ) );
        if any( t2 ),
            if all( t2 ),
                x = exp( log( x ) );
            else
                x = cvx_subsasgn( x, t2, exp( log( cvx_subsref( x, t2 ) ) ) );
            end
            bc( t2 ) = 0;
            tt = tt & ~t2;
            at = any( tt );
        end
    end

    %
    % Replace other multivariable forms with single-variable forms
    %

    if at,
        abc = any( bc( :, tt ) );
        if abc,
            xc = cvx_constant( x );
            x = x - xc;
        end
        forms = cvx___.linforms;
        repls = cvx___.linrepls;
        [ x, forms, repls ] = replcols( x, tt, 'full', forms, repls, isobj );
        cvx___.linforms = forms;
        cvx___.linrepls = repls;
        if abc,
            x = x + xc;
        end
    end

end

%
% Arguments: no constraints on coefficients or constant values
% Objectives:     all replaced with coefficient > 0, constant terms preserved
% Exponentiation: coefficient == 1, constant terms perserved
% Logarithm:      coefficient >  0, constant terms eliminated
%

switch mode,
    case { 'argument', 'constraint' },
        tt = false;
    case 'objective',
        tt = ~tt | sum( x.basis_, 1 ) < x.basis_( 1, : );
        usexc = true;
    case 'logarithm',
        tt = sum( x.basis_, 1 ) < x.basis_( 1, : );
        usexc = false;
    case 'exponential';
        tt = sum( x.basis_, 1 ) ~= x.basis_( 1, : ) + 1;
        usexc = true;
    otherwise,
        error( [ 'Invalid normalization mode: ', mode ] );
end
if any( tt ),
    if usexc,
        abc = any( bc );
        if abc,
            xc = cvx_constant( x );
            x = x - xc;
        end
    else
        abc = false;
    end
    forms = cvx___.uniforms;
    repls = cvx___.unirepls;
    [ x, forms, repls ] = replcols( x, tt, 'none', forms, repls, isobj );
    cvx___.uniforms = forms;
    cvx___.unirepls = repls;
    if abc,
        x = x + xc;
    end
end

function [ x, forms, repls ] = replcols( x, tt, mode, forms, repls, isobj )

%
% Sift through the forms, removing duplicates
%

global cvx___
bN = vec( cvx_subsref( x, tt ) );
nO = length( forms );
nN = length( bN );
if nO ~= 0,
    bN = [ forms ; bN ];
end
[ bNR, bN ] = bcompress( bN, mode, nO );
bNR = bNR( :, nO + 1 : end );
nB = length( bN ) - nO;

%
% Create the replacement variables
%

if nB ~= 0,
    forms   = bN;
    bN      = cvx_subsref( bN, nO + 1 : nO + nB, 1 );
    newrepl = newvar( cvx___.problems( end ).self, '', nB );
    [ ndxs, temp ] = find( newrepl.basis_ ); %#ok
    repls = [ repls ; newrepl ];
    bV = cvx_vexity( bN );
    cvx___.vexity( ndxs ) = bV;
    cvx___.readonly( ndxs ) = vec( cvx_readlevel( bN ) );
    if ~isobj,
        ss = bV == 0;
        if any( ss ),
            temp = cvx_basis( bN );
            temp = any( temp( :, ss ), 2 );
            temp( ndxs( ss ) ) = true;
            cvx___.canslack( temp ) = false;
        end
    end
end

%
% Re-expand the structure
%

x = cvx_subsasgn( x, tt, buncompress( bNR, repls, nN ) );

% 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:
classdef OptimalExponentsComparatorWithPNorm < handle

    properties (Access = private)
        optimalComputer
    end

    properties (Access = private)
        optimalExponentParams
        pNorm
    end

    methods (Access = public)

        function obj = OptimalExponentsComparatorWithPNorm(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            obj.computeOptimalExponents();
            obj.plotStressNormRelationWithQ();
            obj.plotConvergence();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.optimalExponentParams.txi  = pi/4;
            obj.optimalExponentParams.rho  = 0.15;
            obj.optimalExponentParams.phi  = pi/4;
            obj.optimalExponentParams.hMesh = 0.01;
            obj.pNorm = {2,4,16,32,'max'};
        end

        function computeOptimalExponents(obj)
            s = obj.optimalExponentParams;
            for iP = 1:length(obj.pNorm)
                p = obj.pNorm{iP};
                s.pNorm = p;
                s.fileName = ['ExamplePaper',num2str(p)];
                oC = OneOptimalExponentComputerAndFunctionVariation(s);
                oC.compute();
                obj.optimalComputer{iP} = oC;
            end

        end

        function plotStressNormRelationWithQ(obj)
            f = figure();
            hold on
            for iP = 1:length(obj.pNorm)
                h{iP} = plot(obj.optimalComputer{iP}.qValues(2:end),obj.optimalComputer{iP}.fValues(2:end),'-+');
                p = obj.pNorm{iP};
                if isequal(p,'max')
                    leg{iP} = '$p = \infty$';
                else
                    leg{iP} = ['p = ',num2str(p)];
                end
            end
            outPutPath = '/home/alex/git-repos/MicroStructurePaper/';
            fTitle = 'Stress norm';
            xlabel('$q$','interpreter','latex')
            % set(gca, 'XTickLabel',[])
            ylabel(fTitle,'interpreter','latex')

            for iP = 1:length(obj.pNorm)
                x = obj.optimalComputer{iP}.qValues(2:end);
                y = obj.optimalComputer{iP}.fValues(2:end);
                [~,imin] = min(y);
                plot(x(imin),y(imin),'-s','Color',h{iP}.Color,'MarkerSize',10,'MarkerEdgeColor',h{iP}.Color,'MarkerFaceColor',h{iP}.Color)
            end

            legObj = legend(leg);
            set(legObj,'Interpreter','latex','Location','Best');
            %outputName = [outPutPath,'StressNormVsQForPmax'];
            %outputName = [outPutPath,'StressNormVsQForDifferentP'];
            printer = plotPrinter(f,h);
           % printer.print(outputName);
        end

        function plotConvergence(obj)
            f = figure();
            hold on
            for iP = 1:length(obj.pNorm)
                h{iP} = plot(obj.optimalComputer{iP}.fOptIter,'-+');
                p = obj.pNorm{iP};
                if isequal(p,'max')
                    leg{iP} = '$p = \infty$';
                else
                    leg{iP} = ['p = ',num2str(p)];
                end
            end
            outPutPath = '/home/alex/git-repos/MicroStructurePaper/';
            fTitle = 'Stress norm';
            xlabel('iterations','interpreter','latex')
            ylabel(fTitle,'interpreter','latex')

            legObj = legend(leg);
            set(legObj,'Interpreter','latex','Location','Best');
            outputName = [outPutPath,'StressNormMinimizationForPmax'];
            %outputName = [outPutPath,'StressNormMinimizationFineMax'];
            printer = plotPrinter(f,h);
           % printer.print(outputName);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function y = tmat_mxp ( a, x )

%*****************************************************************************80
%
%% TMAT_MXP multiplies a geometric transformation matrix times a point.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics, Principles and Practice,
%    Addison Wesley, Second Edition, 1990.
%
%  Parameters:
%
%    Input, real A(4,4), the geometric transformation matrix.
%
%    Input, real X(3), the point to be multiplied.  The fourth
%    component of X is implicitly assigned the value of 1.
%
%    Output, real Y(3), the result of A*X.  The product is
%    accumulated in a temporary vector, and then assigned to the result.
%    Therefore, it is legal for X and Y to share memory.
%
  y(1:3) = a(1:3,4) + a(1:3,1:3) * x(1:3);

  return
end


text:
function area = polygon_area_2 ( n, v )

%*****************************************************************************80
%
%% POLYGON_AREA_2 computes the area of a polygon.
%
%  Discussion:
%
%    The area is the sum of the areas of the triangles formed by
%    node N with consecutive pairs of nodes.
%
%    If the vertices are given in counterclockwise order, the area
%    will be positive.  If the vertices are given in clockwise order,
%    the area will be negative.
%
%    Thanks to Martin Pineault for noticing that an earlier version
%    of this routine would not correctly compute the area of a nonconvex
%    polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer, John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real V(2,N), the coordinates of the vertices.
%
%    Output, real AREA, the area of the polygon.
%
  area = 0.0;

  for i = 1 : n - 2

    area_triangle = triangle_area ( ...
      v(1,i),   v(2,i),   ...
      v(1,i+1), v(2,i+1), ...
      v(1,n),   v(2,n) );

    area = area + area_triangle;

  end

  return
end


text:
% THE DEMO PRESENTS AN ERROR PROPAGATION USING A MONTE-CARLO METHOD
% The error in the positioning of the joints is propagated to an error
% in the position of the end effector

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

close all

M=1500; %number of particles

%load arm parameters
robot=load_robot('ABB', 'IRB140');

%standard deviation AT EACH JOINT
sigmaq=0.017;%rad

%find errors around this pose
q=[pi/2 -pi/2 0 0 0 0]';

puntos=[];
for i=1:M,
    qi = q + [normrnd(0, sigmaq, robot.DOF, 1)];
    T=directkinematic(robot, qi);
    puntos=[puntos; T(1,4) T(2,4) T(3,4)];
end

adjust_view(robot)
drawrobot3d(robot,q), hold on
plot3(puntos(:,1),puntos(:,2), puntos(:,3),'r.')

repo: ARTE


text:
function [fmat ortirfcell]=tvcfsim1(beta,D,ybarT,data_exo_p,Fperiods,n,m,p,k)

% first recover the set of At matrices, Cbart matrices, and xbart vectors
% initiate At, Cbart and xt
At=cell(Fperiods,1);
Cbart=cell(Fperiods,1);
xbart=cell(Fperiods,1);
% reshape beta to obtain one equation per column
beta=reshape(beta,k,Fperiods*n);
% separate the coefficients related to exogenous variables (the last m rows of each equation)
% from the ones related to endogenous variables (the other entries) in order to be able to obtain At and Cbart
beta_endo=reshape(beta(1:end-m,:),[k-m n Fperiods]);
beta_exo=reshape(beta(end-m+1:end,:),[m n Fperiods]);
% then loop over forecast periods to recover the series of matrices
for ii=1:Fperiods
At{ii,1}=sparse([beta_endo(:,:,ii)';speye(n*(p-1)) sparse(n*(p-1),n)]);
Cbart{ii,1}=sparse([beta_exo(:,:,ii)';sparse(n*(p-1),m)]);
xbart{ii,1}=data_exo_p(ii,:)';
end

% then produce the series of forecasts and impulse response functions
% preliminary element: generate the selection matrix J
J=[speye(n) sparse(n,n*(p-1))];
% then create the cells storing the IRFs and the forecasts
ortirfcell=cell(Fperiods,Fperiods);
fmat=[];

% estimate first the forecasts
% loop over forecasts periods
for hh=1:Fperiods
% compute the first product term
prod1=speye(n*p);
   for ii=1:hh
   prod1=prod1*At{1+hh-ii,1};
   end
% multiply by ybarT to obtain the first term
term1=prod1*ybarT;
% obtain the second term
% initiate the summation
summ=sparse(n*p,1);
   % loop over summation periods
   for ii=1:hh
   % initiate the product
   prod2=speye(n*p);
      % loop over product periods
      for jj=ii:hh-1
      prod2=prod2*At{hh+ii-jj,1};
      end
   % multiply by Cbart*xt and add to summation
   summ=summ+prod2*Cbart{ii,1}*xbart{ii,1};
   end
term2=summ;
% obtain the forecast for the period
fmat(:,hh)=J*(term1+term2);
end
fmat=fmat';

% then compute the period-specific IRFs
% loop over forecast periods
for hh=1:Fperiods
% the first row of the cell represents the term ii=hh in the summation in 273: it is always equal to identity ( and then becomes D once multiplied by the structural matrix)
ortirfcell{1,hh}=D;
   % loop over IRF periods (the summation term in 272)
   for ii=1:hh-1
   % initiate the product
   prod3=speye(n*p);
      % loop over the periods involved into the product and calculate it
      for jj=ii:hh-1
      prod3=prod3*At{hh+ii-jj,1};
      end
   % recover the matrix of interest from the selection matrix J
   irfmat=full(J*prod3*J');
   % obtain orthogonalised IRFs
   ortirfmat=irfmat*D;
   % record in IRFcell
   ortirfcell{hh-ii+1,hh}=ortirfmat;
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
%%*******************************************************************
%%  Read in a problem in SDPA sparse format.
%%
%%  [blk,At,C,b] = read_sdpa(fname)
%%
%%  Input: fname = name of the file containing SDP data in
%%                 SDPA foramt.
%%  Important: the data is assumed to contain only
%%             semidefinite and linear blocks.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%******************************************************************

   function [blk,At,C,b] = read_sdpa(fname);

%%
%%  Open the file for input
%%
   compressed = 0;
   if exist(fname)
      fid = fopen(fname,'r');
   elseif exist([fname,'.Z']);
      compressed = 1;
      unix(['uncompress ',fname,'.Z']);
      fid = fopen(fname,'r');
   elseif exist([fname,'.gz']);
      compressed = 2;
      unix(['gunzip ',fname,'.gz']);
      fid = fopen(fname,'r');
   else
      fprintf('*** Problem not found, please specify the correct path or problem name. \n');
      blk = []; At = []; C = []; b = [];
      return;
   end
%%
%%  Clean up special characters and comments from the file
%%
   [datavec,count] = fscanf(fid,'%c');
   linefeeds = findstr(datavec,char(10));
   comment_chars = '*"=';
   cumidx = [];
   for i=1:length(comment_chars)
      idx = findstr(datavec,comment_chars(i));
      cumidx = [cumidx,idx];
   end
   for j=length(cumidx):-1:1
      if (cumidx(j)==1) | (strcmp(datavec(cumidx(j)-1),char(10)))
         datavec(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds))))='';
      else
         datavec(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds)))-1)='';
      end
   end
   special_chars = ',{}()';
   cumidx=[];
   for i=1:length(special_chars)
      idx = findstr(datavec,special_chars(i));
      cumidx = [cumidx,idx];
   end
   datavec(cumidx) = blanks(length(cumidx));
   clear linefeeds;
%%
%% Close the file
%%
   fclose('all');
   if compressed==1; unix(['compress ',fname]); end;
   if compressed==2; unix(['gzip ',fname]); end;
%%
%%  Next, read in basic problem size parameters.
%%
   datavec = sscanf(datavec,'%f');
   if size(datavec,1) < size(datavec,2); datavec = datavec'; end;
   m = datavec(1);
   numblk  = datavec(2);
   blksize = datavec(2+[1:numblk]);
   if size(blksize,1) > size(blksize,2); blksize = blksize'; end
%%
%% Get input  b.
%%
   idxstrt = 2+numblk;
   b = datavec(idxstrt+[1:m]);
   idxstrt = idxstrt+m;
   b = -b;
%%
%% Construct blk
%%
   deblksize = 100;
   spblkidxtmp = find( (blksize>1) & (blksize < deblksize) );
   spblkidxtmp = sort(spblkidxtmp);
   deblkidx = find( (blksize<=1) | (blksize >= deblksize) );
   denumblk = length(deblkidx);
   linblkidx = zeros(1,denumblk);
   for p = 1:denumblk
      n = blksize(deblkidx(p));
      if (n > 1);
         blk{p,1} = 's'; blk{p,2} = n;
         n2 = n*(n+1)/2;
         At{p,1} = sparse(n2,m);
         C{p,1} = sparse(n,n);
      else
         linblkidx(p) = p;
         blk{p,1} = 'l'; blk{p,2} = abs(n);
         At{p,1} = sparse(abs(n),m);
         C{p,1} = sparse(abs(n),1);
      end
   end
   if ~isempty(spblkidxtmp)
      maxnumblk = 200;
      spnumblk = ceil(length(spblkidxtmp)/maxnumblk);
      for q = 1:spnumblk
         if (q < spnumblk)
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: q*maxnumblk]);
         else
            spblkidxall{q} = spblkidxtmp([(q-1)*maxnumblk+1: length(spblkidxtmp)]);
         end
         tmp = blksize(spblkidxall{q});
         blk{denumblk+q,1} = 's';
         blk{denumblk+q,2} = tmp;
         n2 = sum(tmp.*(tmp+1))/2;
         At{denumblk+q,1} = sparse(n2,m);
         C{denumblk+q,1} = sparse(sum(tmp),sum(tmp));
      end
   else
      spnumblk = 0;
   end
   linblkidx(denumblk+[1:spnumblk]) = zeros(1,spnumblk);
%%
%% Construct single blocks of A,C
%%
   len = length(datavec);
   Y = reshape(datavec(idxstrt+1:len),5,(len-idxstrt)/5)';
   clear datavec;
   Y = sortrows(Y,[1 2]);
   matidx = [0; find(diff(Y(:,1)) ~= 0); size(Y,1)];
%%
   for k = 1:length(matidx)-1
      idx = [matidx(k)+1 : matidx(k+1)];
      Ytmp  = Y(idx,1:5);
      matno = Ytmp(1,1);
      Ytmp2 = Ytmp(:,2);
      for p = 1:denumblk
         n  = blksize(deblkidx(p));
         idx = find(Ytmp2 == deblkidx(p));
         ii = Ytmp(idx,3); jj = Ytmp(idx,4); vv =Ytmp(idx,5);
         len = length(idx);
         if (n > 1)
            idxtmp = find(ii > jj);
            if ~isempty(idxtmp);
               tmp = jj(idxtmp);
               jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
            end
            tmp = -sparse(ii,jj,vv,n,n);
            tmp = tmp + triu(tmp,1)';
         else
            tmp = -sparse(ii,ones(len,1),vv,abs(n),1);
         end
         if (matno == 0)
            C{p,1} = tmp;
         else
            if (n > 1)
               At{p,1}(:,matno) = svec(blk(p,:),tmp,1);
            else
               At{p,1}(:,matno) = tmp;
            end
         end
      end
   end
%%
%% Construct big sparse block of A,C
%%
if (spnumblk > 0)
   Y1 = Y(:,1);
   diffY1 = find(diff([-1; Y1; inf]));
   for kk = 1:length(diffY1)-1
      idx = [diffY1(kk) : diffY1(kk+1)-1];
      matno = Y1(diffY1(kk));
      Ytmp  = Y(idx,1:5);
      Ytmp2 = Ytmp(:,2);
      maxYtmp2  = Ytmp2(length(Ytmp2));
      minYtmp2  = Ytmp2(1);
      diffYtmp2 = Ytmp2(find(diff([-1; Ytmp2])));
      for q = 1:spnumblk
         spblkidx    = spblkidxall{q};
         maxspblkidx = spblkidx(length(spblkidx));
         minspblkidx = spblkidx(1);
         count = 0;
         if (minYtmp2 <= maxspblkidx) & (maxYtmp2 >= minspblkidx)
            tmpblksize = blksize(spblkidx);
            n = sum(tmpblksize);
            cumspblksize = [0 cumsum(tmpblksize)];
            n2 = sum(tmpblksize.*(tmpblksize+1))/2;
            idx = zeros(n2,1); offset = zeros(n2,1);
            for t = [1:length(diffYtmp2)]
  	       p = find(spblkidx == diffYtmp2(t));
               if ~isempty(p)
	          idxtmp = find(Ytmp2 == spblkidx(p));
                  len = length(idxtmp);
     	          idx(count+[1:len]) = idxtmp;
                  offset(count+[1:len]) = cumspblksize(p);
                  count = count + len;
               end
            end
            idx = idx(1:count); offset = offset(1:count);
            ii = Ytmp(idx,3)+offset; jj = Ytmp(idx,4)+offset; vv = Ytmp(idx,5);
            idxtmp = find(ii > jj);
            if ~isempty(idxtmp);
               tmp = jj(idxtmp);
               jj(idxtmp) = ii(idxtmp); ii(idxtmp) = tmp;
            end
	    idxeq = find(ii==jj);
            tmp = spconvert([ii jj -vv; jj ii -vv; n n 0]) ...
                  + spconvert([ii(idxeq) jj(idxeq) vv(idxeq); n n 0]);
            if (matno == 0)
               C{denumblk+q,1} = tmp;
            else
               At{denumblk+q,1}(:,matno) = svec(blk(denumblk+q,:),tmp,1);
            end
         end
      end
   end
end
%%
%% put all linear blocks together as a single linear block
%%
   idx = find(linblkidx);
   if (length(idx) > 1)
      sdpidx = find(linblkidx==0);
      blktmp = 0; Atmp = []; Ctmp = [];
      for k = 1:length(idx)
          tmp = linblkidx(idx(k));
          blktmp = blktmp+blk{tmp,2};
          Atmp = [Atmp; At{tmp}];
          Ctmp = [Ctmp; C{tmp}];
      end
      At = At(sdpidx); C = C(sdpidx); blk = blk(sdpidx,:);
      len = length(sdpidx);
      blk(2:len+1,:) = blk;
      blk{1,1} = 'l'; blk{1,2} = blktmp;
      At(2:len+1,1) = At; C(2:len+1,1) = C;
      At{1,1} = Atmp; C{1,1} = Ctmp;
   end
%%******************************************************************

repo: REST


text:
function [x_lambda,rho,eta,data,X] = maxent(A,b,lambda,w,x0)
%MAXENT Maximum entropy regularization.
%
% [x_lambda,rho,eta] = maxent(A,b,lambda,w,x0)
%
% Maximum entropy regularization:
%    min { || A x - b ||^2 + lambda^2*x'*log(diag(w)*x) } ,
% where -x'*log(diag(w)*x) is the entropy of the solution x.
% If no weights w are specified, unit weights are used.
%
% If lambda is a vector, then x_lambda is a matrix such that
%    x_lambda = [x_lambda(1), x_lambda(2), ... ] .
%
% This routine uses a nonlinear conjugate gradient algorithm with "soft"
% line search and a step-length control that insures a positive solution.
% If the starting vector x0 is not specified, then the default is
%    x0 = norm(b)/norm(A,1)*ones(n,1) .

% Per Christian Hansen, IMM and Tommy Elfving, Dept. of Mathematics,
% Linkoping University, 06/10/92.

% Reference: R. Fletcher, "Practical Methods for Optimization",
% Second Edition, Wiley, Chichester, 1987.

% Set defaults.
flat = 1e-3;     % Measures a flat minimum.
flatrange = 10;  % How many iterations before a minimum is considered flat.
maxit = 150;     % Maximum number of CG iterations.
minstep = 1e-12; % Determines the accuracy of x_lambda.
sigma = 0.5;     % Threshold used in descent test.
tau0 = 1e-3;     % Initial threshold used in secant root finder.

% Initialization.
[m,n] = size(A); x_lambda = zeros(n,length(lambda)); F = zeros(maxit,1);
if (min(lambda) <= 0)
  error('Regularization parameter lambda must be positive')
end
if (nargin ==3), w  = ones(n,1); end
if (nargin < 5), x0 = ones(n,1); end

% Treat each lambda separately.
for j=1:length(lambda);

  % Prepare for nonlinear CG iteration.
  l2 = lambda(j)^2;
  x  = x0; Ax = A*x;
  g  = 2*A'*(Ax - b) + l2*(1 + log(w.*x));
  p  = -g;
  r  = Ax - b;

  % Start the nonlinear CG iteration here.
  delta_x = x; dF = 1; it = 0; phi0 = p'*g;
  while (norm(delta_x) > minstep*norm(x) & dF > flat & it < maxit & phi0 < 0)
    it = it + 1;

    % Compute some CG quantities.
    Ap = A*p; gamma = Ap'*Ap; v = A'*Ap;

    % Determine the steplength alpha by "soft" line search in which
    % the minimum of phi(alpha) = p'*g(x + alpha*p) is determined to
    % a certain "soft" tolerance.
    % First compute initial parameters for the root finder.
    alpha_left = 0; phi_left = phi0;
    if (min(p) >= 0)
      alpha_right = -phi0/(2*gamma);
      h = 1 + alpha_right*p./x;
    else
      % Step-length control to insure a positive x + alpha*p.
      I = find(p < 0);
      alpha_right = min(-x(I)./p(I));
      h = 1 + alpha_right*p./x; delta = eps;
      while (min(h) <= 0)
        alpha_right = alpha_right*(1 - delta);
        h = 1 + alpha_right*p./x;
        delta = delta*2;
      end
    end
    z = log(h);
    phi_right = phi0 + 2*alpha_right*gamma + l2*p'*z;
    alpha = alpha_right; phi = phi_right;

    if (phi_right <= 0)

      % Special treatment of the case when phi(alpha_right) = 0.
      z = log(1 + alpha*p./x);
      g_new = g + l2*z + 2*alpha*v; t = g_new'*g_new;
      beta = (t - g'*g_new)/(phi - phi0);

    else

      % The regular case: improve the steplength alpha iteratively
      % until the new step is a descent step.
      t = 1; u = 1; tau = tau0;
      while (u > -sigma*t)

        % Use the secant method to improve the root of phi(alpha) = 0
        % to within an accuracy determined by tau.
        while (abs(phi/phi0) > tau)
          alpha = (alpha_left*phi_right - alpha_right*phi_left)/...
                  (phi_right - phi_left);
          z = log(1 + alpha*p./x);
          phi = phi0 + 2*alpha*gamma + l2*p'*z;
          if (phi > 0)
            alpha_right = alpha; phi_right = phi;
          else
            alpha_left  = alpha; phi_left  = phi;
          end
        end

        % To check the descent step, compute u = p'*g_new and
        % t = norm(g_new)^2, where g_new is the gradient at x + alpha*p.
        g_new = g + l2*z + 2*alpha*v; t = g_new'*g_new;
        beta = (t - g'*g_new)/(phi - phi0);
        u = -t + beta*phi;
        tau = tau/10;

      end  % End of improvement iteration.

    end  % End of regular case.

    % Update the iteration vectors.
    g = g_new; delta_x = alpha*p;
    x = x + delta_x;
    p = -g + beta*p;
    r = r + alpha*Ap;
    phi0 = p'*g;

    % Compute some norms and check for flat minimum.
    rho(j,1) = norm(r); eta(j,1) = x'*log(w.*x);
    F(it) = rho(j,1)^2 + l2*eta(j,1);
    if (it <= flatrange)
      dF = 1;
    else
      dF = abs(F(it) - F(it-flatrange))/abs(F(it));
    end

    data(it,:) = [F(it),norm(delta_x),norm(g)];
    X(:,it) = x;

  end  % End of iteration for x_lambda(j).

  x_lambda(:,j) = x;

end


text:
function f = tt_loglikelihood(X,M)
%TT_LOGLIKELIHOOD Compute log-likelihood of data X with model M.
%
%   F = TT_LOGLIKELIHOOD(X,M) computes the log-likelihood of model M given
%   data X, where M is a ktensor and X is a tensor or sptensor.
%   Specifically, F = - (sum_i m_i - x_i * log_i) where i is a multiindex
%   across all tensor dimensions.
%
%   See also cp_apr, tensor, sptensor, ktensor.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

N = ndims(X);

if ~isa(M, 'ktensor')
    error('M must be a ktensor');
end

M = normalize(M,1,1);

if isa(X, 'sptensor')
    xsubs = X.subs;
    A = M.U{1}(xsubs(:,1),:);
    for n = 2:N
       A = A .* M.U{n}(xsubs(:,n),:);
    end
    f = sum(X.vals .* log(sum(A,2))) - sum(sum(M.U{1}));
else
    f = sum(sum(double(tenmat(X,1)) .* log(double(tenmat(M,1))))) - sum(sum(M.U{1}));
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt