text:
%* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%*                                                                           *
%*                  This file is part of the program and library             *
%*         SCIP --- Solving Constraint Integer Programs                      *
%*                                                                           *
%*    Copyright (C) 2002-2012 Konrad-Zuse-Zentrum                            *
%*                            fuer Informationstechnik Berlin                *
%*                                                                           *
%*  SCIP is distributed under the terms of the ZIB Academic Licence.         *
%*                                                                           *
%*  You should have received a copy of the ZIB Academic License              *
%*  along with SCIP; see the file COPYING. If not email to scip@zib.de.      *
%*                                                                           *
%* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

%@file    scip.m
%@brief   SCIP installation file for MATLAB interface
%@author  Ambros Gleixner
%@author  Timo Berthold

function[bestsol, objval] = scip(matrix, lhs, rhs, obj, lb, ub, vartype, objsense)

% copyright %
fprintf('SCIP-MATLAB interface\n');
fprintf('Copyright (c) 2002-2012 Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)\n\n');

% syntax description %
usage = '\n';
usage = [usage 'syntax: [bestsol, objval] = scip(matrix, lhs, rhs, obj, lb, ub, vartype, objsense)\n'];
usage = [usage 'input parameters are\n'];
usage = [usage '   matrix   : real matrix with constraint nonzeros\n'];
usage = [usage '   lhs      : real vector with left-hand sides of constraints\n'];
usage = [usage '   rhs      : real vector with right-hand sides of constraints\n'];
usage = [usage '   obj      : real vector with objective coefficient of variables\n'];
usage = [usage '   lb       : real vector with lower bounds of variables\n'];
usage = [usage '   ub       : real vector with upper bounds of variables\n'];
usage = [usage '   vartype  : character vector with types of variables (b: binary, i: general integer, c: continuous)\n'];
usage = [usage '   objsense : value indicating optimization sense (+1: minimize, -1: maximize)\n'];
usage = [usage 'output parameters are\n'];
usage = [usage '   bestsol  : real vector containing best solution found\n'];
usage = [usage '   objval   : objective function value of best solution found\n'];
usage = [usage '\n'];

% output is undefined for an invalid call %
bestsol = NaN;
objval = NaN;

% we do not define default values for input parameters %
if( nargin ~= 8 )
   fprintf('\ninvalid call: provide all input parameters\n');
   fprintf(usage);
   return;
end

% display input parameters %
matrix
lhs
rhs
obj
lb
ub
vartype
objsense

% matrix must be two-dimensional %
if( ndims(matrix) ~= 2 )
   fprintf('\ninvalid call: ndims(matrix) ~= 2\n');
   fprintf(usage);
   return;
end

% get number of variables and constraints %
[nconss,nvars] = size(matrix);

% check that vectors are one-dimensional %
if( ~isvector(lhs) )
   fprintf('\ninvalid call: ~isvector(lhs)\n');
   fprintf(usage);
   return;
elseif( ~isvector(rhs) )
   fprintf('\ninvalid call: ~isvector(rhs)\n');
   fprintf(usage);
   return;
elseif( ~isvector(obj) )
   fprintf('\ninvalid call: ~isvector(obj)\n');
   fprintf(usage);
   return;
elseif( ~isvector(lb) )
   fprintf('\ninvalid call: ~isvector(lb)\n');
   fprintf(usage);
   return;
elseif( ~isvector(ub) )
   fprintf('\ninvalid call: ~isvector(ub)\n');
   fprintf(usage);
   return;
elseif( ~isvector(vartype) )
   fprintf('\ninvalid call: ~isvector(vartype)\n');
   fprintf(usage);
   return;
end

% check that vectors have consistent length %
if( length(lhs) ~= nconss )
   fprintf('\ninvalid call: length(lhs) ~= size(matrix,1)\n');
   fprintf(usage);
   return;
elseif( length(rhs) ~= nconss )
   fprintf('\ninvalid call: length(rhs) ~= size(matrix,1)\n');
   fprintf(usage);
   return;
elseif( length(obj) ~= nvars )
   fprintf('\ninvalid call: length(obj) ~= size(matrix,2)\n');
   fprintf(usage);
   return;
elseif( length(lb) ~= nvars )
   fprintf('\ninvalid call: length(lb) ~= size(matrix,2)\n');
   fprintf(usage);
   return;
elseif( length(ub) ~= nvars )
   fprintf('\ninvalid call: length(ub) ~= size(matrix,2)\n');
   fprintf(usage);
   return;
elseif( length(vartype) ~= nvars )
   fprintf('\ninvalid call: length(vartype) ~= size(matrix,2)\n');
   fprintf(usage);
   return;
end

% check that input values have correct type %
if( ~isreal(matrix) )
   fprintf('\ninvalid call: ~isreal(matrix)\n');
   fprintf(usage);
   return;
elseif( ~isreal(lhs) )
   fprintf('\ninvalid call: ~isreal(lhs)\n');
   fprintf(usage);
   return;
elseif( ~isreal(rhs) )
   fprintf('\ninvalid call: ~isreal(rhs)\n');
   fprintf(usage);
   return;
elseif( ~isreal(obj) )
   fprintf('\ninvalid call: ~isreal(obj)\n');
   fprintf(usage);
   return;
elseif( ~isreal(lb) )
   fprintf('\ninvalid call: ~isreal(lb)\n');
   fprintf(usage);
   return;
elseif( ~isreal(ub) )
   fprintf('\ninvalid call: ~isreal(ub)\n');
   fprintf(usage);
   return;
elseif( ~ischar(vartype) )
   fprintf('\ninvalid call: ~ischar(vartype)\n');
   fprintf(usage);
   return;
elseif( ~strcmp(objsense,'min') && ~strcmp(objsense,'max') )
   fprintf('\ninvalid call: objsense must be min or max\n');
   fprintf(usage);
   return;
else
   for v = 1:length(vartype)
      if( vartype(v) ~= 'b' && vartype(v) ~= 'i' && vartype(v) ~= 'c' )
         fprintf('\ninvalid call: vartype(%d) must be one of {b,i,c}\n', v);
         fprintf(usage);
         return;
      end
   end
end

% call SCIP %
[bestsol, objval] = matscip(matrix, lhs, rhs, obj, lb, ub, vartype, objsense);

bestsol
objval

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% motorseed comparison
% Perform regression of motor traces (seed or raw) with all cells.
% Histograms of regression coefficients (red: motorseed; gray: raw traces)
% show that for each of all fish, regression with motor seeds yields higher
% coefficients for the population.
% not surprising, but shows that (even if a bit circular), the motorseeds
% have a better chance of capturing motor related neurons.

%%
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%%
binEdges = -1:0.05:1;
nbins = length(binEdges)-1;
N_hist_notseed = zeros(18,nbins,2);
N_hist_seed = zeros(18,nbins,2);
for i_fish = 1:18
    ClusterIDs = [2,1];
    [cIX,gIX,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);

    %%
    isMotorseed = 0;
    setappdata(hfig,'isMotorseed',isMotorseed);
    [~,~,behavior] = UpdateTimeIndex(hfig);

    [~,~,regressor_m,names_m] = GetMotorRegressor(behavior,i_fish);
    if isMotorseed
        Reg = regressor_m;
    else
        Reg = regressor_m([1,3],:);
    end
    % regression
    Corr = corr(Reg',M_0');
    for i_reg = 1:2
        R = Corr(i_reg,:);
        [N,~] = histcounts(R,binEdges);
        N_hist_notseed(i_fish,:,i_reg) = N;
    end
    %%
    isMotorseed = 1;
    setappdata(hfig,'isMotorseed',isMotorseed);
    [~,~,behavior] = UpdateTimeIndex(hfig);

    [~,~,regressor_m,names_m] = GetMotorRegressor(behavior,i_fish);
    if isMotorseed
        Reg = regressor_m;
    else
        Reg = regressor_m([1,3],:);
    end
    % regression
    Corr = corr(Reg',M_0');
    for i_reg = 1:2
        R = Corr(i_reg,:);
        [N,~] = histcounts(R,binEdges);
        N_hist_seed(i_fish,:,i_reg) = N;
    end
end

%%
figure;
cmap = [1 0.4 0.4;0.4 0.4 0.4];
range_fish = 1:18;
m = length(range_fish);
N_type = cell(1,2);
N_type{1} = N_hist_seed;
N_type{2} = N_hist_notseed;
% i_plot = 0;

% left half
i_reg = 1;
for i_fish = range_fish
    for i_type = 1:2
        i_plot = 1+2*(i_fish-1);%i_plot+1;
        subplot(m,2,i_plot);
        hold on;

        counts = squeeze(N_type{i_type}(i_fish,:,i_reg));
        histogram('BinEdges',binEdges,'BinCounts',counts,'FaceColor',cmap(i_type,:))

        %     ymax = max(max(N),max(N_shf));
        %     plot([0,0],[0,ymax],'r--');
        ymax = counts(24)+eps;
        xlim([-1,1]);
        ylim([0,ymax]);
        set(gca,'YTick',[])
        ylabel('a.u.')
        title(['Fish ' num2str(i_fish)])
    end
    end

% right half
i_reg = 2;
for i_fish = range_fish
    for i_type = 1:2
        i_plot = 2*i_fish;%i_plot+1;
        subplot(m,2,i_plot);
        hold on;

        counts = squeeze(N_type{i_type}(i_fish,:,i_reg));
        histogram('BinEdges',binEdges,'BinCounts',counts,'FaceColor',cmap(i_type,:))

        %     ymax = max(max(N),max(N_shf));
        %     plot([0,0],[0,ymax],'r--');
        ymax = counts(24)+eps;
        xlim([-1,1]);
        ylim([0,ymax]);
        set(gca,'YTick',[])
        ylabel('a.u.')
        title(['Fish ' num2str(i_fish)])
    end
end

author: xiuyechen

repo: FishExplorer


text:
function [ r ] = EffRank(s, thresh)
%[ r ] = EffRank(s, thresh)
% get the effective rank of an singular value array
% s: the singular values in descending order
% thresh: the total variance to preserve
% author: Liang Xiong (lxiong@cs.cmu.edu)

if nargin < 2; thresh = 0.99; end
assert(all(s > 0), 'not a proper singular value array');

s = sort(s.^2,'descend');
cs = cumsum(s);
if cs(end) < 1e-10
    error('EffRank has encounted a zero matrix');
else
    r = sum(cs./cs(end) <= thresh) + 1;
    r = min(length(s), r);
