text:
function [costs,paths] = dijkstra(AorV,xyCorE,SID,FID,iswaitbar)
%DIJKSTRA Calculate Minimum Costs and Paths using Dijkstra's Algorithm
%
%   Inputs:
%     [AorV] Either A or V where
%         A   is a NxN adjacency matrix, where A(I,J) is nonzero (=1)
%               if and only if an edge connects point I to point J
%               NOTE: Works for both symmetric and asymmetric A
%         V   is a Nx2 (or Nx3) matrix of x,y,(z) coordinates
%     [xyCorE] Either xy or C or E (or E3) where
%         xy  is a Nx2 (or Nx3) matrix of x,y,(z) coordinates (equivalent to V)
%               NOTE: only valid with A as the first input
%         C   is a NxN cost (perhaps distance) matrix, where C(I,J) contains
%               the value of the cost to move from point I to point J
%               NOTE: only valid with A as the first input
%         E   is a Px2 matrix containing a list of edge connections
%               NOTE: only valid with V as the first input
%         E3  is a Px3 matrix containing a list of edge connections in the
%               first two columns and edge weights in the third column
%               NOTE: only valid with V as the first input
%     [SID] (optional) 1xL vector of starting points
%         if unspecified, the algorithm will calculate the minimal path from
%         all N points to the finish point(s) (automatically sets SID = 1:N)
%     [FID] (optional) 1xM vector of finish points
%         if unspecified, the algorithm will calculate the minimal path from
%         the starting point(s) to all N points (automatically sets FID = 1:N)
%     [iswaitbar] (optional) a scalar logical that initializes a waitbar if nonzero
%
%   Outputs:
%     [costs] is an LxM matrix of minimum cost values for the minimal paths
%     [paths] is an LxM cell array containing the shortest path arrays
%
%   Revision Notes:
%     (4/29/09) Previously, this code ignored edges that have a cost of zero,
%     potentially producing an incorrect result when such a condition exists.
%     I have solved this issue by using NaNs in the table rather than a
%     sparse matrix of zeros. However, storing all of the NaNs requires more
%     memory than a sparse matrix. This may be an issue for massive data
%     sets, but only if there are one or more 0-cost edges, because a sparse
%     matrix is still used if all of the costs are positive.
%
%   Note:
%     If the inputs are [A,xy] or [V,E], the cost is assumed to be (and is
%       calculated as) the point-to-point Euclidean distance
%     If the inputs are [A,C] or [V,E3], the cost is obtained from either
%       the C matrix or from the edge weights in the 3rd column of E3
%
%   Example:
%       % Calculate the (all pairs) shortest distances and paths using [A,xy] inputs
%       n = 7; A = zeros(n); xy = 10*rand(n,2)
%       tri = delaunay(xy(:,1),xy(:,2));
%       I = tri(:); J = tri(:,[2 3 1]); J = J(:);
%       IJ = I + n*(J-1); A(IJ) = 1
%       [costs,paths] = dijkstra(A,xy)
%
%   Example:
%       % Calculate the (all pairs) shortest distances and paths using [A,C] inputs
%       n = 7; A = zeros(n); xy = 10*rand(n,2)
%       tri = delaunay(xy(:,1),xy(:,2));
%       I = tri(:); J = tri(:,[2 3 1]); J = J(:);
%       IJ = I + n*(J-1); A(IJ) = 1
%       a = (1:n); b = a(ones(n,1),:);
%       C = round(reshape(sqrt(sum((xy(b,:) - xy(b',:)).^2,2)),n,n))
%       [costs,paths] = dijkstra(A,C)
%
%   Example:
%       % Calculate the (all pairs) shortest distances and paths using [V,E] inputs
%       n = 7; V = 10*rand(n,2)
%       I = delaunay(V(:,1),V(:,2));
%       J = I(:,[2 3 1]); E = [I(:) J(:)]
%       [costs,paths] = dijkstra(V,E)
%
%   Example:
%       % Calculate the (all pairs) shortest distances and paths using [V,E3] inputs
%       n = 7; V = 10*rand(n,2)
%       I = delaunay(V(:,1),V(:,2));
%       J = I(:,[2 3 1]);
%       D = sqrt(sum((V(I(:),:) - V(J(:),:)).^2,2));
%       E3 = [I(:) J(:) D]
%       [costs,paths] = dijkstra(V,E3)
%
%   Example:
%       % Calculate the shortest distances and paths from the 3rd point to all the rest
%       n = 7; V = 10*rand(n,2)
%       I = delaunay(V(:,1),V(:,2));
%       J = I(:,[2 3 1]); E = [I(:) J(:)]
%       [costs,paths] = dijkstra(V,E,3)
%
%   Example:
%       % Calculate the shortest distances and paths from all points to the 2nd
%       n = 7; A = zeros(n); xy = 10*rand(n,2)
%       tri = delaunay(xy(:,1),xy(:,2));
%       I = tri(:); J = tri(:,[2 3 1]); J = J(:);
%       IJ = I + n*(J-1); A(IJ) = 1
%       [costs,paths] = dijkstra(A,xy,1:n,2)
%
%   Example:
%       % Calculate the shortest distance and path from points [1 3 4] to [2 3 5 7]
%       n = 7; V = 10*rand(n,2)
%       I = delaunay(V(:,1),V(:,2));
%       J = I(:,[2 3 1]); E = [I(:) J(:)]
%       [costs,paths] = dijkstra(V,E,[1 3 4],[2 3 5 7])
%
%   Example:
%       % Calculate the shortest distance and path between two points
%       n = 1000; A = zeros(n); xy = 10*rand(n,2);
%       tri = delaunay(xy(:,1),xy(:,2));
%       I = tri(:); J = tri(:,[2 3 1]); J = J(:);
%       D = sqrt(sum((xy(I,:)-xy(J,:)).^2,2));
%       I(D > 0.75,:) = []; J(D > 0.75,:) = [];
%       IJ = I + n*(J-1); A(IJ) = 1;
%       [cost,path] = dijkstra(A,xy,1,n)
%       gplot(A,xy,'k.:'); hold on;
%       plot(xy(path,1),xy(path,2),'ro-','LineWidth',2); hold off
%       title(sprintf('Distance from 1 to 1000 = %1.3f',cost))
%
% Web Resources:
%   Dijkstra's Algorithm
%   Graphs
%   Adjacency Matrix
%
% See also: gplot, gplotd, gplotdc, distmat, ve2axy, axy2ve
%
% Author: Joseph Kirk
% Email: jdkirk630@gmail.com
% Release: 1.1
% Date: 4/29/09

% Process Inputs
error(nargchk(2,5,nargin));
all_positive = 1;
[n,nc] = size(AorV);
[m,mc] = size(xyCorE);
[E,cost] = processInputs(AorV,xyCorE);
if nargin < 5
    iswaitbar = 0;
end
if nargin < 4
    FID = (1:n);
end
if nargin < 3
    SID = (1:n);
end
if max(SID) > n || min(SID) < 1
    eval(['help ' mfilename]);
    error('Invalid [SID] input. See help notes above.');
end
if max(FID) > n || min(FID) < 1
    eval(['help ' mfilename]);
    error('Invalid [FID] input. See help notes above.');
end

isreversed = 0;
if length(FID) < length(SID)
    E = E(:,[2 1]);
    cost = cost';
    tmp = SID;
    SID = FID;
    FID = tmp;
    isreversed = 1;
end

L = length(SID);
M = length(FID);
costs = zeros(L,M);
paths = num2cell(nan(L,M));

% Find the Minimum Costs and Paths using Dijkstra's Algorithm
if iswaitbar, wbh = waitbar(0,'Please Wait ... '); end
for k = 1:L
    % Initializations
    if all_positive, TBL = sparse(1,n); else TBL = NaN(1,n); end
    min_cost = Inf(1,n);
    settled = zeros(1,n);
    path = num2cell(nan(1,n));
    I = SID(k);
    min_cost(I) = 0;
    TBL(I) = 0;
    settled(I) = 1;
    path(I) = {I};

    while any(~settled(FID))
        % Update the Table
        TAB = TBL;
        if all_positive, TBL(I) = 0; else TBL(I) = NaN; end
        nids = find(E(:,1) == I);
        % Calculate the Costs to the Neighbor Points and Record Paths
        for kk = 1:length(nids)
            J = E(nids(kk),2);
            if ~settled(J)
                c = cost(I,J);
                if all_positive, empty = ~TAB(J); else empty = isnan(TAB(J)); end
                if empty || (TAB(J) > (TAB(I) + c))
                    TBL(J) = TAB(I) + c;
                    if isreversed
                        path{J} = [J path{I}];
                    else
                        path{J} = [path{I} J];
                    end
                else
                    TBL(J) = TAB(J);
                end
            end
        end

        if all_positive, K = find(TBL); else K = find(~isnan(TBL)); end
        % Find the Minimum Value in the Table
        N = find(TBL(K) == min(TBL(K)));
        if isempty(N)
            break
        else
            % Settle the Minimum Value
            I = K(N(1));
            min_cost(I) = TBL(I);
            settled(I) = 1;
        end
    end
    % Store Costs and Paths
    costs(k,:) = min_cost(FID);
    paths(k,:) = path(FID);
    if iswaitbar, waitbar(k/L,wbh); end
end
if iswaitbar, close(wbh); end

if isreversed
    costs = costs';
    paths = paths';
end

if L == 1 && M == 1
    paths = paths{1};
end

% -------------------------------------------------------------------
    function [E,C] = processInputs(AorV,xyCorE)
        C = sparse(n,n);
        if n == nc
            if m == n
                if m == mc % Inputs: A,cost
                    A = AorV;
                    A = A - diag(diag(A));
                    C = xyCorE;
                    all_positive = all(C(logical(A)) > 0);
                    E = a2e(A);
                else % Inputs: A,xy
                    A = AorV;
                    A = A - diag(diag(A));
                    xy = xyCorE;
                    E = a2e(A);
                    D = ve2d(xy,E);
                    all_positive = all(D > 0);
                    for row = 1:length(D)
                        C(E(row,1),E(row,2)) = D(row);
                    end
                end
            else
                eval(['help ' mfilename]);
                error('Invalid [A,xy] or [A,cost] inputs. See help notes above.');
            end
        else
            if mc == 2 % Inputs: V,E
                V = AorV;
                E = xyCorE;
                D = ve2d(V,E);
                all_positive = all(D > 0);
                for row = 1:m
                    C(E(row,1),E(row,2)) = D(row);
                end
            elseif mc == 3 % Inputs: V,E3
                E3 = xyCorE;
                all_positive = all(E3 > 0);
                E = E3(:,1:2);
                for row = 1:m
                    C(E3(row,1),E3(row,2)) = E3(row,3);
                end
            else
                eval(['help ' mfilename]);
                error('Invalid [V,E] inputs. See help notes above.');
            end
        end
    end

    % Convert Adjacency Matrix to Edge List
    function E = a2e(A)
        [I,J] = find(A);
        E = [I J];
    end

    % Compute Euclidean Distance for Edges
    function D = ve2d(V,E)
        VI = V(E(:,1),:);
        VJ = V(E(:,2),:);
        D = sqrt(sum((VI - VJ).^2,2));
    end
end


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function bb = bb_hull(bb0)

bb=[min(bb0(1,:)); ...
    min(bb0(2,:)); ...
    max(bb0(3,:)); ...
    max(bb0(4,:))];

repo: OpenTLD


text:
function b = r8vec_bracket5 ( nd, xd, xi )

%*****************************************************************************80
%
%% R8VEC_BRACKET5 brackets data between successive entries of a sorted R8VEC.
%
%  Discussion:
%
%    We assume XD is sorted.
%
%    If XI is contained in the interval [XD(1),XD(N)], then the returned
%    value B indicates that XI is contained in [ XD(B), XD(B+1) ].
%
%    If XI is not contained in the interval [XD(1),XD(N)], then B = -1.
%
%    This code implements a version of binary search which is perhaps more
%    understandable than the usual ones.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ND, the number of data values.
%
%    Input, real XD(N), the sorted data.
%
%    Input, real XD, the query value.
%
%    Output, integer B, the bracket information.
%
  if ( xi < xd(1) || xd(nd) < xi )

    b = -1;

  else

    l = 1;
    r = nd;

    while ( l + 1 < r )
      m = floor ( ( l + r ) / 2 );
      if ( xi < xd(m) )
        r = m;
      else
        l = m;
      end
    end

    b = l;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Demo of Hessenberg QR
clc;
close all;
clearvars;
n = 8;
A = gallery('binomial',n);
H = hess(A);
[Q1, R1] = spx.la.hessenberg.qr(H);
H1 = Q1 * R1;
H
H1
fprintf('::: %f\n', max(max(abs(H-H1))));

[Q2, R2] = qr(H);
R1
R2
R3 = abs(R1) - abs(R2)

Q1
Q2
Q3 = abs(abs(Q1) - abs(Q2))

HDIFF  = Q1 * R1  - Q2 * R2

author: indigits

repo: sparse-plex


text:
% Script demonstrating usage of the olcdl_sgd_msk function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>
%         Jialin Liu <danny19921123@gmail.com>
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% NB: This script uses only 5 training images to avoid running for
% longer than is appropriate for a demo script. Note, however, that
% the advantages of the online learning algorithms relative to the
% batch algorithms are not apparent for such a small training set.

% Training images
S0 = zeros(512, 512, 5, 'single');
S0(:,:,1) = single(stdimage('lena.grey')) / 255;
S0(:,:,2) = single(stdimage('barbara.grey')) / 255;
S0(:,:,3) = single(stdimage('kiel.grey')) / 255;
S0(:,:,4) = single(rgb2gray(stdimage('mandrill'))) / 255;
tmp = single(stdimage('man.grey')) / 255;
S0(:,:,5) = tmp(101:612, 101:612);

%Reduce images size to speed up demo script
tmp = zeros(256, 256, 5, 'single');
for k = 1:size(S0,3),
  tmp(:,:,k) = imresize(S0(:,:,k), 0.5);
end
S0 = tmp;

% Apply 10% salt & pepper noise to images
[Sn, wn] = spnoise(S0, 0.1);

% Filter input images and compute highpass images
Sl = zeros(size(Sn), 'single');
Sh = zeros(size(Sn), 'single');
opt_filt = [];
opt_filt.Verbose = 0;
opt_filt.MaxMainIter = 500;
opt_filt.RelStopTol = 1e-3;
for index = 1:size(Sn,3),
    sn = Sn(:,:,index);
    opt_filt.DatFidWeight = ~wn(:,:,index);
    sl = l2tvdenoise(sn, 1, opt_filt);
    sh = sn - sl;
    Sl(:,:,index) = sl;
    Sh(:,:,index) = sh;
end

% Construct initial dictionary
D0 = zeros(8,8,32, 'single');
D0(3:6,3:6,:) = single(randn(4,4,32));

% Set up olcdl parameters
lambda = 0.2;
opt = [];
opt.MaxMainIter = 5*size(S0,3);

% Do dictionary learning
W = ~wn;
[D, optinf] = olcdl_sgd_msk(D0, Sh, lambda, W, opt);

% Display learned dictionary
figure;
imdisp(tiledict(D));

% Plot functional value evolution
figure;
plot(optinf.itstat(:,2));
xlabel('Iterations');
ylabel('Functional value');

author: thfylsty


text:
function legcoeffs = legvals2legcoeffs(legvals)
%LEGVALS2LEGCOEFFS  Convert Legendre values to Legendre coefficients.
% 	LEGCOEFFS = LEGVALS2LEGCOEFFS(LEGVALS) converts the column vector
%   LEGVALS representing values on a Legendre grid (i.e, F(LEGPTS)) to a
%   vector LEGCOEFFS representing the coefficients of the series
%       F(X) = C_LEG(1)*P0(X) + ... + C_LEG(N)*P{N-1}(X).
%
% See also CHEBFUN.IDLT, LEGPTS.

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

% This command is a wrapper for chebfun/idlt.
legcoeffs = chebfun.idlt(legvals);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [f,relres,iter]=frsyniter(F,c,varargin)
%FRSYNITER  Iterative synthesis
%   Usage:  f=frsyniter(F,c);
%           f=frsyniter(F,c,Ls);
%           [f,relres,iter]=frsyniter(F,c,...);
%
%   Input parameters:
%         F       : Frame
%         c       : Array of coefficients.
%         Ls      : length of signal.
%   Output parameters:
%         f       : Signal.
%         relres  : Vector of residuals.
%         iter    : Number of iterations done.
%
%   `f=frsyniter(F,c)` iteratively inverts the analysis operator of *F*, so
%   `frsyniter` always performs the inverse operation of |frana|, even
%   when a perfect reconstruction is not possible by using |frsyn|.
%
%   `[f,relres,iter]=frsyniter(...)` additionally returns the relative
%   residuals in a vector *relres* and the number of iteration steps *iter*.
%
%   **Note:** If it is possible to explicitly calculate the canonical dual
%   frame then this is usually a much faster method than invoking
%   `frsyniter`.
%
%   `frsyniter` takes the following parameters at the end of the line of
%   input arguments:
%
%     'tol',t      Stop if relative residual error is less than the
%                  specified tolerance. Default is 1e-9 (1e-5 for single precision)
%
%     'maxit',n    Do at most n iterations.
%
%     'cg'         Solve the problem using the Conjugate Gradient
%                  algorithm. This is the default.
%
%     'pcg'        Solve the problem using the Preconditioned Conjugate Gradient
%                  algorithm. Please note that preconditioning is not supported
%                  for all frame types.
%
%     'print'      Display the progress.
%
%     'quiet'      Don't print anything, this is the default.
%
%   Algorithms
%   ----------
%
%   The function uses the (Preconditioned) Conjugate Gradient algorithm
%   to solve the following problem::
%
%   ..   FF*f=Fc
%
%   .. math:: FF* f = Fc
%
%   The preconditioning alters the equations such that
%
%   ..   inv(M)FF*f=inv(M)Fc
%
%   .. math:: M^{-1}FF* f = M^{-1}Fc
%
%   Examples
%   --------
%
%   The following example shows how to rectruct a signal without ever
%   using the dual frame:::
%
%      F=frame('dgtreal','gauss',10,20);
%      c=frana(F,bat);
%      [r,relres]=frsyniter(F,c,'tol',1e-14);
%      norm(bat-r)/norm(bat)
%      semilogy(relres);
%      title('Conversion rate of the CG algorithm');
%      xlabel('No. of iterations');
%      ylabel('Relative residual');
%
%   See also: frame, frana, frsyn, franaiter

% AUTHORS: Nicki Holighaus & Peter L. Søndergaard

complainif_notenoughargs(nargin,2,'FRSYNITER');
complainif_notvalidframeobj(F,'FRSYNITER');

tolchooser.double=1e-9;
tolchooser.single=1e-5;

definput.keyvals.Ls=[];
definput.keyvals.tol=tolchooser.(class(c));
definput.keyvals.maxit=100;
definput.keyvals.Fd = [];
definput.flags.alg={'cg','pcg'};
definput.keyvals.printstep=10;
definput.flags.print={'quiet','print'};

[flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);

% if flags.do_auto
%     varargin2 = varargin;
%     varargin2(strcmpi(varargin2,'auto')) = [];
%
%     try
%         varargin2{end+1} = 'pcg';
%         [f,relres,iter]=frsyniter(F,c,varargin2{:});
%     catch
%         if ~flags.do_quiet
%             warning(sprintf('%s: Falling back to regular CG.',upper(mfilename)));
%         end
%         varargin2{end+1} = 'cg';
%         [f,relres,iter]=frsyniter(F,c,varargin2{:});
%     end
%     return;
% end

L=framelengthcoef(F,size(c,1));

Fd = kv.Fd;
% Compute the preconditioner
if flags.do_pcg && isempty(Fd)
    try
       d = cast(1./framediag(F,L),class(c));
    catch
       switch F.type
            case {'filterbank','ufilterbank'}
                Fd = frame(F.type,{'dual',F.g,'forcepainless'},F.a,numel(F.g));
            case {'filterbankreal','ufilterbankreal'}
                Fd = frame(F.type,{'realdual',F.g,'forcepainless'},F.a,numel(F.g));
            otherwise
                error('%s: No preconditioning method available for given frame type.',...
                upper(mfilename));
       end
    end
end

F=frameaccel(F,L);

A=@(x) F.frsyn(F.frana(x));

% It is possible to specify the initial guess, but this is not
% currently done

if flags.do_pcg && isempty(Fd)

      [f,flag,~,iter,relres]=pcg(A,F.frsyn(c),kv.tol,kv.maxit,@(x)d.*x);
elseif flags.do_pcg

      Fd = frameaccel(Fd,L);
      A=@(x) Fd.frsyn(F.frana(x));
      [f,flag,~,iter,relres]=pcg(A,Fd.frsyn(c),kv.tol,kv.maxit);
else

      [f,flag,~,iter,relres]=pcg(A,F.frsyn(c),kv.tol,kv.maxit);
end

if nargout>1
      relres=relres/norm(c(:));
end

% Cut or extend f to the correct length, if desired.
if ~isempty(Ls)
    f=postpad(f,Ls);
else
    Ls=L;
end

if 0
      % This code has been disabled, as the PCG algorithm is so much faster.
    if flags.do_unlocbox

          % Get the upper frame bound (Or an estimation bigger than the bound)
          [~,B]=framebounds(F,L,'a');

          % Set the parameter for the fast projection on a B2 ball
          param.At=@(x) frsyn(F,x);     % adjoint operator
          param.A=@(x)  frana(F,x);     % direct operator
          param.y=c;                    % coefficient
          param.tight=0;                % It's not a tight frame
          param.max_iter=kv.maxit;
          param.tol=kv.tol;
          param.nu=B;

          % Display parameter 0 nothing, 1 summary at convergence, 2 all
          % steps
      if flags.do_print
          param.verbose=1;
      else
          param.verbose=0;
      end

      % Make the projection. Requires UNLocBOX
      [f, ~] = fast_proj_B2(zeros(L,1), 0, param);

      % compute the residue
      res = param.A(f) - param.y; norm_res = norm(res(:), 2);
      relres=norm_res/norm(c(:), 2);

      iter=0; % The code of the fast_proj_B2 is not yet compatible with this
    end

end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% about
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           setup.m
% This file contains all the simulation parameters that can be tuned by the
% user for different test cases. The code in any of the other files
% shouldn't be varied for trying differennt test cases.
%
% Author : Ajinkya Khade, askhade@ncsu.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~license('test','optimization_toolbox')
    error('Optimization Toolbox License is needed for running this program.')
end
%% Time parameters

tau  = 0.5;                 % in s, time-lag of lower-level controller
ts   = 0.1;                 % in s, sampling time of system
tsim = 10;                  % in s, simulation time

%% Constraints

% Input
g = 9.81;

umin = -0.5*g;              % in m/s^2, maximum deceleration
umax = 0.25*g;              % in m/s^2, maximum acceleration

% Constraint Modes
% 0 = no constraints; 1 = all;
ineqMode = 1;               % 2 = input constraint only
eqMode   = 1;

chkInputBoundFlag = true;	% input saturation
%% MPC Parameters

% Horizons
Np = 100;                    % Prediction Horizon
Nc = 100;                    % Control Horizon

% Weighting Matrices
Q = eye(3);                  % state weighting matrix
R = eye(1);                  % input weighting matrix
S = eye(3);                  % terminal state weights

%% Initial Conditions
% Legend -
% e0 - initial error vector
%
% sivdMode -
%       0 - SIVD varies linearly with reference trajectory
%       1 - SIVD varies dynamically with range rate
%
% sivdPlotMode -
%       0 - Position is plotted with respect to origin of refrence frame only
%       1 - Position is plotted with respect to tracking vehicle also

e0 = [-110; 30; 0];
sivdMode = 0;
sivdPlotMode = false;

% For Stalled Target Vehicle
a = 0;         % m/s^2
v0 = 0;        % m/s
vdes = 0;      % m/s
Xref = getXref(a, v0, vdes, ts, tsim);      % generate reference trajectory for tracking vehicle

%% Sweep Parameters & Plot Specific Parameters
%
% Legend -
% sweepParamName - names of parameters over which a sweep is being
%                  performed for the required analysis
% sweepParamVals - sets of values to sweep over for the sweepParams
% sweepParamDisp - names to display for the different cases in the sweep
% legendMode     -
%       1 - concatenates name of sweep parameter to display name
%       2 - uses display name as is
% e0             - initial error vector

% Fig 2
sweepParamName = {'ineqMode','chkInputBoundFlag'};
sweepParamVals = {0,false;...
                  0,true;...
                  2,false};
sweepParamDisp = {'unlimited deceleration', ...
                  'saturated deceleration', ...
                  'constrained deceleration'};
legendMode = 2;
e0 = [-110; 30; 0];

% % Fig 3
% sweepParamName = {'sivdMode'};
% sweepParamVals = {0;1};
% sweepParamDisp = {'constant sivd','dynamic sivd'};
% legendMode = 2;
% sivdPlotMode = true;

% % Fig 4
% sweepParamName = {'tau'};
% sweepParamVals = {0.5;0.4;0.6};
% sweepParamDisp = {'0.5','0.4','0.6'};
% legendMode = 1;

% % Fig 5
% e0 = [-110; 10; 0];
% sweepParamName = {'Nc'};
% sweepParamVals = {25;50;75;100};
% sweepParamDisp = {'25','50','75','100'};
% legendMode = 1;

% % Fig 6
% e0 = [-50; 5; 0];
% sweepParamName = {'Q'};
% sweepParamVals = {diag([1  100 1]); ...
%                   eye(3); ...
%                   diag([100 1  1])};
% sweepParamDisp = {'Q_{Vel} = 100*Q_{Pos}'; ...
%                   'Q_{Vel} = Q_{Pos}'; ...
%                   'Q_{Pos} = 100*Q_{Vel}'};
% legendMode = 2;

% % Fig 7
% a = 2;          %m/s^2
% v0 = 10;        % m/s
% vdes = 29;      % m/s
% Xref = getXref(a, v0, vdes, ts, tsim);
%
% e0 = [-50; 20; 0];
% sweepParamName = {'tau'};
% sweepParamVals = {0.25;0.5;0.75};
% sweepParamDisp = {'0.25';'0.5';'0.75'};
% % sweepParamVals = {0.5};
% % sweepParamDisp = {'0.5'};
% legendMode = 1;

author: ajinkya-khade


text:
clearvars;
close all;
% clc;

rng default;

M = 300;
N = 1000;
K = 60;

A = randn(M, N);
x0 = zeros(N,1);

permutation = randperm(N);
indices = permutation(1:K);
% non-zero entries
x0(indices) = randn(K,1);

b0 = A * x0;
% noise variance
sigma = 0.1;
e = sigma * randn(M, 1);
% add noise
b = b0 + e;

% orthonormalize the rows of A
% perform QR decomposition of rows of A
[Q, R] = qr(A', 0);
% keep the orthogonal rows
A = Q';
% change b accordingly to the new coordinates
b = R'\b;

options.verbose = 0;
options.max_iterations = 200;
options.tolerance = 5e-3;
solver = spx.pursuit.single.L1_ADMM_YZ(A, options);
% an upper bound on residual norm
delta = norm(e);
x = solver.solve_bpic(b, delta);

r = x - x0;
max_diff = max(abs(r));
rel_error = norm(x-x0) /norm(x0);
iterations = solver.details.iterations(1);
elapsed_time = solver.details.elapsed_times(1);
fprintf('Iterations: %d, Relative error: %e, max diff: %.4f, time: %.4f seconds\n', ...
    iterations, rel_error, max_diff, elapsed_time);

print = 1;
if print
    subplot(211);
    stem(x0, '.');
    subplot(212);
    stem(x, '.');
    figure;
    iterations = solver.details.iterations(1);
    iters = 1:iterations;
    primal_objectives = solver.details.primal_objectives(iters, 1);
    dual_objectives = solver.details.dual_objectives(iters, 1);
    plot(iters, primal_objectives);
    hold on;
    plot(iters, dual_objectives);
    xlabel('Iterations');
    ylabel('Objective Value');
    legend({'Primal Objective', 'Dual Objective'});
    grid on;
end


text:
function res = formatMeshOutput(nbArgs, vertices, edges, faces)
%FORMATMESHOUTPUT Format mesh output depending on nargout.
%
%   OUTPUT = formatMeshOutput(NARGOUT, VERTICES, EDGES, FACES)
%   Utilitary function to convert mesh data .
%   If NARGOUT is 0 or 1, return a matlab structure with fields vertices,
%   edges and faces.
%   If NARGOUT is 2, return a cell array with data VERTICES and FACES.
%   If NARGOUT is 3, return a cell array with data VERTICES, EDGES and
%   FACES.
%
%   OUTPUT = formatMeshOutput(NARGOUT, VERTICES, FACES)
%   Same as before, but do not intialize EDGES in output. NARGOUT can not
%   be equal to 3.
%
%   Example
%     % Typical calling sequence (for a very basic mesh of only one face)
%     v = [0 0; 0 1;1 0;1 1];
%     e = [1 2;1 3;2 4;3 4];
%     f = [1 2 3 4];
%
%     varargout = formatMeshOutput(nargout, v, e, f);
%
%   See also
%   meshes3d, parseMeshData
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-12-06,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

if nargin < 4
    faces = edges;
    edges = [];
end

switch nbArgs
    case {0, 1}
        % output is a data structure with fields vertices, edges and faces
        mesh.vertices = vertices;
        if ~isempty(edges)
            mesh.edges = edges;
        end
        mesh.faces = faces;
        res = {mesh};

    case 2
        % keep only vertices and faces
        res = cell(nbArgs, 1);
        res{1} = vertices;
        res{2} = faces;

    case 3
        % return vertices, edges and faces as 3 separate outputs
        res = cell(nbArgs, 1);
        res{1} = vertices;
        res{2} = edges;
        res{3} = faces;

    otherwise
        error('Can not manage more than 3 outputs');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%

for i = 1:length(lambda_vec)
    lambda = lambda_vec(i);

    theta = trainLinearReg(X, y, lambda);

    error_train(i) = linearRegCostFunction(X, y, theta, 0);
    error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end

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

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 6 Assignments/Regularized Linear Regression and Bias,Variance/mlclass-ex5-005/mlclass-ex5/validationCurve.m


text:
function [node,elem,bdFlag,HB] = uniformrefine3(node,elem,bdFlag,HB)
%% UNIFORMREFINE3 uniformly refine a 3-D triangulation.
%
% [node,elem] = uniformrefine3(node,elem) divides each tetrahedra into
% eight small similar sub-tetrahedrons.
%
% [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag) also update boundary
% conditions represented by bdFlag.
%
% [node,elem,~,HB] = uniformrefine3(node,elem) outpus HB array which is
% useful for nodal interpolation. Unlike bisect3, HB is not useful for the
% coarsening. See uniformcoarsen3red.
%
% Warning: uniform refine in 3D is not orientation presereved. Some volumes
% of the sub-tetrahedron are negative. The ordering is special such that
% the refinement can be easy. See uniformrefine3doc for details.
%
% Example
%
%     [node,elem] = cubemesh([-1,1,-1,1,-1,1]);
%     figure(1); subplot(1,3,1);
%     set(gcf,'Units','normal'); set(gcf,'Position',[0.25,0.25,0.5,0.3]);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%     [node,elem] = uniformrefine3(node,elem);
%     figure(1); subplot(1,3,2);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%     bdFlag = setboundary3(node,elem,'Dirichlet');
%     [node,elem,~,bdFlag] = uniformrefine3(node,elem,[],bdFlag);
%     figure(1); subplot(1,3,3);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%
% See also uniformbisect, uniformrefine, bisect, bisect3, uniformcoarsen3red
%
% Reference page in Help browser
%  ifem uniformrefine3doc
%
% Reference: J. Bey. Simplicial grid refinement: on Freudenthal's algorithm
% and the optimal number of congruence classes. Numer. Math.. 85(1):1--29,
% 2000. p11 Algorithm: RedRefinement3D.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('bdFlag','var'), bdFlag =[]; end
if ~exist('HB','var'),    HB = [];     end

%% Construct data structure
[elem2dof,edge] = dof3P2(elem);
N = size(node,1); NT = size(elem,1); NE = size(edge,1);

%% Add new nodes
node(N+1:N+NE,:) = (node(edge(:,1),:)+node(edge(:,2),:))/2;
if ~isempty(HB)
    maxgeneration = max(HB(:,4));
    HB(N+1:N+NE,[1 2 3]) = [(N+1:N+NE)', edge(:,1:2)];
    HB(N+1:N+NE,4) = maxgeneration + 1;
end

%% Refine each tetrahedron into 8 tetrahedrons
t = 1:NT;
p(t,1:10) = elem2dof;
elem(8*NT,:) = [0 0 0 0]; % enlarge the elem array
elem(t,:) = [p(t,1), p(t,5), p(t,6), p(t,7)];
elem(NT+1:2*NT,:) = [p(t,5), p(t,2), p(t,8), p(t,9)];
elem(2*NT+1:3*NT,:) = [p(t,6), p(t,8), p(t,3), p(t,10)];
elem(3*NT+1:4*NT,:) = [p(t,7), p(t,9), p(t,10), p(t,4)];
% always use diagonal 6-9. The ordering is important. See the reference.
elem(4*NT+1:5*NT,:) = [p(t,5), p(t,6), p(t,7), p(t,9)];
elem(5*NT+1:6*NT,:) = [p(t,5), p(t,6), p(t,8), p(t,9)];
elem(6*NT+1:7*NT,:) = [p(t,6), p(t,7), p(t,9), p(t,10)];
elem(7*NT+1:8*NT,:) = [p(t,6), p(t,8), p(t,9), p(t,10)];

%% Update boundary edges
if ~isempty(bdFlag)
    bdFlag(8*NT,:) = [0 0 0 0]; % enlarge the bdFlag array
    bdFlag(NT+1:2*NT,[1 3 4]) = bdFlag(t,[1 3 4]);
    bdFlag(2*NT+1:3*NT,[1 2 4]) = bdFlag(t,[1 2 4]);
    bdFlag(3*NT+1:4*NT,[1 2 3]) = bdFlag(t,[1 2 3]);
    % always use diagonal 6-9
    bdFlag(4*NT+1:5*NT,2) = bdFlag(t,3);
    bdFlag(5*NT+1:6*NT,4) = bdFlag(t,4);
    bdFlag(6*NT+1:7*NT,3) = bdFlag(t,2);
    bdFlag(7*NT+1:8*NT,1) = bdFlag(t,1);
    % change t in the last
    bdFlag(t,1) = 0;
else
    bdFlag = [];
end

repo: ifem


text:
% this is an example code for training a PixelNet model using caffe
% here we consider a 224x224 image
% uniform sampling of pixels in an image, used for segmentation
function trainSeg(gpuid, options)

	method = 'seg';
	cachepath = [options.cachepath, method, '/'];
	if(~isdir(cachepath))
		mkdir(cachepath);
	end

	% --
	caffe.reset_all;
	caffe.set_device(gpuid);
	caffe.set_mode_gpu;

	% train the network --
	trainNet(options, cachepath);
	caffe.reset_all;
end

function trainNet(options,  cachepath)

	% check if it has already been trained
	trainfolder = [cachepath, 'TRAIN/'];
	if(~isdir(trainfolder))
		mkdir(trainfolder);
	end

	% check if the dataset is ready
	% this would be available in the download_data.sh script
	datasetfile = [options.datapath,...
			'/train/seg-voc2012-aug_trainval.mat'];
	if ~(exist(datasetfile, 'file'))
    		error('Dataset is not prepared!');
	end
	load(datasetfile, 'imagelist', 'seglist');

	imagelist = strcat(options.datapath, '/',imagelist);
	seglist = strcat(options.datapath, '/',seglist);

	% load the network --
	solverpath = [options.datapath, '/train/solver.prototxt'];
	initmodelpath = [options.datapath,'/train/VGG16_fconv.caffemodel'];
	trainModel(options,trainfolder,imagelist,seglist,solverpath,initmodelpath);
end

function trainModel(options,trainfolder,imagelist,seglist,solverpath,initmodelpath)

	% Though this is not required --
	imagelist = repmat(imagelist, 10,1);
	seglist = repmat(seglist,10,1);
	rand_ids = randperm(length(imagelist));
	imagelist = imagelist(rand_ids);
	seglist = seglist(rand_ids);

	%
	li = length(imagelist);
	caffe.reset_all;

	% initialize network
	solver = caffe.get_solver(solverpath);
	% load the model
	solver.net.copy_from(initmodelpath);

	% --
	maxsize = options.cnn_input_size + 200;
	input_data = zeros(maxsize, maxsize, 3, options.segimbatch, 'single');
	input_index = zeros(3,options.segbatchsize, 'single');
	input_label = zeros(1,options.segbatchsize, 'single');
	solver.net.blobs('data').reshape([maxsize, maxsize, 3, options.segimbatch]);
	solver.net.blobs('pixels').reshape([3,options.segbatchsize]);
	solver.net.blobs('labels').reshape([1,options.segbatchsize]);

	% set up memory
	solver.net.forward_prefilled();
	% then start training
	oldrng = rng;
	rng(options.seed, 'twister');

	for epoch = 1:options.segepoch
   	 index = randperm(li);
   	 for i = 1:options.segimbatch:li
        	j = i+options.segimbatch-1;
        	if j > li
            	continue;
        	end

	        st = 1;
       		ed = options.segsamplesize;
        	im = 0;
        	for k=i:j
            		ik = index(k);

            		if options.trainFlip
		                flip = rand(1) > 0.5;
            		end

	    		im_data = seg_image_provider(options, imagelist{ik}, flip);
	    		[sampled, label] = seg_label_provider(...
				options, seglist{ik}, options.segsamplesize, flip);

            		% notice the zero-index vs. the one-index
            		sampled(1,:) = im;
            		sampled(2,:) = sampled(2,:) + 100;
            		sampled(3,:) = sampled(3,:) + 100;

            		im = im + 1;
            		input_data(101:options.cnn_input_size+100,....
		      	 101:options.cnn_input_size+100, :, im) = im_data;
            		input_index(:, st:ed) = sampled;
            		input_label(st:ed) = label;
           		st = st + options.segsamplesize;
            		ed = ed + options.segsamplesize;
        	end

	        solver.net.blobs('data').set_data(input_data);
        	solver.net.blobs('pixels').set_data(input_index);
        	solver.net.blobs('labels').set_data(input_label);
        	solver.step(1);
        	% clean up everything
        	input_data(:) = 0;
    	end
    	% add another condition
    	if mod(epoch, options.saveEpoch) == 0 && epoch < options.segepoch
        	epochfile = [trainfolder, sprintf('(%02d).caffemodel',epoch)];
        	solver.net.save(epochfile);
    	end

	% --
      end

	targetfile = [trainfolder,'final_model.caffemodel'];
	solver.net.save(targetfile);
	caffe.reset_all;
	rng(oldrng);

end

author: aayushbansal

repo: PixelNet


text:
function Power = actuatorPower(States, Actuators, Phase)
% function Power = actuatorPower(States, Actuators, Phase)
%
% Computer Generated File -- DO NOT EDIT
%
% This function was created by the function Write_Power()
% 13-Dec-2013 11:41:00
%
%  ARGUMENTS:
%   States = [Nstate x Ntime] matrix of states
%   Actuator = [Nactuator x Ntime] matrix of actuator values
%   Phase = {'D','F','S1',S2'} or omit.
%       Since the ankle torques are acting against the ground, they should
%       be set to zero when a given leg is in flight. If Phase is included,
%       then this function will force these inputs to be zero. If it is
%       ommitted, then it is assumed that the user has done that error
%       checking.
%
%  RETURNS:
%   Power = a struct with fields for each actuator:
%       leg one (linear force actuator)
%       leg two (linear force actuator)
%       ankle one (torque actuator)
%       ankle two (torque actuator)
%       hip (torque from leg one on leg two
%
%
% Matthew Kelly
% Cornell University
%

%
% See also DERIVE_EOM
x0 = States(:,1); % (m) Hip horizontal position
y0 = States(:,2); % (m) Hip vertical position
x1 = States(:,3); % (m) Foot One horizontal position
y1 = States(:,4); % (m) Foot One vertical position
x2 = States(:,5); % (m) Foot Two horizontal position
y2 = States(:,6); % (m) Foot Two vertical position
dx0 = States(:,7); % (m/s) Hip horizontal velocity
dy0 = States(:,8); % (m/s) Hip vertical velocity
dx1 = States(:,9); % (m/s) Foot One horizontal velocity
dy1 = States(:,10); % (m/s) Foot One vertical velocity
dx2 = States(:,11); % (m/s) Foot Two horizontal velocity
dy2 = States(:,12); % (m/s) Foot Two vertical velocity

if nargin==2
    Phase = '';
end

F1 = Actuators(:,1); % (N) Compresive axial force in Leg One
F2 = Actuators(:,2); % (N) Compresive axial force in Leg Two
if strcmp(Phase,'F') || strcmp(Phase, 'S2')
    T1 = zeros(size(Actuators,1),1); %Foot One not in contact with ground!
else
    T1 = Actuators(:,3); % (Nm) External torque applied to Leg One
end
if strcmp(Phase,'F') || strcmp(Phase, 'S1')
    T2 = zeros(size(Actuators,1),1); %Foot Two not in contact with ground!
else
    T2 = Actuators(:,4); % (Nm) External torque applied to Leg Two
end
Thip = Actuators(:,5); % (Nm) Torque acting on Leg Two from Leg One

% Intermediate Calculations
L1 = ((x0 - x1).^2 + (y0 - y1).^2).^(1./2);
L2 = ((x0 - x2).^2 + (y0 - y2).^2).^(1./2);
dL1 = ((dx0 - dx1).*(x0 - x1) + (dy0 - dy1).*(y0 - y1))./L1;
dL2 = ((dx0 - dx2).*(x0 - x2) + (dy0 - dy2).*(y0 - y2))./L2;
% th1 = atan2(x1 - x0, y1 - y0);
% th2 = atan2(x2 - x0, y2 - y0);
dth1 = -((dy0 - dy1).*(x0 - x1) - (dx0 - dx1).*(y0 - y1))./L1.^2;
dth2 = -((dy0 - dy2).*(x0 - x2) - (dx0 - dx2).*(y0 - y2))./L2.^2;

Power.legOne = F1.*dL1;
Power.legTwo = F2.*dL2;
Power.ankleOne = T1.*dth1;
Power.ankleTwo = T2.*dth2;
Power.hip = -Thip.*(dth1 - dth2);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% OP_U_OTIMES_N_V_OTIMES_N: assemble the matrix A = [A(i,j)], A(i,j) = (epsilon (u \otimes n)_j, (v \otimes n)_i ).
%
%   A = op_u_otimes_n_v_otimes_n (spu, spv, msh, coeff);
%
% INPUT:
%
%   spu:   structure representing the space of trial functions (see sp_vector/sp_eval_boundary_side)
%   spv:   structure representing the space of test functions (see sp_vector/sp_eval_boundary_side)
%   msh:   structure containing the domain partition and the quadrature rule for the boundary,
%           since it must contain the normal vector (see msh_cartesian/msh_eval_boundary_side)
%   coeff: vector-valued function f, evaluated at the quadrature points
%
% OUTPUT:
%
%   A: assembled matrix
%
% Copyright (C) 2015, Rafael, 2017 Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function varargout = op_u_otimes_n_v_otimes_n (spu, spv, msh, mshv, coeff)

  shpu = reshape (spu.shape_functions, spu.ncomp, msh.nqn, spu.nsh_max, msh.nel);
  shpv = reshape (spv.shape_functions, spv.ncomp, msh.nqn, spv.nsh_max, msh.nel);

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:, iel)))
      shpu_iel = reshape (shpu(:, :, :, iel), spu.ncomp, msh.nqn, spu.nsh_max);
      shpv_iel = reshape (shpv(:, :, :, iel), spv.ncomp, msh.nqn, spv.nsh_max);

      u_otimes_n_iel = zeros (spu.ncomp, spu.ncomp, msh.nqn, spu.nsh_max);
      v_otimes_n_iel = zeros (spv.ncomp, spv.ncomp, msh.nqn, spv.nsh_max);
% I need the normal from both sides
      normalu_iel = msh.normal (:, :, iel);
      normalv_iel = mshv.normal (:, :, iel);
      for ii = 1:spu.ncomp
        for jj = 1:spu.ncomp
          u_otimes_n_iel(ii,jj,:,:) = bsxfun (@times, shpu_iel(ii,:,:), normalu_iel(jj,:));
          v_otimes_n_iel(ii,jj,:,:) = bsxfun (@times, shpv_iel(ii,:,:), normalv_iel(jj,:));
        end
      end
% Should I permute it, before reshaping?
      u_otimes_n_iel = reshape (u_otimes_n_iel, spu.ncomp*spu.ncomp, msh.nqn, 1, spu.nsh_max);
      v_otimes_n_iel = reshape (v_otimes_n_iel, spv.ncomp*spv.ncomp, msh.nqn, spv.nsh_max, 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1, msh.nqn, 1, 1]);
      v_oxn_times_jw = bsxfun (@times, jacdet_iel, v_otimes_n_iel);
      tmp1 = sum (bsxfun (@times, v_oxn_times_jw, u_otimes_n_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);

    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_u_otimes_n_v_otimes_n: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_u_otimes_n_v_otimes_n: wrong number of output arguments')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [ output_meta ] = meta2sicd_s1noise( domnode, meta_product )
%META2SICD_S1NOISE Converts Sentinel-1 noise XML file into SICD format
%
% Written by: Wade Schwartzkopf, NGA Research
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Setup
xp=javax.xml.xpath.XPathFactory.newInstance.newXPath();

% Compute the first line of each burst in overall TIFF file
first_line = zeros(numel(meta_product),1);
for i = 1:(numel(meta_product)-1)
    first_line(i+1) = first_line(i) + meta_product{i}.ImageData.NumCols;
end
% Read noise parameters from XML
n_str = 'noise';
num_noise_vecs=str2double(xp.evaluate(...
    ['count(noise/' n_str 'VectorList/' n_str 'Vector)'],domnode));
if num_noise_vecs == 0  % Data after March 2018 had different format
    n_str = 'noiseRange';
    num_noise_vecs=str2double(xp.evaluate(...
        ['count(noise/' n_str 'VectorList/' n_str 'Vector)'],domnode));
end
noisevals = cell(num_noise_vecs,1);
pixel = cell(num_noise_vecs,1);
noise_line = zeros(num_noise_vecs,1);
for i = 1:num_noise_vecs
    noise_line(i) = str2double(xp.evaluate(...
        ['noise/' n_str 'VectorList/' n_str 'Vector[' num2str(i) ']/line'],...
        domnode));
    pixel{i} = str2num(xp.evaluate(...
        ['noise/' n_str 'VectorList/' n_str 'Vector[' num2str(i) ']/pixel'],...
        domnode));
    noisevals{i} = str2num(xp.evaluate(...
        ['noise/' n_str 'VectorList/' n_str 'Vector[' num2str(i) ']/' n_str 'Lut'],...
        domnode));
end
mode = char(xp.evaluate('noise/adsHeader/mode',domnode));
fit_tolerance = 0.1; % Simple 1D data for IW should fit quite tightly
noisepoly_order = 6; % For IW, this fits quite tightly
if upper(mode(1))=='S'
    % Average all of the noise values across azimuth for a single 1D
    % polynomial that varies in range.  They are roughly similar.  A
    % precise solution would include a 2D fit.
    pixel = {cell2mat(pixel)};
    noisevals = {cell2mat(noisevals)};
    noise_line = 0;
    % Allow for more slop for SM since we have 2D noise values, but aren't
    % doing a full 2D fit.  Because we average across all azimuth, there is
    % no reason to attempt a tight fit, since it won't match across all
    % lines for any fit.
    fit_tolerance = 0.5;
    noisepoly_order = 3;
% Double check our assumptions about the noise data
elseif any(mod(noise_line(1:(end-1)), ...
        double(meta_product{1}.ImageData.NumCols))~=0)
    warning('META2SICD_S1NOISE:UNEXPECTED_LINE','Expected noise values at beginning of each burst.');
    return;
end
% Compute SICD noise polynomials
output_meta = cell(numel(meta_product),1);
for i = 1:numel(meta_product) % First and last are usually outside bursts
    % Which noise values should we use for this burst?  We use the ones
    % that correspond to the first column of the burst.
    line_ind = find(noise_line>=first_line(i) & ...
        noise_line<(first_line(i)+meta_product{i}.ImageData.NumCols),1);
    if isempty(line_ind), continue, end
    coords_rg_m = (double(0:(meta_product{i}.ImageData.NumRows-1)) - ...
        double(meta_product{i}.ImageData.SCPPixel.Row)) * ...
        meta_product{i}.Grid.Row.SS;
    coords_rg_m = coords_rg_m(pixel{line_ind}(:)+1);
    old_state = warning('off','MATLAB:polyfit:RepeatedPointsOrRescale');
    noisepoly = polyfit(coords_rg_m(:), 10*log10(noisevals{line_ind}(:)), noisepoly_order);
    warning(old_state);
    if any(polyval(noisepoly, coords_rg_m(:)) - ...
            10*log10(noisevals{line_ind}(:)) > fit_tolerance)
        warning('META2SICD_S1NOISE:NOISE_TOLERANCE', ...
            'Noise model does not appear to fit data tightly.');
    end
    % These are typically very close across all bursts, but we compute for
    % each bursts independently anyway.
    output_meta{i}.Radiometric.NoiseLevel.NoiseLevelType = 'ABSOLUTE';
    output_meta{i}.Radiometric.NoiseLevel.NoisePoly = noisepoly(end:-1:1).';
end
% Throw warning for old, potentially inaccurate data
last_az_time = xp.evaluate(...
        ['noise/' n_str 'VectorList/' n_str 'Vector[' num2str(num_noise_vecs) ']/azimuthTime'],...
        domnode);
if datenum(char(last_az_time),'yyyy-mm-ddTHH:MM:SS')<datenum(2015, 11, 25)
    warning('META2SICD_S1NOISE:EARLY_DATE','Sentinel-1 radiometric data prior to November 25, 2015 might not be accurate.');
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function y = log_sum_exp( varargin )

%LOG_SUM_EXP    log(sum(exp(x))).
%   LOG_SUM_EXP(X) = LOG(SUM(EXP(X)).
%
%   When used in a CVX model, LOG_SUM_EXP(X) causes CVX's successive
%   approximation method to be invoked, producing results exact to within
%   the tolerance of the solver. This is in contrast to LOGSUMEXP_SDP,
%   which uses a single SDP-representable global approximation.
%
%   If X is a matrix, LOGSUMEXP_SDP(X) will perform its computations
%   along each column of X. If X is an N-D array, LOGSUMEXP_SDP(X)
%   will perform its computations along the first dimension of size
%   other than 1. LOGSUMEXP_SDP(X,DIM) will perform its computations
%   along dimension DIM.
%
%   Disciplined convex programming information:
%       LOG_SUM_EXP(X) is convex and nondecreasing in X; therefore, X
%       must be convex (or affine).

persistent P
if isempty( P ),
    P.map = cvx_remap( { 'real' ; 'convex' } );
    P.funcs = { @lse_1, @lse_2 };
    P.zero = -Inf;
    P.reduce = true;
    P.constant = 1;
    P.reverse = false;
    P.name = 'log_sum_exp';
    P.dimarg = 2;
end
cvx_expert_check( 'log_sum_exp', varargin{1} );
y = cvx_reduce_op( P, varargin{:} );

function x = lse_1( x )
xmid = 0.5 * ( max( x, [], 1 ) + min( x, [], 1 ) );
x = log( sum( exp( bsxfun( @minus, x, xmid ) ), 1 ) ) + xmid;

function y = lse_2( x ) %#ok
[nx,nv] = size(x); %#ok
cvx_begin
    variable w( nx, nv )
    epigraph variable y( 1, nv )
    exp( x - repmat( y, [ nx, 1 ] ) ) <= w; %#ok
    sum( w, 1 ) == 1; %#ok
cvx_end

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

author: yu-jiang

repo: radpbook

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


text:
%
% Extract variable names from any number of polynomial strings
%
% Syntax:  (pvar is the shortened alias for GetVariables)
%          >> z = GetVariables(f)
%          >> z = GetVariables(f,g,h)
%
%  to extract variables from a cell array f of polynomial strings:
%          >> z = GetVariables(f{:})
%
%  Input:   f --- (string) polynomial
% Output:   z --- (cell)   variable names of f
%
% Example:  >> f = '-2*x + 9*x*y^2 - 9*x^2*y^2 + 8*x^3*y^2';
%           >> g = '5 - x^3*z';
%           >> z = GetVariables(f,g)
%           z =
%               'x'    'y'    'z'

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [map2] = changeNodesArea(map, metList, areaHeight, areaWidth)
% Change the area size of a list of metabolites (standardize the map)
%
% USAGE:
%
%    [map2] = changeNodesArea(map,metList,areaHeight,areaWidth)
%
% INPUTS:
%    map:          A parsed model structure generated by
%                  'transformXML2Map' function
%    metList:      list of metabolites which area wants to be changed
%
% OPTIONAL INPUTS:
%    areaHeight:       change heigt
%    areaWidth:        change areaWidth
%
% OUTPUTS:
%    map2:          New parsed file with the changes in the reactions
%
% .. Authors:
%       - J.modamio 18/07/2017. Belval, Luxembourg, 18/07/2017.
%       - N.Sompairac - Institut Curie, Paris, 11/10/2017. (Code checking)

    if nargin < 3
        areaHeight = 25.0;
        if nargin < 4
            areaWidth = 70.0;
        end
    end

    % Rename map
    map2 = map;

    % Find mets in the map
    id = map2.specID(ismember(map2.specName, metList));
    idx = find(ismember(map2.molID, id));

    % Change areaHeight
    map2.molHeight(idx, 1) = {areaHeight};

    % Change areaWidth
    map2.molWidth(idx, 1) = {areaWidth};

end

author: opencobra

repo: cobratoolbox