end

author: andrewssobral

repo: lrslibrary


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

A1_intercept = [ones(1,size(X',2)); X'];

Z2 = Theta1*A1_intercept;
A2 = sigmoid(Z2);
A2_intercept = [ones(1,size(A2,2)); A2];

Z3 = Theta2*A2_intercept;
A3 = sigmoid(Z3);

[M, I] = max(A3);
p = I(:);

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

end

author: zlotus


text:
function F = fliplr(F)
%FLIPLR   Flip the columns of a CHEBMATRIX.

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

F.blocks = fliplr(F.blocks);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function a = i4row_swap ( m, n, a, irow1, irow2 )

%*****************************************************************************80
%
%% I4ROW_SWAP swaps two rows of an I4ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an array of data.
%
%    Input, integer IROW1, IROW2, the two rows to swap.
%
%    Output, integer A(M,N), the array of data, with two rows swapped.
%

%
%  Check.
%
  if ( irow1 < 1 || m < irow1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4ROW_SWAP - Fatal error!\n' );
    fprintf ( 1, '  IROW1 is out of range.\n' );
    error ( 'I4ROW_SWAP - Fatal error!' );
  end

  if ( irow2 < 1 || m < irow2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4ROW_SWAP - Fatal error!\n' );
    fprintf ( 1, '  IROW2 is out of range.\n' );
    error ( 'I4ROW_SWAP - Fatal error!' );
  end

  if ( irow1 == irow2 )
    return
  end

  row(1:n) = a(irow1,1:n);
  a(irow1,1:n) = a(irow2,1:n);
  a(irow2,1:n) = row(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function time_skew = offset_to_timeskew(offset, N, fs)
    time_skew = mod(offset, N) / fs;
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function plugsStruct = readLeksellPlugFile(filename)
%"readLeksellPlugFile"
%   Reads a Leksell Plug File into meaningful data structures.  Some values
%   are still undetermined.  However, none of these values are currently
%   used in CERR.
%
%KRK 05/28/07
%
%Usage:
%   function plugStruct = readLeksellPlugFile(filename)
%
% 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/>.

fid = fopen(filename, 'r', 'b');

data = decodeLeksellData(fid);

fclose(fid);

%If data only has one cell, it's empty since there is always one cell that
%   contains no meaningful data
if length(data) < 2
    plugsStruct = [];
    return;
end

% Get meaningful data
plugsData = data{1};

% Type of plugging
plugsStruct.plugType = plugsData{1};
% Number of plugged holes in the helmet
plugsStruct.numSourcesPlugged = plugsData{4};

% Unsolved variables
plugsStruct.mystery_value_2 = plugsData{2}; % integer value
% mystery_value_3 appears to be a binary sequence showing which holes are
%   plugged, but this may be incorrect
plugsStruct.mystery_value_3 = plugsData{3};

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [idx, iter, obj, H] = symnmf_cluster(X, k, options)
%SYMNMF_CLUSTER
% [idx, iter, obj, H] = symnmf_cluster(X, k, options)
%
% This function performs graph clustering on a data matrix.
%
% Input:
% X - NxP data matrix (N observations with dimension P),
%     where each row is one observation,
%     and each column is one feature/variable.
%     (If an NxN similarity matrix is available,
%      please call 'symnmf_anls' directly.)
% k - Specifying the number of clusters.
% options - A structure of optional parameters for
%           clustering (see details below).
%
% Default configurations:
%   options.graph_type = 'sparse';
%   options.similarity_type = 'gaussian';
%   options.graph_objfun = 'ncut';
%   (no options.kk given)
%   options.nn = 7;
%   options.tol = 1e-3;
%   options.maxiter = 10000;
%   options.rep = 1;
%   (no options.Hinit given)
%   options.computeobj = true;
%   options.alg = 'anls';
%
% Output:
% idx - Nx1 vector, containing the clustering assignment
%       of each observation.
% iter - Number of iterations actually used
% obj - Final value of the objective function
%           f(H) = ||A - HH'||_F^2
% H - Final result of low-rank matrix H
%
% Available options:
% 1. options.graph_type
%    = 'full': to construct a fully-connected graph
%    = 'sparse': to construct a sparse graph where
%                each observation is connected to
%                its KK nearest neighbors. (DEFAULT)
%
%    (Note: Only undirected graph is supported,
%           so the resulting similarity matrix
%           is symmetric.)
%
% 2. options.similarity_type
%    = 'gaussian': to use self-tuning Gaussian similarity,
%                  proposed in [Zelnik-Manor and Perona,
%                  2004] (see more details in functions
%                  'scale_dist3' and 'scale_dist3_knn') (DEFAULT)
%    = 'inner_product': to use inner product similarity
%                       (good choice for text data)
%
% 3. options.graph_objfun
%    = 'ncut': to use normalized cut (DEFAULT)
%    = 'rcut': to use ratio cut
%
% 4. options.kk (DEFAULT is unset, i.e. depending on N)
%    specifies the number of nearest neighbors WHEN
%    options.graph_type = 'sparse'.
%
%    The default value depends on N:
%    option.kk = floor(log2(N)) + 1;
%
% 5. options.nn (DEFAULT is 7)
%    specifies the nn-th neighbor, which is used in
%    the self-tuning Gaussian similarity WHEN
%    options.similarity_type = 'gaussian'.
%
% 6. options.tol (DEFAULT is 1e-3)
%    controls the termination of SymNMF algorithm.
%
% 7. options.maxiter (DEFAULT is 10000)
%    limits the number of iteration allowed in each run.
%
% 8. options.rep (DEFAULT is 1)
%    sets the nubmer of runs of SymNMF algorithm.
%
%    The return values 'idx', 'iter', 'obj', and 'H' are
%    corresponding to a single run, which is the run that
%    yields the lowest objective function value.
%
%    Using multiple runs of SymNMF may help avoid local minima.
%
% 9. options.Hinit (DEFAULT is unset, i.e. random initialization)
%    sets the initialization of matrix H in each run.
%
%    options.Hinit is a 3-dim array. The 3rd dimension may imply
%    the choice of options.rep, and each 2-dim array in the 1st
%    and 2nd dimensions is a NxK nonnegative matrix.
%
%    Random initialization is recommended. If options.Hinit has to
%    be a fixed set of matrices for testing, the instructions on
%    how to initialize H should be followed. Please consult the
%    usage of 'symnmf_newton'.
%
% 10. options.computeobj (DEFAULT is true)
%     specifies whether to compute the objective value
%     f(H) at the final solution H.
%
% 11. options.alg
%     = 'anls' to use ANLS algorithm (DEFAULT)
%     = 'newton' to use Newton-like algorithm
%
% This function is used for experiments in the following paper:
%     Da Kuang, Chris Ding, Haesun Park,
%     Symmetric Nonnegative Matrix Factorization for Graph Clustering,
%     The 12th SIAM International Conference on Data Mining (SDM '12), pp. 106--117.
% Please cite this paper if you find this code useful.
%

[n, p] = size(X);

if ~exist('options', 'var')
    graph_type = 'sparse';
    similarity_type = 'gaussian';
    graph_objfun = 'ncut';
    kk = floor(log2(n)) + 1;
    nn = 7;
    tol = 1e-3;
    maxiter = 10000;
    rep = 1;
    Hinit = [];
    computeobj = true;
    alg = 'anls';
else
    if isfield(options, 'graph_type')
        graph_type_names = {'full', 'sparse'};
        j = find(strcmpi(options.graph_type, graph_type_names));
        if ~isempty(j)
            graph_type = graph_type_names{j};
        else
            error('Invalid options.graph_type value!');
        end
    else
        graph_type = 'sparse';
    end
    if isfield(options, 'similarity_type')
        similarity_type_names = {'gaussian', 'inner_product'};
        j = find(strcmpi(options.similarity_type, similarity_type_names));
        if ~isempty(j)
            similarity_type = similarity_type_names{j};
        else
            error('Invalid options.similarity_type value!');
        end
    else
        similarity_type = 'gaussian';
    end
    if isfield(options, 'graph_objfun')
        graph_objfun_names = {'ncut', 'rcut'};
        j = find(strcmpi(options.graph_objfun, graph_objfun_names));
        if ~isempty(j)
            graph_objfun = graph_objfun_names{j};
        else
            error('Invalid options.graph_objfun value!');
        end
    else
        graph_objfun = 'ncut';
    end
    if isfield(options, 'kk')
        if ~isempty(options.kk) & isnumeric(options.kk) & options.kk < n
            kk = options.kk;
        else
            error('options.kk must be an integer less than N!');
        end
    else
        kk = floor(log2(n)) + 1;
    end
    if isfield(options, 'nn')
        if ~isempty(options.nn) & isnumeric(options.nn) & options.nn < n
            nn = options.nn;
        else
            error('options.nn must be an integer less than N!');
        end
    else
        nn = 7;
    end
    if isfield(options, 'tol')
        if ~isempty(options.tol) & isnumeric(options.tol) & options.tol > 0 & options.tol < 1
            tol = options.tol;
        else
            error('options.tol must be a real number and 0 < options.tol < 1!');
        end
    else
        tol = 1e-3;
    end
    if isfield(options, 'maxiter')
        if ~isempty(options.maxiter) & isnumeric(options.maxiter) & options.maxiter > 0
            maxiter = options.maxiter;
        else
            error('options.maxiter must be a positive integer!');
        end
    else
        maxiter = 10000;
    end
    if isfield(options, 'rep')
        if ~isempty(options.rep) & isnumeric(options.rep) & options.rep > 0
            rep = options.rep;
        else
            error('options.rep must be a positive integer!');
        end
        if isfield(options, 'Hinit') & size(options.Hinit, 3) ~= rep
            error('The third dimension of options.Hinit must match options.rep!')
        end
    else
        if isfield(options, 'Hinit') & ~isempty(options.Hinit)
            rep = size(options.Hinit, 3);
        else
            rep = 1;
        end
    end
    if isfield(options, 'Hinit')
        if ~isempty(options.Hinit) & isnumeric(options.Hinit) & size(options.Hinit, 1) == n & size(options.Hinit, 2) == k
            Hinit = options.Hinit;
        else
            error('The size of each initialization of H must be Nxk!');
        end
    else
        Hinit = [];
    end
    if isfield(options, 'computeobj')
        if ~isempty(options.computeobj) & (isnumeric(options.computeobj) | isboolean(options.computeobj))
            computeobj = options.computeobj;
        else
            error('options.computeobj must be boolean or a real number!');
        end
    else
        computeobj = true;
    end
    if computeobj == false & rep > 1
        error('options.computeobj must be true if options.rep > 1!');
    end
    if isfield(options, 'alg')
        alg_names = {'anls', 'newton'};
        j = find(strcmpi(options.alg, alg_names));
        if ~isempty(j)
            alg = alg_names{j};
        else
            error('Invalid options.alg value!');
        end
    else
        alg = 'anls';
    end
end

D = dist2(X, X);
if strcmp(graph_type, 'full') & strcmp(similarity_type, 'gaussian')
    A = scale_dist3(D, nn);
elseif strcmp(graph_type, 'full') & strcmp(similarity_type, 'inner_product')
    A = X * X';
elseif strcmp(graph_type, 'sparse') & strcmp(similarity_type, 'gaussian')
    A = scale_dist3_knn(D, nn, kk, true);
else % graph_type == 'sparse' & similarity_type == 'inner_product'
    Xnorm = X';
    d = 1./sqrt(sum(Xnorm.^2));
    Xnorm = bsxfun(@times, Xnorm, d);
    A = inner_product_knn(D, Xnorm, knn, true);
    clear Xnorm, d;
end
clear D;

if strcmp(graph_objfun, 'ncut')
    dd = 1 ./ sum(A);
    dd = sqrt(dd);
    A = bsxfun(@times, A, dd);
    A = A';
    A = bsxfun(@times, A, dd);
    clear dd;
end
A = (A + A') / 2;

params.maxiter = maxiter;
params.tol = tol;

obj_best = Inf;
for i = 1 : rep
    if isempty(Hinit)
        if strcmp(alg, 'newton')
            [sol, infos] = symm_newton(A, k, params);
        else % strcmp(alg, 'anls')
            [sol, infos] = symm_anls(A, k, params);
        end
    else
        params.Hinit = Hinit(:, :, i);
        if strcmp(alg, 'newton')
            [sol, infos] = symm_newton(A, k, params);
        else % strcmp(alg, 'anls')
            [sol, infos] = symm_anls(A, k, params);
        end
    end
    H = sol.H;
    [max_val, idx] = max(H, [], 2);
    if infos.cost(end) < obj_best
        idx_best = idx;
        iter_best = infos.iter(end);
        obj_best = infos.cost(end);
        H_best = H;
    end
end

idx = idx_best;
iter = iter_best;
obj = obj_best;
H = H_best;

author: hiroyuki-kasai

repo: NMFLibrary


text:
function varargout = quiver(sVF,varargin)
% quiver spherical vector field
%
% Syntax
%   quiver3(sVF)
%
% Options
%  normalized - normalize vectors
%  arrowSize  - arrow size
%  maxHeadSize - head size

% See also
%   S2VectorField/plot
%

% plot the function values
[varargout{1:nargout}] = plot(sVF,varargin{:});

end

author: mtex-toolbox

repo: mtex

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


text:
% Build a distribution of RTB

!make clean
!make
matlab.addons.toolbox.packageToolbox('RTB.prj', '~/Desktop/RTB')

author: petercorke

repo: robotics-toolbox-matlab


text:
function yn = has_mex_jmh
% for advanced features of this toolbox.

yn = which('jmh_sse_mex'); % seems to force a path check ??

yn =	...
	exist('jmh_sse_mex') == 3;

persistent Mexwarned
if isempty(Mexwarned)
	Mexwarned = 0;
end

if ~yn & ~Mexwarned
  warning('mex files not found: are your paths correct?')
  Mexwarned = 1;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function IMs = slinv2x2(Ms)
%SLINV2X2 Computes inverse matrices for 2 x 2 matrices in a fast way
%
% $ Syntax $
%   - IMs = slinv2x2(Ms)
%
% $ Arguments $
%   - Ms:       the 2 x 2 matrix (matrices)
%   - IMs:      the computed inverse matrix (matrices)
%
% $ Description $
%   - IMs = slinv2x2(Ms) computes the inverse of 2x2 matrices. If Ms is
%     a 2 x 2 matrix, then r is its inverse.
%     Or, if Ms is a 2 x 2 x ... array, then r would be a 2 x 2 x ... array
%     storing the inverse matrices.
%
% $ Remarks $
%   - The function uses the following formula for fast calculation of
%     inverse of 2 x 2 matrices:
%       inv(A) = [a22, -a12; -a21, a11] / det(A)
%
% $ History $
%   - Created by Dahua Lin on Apr 22nd, 2006
%

%% parse and verify input arguments

if size(Ms, 1) ~= 2 || size(Ms, 2) ~= 2
    error('sltoolbox:invalidarg', 'Ms should be set of 2 x 2 matrices');
end

%% compute

if ndims(Ms) == 2       % single matrix
    IMs = [Ms(4), -Ms(2); -Ms(3), Ms(1)] / (Ms(1) * Ms(4) - Ms(2) * Ms(3));

else                    % a set of matrices

    IMs = zeros(size(Ms));

    % compute adjacent matrix
    IMs(1,1,:) = Ms(2,2,:);
    IMs(1,2,:) = -Ms(1,2,:);
    IMs(2,1,:) = -Ms(2,1,:);
    IMs(2,2,:) = Ms(1,1,:);

    % compute the determinants
    d = sldet2x2(Ms);

    % scale
    d = reshape(d, [1, 1, size(d)]);
    IMs = slmul(IMs, 1 ./ d);

end

author: lmthang


text:
function results = run_OTBReadyTracker(seq, res_path, bSaveImage)
 % rename OTBReadyTracker to the name of the tracker you are using in the OTB evaluation
    startup;
    paths = env_paths_tracking();
    tracker_params.visualization = false;
    tracker_params.gpus = 1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    tracker_params.net = 'some-pretrained-network';
    tracker_params.net_gray = 'some-pretrained-network_gray';
    tracker_params.visualization = false;
    tracker_params.join.method = 'corrfilt'; % or 'xcorr'
    root_folder = 'your/path/to/otb/benchmark/root/';
    % % hyperparameters that work well for the tracker version you are using
    % tracker_params.scaleStep =
    % tracker_params.scalePenalty =
    % tracker_params.scaleLR =
    % tracker_params.wInfluence =
    % tracker_params.zLR =
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    tracker_params.imgFiles = vl_imreadjpeg(strcat(root_folder,seq.s_frames),'numThreads', 12);
    [cx, cy, w, h] = get_axis_aligned_BB(seq.init_rect);
    tracker_params.targetPosition = [cy cx];
    tracker_params.targetSize = round([h w]);
    % Call the main tracking function
    [bboxes, ~] = tracker(tracker_params);
    results = struct();
    results.res = bboxes;
    results.type = 'rect';
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
% VL_FISHER    Fisher vector feature encoding
%   ENC = VL_FISHER(X, MEANS, COVARIANCES, PRIORS) computes the Fisher
%   vector encoding of the vectors X relative to the Gaussian mixture
%   model with means MEANS, covariances COVARIANCES, and prior mode
%   probabilities PRIORS.
%
%   X has one column per data vector (e.g. a SIFT descriptor), and
%   MEANS and COVARIANCES one column per GMM component (covariance
%   matrices are assumed diagonal, hence these are simply the variance
%   of each data dimension). PRIORS has size equal to the number of
%   GMM components. All data must be of the same class, either SINGLE
%   or DOUBLE.
%
%   ENC is a vector of the same class of X of size equal to the
%   product of the data dimension and the number of components.
%
%   By default, the standard Fisher vector is computed. VL_FISHER()
%   accepts the following options:
%
%   Normalized::
%     If specified, L2 normalize the Fisher vector.
%
%   SquareRoot::
%     If specified, the signed square root function is applied to
%     ENC before normalization.
%
%   Improved::
%     If specified, compute the improved variant of the Fisher
%     Vector. This is equivalent to specifying the Normalized and
%     SquareRoot options.
%
%   Fast::
%     If specified, uses slightly less accurate computations but
%     significantly increase the speed in some cases (particularly
%     with a large number of Gaussian modes).
%
%   Verbose::
%     Increase the verbosity level (may be specified multiple times).
%
%   See: Fisher vectors, VL_HELP().

% Authors: David Novotny, Andrea Vedaldi

% Copyright (C) 2013 David Novotny and Andrea Vedaldi
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function pass = test_deriv(~)
% TEST_DERIV   Test the chebfun/deriv method

% This test creates some chebfun, evaluates the derivative and compares with the
% results of feval(diff(...)).

% Our CHEBFUN
u = chebfun(@(x) sin(exp(x)), [0 2]);

%% Evaluation at one point, default value of derivative
pass(1) = norm(deriv(u, 1) - feval(diff(u), 1)) == 0;

%% Higher order derivative
pass(2) = norm(deriv(u, .5, 3) - feval(diff(u, 3), .5)) == 0;

%% Default derivative at a vector of points
xx = linspace(0.1, 0.5, 11);
pass(3) = norm(deriv(u, xx) - feval(diff(u), xx)) == 0;

%% Higher derivative, vector of points
pass(4) = norm(deriv(u, xx, 4) - feval(diff(u, 4), xx)) == 0;

%% Array valued CHEBFUN at a vector of points
pass(5) = norm(deriv([u cos(u)], xx) - feval(diff([u cos(u)]), xx)) == 0;
pass(6) = norm(deriv([u cos(u)], xx, 2) - feval(diff([u cos(u)], 2), xx)) == 0;

%% Left- and right-sided evaluation
x = chebfun(@(x) x, [-1 1]);
v = cumsum(sign(x));
pass(7) = norm(feval(diff(v), 0, 'left') - deriv(v, 0, 'left')) == 0;
pass(8) = norm(feval(diff(v, 2), 0, '+') - deriv(v, 0, '+', 2)) == 0;

end


text:
function [Xmat,pobjround,info] = local_search_sra(Zmat,C,rrPar,opt,roundonly)
%% local search for single rotation averaging
%% used as a subroutine for STRIDE
if nargin < 5
    roundonly = false;
end
if nargin < 4
    % default round the first two eigenvectors
    opt = [1,2];
end

if roundonly
    [R,theta] = round_sra(Zmat,[1]);
    xtld      = lift_sra(R(:),theta);
    pobjround = xtld{1}' * C{1} * xtld{1};
    Xmat      = rank_one_lift(xtld);
else
    [R,theta] = round_sra(Zmat,opt);
    Zmatround = {};
    pobjround = zeros(length(opt),1);
    for i = 1:length(opt)
        Ri              = squeeze(R(:,:,i));
        thetai          = theta(:,i);
        [Ztmp,~,~,pobjtmp] = nlp_sra(C,Ri,thetai);
        Zmatround{end+1} = Ztmp;
        pobjround(i) = pobjtmp;
    end
    pobjs            = pobjround;
    [pobjround,idx]  = min(pobjround);
    if pobjround == inf
        fprintf('        NLP fails to find a good solution, return a rounded solution only.\n');
        Ropt         = squeeze(R(:,:,1));
        ropt         = Ropt(:);
        thetaopt     = theta(:,1);
        xopt         = lift_sra(ropt,thetaopt);
        Xmat         = rank_one_lift(xopt);
        pobjround    = xopt{1}'*C{1}*xopt{1};
        nlpsuccess   = false;
    else
        Xmat = Zmatround{idx};
        nlpsuccess   = true;
    end
end

if nargout > 2
    info.minidx     = idx;
    info.nlpsuccess = nlpsuccess;
    info.pobjs      = pobjs;
    info.diffpobj   = pobjs(1) - pobjround;
end

end

function f = sra_cost(x,C)
R       = x.A;
r       = R(:);
theta   = x.B;
theta   = theta(:);
v       = lift_sra(r,theta);
f       = v{1}'*C*v{1};
end

function g = sra_egrad(x,C)
R       = x.A;
r       = R(:);
theta   = x.B;
theta   = theta(:);
N       = length(theta);
v       = lift_sra(r,theta);
v       = v{1};

gv      = 2 * v' * C; % 1 x n

vdr     = [sparse(1,9);...
           speye(9);...
           sparse(N,9);...
           kron(theta,speye(9))]; % n x 9
gr      = (gv * vdr)'; % 9 x 1

vdtheta = [sparse(10,N);...
           speye(N);...
           kron(speye(N),r)]; % n x N
gtheta  = (gv * vdtheta); % 1 x N

g.A     = reshape(gr,3,3);
g.B     = gtheta;
end

function [Xmat,Ropt,thetaopt,fopt] = nlp_sra(C,R,theta)
N           = length(theta);
elements.A  = rotationsfactory(3,1);
elements.B  = obliquefactory(1,N);
manifold    = productmanifold(elements);

problem.M   = manifold;

warning('off', 'manopt:getHessian:approx')
% Define the problem cost function and its Euclidean gradient.
problem.cost  = @(x) sra_cost(x,C{1});
problem.egrad = @(x) sra_egrad(x,C{1});

% Numerically check gradient consistency (optional).
% checkgradient(problem);
% Solve.
x0.A                 = R;
x0.B                 = theta';
options.verbosity    = 0;
options.tolgradnorm  = 1e-6;
[xopt, fopt, output, options] = trustregions(problem,x0,options);
Ropt = xopt.A;
thetaopt = xopt.B';
constraintviolationR          = norm(Ropt*Ropt'-eye(3),'fro');
constraintviolationtheta      = max(abs(thetaopt.^2 - 1));
firstorderopt                 = output(end).gradnorm;
fprintf('        MANOPT: itr: %3d, constraint violation: %3.2e, %3.2e, gradnorm: %3.2e, cost: %3.8e.\n',...
    length(output),constraintviolationR,constraintviolationtheta,firstorderopt,fopt);
if max(constraintviolationR,constraintviolationtheta) < 1e-8 ...
        && firstorderopt < 1e-6
    % Do nothing
else
    fopt    = inf;
end
vnew     = lift_sra(Ropt(:),thetaopt(:));
Xmat     = rank_one_lift(vnew);
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
%% Copyright (C) 2014-2016, 2018-2019 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
%% @deftypemethod  @@sym {} disp (@var{x})
%% @deftypemethodx @@sym {} disp (@var{x}, 'unicode')
%% @deftypemethodx @@sym {} disp (@var{x}, 'ascii')
%% @deftypemethodx @@sym {} disp (@var{x}, 'flat')
%% @deftypemethodx @@sym {@var{s} =} disp (@var{x})
%% Display the value of a symbolic expression.
%%
%% Examples:
%% @example
%% @group
%% syms x a c
%% str = disp(sin(2*sym(pi)*x))
%%   @result{} str =   sin(2⋅π⋅x)
%%
%% A = [sin(x/2) floor(a^(x*c)); acosh(2*x/pi) ceil(sin(x/gamma(x)))];
%% disp(A, 'unicode')
%%   @print{}   ⎡     ⎛x⎞      ⎢ c⋅x⎥   ⎤
%%   @print{}   ⎢  sin⎜─⎟      ⎣a   ⎦   ⎥
%%   @print{}   ⎢     ⎝2⎠               ⎥
%%   @print{}   ⎢                       ⎥
%%   @print{}   ⎢     ⎛2⋅x⎞  ⎡   ⎛ x  ⎞⎤⎥
%%   @print{}   ⎢acosh⎜───⎟  ⎢sin⎜────⎟⎥⎥
%%   @print{}   ⎣     ⎝ π ⎠  ⎢   ⎝Γ(x)⎠⎥⎦
%% @end group
%%
%% @group
%% disp(A, 'ascii')
%%   @print{}   [     /x\              / c*x\      ]
%%   @print{}   [  sin|-|         floor\a   /      ]
%%   @print{}   [     \2/                          ]
%%   @print{}   [                                  ]
%%   @print{}   [     /2*x\         /   /   x    \\]
%%   @print{}   [acosh|---|  ceiling|sin|--------||]
%%   @print{}   [     \ pi/         \   \Gamma(x)//]
%%
%% disp(A, 'flat')
%%   @print{} Matrix([[sin(x/2), floor(a**(c*x))], [acosh(2*x/pi), ceiling(sin(x/gamma(x)))]])
%% @end group
%% @end example
%%
%% @seealso{@@sym/pretty}
%% @end deftypemethod

function varargout = disp(x, wh)

  if (nargin == 1)
    %% read config to see how to display x
    wh = sympref('display');
  end

  switch lower(wh)
    case 'flat'
      s = x.flat;
    case 'ascii'
      s = x.ascii;
    case 'unicode'
      s = x.unicode;
    otherwise
      print_usage ();
  end
  s = make_indented(s);

  if (nargout == 0)
    disp(s)
  else
    varargout = {[s sprintf('\n')]};  % add a newline
  end
end

function s = make_indented(s, n)
  if (nargin == 1)
    n = 2;
  end
  pad = char (double (' ')*ones (1,n));
  newl = sprintf('\n');
  s = strrep (s, newl, [newl pad]);
  s = [pad s];  % first line
end

%!test
%! syms x
%! s = disp(sin(x));
%! assert(strcmp(s, sprintf('  sin(x)\n')))

%!test
%! syms x
%! s = disp(sin(x/2), 'flat');
%! assert(strcmp(s, sprintf('  sin(x/2)\n')))

%!test
%! % Examples of 2x0 and 0x2 empty matrices:
%! a = sym([1 2; 3 4]);
%! b2x0 = a([true true], [false false]);
%! b0x2 = a([false false], [true true]);
%! assert (isequal (size (b2x0), [2 0]))
%! assert (isequal (size (b0x2), [0 2]))
%! s = disp(b2x0);
%! assert(strcmp(s, sprintf('  []\n')))
%! s = disp(b0x2);
%! assert(strcmp(s, sprintf('  []\n')))

repo: octsympy


text:
function [dat, state] = ft_preproc_standardize(dat, begsample, endsample, state)

% FT_PREPROC_STANDARDIZE performs a z-transformation or standardization
% of the data. The standardized data will have a zero-mean and a unit
% standard deviation.
%
% Use as
%   [dat] = ft_preproc_standardize(dat, begsample, endsample)
% where
%   dat        data matrix (Nchans dat Ntime)
%   begsample  index of the begin sample for the mean and stdev estimate
%   endsample  index of the end sample for the mean and stdev estimate
%
% If no begin and end sample are specified, it will be estimated on the
% complete data.
%
% See also PREPROC

% Copyright (C) 2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin<2 || isempty(begsample)
  begsample = 1;
end

if nargin<3 || isempty(endsample)
  endsample = size(dat,2);
end

if nargin<4
  state = [];
end

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

% get the data selection
y = dat(:,begsample:endsample);

% determine the size of the selected data: nChans dat nSamples
[m, n] = size(y);

% compute the sum and sum of squares
s  = sum(y,2);
ss = sum(y.^2,2);

% include the state information from the previous calls
if ~isempty(state)
  s  = s  + state.s;
  ss = ss + state.ss;
  n  = n  + state.n;
end

% compute the mean and standard deviation
my = s ./ n;
sy = sqrt((ss - (s.^2)./n) ./ (n-1));

% standardize the complete input data
dat = (dat - repmat(my, 1, size(dat, 2))) ./ repmat(sy, 1, size(dat, 2));

% remember the state
state.s  = s;  % sum
state.ss = ss; % sum of sqares
state.n  = n;  % number of samples

author: spm


text:
function [G,H] = mass_spring_gradient_hessian(V,E,k,U)

  % Compute the gradient and hessian of a mass-spring system defined for (V,E)
  % with the stiffness k and displacement U
  % Inputs:
  %   V  #V by 2 list of rest curve vertex positions
  %   E  #E by 2 list of edge indices
  %   k  scalar spring stiffness parameter
  %   U  #V by 2 list of deformed curve vertex displacements from the rest pose
  % Outputs:
  %   G  #V*2 list of gradient vectors per vertex defined at the current deformed position
  %   H  #V*2 by #V*2 hessian matrix defined at the current deformed position

  assert(size(E,2) == 2);
  U = reshape(U,size(V));

  % rest edge lengths
  LV = edge_lengths(V,E);
  LU = edge_lengths(U,E);
  f = k*0.5*sum((LV-LU).^2);

  if nargout==0
    return;
  end

  I = E(:,2);
  J = E(:,1);
  T0 = U(I,:) - U(J,:);
  t1 = normrow(T0);

  GI = ((t1-LV)./t1).*T0;
  GJ = -GI;

  dim = size(V,2);
  G = k*full(sparse([repmat(I,1,dim) repmat(J,1,dim)],repmat(1:dim,size(E,1),2),[GI GJ],size(V,1),dim));
  G = G(:); % to match the dimension: not sure if correct

  if nargout==1
    return;
  end

  % Outer product
  vec = @(X) X(:);
  T2 = T0(:,vec(repmat(1:dim,dim,1))').*T0(:,vec(repmat(1:dim,1,dim))');
  t3 = t1 - LV;

  HIJ = (t3./(t1.^3) - 1./(t1.^2)).*T2 - t3./t1.*(vec(eye(dim))');
  HII = -HIJ;
  HJJ = -HIJ;
  HJI = HIJ;
  n = size(V,1);
  H = sparse(n*dim,n*dim);
  for ii = 1:dim
    for jj = 1:dim
      oi = (ii-1)*dim+jj;
      Hiijj = sparse( ...
        [I I J J],[I J I J],[HII(:,oi) HIJ(:,oi) HJI(:,oi) HJJ(:,oi)],n,n);
      H(((ii-1)*n)+(1:n),((jj-1)*n)+(1:n)) = Hiijj;
    end
  end
  H = k*H;

end

author: odedstein

repo: sgi-introduction-course


text:
function star = getStarVerticies(n,r)
% star = getStarVerticies(n,r)
%
% This function draws a star with n points centered on the origin.
%
% INPUTS:
%   n = number of points
%   r = inner radius (outer radius == 1)
%
% OUTPUTS:
%   star = [2, 2*n+1] = [x;y] coordinates of verticies
%

th = linspace(0,2*pi,2*n+1);
star = [...
    sin(th);
    cos(th)];
idx = 2:2:(2*n);
star(:,idx) = r*star(:,idx);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [embedding_layer_state, hidden_layer_state, output_layer_state] = ...
  fprop(input_batch, word_embedding_weights, embed_to_hid_weights,...
  hid_to_output_weights, hid_bias, output_bias)
% This method forward propagates through a neural network.
% Inputs:
%   input_batch: The input data as a matrix of size numwords X batchsize where,
%     numwords is the number of words, batchsize is the number of data points.
%     So, if input_batch(i, j) = k then the ith word in data point j is word
%     index k of the vocabulary.
%
%   word_embedding_weights: Word embedding as a matrix of size
%     vocab_size X numhid1, where vocab_size is the size of the vocabulary
%     numhid1 is the dimensionality of the embedding space.
%
%   embed_to_hid_weights: Weights between the word embedding layer and hidden
%     layer as a matrix of soze numhid1*numwords X numhid2, numhid2 is the
%     number of hidden units.
%
%   hid_to_output_weights: Weights between the hidden layer and output softmax
%               unit as a matrix of size numhid2 X vocab_size
%
%   hid_bias: Bias of the hidden layer as a matrix of size numhid2 X 1.
%
%   output_bias: Bias of the output layer as a matrix of size vocab_size X 1.
%
% Outputs:
%   embedding_layer_state: State of units in the embedding layer as a matrix of
%     size numhid1*numwords X batchsize
%
%   hidden_layer_state: State of units in the hidden layer as a matrix of size
%     numhid2 X batchsize
%
%   output_layer_state: State of units in the output layer as a matrix of size
%     vocab_size X batchsize
%

[numwords, batchsize] = size(input_batch);
[vocab_size, numhid1] = size(word_embedding_weights);
numhid2 = size(embed_to_hid_weights, 2);

%% COMPUTE STATE OF WORD EMBEDDING LAYER.
% Look up the inputs word indices in the word_embedding_weights matrix.
embedding_layer_state = reshape(...
  word_embedding_weights(reshape(input_batch, 1, []),:)',...
  numhid1 * numwords, []);

%% COMPUTE STATE OF HIDDEN LAYER.
% Compute inputs to hidden units.
inputs_to_hidden_units = embed_to_hid_weights' * embedding_layer_state + ...
  repmat(hid_bias, 1, batchsize);

% Apply logistic activation function.
% FILL IN CODE. Replace the line below by one of the options.
% hidden_layer_state = zeros(numhid2, batchsize);
% Options
% (a) hidden_layer_state = 1 ./ (1 + exp(inputs_to_hidden_units));
% (b) hidden_layer_state = 1 ./ (1 - exp(-inputs_to_hidden_units));
hidden_layer_state = 1 ./ (1 + exp(-inputs_to_hidden_units));
% (d) hidden_layer_state = -1 ./ (1 + exp(-inputs_to_hidden_units));

%% COMPUTE STATE OF OUTPUT LAYER.
% Compute inputs to softmax.
% FILL IN CODE. Replace the line below by one of the options.
% inputs_to_softmax = zeros(vocab_size, batchsize);
% Options
inputs_to_softmax = hid_to_output_weights' * hidden_layer_state +  repmat(output_bias, 1, batchsize);
% (b) inputs_to_softmax = hid_to_output_weights' * hidden_layer_state +  repmat(output_bias, batchsize, 1);
% (c) inputs_to_softmax = hidden_layer_state * hid_to_output_weights' +  repmat(output_bias, 1, batchsize);
% (d) inputs_to_softmax = hid_to_output_weights * hidden_layer_state +  repmat(output_bias, batchsize, 1);

% Subtract maximum.
% Remember that adding or subtracting the same constant from each input to a
% softmax unit does not affect the outputs. Here we are subtracting maximum to
% make all inputs <= 0. This prevents overflows when computing their
% exponents.
inputs_to_softmax = inputs_to_softmax...
  - repmat(max(inputs_to_softmax), vocab_size, 1);

% Compute exp.
output_layer_state = exp(inputs_to_softmax);

% Normalize to get probability distribution.
output_layer_state = output_layer_state ./ repmat(...
  sum(output_layer_state, 1), vocab_size, 1);

repo: neural-nets


text:
% At_fhp.m
%
% Adjoint of At_fhp (2D Fourier half plane measurements).
%
% Usage: x = At_fhp(b, OMEGA, n)
%
% b - K vector = [mean; real part(OMEGA); imag part(OMEGA)]
%
% OMEGA - K/2-1 vector denoting which Fourier coefficients to use
%         (the real and imag parts of each freq are kept).
%
% n - Image is nxn pixels
%
% x - N vector
%
% Written by: Justin Romberg, Caltech
% Created: October 2005
% Email: jrom@acm.caltech.edu
%

function x = At_fhp(y, OMEGA, n)

K = length(y);

fx = zeros(n,n);
fx(1,1) = y(1);
fx(OMEGA) = sqrt(2)*(y(2:(K+1)/2) + i*y((K+3)/2:K));
x = reshape(real(n*ifft2(fx)), n*n, 1);

repo: REST


text:
function [Model, Info] = linear_stepwise_reg(X,Y,Model,parm)
%  Estimate linear weight matrix for input-output mapping
%     regularization version
%
% Notice !!
%     Delay embedding is done in this module,
%     then input vector should be original input data without embedding
%
%   [Model, Info] = linear_stepwise_reg(X,Y,Model,parm)
%     Scalar iteration version
% --- Input
%  Y  : Output data ( N x T x Ntrial )
%  X  : Input data  ( M x (T + (Dtau-1)*Tau)  x Ntrial)
%  N  =  # of output
%  M  =  # of input (original input space dimension)
%  T  =  # of time sample
%
%  Estimate the following model
%    Y(t) = W * [X(:, t + (Dtau-1)*Tau ); ...; X(:, t)]
%
%  Model : Structure for estimated model
%  if Model is empty, initialization is done before training
%  if Model is previous training result, re-training us done
%
%  parm  : Structure for learning parameter
%  parm.Npre_train :  # of VB-update in initial training
%  parm.Ntrain :  # of training
%  parm.Nskip  :  skip # for print
%  parm.a_min  :  Min value for pruning small variance component
%  parm.Prune  :  = 1 : Prune small variance & irrelevant input dimension
%
%  parm.Tau       = Lag time
%  parm.Dtau      = Number of embedding dimension
%    Total input dimension in embedding space is (M * parm.Dtau)
% --- Output
%  Model : Structure for estimated model
%  Model.SY  :  Noise variance         ( 1 x 1 )
%  Model.W   :  Weight matrix          ( N x M*D ) , D = parm.Dtau
%  Model.A   :  Prior weight variance  ( 1 x M*D )
%  Model.ix_act : Active index for W after pruning
%
%  Info  : Structure for learning process history
%  Info.FE  = LP + H : Free energy
%  Info.LP  = Log likelihood
%  Info.H   = - Model entropy
%
% 2009-10-18 Made by M. Sato

MINVAL = 1.0e-15;

fprintf('linear_stepwise (regularization) start\n')

% Dimension
[N ,T  ,Ntrial]  = size(Y); % N = # of output
[M ,Tx ,Ntrialx] = size(X); % M = # of input without embedding

Tall = T * Ntrial;

if Ntrial~=Ntrialx, error('# of trial is different in X and Y'); end;

Nskip  = 100;   % skip steps for display info
a_min = 1e-14;	% Minimum value for weight pruning
Fdiff = 1e-12; % Threshold for convergence
Ncheck = 100;   % Minimum number of training iteration
Fstep  = 5;     % Free energy convergence check step

if isfield(parm,'Nskip'), Nskip  = parm.Nskip; end;
if isfield(parm,'Fdiff'), Fdiff   = parm.Fdiff; end;
if isfield(parm,'a_min'), a_min   = parm.a_min ; end;
if isfield(parm,'Ncheck'), Ncheck = parm.Ncheck; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end

% # of embedding dimension
if isfield(parm,'Dtau')
	D    = parm.Dtau;
	tau  = parm.Tau;
else
	D    = 1;
	tau  = 1;
end
% Parameters
Ntrain = parm.Ntrain;

if Ntrain < 1, Info = []; return; end;

fprintf('--- Output Dimension  = %d\n',N)
fprintf('--- Input  Dimension  = %d\n',M)
fprintf('--- Embedding  Dimension  = %d\n',D)
fprintf('--- Number of trials  = %d\n',Ntrial)
fprintf('--- Number of training sample = %d\n',T)
fprintf('--- Total update iteration    = %d \n',Ntrain)

% original input dimension
Xdim  = M;
M = Xdim*D;
M_ALL = M;

%
% --- Initialization
%

% Input/Output variance
sx = mean((X(:) - mean(X(:))).^2);
sy = mean((Y(:) - mean(Y(:))).^2);

SY0 = mean(sy);
A0  = mean(sx);

% Input covariance
XX  = sum(sum(X.^2,3),2)/(Tx*Ntrial);
XX  = repmat(XX', [1 D]);% 1 x M

if isfield(Model,'A')
	A = Model.A;
	A = mean(A);
else
	A = A0;
	W = zeros(N,M_ALL);
end

if size(A,1)==Xdim && size(A,2)==1
	A   = repmat(A',[1 D]);
elseif size(A,1)==1 && size(A,2)==Xdim
	A   = repmat(A,[1 D]);
elseif size(A,1)==1 && size(A,2)==1
	A   = repmat(A,[1 M_ALL]);
end

W = zeros(N,M_ALL);
SY  = SY0;

% Save original input for pruning
IX_act  = 1:M_ALL;
IX_dim  = 1:Xdim;
ix_act  = 1:M;
ix_dim  = 1:Xdim;
M_all   = M;

fprintf('SY0 = %g\n', SY0)
fprintf('SY  = %g\n', SY)

% Free energy histry
FE  = zeros(Ntrain,1);
LP  = zeros(Ntrain,1);
H   = zeros(Ntrain,1);
Err = zeros(Ntrain,1);
Mhist = zeros(Ntrain,1);

[dY] = error_delay_time(X,Y,W,tau);

k_save  = 0;
% Ainv = alpha , A = 1/alpha
Ainv    = A;

%%%%%% Learning Loop %%%%%%
for k=1:Ntrain

	% E = (Y-W*X)^2/SY +  W^2/A
	%   = ( (Y-W*X)^2 +  W^2 * (SY/A) )/SY
    [W]  = weight_update_embed(X, dY, W, Tall*XX, Ainv, tau);
	[dY] = error_delay_time(X, Y, W, tau);

    dYY = sum(dY(:).^2)/(N*Tall);
	WW  = sum(W.^2,1);

    % Noise variance update
    SY  = dYY + sum(WW .* Ainv)/(N*Tall);
    % Prevent zero variance
    SY  = max( SY, MINVAL);

    % Log variance
    log_sy  = N*( log(SY) );

    % Free energy
    LP(k)  = - 0.5 * log_sy ;
    FE(k)  = LP(k) ;
    Err(k) = sum(dYY)/(SY0);

    if mod(k, Nskip)==0
        fprintf('Iter = %4d, err = %g, WW=%g, \n', ...
               k, Err(k), mean(WW));
    end

    % Convergence check
	if k > Ncheck,
		Fdif = (FE(k) - FE(k-Fstep))/(abs(FE(k))+eps);
	else
		Fdif = Fdiff + 1;
	end

	if (Fdiff > abs(Fdif)),
		fprintf('Converged : Free energy change = %g\n',Fdif)
		break;
	end;
end

ix_act = IX_act(ix_act);

% Active index
%Model.ix_act = ix_act;
Model.M_all  = M_ALL ;

% Save output variable
Model.A    = A ;
Model.W    = W ;
Model.SY   = SY;
Model.method = 'linear_stepwise';
Model.mode   = 'scalar';
Model.sparse = 'reg';

% Save history
Info.FE  = FE(1:k);
Info.LP  = LP(1:k);
Info.H   = H(1:k) ;
Info.Err = Err(1:k);

author: KamitaniLab

repo: GenericObjectDecoding


text:
%DATE:21/04/2005
%COFDM SIMULATOR
%BY:Apurva Gupta
%B.ENG,ELECTRONICs & COMMUNICATION ENGINEERING,
%UNIVERSITY OF LEEDS,UK
%SUPERVISOR:Prof.G.Markarian

function varargout = final_simulator(varargin)
% FINAL_SIMULATOR M-file for final_simulator.fig
%      FINAL_SIMULATOR, by itself, creates a new FINAL_SIMULATOR or raises the existing
%      singleton*.
%
%      H = FINAL_SIMULATOR returns the handle to a new FINAL_SIMULATOR or the handle to
%      the existing singleton*.
%
%      FINAL_SIMULATOR('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in FINAL_SIMULATOR.M with the given input arguments.
%
%      FINAL_SIMULATOR('Property','Value',...) creates a new FINAL_SIMULATOR or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before final_simulator_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to final_simulator_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help final_simulator

% Last Modified by GUIDE v2.5 24-Apr-2005 12:30:13

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @final_simulator_OpeningFcn, ...
                   'gui_OutputFcn',  @final_simulator_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before final_simulator is made visible.
function final_simulator_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to final_simulator (see VARARGIN)

% Choose default command line output for final_simulator
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes final_simulator wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = final_simulator_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function [varargout] = code1(varargin)
M=2;
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 10; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo; % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      msg= zeros(1,1e6);%message is a string of a thousand zeros.
      %t = poly2trellis([4 3],[4 5 17;7 4 2]); % Define trellis.
      %code = convenc(msg,t); % Encode a string of ones.
      txsig = pskmod(msg,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      %qcode = quantiz(decodmsg,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      %tblen = 48; delay = tblen; % Traceback length
      %decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decodmsg); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * 1e6),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
%hold off;
if nargout == 0
  builtin('code1', varargin{:});
else
  [varargout{1:nargout}] = builtin('code1', varargin{:});
end

%--------------------------------------------------------------
function [varargout] = code2(varargin)
M=2;
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 10; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo; % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      msg= zeros(1,1e6);%message is a string of a thousand zeros.
      %t = poly2trellis([4 3],[4 5 17;7 4 2]); % Define trellis.
      %code = convenc(msg,t); % Encode a string of ones.
      txsig = pskmod(msg,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      %qcode = quantiz(decodmsg,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      %tblen = 48; delay = tblen; % Traceback length
      %decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decodmsg); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * 1e6),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
%hold off;
if nargout == 0
  builtin('code2', varargin{:});
else
  [varargout{1:nargout}] = builtin('code2', varargin{:});
end

%-------16QAM-------------------------------------------
function [varargout] = code3(varargin)
bits=1e5;
M=16;
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 10; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo;%*0.5*log2(M);% Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      msg= zeros(bits,1);%message is a string of a thousand zeros.
      txsig = qammod(msg,M); % Modulate.
      noise = awgn(txsig, snr,'measured'); % Add noise.
      decodmsg = qamdemod(noise,M); % Demodulate.
      newerrs = biterr(msg,decodmsg); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * bits),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%set(gcf,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1])
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
if nargout == 0
  builtin('code3', varargin{:});
else
  [varargout{1:nargout}] = builtin('code3', varargin{:});
end

%------64QAM------------------------------------------------
function [varargout] = code4(varargin)
bits=1e5;
M=64;
EbNomin = 5; EbNomax = 15; % EbNo range, in dB
numerrmin = 5; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo;%*0.5*log2(M);% Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      msg= zeros(bits,1);%message is a string of a thousand zeros.
      txsig = qammod(msg,M); % Modulate.
      noise = awgn(txsig, snr,'measured'); % Add noise.
      decodmsg = qamdemod(noise,M); % Demodulate.
      newerrs = biterr(msg,decodmsg); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * bits),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%set(gcf,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1])
%axes(handles.axes2);
for jj=1:numEbNos
    %legend('Emperical BER','64QAM','Poly Ratio fit');
    semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
if nargout == 0
  builtin('code4', varargin{:});
else
  [varargout{1:nargout}] = builtin('code4', varargin{:});
end

function [varargout] = conv_code1(varargin)
M=2;
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 5; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo; % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      msg= zeros(1,1e6);%message is a string of a thousand zeros.
      t = poly2trellis([4 3],[4 5 17;7 4 2]); % Define trellis.
      code = convenc(msg,t); % Encode a string of ones.
      txsig = pskmod(code,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      qcode = quantiz(decodmsg,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      tblen = 48; delay = tblen; % Traceback length
      decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decoded); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * 1e6),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
%hold off;
if nargout == 0
  builtin('conv_code1', varargin{:});
else
  [varargout{1:nargout}] = builtin('conv_code1', varargin{:});
end
%------------------------------------------------------------
function [varargout] = conv_code2(varargin)
M=4;
EbNomin = 1; EbNomax = 10; % EbNo range, in dB
numerrmin = 5; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo; % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      msg= zeros(1,1e6);%message is a string of a thousand zeros.
      t = poly2trellis([4 3],[4 5 17;7 4 2]); % Define trellis.
      code = convenc(msg,t); % Encode a string of ones.
      txsig = pskmod(code,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      qcode = quantiz(decodmsg,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      tblen = 48; delay = tblen; % Traceback length
      decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decoded); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * 1e6),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
%hold off;
if nargout == 0
  builtin('conv_code2', varargin{:});
else
  [varargout{1:nargout}] = builtin('conv_code2', varargin{:});
end
%-------------------------------------------------------------
function [varargout] = conv_code3(varargin)
M=2;
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 10; % Compute BER only after 10 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo;%*0.5*log2(M); % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)

      msg = zeros(1e6,1);%message is a string of a thousand zeros.
      t = poly2trellis(7,[171 133]); % Define trellis.
      code = convenc(msg,t); % Encode a string of ones.
      txsig = pskmod(code,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      qcode = quantiz(decodmsg,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      tblen = 48; delay = tblen; % Traceback length
      decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decoded); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * 1e6),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
%hold off;
if nargout == 0
  builtin('conv_code3', varargin{:});
else
  [varargout{1:nargout}] = builtin('conv_code3', varargin{:});
end

%-------------------------------------------------------
function [varargout] = conv_code4(varargin)
M=4;
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 10; % Compute BER only after 10 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo*0.5*log2(M); % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)

      msg = zeros(1e6,1);%message is a string of a thousand zeros.
      t = poly2trellis(7,[171 133]); % Define trellis.
      code = convenc(msg,t); % Encode a string of ones.
      txsig = pskmod(code,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      qcode = quantiz(decodmsg,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      tblen = 48; delay = tblen; % Traceback length
      decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decoded); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials * 1e6),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
%hold off;
if nargout == 0
  builtin('conv_code4', varargin{:});
else
  [varargout{1:nargout}] = builtin('conv_code4', varargin{:});
end

% --- Executes on selection change in modulation.
function modulation_Callback(hObject, eventdata, handles)
% hObject    handle to modulation (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns modulation contents as cell array
%        contents{get(hObject,'Value')} returns selected item from modulation

% --- Executes during object creation, after setting all properties.
function modulation_CreateFcn(hObject, eventdata, handles)
% hObject    handle to modulation (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popupmenu2.
function popupmenu2_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popupmenu2 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu2

% --- Executes during object creation, after setting all properties.
function popupmenu2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in coding.
function coding_Callback(hObject, eventdata, handles)
% hObject    handle to coding (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns coding contents as cell array
%        contents{get(hObject,'Value')} returns selected item from coding

% --- Executes during object creation, after setting all properties.
function coding_CreateFcn(hObject, eventdata, handles)
% hObject    handle to coding (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in interleaving.
function interleaving_Callback(hObject, eventdata, handles)
% hObject    handle to interleaving (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns interleaving contents as cell array
%        contents{get(hObject,'Value')} returns selected item from interleaving

% --- Executes during object creation, after setting all properties.
function interleaving_CreateFcn(hObject, eventdata, handles)
% hObject    handle to interleaving (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on button press in int.
function int_Callback(hObject, eventdata, handles)
% hObject    handle to int (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
popup_sel_index = get(handles.interleaving, 'Value');
switch popup_sel_index
    case 1

bits =1e6
M=2
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 10; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo; % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      % st1 = 27221; st2 = 4831; % States for random number generator
      msg= zeros(bits,1);%message is a string of a thousand zeros.
      t = poly2trellis(7,[171 133]); % Define trellis.
      %t = poly2trellis([4 3],[4 5 17;7 4 2]); % Define trellis.
      code = convenc(msg,t); % Encode a string of ones.
      nrows = 5; % Use 5 shift registers
      slope = 3; % Delays are 0, 3, 6, 9, and 12.
      inter = convintrlv(code,nrows,slope); % Interleave.
      txsig = pskmod(inter,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      deinter = convdeintrlv(decodmsg,nrows,slope); % Deinterleave.
      qcode = quantiz(deinter,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      tblen = 48; delay = tblen; % Traceback length
      decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decoded); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials *bits ),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%set(gcf,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1])
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end

    case 2

bits =1e6
M=4
EbNomin = 1; EbNomax = 5; % EbNo range, in dB
numerrmin = 10; % Compute BER only after 5 errors occur.
EbNovec = EbNomin:1:EbNomax; % Vector of EbNo values
numEbNos = length(EbNovec); % Number of EbNo values
% Preallocate space for certain data.
ber = zeros(1,numEbNos); % BER values
intv = cell(1,numEbNos); % Cell array of confidence intervals
% Loop over the vector of EbNo values.
for jj = 1:numEbNos
   EbNo = EbNovec(jj);
   snr = EbNo*0.5*log2(M); % Because of binary modulation
   ntrials = 0; % Number of passes through the while loop below
   numerr = 0; % Number of errors for this EbNo value
   % Simulate until numerrmin errors occur.
   while (numerr < numerrmin)
      % st1 = 27221; st2 = 4831; % States for random number generator
      msg= zeros(bits,1);%message is a string of a thousand zeros.
      t = poly2trellis(7,[171 133]); % Define trellis.
      %t = poly2trellis([4 3],[4 5 17;7 4 2]); % Define trellis.
      code = convenc(msg,t); % Encode a string of ones.
      nrows = 5; % Use 5 shift registers
      slope = 3; % Delays are 0, 3, 6, 9, and 12.
      inter = convintrlv(code,nrows,slope); % Interleave.
      txsig = pskmod(inter,M); % Modulate.
      rxsig = awgn(txsig, snr, 'measured'); % Add noise.
      decodmsg = pskdemod(rxsig,M); % Demodulate.
      deinter = convdeintrlv(decodmsg,nrows,slope); % Deinterleave.
      qcode = quantiz(deinter,[0.001,.1,.3,.5,.7,.9,.999]);% Quantize to prepare for soft-decision decoding.
      tblen = 48; delay = tblen; % Traceback length
      decoded = vitdec(qcode,t,tblen,'cont','soft',3); % Decode.
      newerrs = biterr(msg,decoded); % Errors in this trial
      numerr = numerr + newerrs; % Total errors for this EbNo value
      ntrials = ntrials + 1; % Update trial index.
   end
   % Error rate and 98% confidence interval for this EbNo value
   [ber(jj), intv1] = berconfint(numerr,(ntrials *bits ),.98);
   intv{jj} = intv1; % Store in cell array for later use.
   disp(['EbNo = ' num2str(EbNo) ' dB, ' num2str(numerr) ...
         ' errors, BER = ' num2str(ber(jj))])
end
% Use BERFIT to plot the best fitted curve,
% interpolating to get a smooth plot.
fitEbNo = EbNomin:0.25:EbNomax; % Interpolation values
berfit(EbNovec,ber,fitEbNo);

% Also plot confidence intervals.
hold on;
grid on;
%set(gcf,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1])
%axes(handles.axes2);
for jj=1:numEbNos
   semilogy([EbNovec(jj) EbNovec(jj)],intv{jj},'g-+');
end
end

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
popup_sel_index = get(handles.modulation, 'Value');
switch popup_sel_index
    case 1
popup_sel_index = get(handles.coding, 'Value');
switch popup_sel_index
    case 1
        code1;
    case 2
        conv_code3;
    case 3
        conv_code1;

end
    case 2
popup_sel_index = get(handles.coding, 'Value');
switch popup_sel_index
    case 1
        code2;
    case 2
        conv_code4;
    case 3
        conv_code2;
end
    case 3
        code3;
    case 4
        code4;
end

% --- Executes on button press in clr.
function clr_Callback(hObject, eventdata, handles)
% hObject    handle to clr (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
cla

% --- Executes on button press in hlp.
function hlp_Callback(hObject, eventdata, handles)
% hObject    handle to hlp (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
helpPath = which('help.txt');
web(helpPath);

% --------------------------------------------------------------------
function close_Callback(hObject, eventdata, handles)
% hObject    handle to close (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
selection = questdlg(['Close ' get(handles.figure1,'Name') '?'],...
                     ['Close ' get(handles.figure1,'Name') '...'],...
                     'Yes','No','Yes');
if strcmp(selection,'No')
    return;
end

delete(handles.figure1)

% --------------------------------------------------------------------
function print_Callback(hObject, eventdata, handles)
% hObject    handle to print (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
printdlg(handles.figure1)

% --------------------------------------------------------------------
function file_Callback(hObject, eventdata, handles)
% hObject    handle to file (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

author: Sable

repo: mcbench-benchmarks


text:
addpath('../src/')
addpath('../src/utils/')

% specify your own discrete linear system
A = [1 1; 0 1];
B = [0.5; 1];
Q = diag([1, 1]);
R = 0.1;

% construct a convex set of system noise (2dim here)
W_vertex = [0.15, 0.15; 0.15, -0.15; -0.15, -0.15; -0.15, 0.15];
W = Polyhedron(W_vertex);

% construct disturbance Linear system
% note that disturbance invariant set Z is computed and stored as member variable in the constructor.
disturbance_system = DisturbanceLinearSystem(A, B, Q, R, W);

% you can see that with any disturbance bounded by W, the state is guaranteed to inside Z
x = zeros(2); % initial state
for i = 1:50
    u = disturbance_system.K * (x - 0);
    x = disturbance_system.propagate(x, u); % disturbance is considered inside the method
    clf;
    Graphics.show_convex(disturbance_system.Z, 'g', 'FaceAlpha', .3); % show Z
    scatter(x(1, :), x(2, :)); % show particle
    pause(0.01);
end

author: HiroIshida

repo: robust-tube-mpc


text:
function prob_test101 ( )

%*****************************************************************************80
%
%% TEST101 tests LOG_SERIES_CDF, LOG_SERIES_CDF_INV, LOG_SERIES_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST101\n' );
  fprintf ( 1, '  For the Logseries PDF,\n' );
  fprintf ( 1, '  LOG_SERIES_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  LOG_SERIES_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  LOG_SERIES_PDF evaluates the PDF;\n' );

  a = 0.25;

  check = log_series_check ( a );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST101 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =  %14f\n', a );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = log_series_sample ( a, seed );

    pdf = log_series_pdf ( x, a );

    cdf = log_series_cdf ( x, a );

    x2 = log_series_cdf_inv ( cdf, a );

    fprintf ( 1, '  %14d  %14f  %14f  %14d\n', x, pdf, cdf, x2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y]=funcrs2(tt,fun,eps,y,nswp)
%Cross approximation of a function of a TT-tensor, Method 2
%   [Y]=FUNCRS2(TT,FUN,EPS,Y,NSWP)
%   Computes approximation to the function FUN(TT) with accuracy EPS
%   Auxiliary parameters:  Y (initial approximation), NSWP
%   (number of sweeps in the method
%   Much faster then usual cross by vectorized computation of subtensors
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

%PARAMETERS SECTION
rmin=1;
verb=true;
kick_rank=5;
if (~isempty(y))
   yold=y;
end
y1=y;
%For this procedure we need only local (!) indices, since it
%is more or less equivalent to orthogonalization;
d=tt.d;
ps=tt.ps;
core0=tt.core;
n=tt.n;
r=tt.r;

ry=y.r;
psy=y.ps;
cry=y.core;
phi=cell(d+1,1);
phi{d+1}=1;
phi{1}=1;
phx=cell(d+1,1);
phx{d+1}=1; %For storing submatrices in U & V
phx{1}=1;
%Warmup procedure: orthogonalize from right-to-left & maxvol
   cry_old=cry; %This is for checking the accuracy
   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
   pos1=psy(d+1);
   %The first is the the o
   for i=d-1:-1:1
      %Do right-to-left SVD + maxvol (i.e., no fun() is employed, just the
      %current approximation)
      cr=cry(pos1-ry(i+1)*n(i+1)*ry(i+2):pos1-1);
      cr2=cry_old(psy(i):psy(i+1)-1);
      cr2=reshape(cr2,[ry(i)*n(i),ry(i+1)]);
      cr=reshape(cr,[ry(i+1),n(i+1)*ry(i+2)]);
      %cr=cr.';
      %[cr,rm]=qr(cr,0);
      %rm=2*eye(ry(i+1));
      %cr=cr/2;
      [u,s,v]=svd(cr,'econ'); s=diag(s);
      ry(i+1) = my_chop2(s,norm(s)*eps/sqrt(d-1));
      %cr = u * s * v', I think we should leave v orthogonal
      u=u(:,1:ry(i+1)); v=v(:,1:ry(i+1));
      s=s(1:ry(i+1)); u=u*diag(s);
      cr=conj(v); %cr is n(i+1),ry(i+2),ry(i+1) ---  No. it is conj(v)
      rm=u.'; %This is discussable --- maybe u' (or u.')?
      ry(i+1)=size(cr,2);
      %Maxvol should be computed in a different matrix
      cr0=reshape(cr,[n(i+1),ry(i+2),ry(i+1)]);
      cr0=permute(cr0,[3,1,2]);
      cr0=reshape(cr0,[ry(i+1)*n(i+1),ry(i+2)]);
      cr0=cr0*phx{i+2}.';
      cr0=reshape(cr0,[ry(i+1),n(i+1)*ry(i+2)]);
      cr0=cr0.';
      ind=maxvol2(cr0);
      r1=cr0(ind,:);
      phx{i+1}=r1;

      cr=cr.';

      cry(pos1-ry(i+1)*n(i+1)*ry(i+2):pos1-1)=cr(:);
      pos1=pos1-ry(i+1)*n(i+1)*ry(i+2);
      cr2=cr2*(rm).';
      cry(pos1-ry(i)*n(i)*ry(i+1):pos1-1)=cr2(:);
      %Take phi matrix; convolve from right with current cors of V
      cr0=core0(ps(i+1):ps(i+2)-1);
      cr0=reshape(cr0,[r(i+1)*n(i+1),r(i+2)]);
      cr0=cr0*phi{i+2}; %cr0 is now r(i)*n(i)*ry(i+1);
      cr0=reshape(cr0,[r(i+1),n(i+1)*ry(i+2)]);
      phi{i+1}=cr0(:,ind);
   end
   pos1=pos1-ry(1)*n(1)*ry(2);
   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
   cry=cry(pos1:numel(cry));
   y.core=cry;
   y.r=ry;
   y.ps=psy;

   swp=1;
yold=[];
not_converged = true;
while ( swp < nswp && not_converged )
    max_er=0;
cry_old=cry; %This is for checking the accuracy
psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
pos1=1;
 for i=1:d-1
     %fprintf('i=%d \n',i);
     %We care for two cores, with number i & number i+1, and use
     %psi(i) and psi(i+2) as a basis; also we will need to recompute
     %psi(i+1)
     ps1=phi{i}; ps2=phi{i+2}; px1=phx{i}; px2=phx{i+2};
     %Compute (!) superblock and function (!) of it
     cr1=core0(ps(i):ps(i+1)-1);
     cr2=core0(ps(i+1):ps(i+2)-1);
     cr1=reshape(cr1,[r(i),n(i)*r(i+1)]);
     cr1=ps1*cr1;
     cr2=reshape(cr2,[r(i+1)*n(i+1),r(i+2)]);
     cr2=cr2*ps2;
     cr1=reshape(cr1,[ry(i)*n(i),r(i+1)]);
     cr2=reshape(cr2,[r(i+1),n(i+1)*ry(i+2)]);
     cr=cr1*cr2;
     %cr=reshape(cr,[ry(i)*n(i)*n(i+1),ry(i+2)]);
     cr=fun(cr); %Elements are evaluated here!
     cr=reshape(cr,[ry(i)*n(i)*n(i+1),ry(i+2)]);
     cr = cr/(px2.');
     cr=reshape(cr,[ry(i),n(i)*n(i+1)*ry(i+2)]);
     cr=px1 \cr;
     cr=reshape(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
     %Check for local approximation of cr for the error
     cry1=cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1);
     %cry2=cry(pos1+ry(i)*n(i)*ry(i+1):pos1+ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2)-1);
     %if ( swp == 2 )
     % if ( i == 5 )
     %   keyboard;
     % end
     %end
     cry2=cry_old(psy(i+1):psy(i+2)-1);

     cry1=reshape(cry1,[ry(i)*n(i),ry(i+1)]);
     cry2=reshape(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
     appr=cry1*cry2;
     er=norm(appr-cr,'fro')/norm(cr,'fro');
     %keyboard;
     max_er=max(er,max_er);

     %Compute SVD of cr

     [u,s,v]=svd(cr,'econ');
     s=diag(s);
     r2=my_chop2(s,eps*norm(s)/sqrt(d-1));
     s=s(1:r2); u=u(:,1:r2); v=conj(v(:,1:r2));
     v=v*diag(s);

     %Kick rank of u
     ur=randn(size(u,1),kick_rank);
     %Orthogonalize ur to u by Golub-Kahan reorth <- it sucks!!
     [u,rv]=qr([u,ur], 0);
%      u=reort(u,ur);
%      radd=size(u,2)-r2;
%      if ( radd > 0 )
       vr=zeros(size(v,1),kick_rank);
       v=[v,vr];
       v = v*(rv.');
%      end
     r2=size(u,2);

     %ind=maxvol2(u);
     %r1=u(ind,:);
     %u=u/r1; v=v*r1'; v=v.';
     v=v.';
     ry(i+1)=r2;
     cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);
     pos1=pos1+ry(i)*n(i)*ry(i+1);
     cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:);

     %Compute new maxvol and (new) submatrices
     u=reshape(u,[ry(i),n(i)*ry(i+1)]);
     u=px1*u;
     u=reshape(u,[ry(i)*n(i),ry(i+1)]);
     ind=maxvol2(u);
     phx{i+1}=u(ind,:);
     %Now we have to: cr1 with phi from the left
     phi{i+1}=cr1(ind,:); %phi{i+1}=phi{i+1};

 end
  %Truncate local memory
  cry=cry(1:pos1+ry(d)*n(d)*ry(d+1)-1);
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
y.core=cry;
y.r=ry;
y.ps=psy;
cry_old=cry;
%m=numel(cry);
%cry=[zeros(size(cry)),cry];
%pos1=pos1+m;
%cry2=cry_old(psy(i+1):psy(i+2)-1);
cry=cry(psy(d):psy(d+1)-1); %Start--only two cores
 for i=d-1:-1:1
     %fprintf('i=%d \n',i);
     %We care for two cores, with number i & number i+1, and use
     %psi(i) and psi(i+2) as a basis; also we will need to recompute
     %psi(i+1)
     ps1=phi{i}; ps2=phi{i+2}; px1=phx{i}; px2=phx{i+2};
     %Take current core; convolve it with
     %core=core(ps
     %Compute (!) superblock and function (!) of it
     cr1=core0(ps(i):ps(i+1)-1);
     cr2=core0(ps(i+1):ps(i+2)-1);
     cr1=reshape(cr1,[r(i),n(i)*r(i+1)]);
     cr1=ps1*cr1;
     cr2=reshape(cr2,[r(i+1)*n(i+1),r(i+2)]);
     cr2=cr2*ps2;
     cr1=reshape(cr1,[ry(i)*n(i),r(i+1)]);
     cr2=reshape(cr2,[r(i+1),n(i+1)*ry(i+2)]);
     cr=cr1*cr2;
     cr=reshape(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
     cr=fun(cr); %Elements are evaluated here!
     cr=reshape(cr,[ry(i)*n(i)*n(i+1),ry(i+2)]);
     cr = cr/(px2.');
     cr=reshape(cr,[ry(i),n(i)*n(i+1)*ry(i+2)]);
     cr=px1 \cr;
     cr=reshape(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
     %Check for local approximation of cr for the error
     %cry1=cry(pos1-ry(i)*n(i)*ry(i+1):pos1-1);
     cry1=cry_old(psy(i):psy(i+1)-1);
     %cry2=cry(ry(:ry(i+1)*n(i+1)*ry(i+2));
     %cry1=cry(1:ry(i)*n(i)*ry(i+1));
     %cry2=cry(ry(i)*n(i)*ry(i+1)+1:ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2));
     cry2=cry(1:ry(i+1)*n(i+1)*ry(i+2));
     cry(1:ry(i+1)*n(i+1)*ry(i+2))=[];
     %cry2=cry_old(psy(i+1):psy(i+2)-1);
     %cry(1:(ry(i+1)*n(i+1)*ry(i+2)))=[]; %Delete both of first cores
     %cry(1:ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2))=[];
     cry1=reshape(cry1,[ry(i)*n(i),ry(i+1)]);
     cry2=reshape(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
     appr=cry1*cry2;
     er=norm(appr-cr,'fro')/norm(cr,'fro');
     %er
          max_er=max(er,max_er);

     %Compute SVD of cr
     [u,s,v]=svd(cr,'econ');
     s=diag(s);
     r2=my_chop2(s,eps*norm(s)/sqrt(d-1));
     s=s(1:r2); u=u(:,1:r2); v=v(:,1:r2);

     v = conj(v); % <- bug was here

     %Make it standard
     u=u*diag(s);

     %Kick rank

      vr=randn(size(v,1),kick_rank);
      [v,rv]=qr([v,vr], 0);
%       v=reort(v,vr);
%       radd=size(v,2)-r2;
%       if ( radd > 0 )
         ur=zeros(size(u,1),kick_rank);
         u=[u,ur];
         u = u*(rv.');
%       end
      r2=size(v,2);

     v=v.';v0=v;
     %Compute new phi;
     ry(i+1)=r2;
     u=u(:); u=u.'; v=v(:); v=v.';
     cry=[u,v,cry];
     %keyboard;
     %We need new memory for
     %cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:); %Here memory has to be (?) enlarged
     %if ( pos1 <= ry(i)*n(i)*ry(i+1) ) %Have to enlarge memory
     %  cry=cry(pos1:numel(cry));
     %  cry=[u(:)',cry];
     %  pos1=ry(i)*n(i)*ry(i+1)+1;
     %else
     %   pos1=pos1-ry(i)*n(i)*ry(i+1);
      %  cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);
     %end
     %Compute new maxvol and (new) submatrices
     v0=reshape(v0,[ry(i+1)*n(i+1),ry(i+2)]);
     v0=v0*px2.';
     v0=reshape(v0,[ry(i+1),n(i+1)*ry(i+2)]);
     v0=v0.';
     ind=maxvol2(v0);
     phx{i+1}=v0(ind,:);
     %Now we have to: cr1 with phi from the left
     phi{i+1}=cr2(:,ind);

 end

  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);

y.core=cry;
y.r=ry;
y.ps=psy;
if ( isempty(yold) )
 yold=y;
 er_nrm=1;
else

   er_nrm=norm(yold-y)/norm(y);
   yold=y;
end
 fprintf('sweep=%d, er=%3.2e er_nrm=%3.2e \n',swp,max_er,er_nrm);

swp=swp+1;
end
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);

y.core=cry;
y.r=ry;
y.ps=psy;

end

author: oseledets

repo: TT-Toolbox