text:
function c = times(a,b,dummy)
%TIMES        Implements  a .* b  for intervals
%

% written  10/16/98     S.M. Rump
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    remove check for 'double'
%                                    extra argument for non-interval output (internal use only)
%                                    take care of Matlab sparse Inf/NaN bug
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 02/11/06     S.M. Rump  SparseInfNanFlag removed
%                                    sparse radius for complex data
% modified 05/23/06     S.M. Rump  sparse Inf/NaN bug corrected in Version 7.2+
% modified 12/03/06     S.M. Rump  Sparse Bug global flag (thanks to Arnold)
% modified 12/05/06     S.M. Rump  0*infsup(0,inf) (thanks to Arnold)
% modified 01/19/07     S.M. Rump  zero radius for huge arrays
% modified 05/22/07     S.M. Rump  check for sparse zero factor
% modified 10/21/07     S.M. Rump  Matlab bug for sparse(0)*inf
% modified 10/18/08     S.M. Rump  abss replaced by mag, improved performance
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
  end
  huge = ( max(prod(size(a)),prod(size(b)))>1e9 );

  if ~isa(a,'intval')                     % a is double
    if ~isreal(a) | b.complex             % complex case
      if ~b.complex
        setround(1)
        b.mid = b.inf + 0.5*(b.sup-b.inf);
        b.rad = b.mid - b.inf;
      end
      c.complex = 1;
      c.inf = [];
      c.sup = [];
      if isreal(a) | ~b.complex           % R .* IC  or  C .* IC
        setround(-1)
        c1 = a .* b.mid;
        setround(1)
        c2 = a .* b.mid;
      else                                % C .* IC
        setround(-1)
        c1 = real(a) .* real(b.mid) + (-imag(a)) .* imag(b.mid) + ...
             ( real(a) .* imag(b.mid) + imag(a) .* real(b.mid) ) * j;
        setround(1)
        c2 = real(a) .* real(b.mid) + (-imag(a)) .* imag(b.mid) + ...
             ( real(a) .* imag(b.mid) + imag(a) .* real(b.mid) ) * j;
      end
      c.mid = c1 + 0.5*(c2-c1);           % R .* IC  or  C .* IC
      c.rad = abs(c.mid-c1);
      if ~isequal(b.rad,0)                % make sure c.rad remains sparse
        c.rad = c.rad + abs(a) .* b.rad;
      end
      % too rare, improvement doesn't pay
%       if huge
%         [cradi,cradj] = find(c.rad);
%         if ~any(find(cradi))              % take care of huge arrays
%           c.rad = 0;
%         end
%       else
%         if ~any(find(c.rad))
%           c.rad = 0;
%         end
%       end
    else                                  % real case  R .* IR
      c.complex = 0;
      if ( prod(size(a))==1 ) & isequal(a(1,1),0)    % careful with 0*inf and sparse a
        if isinf(b.inf) | isinf(b.sup)
          c.inf = repmat(NaN,size(b.inf));
        else
          if issparse(b.inf)
            [mb,nb] = size(b.inf);
            c.inf = sparse([],[],[],mb,nb);
          else
            c.inf = zeros(size(b.inf));
          end
        end
        c.sup = c.inf;
      else
        setround(-1)
        c.inf = min( a .* b.inf , a .* b.sup );
        setround(1)
        c.sup = max( a .* b.inf , a .* b.sup );
      end
      c.mid = [];
      c.rad = [];
    end
  elseif ~isa(b,'intval')                 % b is double
    if a.complex | ~isreal(b)             % complex case
      if ~a.complex
        setround(1)
        a.mid = a.inf + 0.5*(a.sup-a.inf);
        a.rad = a.mid - a.inf;
      end
      c.complex = 1;
      c.inf = [];
      c.sup = [];
      if ~a.complex | isreal(b)           % IC .* R  or  IC .* C
        setround(-1)
        c1 = a.mid .* b;
        setround(1)
        c2 = a.mid .* b;
      else                                % IC .* C
        setround(-1)
        c1 = real(a.mid) .* real(b) + (-imag(a.mid)) .* imag(b) + ...
             ( real(a.mid) .* imag(b) + imag(a.mid) .* real(b) ) * j;
        setround(1)
        c2 = real(a.mid) .* real(b) + (-imag(a.mid)) .* imag(b) + ...
             ( real(a.mid) .* imag(b) + imag(a.mid) .* real(b) ) * j;
      end
      c.mid = c1 + 0.5*(c2-c1);           % R .* IC  or  C .* IC
      c.rad = abs(c.mid-c1);
      if ~isequal(a.rad,0)                % make sure c.rad remains sparse
        c.rad = c.rad + a.rad .* abs(b) ;
      end
      % too rare, improvement doesn't pay
%       if huge
%         [cradi,cradj] = find(c.rad);
%         if ~any(find(cradi))              % take care of huge arrays
%           c.rad = 0;
%         end
%       else
%         if ~any(find(c.rad))
%           c.rad = 0;
%         end
%       end
    else                                  % real case  IR .* R
      c.complex = 0;
      setround(-1)
      c.inf = min( a.inf .* b , a.sup .* b );
      setround(1)
      c.sup = max( a.inf .* b , a.sup .* b );
      c.mid = [];
      c.rad = [];
    end
  else                                    % both a and b interval
    if a.complex | b.complex              % complex case
      if ~a.complex
        setround(1)
        a.mid = a.inf + 0.5*(a.sup-a.inf);
        a.rad = a.mid - a.inf;
      end
      if ~b.complex
        setround(1)
        b.mid = b.inf + 0.5*(b.sup-b.inf);
        b.rad = b.mid - b.inf;
      end
      c.complex = 1;
      c.inf = [];
      c.sup = [];
      if ~a.complex | ~b.complex          % one real factor
        setround(-1)
        c1 = a.mid .* b.mid;
        setround(1)
        c2 = a.mid .* b.mid;
      else
        setround(-1)
        c1 = real(a.mid) .* real(b.mid) + (-imag(a.mid)) .* imag(b.mid) + ...
             ( real(a.mid) .* imag(b.mid) + imag(a.mid) .* real(b.mid) ) * j;
        setround(1)
        c2 = real(a.mid) .* real(b.mid) + (-imag(a.mid)) .* imag(b.mid) + ...
             ( real(a.mid) .* imag(b.mid) + imag(a.mid) .* real(b.mid) ) * j;
      end
      c.mid = c1 + 0.5*(c2-c1);           % IR .* IC  or  IC .* IC
      c.rad = abs(c.mid-c1);
      if ~isequal(a.rad,0)                % make sure c.rad remains sparse
        if ~isequal(b.rad,0)
          c.rad = c.rad + a.rad .* ( abs(b.mid) + b.rad );
        else
          c.rad = c.rad + a.rad .* abs(b.mid);
        end
      end
      if ~isequal(b.rad,0)
        c.rad = c.rad + abs(a.mid) .* b.rad;
      end
      % too rare, improvement doesn't pay
%       if huge
%         [cradi,cradj] = find(c.rad);
%         if ~any(find(cradi))              % take care of huge arrays
%           c.rad = 0;
%         end
%       else
%         if ~any(find(c.rad))
%           c.rad = 0;
%         end
%       end
    else                                  % real case  IR .* IR
      c.complex = 0;
      %setround(-1)
      c.inf = min( a.inf .* b.inf , a.inf .* b.sup );
      c.inf = min( c.inf , a.sup .* b.inf );
      c.inf = min( c.inf , a.sup .* b.sup );
      %setround(1)
      c.sup = max( a.inf .* b.inf , a.inf .* b.sup );
      c.sup = max( c.sup , a.sup .* b.inf );
      c.sup = max( c.sup , a.sup .* b.sup );
      c.mid = [];
      c.rad = [];
    end
  end

  INTLAB_SPARSE_BUG = getappdata(0,'INTLAB_SPARSE_BUG');
  if INTLAB_SPARSE_BUG
    % take care of Matlab sparse NaN bug
    if huge                                       % huge linear indices, for simplicity inf*x=NaN
      if c.complex
        [m,n] = size(c.mid);
      else
        [m,n] = size(c.inf);
      end
      if issparse(a)
        index = any( isnan(b) | isinf(b) );
        if any(index(:))                          % keep linear index small
          [indexi,indexj] = find(isnan(b));
          cNaN = sparse(indexi,indexj,NaN,m,n);
          [ia,ja,aa] = find(mag(a));
          [ib,jb,bb] = find(mag(b));
          ab = sparse([ia;ib],[ja;jb],[complex(aa(:),0);complex(0,bb(:))],m,n);
          [indexi,indexj,ab] = find(ab);
          index = isinf(imag(ab));
          ab = real(ab);
          ab(index & (ab==0)) = NaN;
          ab(~isnan(ab)) = 0;
          cNaN = cNaN + sparse(indexi,indexj,ab,m,n);
          if c.complex
            c.mid = c.mid + cNaN;
            if isequal(c.rad,0)
              c.rad = cNaN;                       % scalar + sparse = full also for scalar=0
            else
              c.rad = c.rad + cNaN;
            end
          else
            c.inf = c.inf + cNaN;
            c.sup = c.sup + cNaN;
          end
        end
      end
      if issparse(b)
        index = any( isnan(a) | isinf(a) );
        if any(index(:))                          % keep linear index small
          [indexi,indexj] = find(isnan(a));
          cNaN = sparse(indexi,indexj,NaN,m,n);
          [ia,ja,aa] = find(mag(a));
          [ib,jb,bb] = find(mag(b));
          ab = sparse([ia;ib],[ja;jb],[complex(0,aa(:));complex(bb(:),0)],m,n);
          [indexi,indexj,ab] = find(ab);
          index = isinf(imag(ab));
          ab = real(ab);
          ab(index & (ab==0)) = NaN;
          ab(~isnan(ab)) = 0;
          cNaN = cNaN + sparse(indexi,indexj,ab,m,n);
          if c.complex
            c.mid = c.mid + cNaN;
            if isequal(c.rad,0)
              c.rad = cNaN;                       % scalar + sparse = full also for scalar=0
            else
              c.rad = c.rad + cNaN;
            end
          else
            c.inf = c.inf + cNaN;
            c.sup = c.sup + cNaN;
          end
        end
      end
    else
      Index = logical(0);
      if issparse(a) & ( prod(size(a))~=1 )     % sparse scalar is handled correctly by Matlab
        Index = isnan(b);
        if any(Index(:))
          if prod(size(b))==1
            Index = find( a==0 );					      % may be almost full
          end
        end
        index = find(isinf(b));
        if any(index)
          if prod(size(b))==1                   % a is scalar
            index = find( a==0 );               % may be almost full
          else                                  % a and b not scalar, same size
            if isa(a,'intval')
              if a.complex
                if isequal(a.rad,0)
                  index( a.mid(index)~=0 ) = [];
                else
                  index( ( a.mid(index)~=0 ) | ( a.rad(index)~=0 ) ) = [];
                end
              else
                index( ( a.inf(index)~=0 ) | ( a.sup(index)~=0 ) ) = [];
              end
            else
              index( a(index)~=0 ) = [];
            end
          end
          Index(index) = 1;
        end
      end
      if issparse(b) & ( prod(size(b))~=1 )     % sparse scalar is handled correctly by Matlab
        Index = Index | isnan(a);
        if any(Index(:))
          if prod(size(a))==1
            Index = find( b==0 );					      % may be almost full
          end
        end
        index = find(isinf(a));
        if any(index)
          if prod(size(a))==1                   % a is scalar
            index = find( b==0 );               % may be almost full
          else                                  % a and b not scalar, same size
            if isa(b,'intval')
              if b.complex
                if isequal(b.rad,0)
                  index( b.mid(index)==0 ) = [];
                else
                  index( ( b.mid(index)~=0 ) | ( b.rad(index)~=0 ) ) = [];
                end
              else
                index( ( b.inf(index)~=0 ) | ( b.sup(index)~=0 ) ) = [];
              end
            else
              index( b(index)~=0 ) = [];
            end
          end
          Index(index) = 1;
        end
      end
      if any(Index(:))
        if c.complex
          c.mid(Index) = NaN;
          if isequal(c.rad,0)
            [m,n] = size(c.mid);
            c.rad = sparse([],[],[],m,n);
          end
          c.rad(Index) = NaN;
        else
          c.inf(Index) = NaN;
          c.sup(Index) = NaN;
        end
      end
    end
  end

  % non-interval output for performance improvement for hessians
  if nargin==2
    c = class(c,'intval');
  end

  setround(rndold)

repo: OpenMAS


text:
%%
dt = 0.2;
b = Bicycle('L', 2.6, 'speedmax', 30, 'steermax', 0.6, 'dt', dt, 'x0', [0 0 0])

b.add_driver( RandomPath(10) );

p = b.run(1000);

%%

clf
[car.image,~,car.alpha] = imread('car2.png');
car.rotation = 180;
car.centre = [81,110];
car.centre = [648; 173];
car.length = 4.2;

plotvol(10)

a = gca

x = [0 0 0]';
h=plot_vehicle(p, 'model', car, 'trail', 'r:')
a.XLimMode = 'manual';
a.YLimMode = 'manual';
set(gcf, 'Color', 'w')
grid on
a = gca;
xyzlabel

author: petercorke

repo: robotics-toolbox-matlab


text:
%%*****************************************************************************
%% HSDsqlp: solve an semidefinite-quadratic-linear program
%%   by infeasible path-following method on the homogeneous self-dual model.
%%
%%  [obj,X,y,Z,info,runhist] =
%%      HSDsqlp(blk,At,C,b,OPTIONS,X0,y0,Z0);
%%
%%  Input: blk: a cell array describing the block diagonal structure of SQL data.
%%          At: a cell array with At{p} = [svec(Ap1) ... svec(Apm)]
%%         b,C: data for the SQL instance.
%%     OPTIONS: a structure that specifies parameters required in HSDsqlp.m,
%%              (if it is not given, the default in sqlparameters.m is used).
%%
%%  (X0,y0,Z0): an initial iterate (if it is not given, the default is used).
%%
%%  Output: obj  = [<C,X> <b,y>].
%%          (X,y,Z): an approximately optimal solution or a primal or dual
%%                   infeasibility certificate.
%%          info.termcode = termination-code
%%          info.iter     = number of iterations
%%          info.obj      = [primal-obj, dual-obj]
%%          info.cputime  = total-time
%%          info.gap      = gap
%%          info.pinfeas  = primal_infeas
%%          info.dinfeas  = dual_infeas
%%          runhist.pobj    = history of primal objective value.
%%          runhist.dobj    = history of dual   objective value.
%%          runhist.gap     = history of <X,Z>.
%%          runhist.pinfeas = history of primal infeasibility.
%%          runhist.dinfeas = history of dual   infeasibility.
%%          runhist.cputime = history of cputime spent.
%%----------------------------------------------------------------------------
%%  The OPTIONS structure specifies the required parameters:
%%      vers  gam  predcorr  expon  gaptol  inftol  steptol
%%      maxit  printlevel  ...
%%      (all have default values set in sqlparameters.m).
%%
%%*************************************************************************
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

function [obj,X,y,Z,info,runhist] = ...
    HSDsqlp(blk,At,C,b,OPTIONS,X0,y0,Z0,kap0,tau0,theta0)

if (nargin < 5); OPTIONS = []; end

isemptyAtb = 0;
if isempty(At) && isempty(b);
    %% Add redundant constraint: <-I,X> <= 0
    b = 0;
    At = ops(ops(blk,'identity'),'*',-1);
    numblk = size(blk,1);
    blk{numblk+1,1} = 'l'; blk{numblk+1,2} = 1;
    At{numblk+1,1} = 1; C{numblk+1,1} = 0;
    isemptyAtb = 1;
end
%%
%%-----------------------------------------
%% get parameters from the OPTIONS structure.
%%-----------------------------------------
%%
% warning off;

matlabversion = sscanf(version,'%f');
if strcmp(computer,'PCWIN64') || strcmp(computer,'GLNXA64')
    par.computer = 64;
else
    par.computer = 32;
end
par.matlabversion = matlabversion(1);
par.vers        = 0;
par.predcorr    = 1;
par.gam         = 0;
par.expon       = 1;
par.gaptol      = 1e-8;
par.inftol      = 1e-8;
par.steptol     = 1e-6;
par.maxit       = 100;
par.printlevel  = 3;
par.stoplevel   = 1;
par.scale_data  = 0;
par.spdensity   = 0.4;
par.rmdepconstr = 0;
par.smallblkdim = 50;
par.schurfun     = cell(size(blk,1),1);
par.schurfun_par = cell(size(blk,1),1);
%%
parbarrier  = cell(size(blk,1),1);
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if strcmp(pblk{1},'s') || strcmp(pblk{1},'q')
        parbarrier{p} = zeros(1,length(pblk{2}));
    elseif strcmp(pblk{1},'l') || strcmp(pblk{1},'u' )
        parbarrier{p} = zeros(1,sum(pblk{2}));
    end
end
parbarrier_0 = parbarrier;
%%
if nargin > 4,
    if isfield(OPTIONS,'vers');        par.vers     = OPTIONS.vers; end
    if isfield(OPTIONS,'predcorr');    par.predcorr = OPTIONS.predcorr; end
    if isfield(OPTIONS,'gam');         par.gam      = OPTIONS.gam; end
    if isfield(OPTIONS,'expon');       par.expon    = OPTIONS.expon; end
    if isfield(OPTIONS,'gaptol');      par.gaptol   = OPTIONS.gaptol; end
    if isfield(OPTIONS,'inftol');      par.inftol   = OPTIONS.inftol; end
    if isfield(OPTIONS,'steptol');     par.steptol  = OPTIONS.steptol; end
    if isfield(OPTIONS,'maxit');       par.maxit    = OPTIONS.maxit; end
    if isfield(OPTIONS,'printlevel');  par.printlevel  = OPTIONS.printlevel; end
    if isfield(OPTIONS,'stoplevel');   par.stoplevel   = OPTIONS.stoplevel; end
    if isfield(OPTIONS,'scale_data');  par.scale_data  = OPTIONS.scale_data; end
    if isfield(OPTIONS,'spdensity');   par.spdensity   = OPTIONS.spdensity; end
    if isfield(OPTIONS,'rmdepconstr'); par.rmdepconstr = OPTIONS.rmdepconstr; end
    if isfield(OPTIONS,'smallblkdim'); par.smallblkdim = OPTIONS.smallblkdim; end
    if isfield(OPTIONS,'parbarrier');
        parbarrier = OPTIONS.parbarrier;
        if isempty(parbarrier); parbarrier = parbarrier_0; end
        if ~iscell(parbarrier);
            tmp = parbarrier; clear parbarrier; parbarrier{1} = tmp;
        end
        if (length(parbarrier) < size(blk,1))
            len = length(parbarrier);
            parbarrier(len+1:size(blk,1)) = parbarrier_0(len+1:size(blk,1));
        end
    end
    if isfield(OPTIONS,'schurfun');
        par.schurfun = OPTIONS.schurfun;
        if ~isempty(par.schurfun); par.scale_data = 0; end
    end
    if isfield(OPTIONS,'schurfun_par'); par.schurfun_par = OPTIONS.schurfun_par; end
    if isempty(par.schurfun);     par.schurfun = cell(size(blk,1),1); end
    if isempty(par.schurfun_par); par.schurfun_par = cell(size(blk,1),1); end
end
if (size(blk,2) > 2); par.smallblkdim = 0; end
%%
%%-----------------------------------------
%% convert matrices to cell arrays.
%%-----------------------------------------
%%
if ~iscell(At); At = {At}; end;
if ~iscell(C);  C = {C}; end;
if all(size(At) == [size(blk,1), length(b)]);
    convertyes = zeros(size(blk,1),1);
    for p = 1:size(blk,1)
        if strcmp(blk{p,1},'s') && all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
        end
    end
    if any(convertyes)
        if (par.printlevel);
            fprintf('\n sqlp: converting At into required format');
        end
        At = svec(blk,At,ones(size(blk,1),1));
    end
end
%%
%%-----------------------------------------
%% validate SQLP data.
%%-----------------------------------------
%%
% tstart = cputime;
[blk,At,C,b,blkdim,numblk] = validate(blk,At,C,b,par);
[blk,At,C,b,iscmp] = convertcmpsdp(blk,At,C,b);
if (iscmp) && (par.printlevel>=2)
    fprintf('\n SQLP has complex data');
end
if (nargin <= 5) || (isempty(X0) || isempty(y0) || isempty(Z0));
    if (max([ops(At,'norm'),ops(C,'norm'),norm(b)]) > 1e2)
        [X0,y0,Z0] = infeaspt(blk,At,C,b,1);
    else
        [X0,y0,Z0] = infeaspt(blk,At,C,b,2,1);
    end
end
if ~iscell(X0);  X0 = {X0}; end;
if ~iscell(Z0);  Z0 = {Z0}; end;
if (length(y0) ~= length(b))
    error('HSDsqlp: length of b and y0 not compatible');
end
[X0,Z0] = validate_startpoint(blk,X0,Z0,par.spdensity);
%%
if (nargin <= 8) || (isempty(kap0) || isempty(tau0) || isempty(theta0))
    if (max([ops(At,'norm'),ops(C,'norm'),norm(b)]) > 1e6)
        kap0 = 10*blktrace(blk,X0,Z0);
    else
        kap0 = blktrace(blk,X0,Z0);
    end
    tau0 = 1; theta0 = 1;
end
%%
if (par.printlevel>=2)
    fprintf('\n num. of constraints = %2.0d',length(b));
    if blkdim(1);
        fprintf('\n dim. of sdp    var  = %2.0d,',blkdim(1));
        fprintf('   num. of sdp  blk  = %2.0d',numblk(1));
    end
    if blkdim(2);
        fprintf('\n dim. of socp   var  = %2.0d,',blkdim(2));
        fprintf('   num. of socp blk  = %2.0d',numblk(2));
    end
    if blkdim(3); fprintf('\n dim. of linear var  = %2.0d',blkdim(3)); end
    if blkdim(4); fprintf('\n dim. of free   var  = %2.0d',blkdim(4)); end
end
%%
%%-----------------------------------------
%% detect unrestricted blocks in linear blocks
%%-----------------------------------------
%%
user_supplied_schurfun = 0;
for p = 1:size(blk,1)
    if ~isempty(par.schurfun{p}); user_supplied_schurfun = 1; end
end
if (user_supplied_schurfun == 0)
    [blk2,At2,C2,ublkinfo,parbarrier2,X02,Z02] = ...
        detect_ublk(blk,At,C,parbarrier,X0,Z0,par.printlevel);
else
    blk2 = blk; At2 = At; C2 = C;
    parbarrier2 = parbarrier; X02 = X0; Z02 = Z0;
    ublkinfo = cell(size(blk2,1),1);
end
ublksize = blkdim(4);
for p = 1:size(ublkinfo,1)
    ublksize = ublksize + length(ublkinfo{p});
end
%%
%%-----------------------------------------
%% detect diagonal blocks in semidefinite blocks
%%-----------------------------------------
%%
if (user_supplied_schurfun==0)
    [blk3,At3,C3,diagblkinfo,diagblkchange,parbarrier3,X03,Z03] = ...
        detect_lblk(blk2,At2,C2,b,parbarrier2,X02,Z02,par.printlevel); %#ok
else
    blk3 = blk2; At3 = At2; C3 = C2;
    % parbarrier3 = parbarrier2;
    X03 = X02; Z03 = Z02;
    diagblkchange = 0;
    diagblkinfo = cell(size(blk3,1),1);
end
%%
%%-----------------------------------------
%% main solver
%%-----------------------------------------
%%
%exist_analytic_term = 0;
%for p = 1:size(blk3,1);
%    idx = find(parbarrier3{p} > 0);
%    if ~isempty(idx); exist_analytic_term = 1; end
%end
%%
if (par.vers == 0);
    if blkdim(1); par.vers = 1; else par.vers = 2; end
end
par.blkdim = blkdim;
par.ublksize = ublksize;
[obj,X3,y,Z3,info,runhist] = ...
    HSDsqlpmain(blk3,At3,C3,b,par,X03,y0,Z03,kap0,tau0,theta0);
%%
%%-----------------------------------------
%% recover semidefinite blocks from linear blocks
%%-----------------------------------------
%%
if any(diagblkchange)
    X2 = cell(size(blk2,1),1); Z2 = cell(size(blk2,1),1);
    count = 0;
    for p = 1:size(blk2,1)
        pblk = blk2(p,:);
        n = sum(pblk{2});
        blkno      = diagblkinfo{p,1};
        idxdiag    = diagblkinfo{p,2};
        idxnondiag = diagblkinfo{p,3};
        if ~isempty(idxdiag)
            len = length(idxdiag);
            Xtmp = [idxdiag,idxdiag,X3{end}(count+1:count+len); n, n, 0];
            Ztmp = [idxdiag,idxdiag,Z3{end}(count+1:count+len); n, n, 0];
            if ~isempty(idxnondiag)
                [ii,jj,vv] = find(X3{blkno});
                Xtmp = [Xtmp; idxnondiag(ii),idxnondiag(jj),vv]; %#ok
                [ii,jj,vv] = find(Z3{blkno});
                Ztmp = [Ztmp; idxnondiag(ii),idxnondiag(jj),vv]; %#ok
            end
            X2{p} = spconvert(Xtmp);
            Z2{p} = spconvert(Ztmp);
            count = count + len;
        else
            X2(p) = X3(blkno); Z2(p) = Z3(blkno);
        end
    end
else
    X2 = X3; Z2 = Z3;
end
%%
%%-----------------------------------------
%% recover linear block from unrestricted block
%%-----------------------------------------
%%
numblk = size(blk,1);
numblknew = numblk;
X = cell(numblk,1); Z = cell(numblk,1);
for p = 1:numblk
    n = blk{p,2};
    if isempty(ublkinfo{p,1})
        X{p} = X2{p}; Z{p} = Z2{p};
    else
        Xtmp = zeros(n,1); Ztmp = zeros(n,1);
        Xtmp(ublkinfo{p,1}) = max(0,X2{p});
        Xtmp(ublkinfo{p,2}) = max(0,-X2{p});
        Ztmp(ublkinfo{p,1}) = max(0,Z2{p});
        Ztmp(ublkinfo{p,2}) = max(0,-Z2{p});
        if ~isempty(ublkinfo{p,3})
            numblknew = numblknew + 1;
            Xtmp(ublkinfo{p,3}) = X2{numblknew};
            Ztmp(ublkinfo{p,3}) = Z2{numblknew};
        end
        X{p} = Xtmp; Z{p} = Ztmp;
    end
end
%%
%%-----------------------------------------
%% recover complex solution
%%-----------------------------------------
%%
if (iscmp)
    for p = 1:numblk
        pblk = blk(p,:);
        n = sum(pblk{2})/2;
        if strcmp(pblk{1},'s');
            X{p} = X{p}(1:n,1:n) + sqrt(-1)*X{p}(n+1:2*n,1:n);
            Z{p} = Z{p}(1:n,1:n) + sqrt(-1)*Z{p}(n+1:2*n,1:n);
            X{p} = 0.5*(X{p}+X{p}');
            Z{p} = 0.5*(Z{p}+Z{p}');
        end
    end
end
if (isemptyAtb)
    X = X(1:end-1); Z = Z(1:end-1);
end
%%*****************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% GLG
%
% Files
%   compute_moments1D       - Compute moments of noise free wavelet coefficients
%   compute_moments2D       - Compute moments of noiseless wavelet coefficients
%   EM_root                 - The EM algorithm for the top level of the GLG model for 1D signals
%   EM_root_exact           - The EM algorithm for the top level of the GLG model for 1D signals
%   EM_tree                 - The EM algorithm for the tree of the GLG model
%   GLG_EM_wrapper          - Fit GLG model to 1D wavelet tree with composite EM algorithms
%   GLG_kld                 - Kullback-Leibler divergence between GLG models
%   GLG_llh                 - Evaluate the (composite) likelihood of data under the estimated GLG model
%   GLG_marginal_density    - Compute the marginal density of wavelet coefs in the GLG model
%   GLG_noise_removal       - Noise removal in the GLG model
%   GLG_noise_removal_exact - Noise removal in the GLG model
%   GLG_simulation          - Simulation from a GLG model
%   hermite_quad            - Find the Gauss-Hermite abscissae and weights.
%   joint_density           - Evaluate the joint density of hidden and observed variables
%   moment_estimation       - Parameter estimates from moment equations
%   noisy_joint_density     - Evaluate the joint density of hidden and observed variables
%   noisy_moment_estimation - Parameter estimates from moment equations

author: Sable

repo: mcbench-benchmarks


text:
function c8_sinh_test ( )

%*****************************************************************************80
%
%% C8_SINH_TEST tests C8_SINH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2015
%
%  Author:
%
%    John Burkardt
%
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_SINH_TEST\n' );
  fprintf ( 1, '  C8_SINH computes the hyperbolic sine of a C8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '       C1=C8_UNIFORM_01          C2 = C8_SINH(C1)           C3 = C8_ASINH(C1)\n' );
  fprintf ( 1, '     ---------------------     ---------------------     ---------------------\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );
    c2 = c8_sinh ( c1 );
    c3 = c8_asinh ( c2);

    fprintf ( 1, '  (%12f  %12f)  (%12f  %12f)  (%12f  %12f)\n', ...
      real ( c1 ), imag ( c1 ), real ( c2 ), imag ( c2 ), real ( c3 ), imag ( c3 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear all
close all

disp('Temporal discounting model');

% Parameters for first level
Nobs=200;

[M,U] = mci_discount_struct (Nobs);

% Generate Data
%w_true = spm_normrnd(M.pE,M.pC,1);
w_true = [0.25, 1.5]';
w_true = log(w_true);
[g,Y] = mci_discount_gen (w_true,M,U);

[a,v1,v2,k] = mci_discount_act (w_true,M,U);

t=sort(U.t1);
r=mean(U.r1);
figure;plot(t,r./(1+k*t),'k.');
set(gca,'FontSize',16);
xlabel('Delay, Weeks');

t=sort(U.t2);
r=mean(U.r2);
hold on
plot(t,r./(1+k*t),'r.');
ylabel('Mean Reward, £');
grid on

disp(sprintf('Average prob of choosing first option = %1.2f',mean(g)));

%mcmc.inference='amc';
%mcmc.inference='vl';
mcmc.inference='langevin';

mcmc.verbose=0;
mcmc.maxits=1024;

post = spm_mci_post (mcmc,M,U,Y,w_true);
disp('Posterior Mean from MCI');
disp(post.Ep)

plot_post=1;
if plot_post
    % Plot posterior surface
    S.Nbins=50;
    r=2;
    S.pxy(1,:)=linspace(-r*w_true(1),r*w_true(1),S.Nbins);
    S.pxy(2,:)=linspace(-0.5*r*w_true(2),1.5*r*w_true(2),S.Nbins);
    S.param{1}='P(1)';
    S.param{2}='P(2)';
    S.name={'log k','log \beta'};
    [L,S] = mci_plot_surface (w_true,M,U,Y,S,'post');
    hold on
    ms=10;
    plot(M.pE(1),M.pE(2),'wo','MarkerSize',ms);
    plot(post.Ep(1),post.Ep(2),'wx','MarkerSize',ms);
    plot(w_true(1),w_true(2),'w+','MarkerSize',ms);
end

stats = spm_mci_mvnpost (post,'ESS')
stats = spm_mci_mvnpost (post,'thinning')

author: spm


text:
function percent = parfor_progress(N,filenameParfor)
%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor.
%   PARFOR_PROGRESS works by creating a file called parfor_progress.txt in
%   your working directory, and then keeping track of the parfor loop's
%   progress within that file. This workaround is necessary because parfor
%   workers cannot communicate with one another so there is no simple way
%   to know which iterations have finished and which haven't.
%
%   PARFOR_PROGRESS(N) initializes the progress monitor for a set of N
%   upcoming calculations.
%
%   PARFOR_PROGRESS updates the progress inside your parfor loop and
%   displays an updated progress bar.
%
%   PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress
%   bar.
%
%   To suppress output from any of these functions, just ask for a return
%   variable from the function calls, like PERCENT = PARFOR_PROGRESS which
%   returns the percentage of completion.
%
%   Example:
%
%      N = 100;
%      parfor_progress(N);
%      parfor i=1:N
%         pause(rand); % Replace with real code
%         parfor_progress;
%      end
%      parfor_progress(0);
%
%   See also PARFOR.

% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/

%error(narginchk(0, 1, nargin, 'struct'));

if nargin < 2
    N = -1;
end

percent = 0;
w = 50; % Width of progress bar

if N > 0
    f = fopen(filenameParfor, 'w');
    if f<0
        error('Do you have write permissions for %s?', pwd);
    end
    fprintf(f, '%d\n', N); % Save N at the top of progress.txt
    fclose(f);

    if nargout == 0
        disp(['  0%[>', repmat(' ', 1, w), ']']);
    end
elseif N == 0
    delete(filenameParfor);
    percent = 100;

    if nargout == 0
        disp([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']);
    end
else
    if ~exist(filenameParfor, 'file')
        error([filenameParfor 'not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize' filenameParfor '.']);
    end

    f = fopen(filenameParfor, 'a');
    fprintf(f, '1\n');
    fclose(f);

    f = fopen(filenameParfor, 'r');
    progress = fscanf(f, '%d');
    fclose(f);
    percent = (length(progress)-1)/progress(1)*100;

    if nargout == 0
        perc = sprintf('%3.0f%%', percent); % 4 characters wide, percentage
        disp([repmat(char(8), 1, (w+9)), char(10), perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']']);
    end
end

author: opencobra

repo: cobratoolbox


text:
function [w] = hann_p(npt)
%
% Type: [w] = hann_p(npt);
%
% Inputs:
%
% npt   := number of points in the window vector
%
% Outputs:
%
% w     := periodic hanning window vector, npt x 1
%
% Compute the periodic hanning window.

% Scot McNeill, University of Houston, Fall 2007.
%
w=0.5*(1-cos(2*pi*[1:npt].'/(npt)));
%
%%


text:
%% MULTIGRID FOR STOKES EQUATIONS IN 2D
%
% This example is to show the convergence of multigrid methods for various
% finite element approximation of the Stokes equation on the unit square:
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0  in \Omega,
%                              u = g_D  on \Gamma.
%
% with the pure Dirichlet boundary condition. The solver is based on a DGS
% type smoother and summarized in the following reference.
%
% Reference
%
% M. Wang and L. Chen. Multigrid Methods for the Stokes equations
% using Distributive Gauss-Seidel Relaxations based on the Least Squares
% Commutator. Journal of Scientific Computing. 56(2): 409-431, 2013.

clear variables;
close all;

%% Setting
% mesh
[node,elem] = squaremesh([0,1,0,1],0.25);
% [node,elem] = circlemesh(0,0,1,0.25);
showmesh(node,elem);
[node,elem] = uniformrefine(node,elem);
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);
figure; showmesh(node,elem); pause(0.5);
% pde
pde = Stokesdata2;
%  pde = StokesZulehnerdata;
% options
option.L0 = 0;
option.maxIt = 4;
option.printlevel = 1;
option.plotflag = 0;
option.rateflag = 0;

%% MG options
option.solver = 'mg';
option.smoothingstep = 2;
option.smootherbarSp = 'SGS';

%% RT0-P0
disp('RT0-P0')
option.elemType = 'RT0-P0';
% option.refType = 'bisect';
femStokesHdiv(mesh,pde,option);

%% BDM1B-P0
disp('BDM1B-P0')
option.elemType = 'BDM1B-P0';
femStokesHdiv(mesh,pde,option);

%% CR-P0 element
disp('CR-P0')
option.elemType = 'CRP0';
femStokes(mesh,pde,option);

%% P2-P0 element
disp('P2-P0')
option.elemType = 'P2P0';
femStokes(mesh,pde,option);

%% isoP2-P0 element
disp('isoP2-P0')
option.elemType = 'isoP2P0';
femStokes(mesh,pde,option);

%% isoP2-P1 element
disp('isoP2-P1')
option.elemType = 'isoP2P1';
femStokes(mesh,pde,option);

%% P1b-P1 element
disp('P1b-P0')
option.elemType = 'P1bP1';
option.solver = 'asmg';
femStokes(mesh,pde,option);

%% P2-P1 element
disp('P2-P1')
option.elemType = 'P2P1';
% option.smoothingStep = 3;
% option.smootherbarSp   = 'VCYCLE';
% option.smootherbarSpPara = 0.75;
femStokes(mesh,pde,option);

repo: ifem


text:
% module = factor_module_gp_factorized(N, covfuncs, theta, varargin)

function module = factor_module_gp_factorized(N, covfuncs, theta, varargin)

D = length(covfuncs);

%
% Parse options
%
options = struct('is_pseudo', false(D,1), ...
                 'noise_prior', [], ...
                 'update_hyperparameters', true, ...
                 'maxiter_hyperparameters', 30, ...
                 'init', [], ...
                 'update_schedule', []);
[options, errmsg] = argparse(options, varargin{:});
error(errmsg);

if isempty(options.update_schedule)
  options.update_schedule = {1:D};
end

K_xx = cell(D,1);
K_px = cell(D,1);
K_pp = cell(D,1);

linsolve_K = cell(D,1);

X = zeros(D,N);
CovX = zeros(D,N);
for d=1:D
  if ~options.is_pseudo(d)
    K_xx{d} = covfuncs{d}(theta{d});
  else
    [K_pp{d}, K_px{d}, K_xx{d}] = covfuncs{d}(theta{d});
  end
end
if isfield(options.init, 'X')
  X(:,:) = options.init.X;
else
  % Initialize from the prior
  for d=1:D
    if ~options.is_pseudo(d)
      L = chol(K_xx{d}, 'lower');
      X(d,:) = L * randn(N,1);
    else
      L = chol(K_pp{d}, 'lower');
      Xp = L * randn(length(L),1);
      X(d,:) = K_px{d}' * linsolve_lchol(L, Xp);
    end
  end
end
if isfield(options.init, 'CovX')
  CovX(:,:) = options.init.CovX;
else
  % Initialize with prior
  for d=1:D
    if ~options.is_pseudo(d)
      CovX(d,:) = 0.1*diag(K_xx{d});
    else
      CovX(d,:) = 0.1*K_xx{d};
    end
  end
end
rho = 1;

module.get_struct = @get_struct;
%module.initialize = @initialize;
module.update = @update;
module.rotate = @rotate;
module.rotation_cost = @rotation_cost;

  function S = get_struct()
  % TODO: not ready yet..
  S.module = 'factor_module_gp_factorized';
  S.X = X;
  S.CovX = CovX;
  S.rho = rho;
  S.covfuncs = covfuncs;
  S.theta = theta;
  S.init = options.init;
  S.options = options;
  end

% $$$   function [Z, CovZ, rhoz] = initialize()
% $$$   end

  function [Z, CovZ, rhoz, logrhoz, KL] = update(iter, Y, Obs, W, CovW, ...
                                                 rhow, Tau)

% $$$   whos
% $$$   for d=1:D
% $$$     nonzeros = nnz(K_xx{d}) / numel(K_xx{d});
% $$$     fprintf('Nonzeros in %d=%.3f\n', d, nonzeros);
% $$$   end

  linsolve_K(:) = {[]};

  M = size(Tau,1);

  % Update each component (component-wise factorization)
  KL = 0;
  Tau(~Obs) = 0;

  if any(rhow~=1)
    error('rhow ~= 1 not yet implemented')
  else
    v_W = ones(M,1);
  end

  updateable = options.update_schedule{min(iter,length(options.update_schedule))};

  for d=1:D

    if any(updateable==d)

      % Indeces of the other components
      others = [1:(d-1), (d+1):D];

      %
      % Pre-evaluations
      %

      % TODO: What if there are totally missing rows in Y??
      % Then u will contain zeros..

      if ndims(CovW) == 2 && all(size(CovW)==[D,M])
        Yh = W(others,:)'*X(others,:);
        u = Tau' * (W(d,:).^2+CovW(d,:))';
        z = (Tau.*(Y-Yh))' * W(d,:)';
      else
        % TODO: Optimize this.. Do not use loops..
        %
        % See help from ARD-module..
        %warning('Optimize this..');
        u = zeros(N,1);
        z = zeros(N,1);
        for n=1:N
          obs = Obs(:,n);

          % Is this correct
          covw = reshape(CovW(d,others,obs),[length(others), sum(obs)]);
% $$$           z(n) = (Tau(obs,n).*v_W(obs).*W(d,obs)')' * (Y(obs,n)-W(others,obs)) - ...
% $$$                  (Tau(obs,n)'*covw')*X(others,n);

          z(n) = (Tau(obs,n).*v_W(obs).*W(d,obs)')'*Y(obs,n) - ...
                 (Tau(obs,n)'*covw')*X(others,n);

          varw = vec(CovW(d,d,obs));
          ww = v_W(obs).*(W(d,obs)'.^2) + varw;
          u(n) = vec(Tau(obs,n))' * ww;
        end
      end
      % Remove totally missing observations
      ind_obs = (u~=0);
      z = z(ind_obs);
      u = u(ind_obs);
      % Helpful variables
      invU_z = z./u;
      invU = spdiag(1./u);
      % Noise covariance function
      covfunc_noise = gp_cov_wrap(invU);

      %
      % Update variables
      %

      if ~options.is_pseudo(d)

        % Standard full GP

        % Form the joint covariance function
        covfunc = gp_cov_sum(gp_cov_select(covfuncs{d}, ind_obs), covfunc_noise);

        % Learn the hyperparameters
        if index_selected(iter, options.update_hyperparameters) && ...
              any(d==updateable)
          [theta{d}, logpdf] = gp_learn(invU_z, covfunc, theta{d}, ...
                                        'maxiter', options.maxiter_hyperparameters, ...
                                        'checkgrad', false);
          K_xx{d} = covfuncs{d}(theta{d});
        else
          if isempty(K_xx{d})
            K_xx{d} = covfuncs{d}(theta{d});
          end
          logpdf = gp_loglikelihood(invU_z, covfunc, theta{d});
        end

        % Evaluate the mean and the variance
        y = invU_z;
        K_y = K_xx{d}(ind_obs,ind_obs) + invU;
        K_yf = K_xx{d}(ind_obs,:);
        k_f = diag(K_xx{d});
        [X(d,:),CovX(d,:)] = gp_predict(y, K_y, K_yf, k_f);
% $$$       [X(d,:),CovX(d,:)] = gp_predict(invU_z, K_xx{d}+invU, K_xx{d}, ...
% $$$                                  diag(K_xx{d}));

      else

        % Sparse approximation GP using pseudo inputs

        % Remove totally missing observations
        covfunc = gp_cov_select_pseudo(covfuncs{d}, ind_obs);

        % Learn the hyperparameters
        if index_selected(iter, options.update_hyperparameters)

          [theta{d}, ~, logpdf] = gp_learn_pseudo(invU_z, ...
                                                  covfunc, ...
                                                  theta{d}, ...
                                                  covfunc_noise, ...
                                                  [], ...
                                                  'maxiter', ...
                                                  options.maxiter_hyperparameters, ...
                                                  'checkgrad', false);
          [K_pp{d}, K_px{d}, K_xx{d}] = covfuncs{d}(theta{d});
        else
          logpdf = gp_loglikelihood_pseudo(invU_z, covfunc, theta{d}, ...
                                           covfunc_noise, []);
          if isempty(K_pp{d}) || isempty(K_px{d}) || isempty(K_xx{d})
            [K_pp{d}, K_px{d}, K_xx{d}] = covfuncs{d}(theta{d});
          end
        end

        % Evaluate the mean and the variance
        [X(d,:),CovX(d,:)] = gp_predict_pseudo(invU_z, ...
                                               @(A) bsxfun(@times, u, A), ...
                                               K_px{d}(:,ind_obs), ...
                                               K_pp{d}, ...
                                               K_px{d}, ...
                                               K_xx{d});

      end

      %
      % Compute Kullback-Leibler divergence KL(q(X)||p(X))
      %
      logpdf_bound = gaussian_logpdf(z'*invU_z, ...
                                     z'*X(d,ind_obs)', ...
                                     (X(d,ind_obs).^2+CovX(d,ind_obs))*u, ...
                                     -sum(log(u)), ...
                                     N);

      KL = KL + logpdf_bound - logpdf;
    else

      % Do not update the variable
      if ~options.is_pseudo(d)
        if isempty(K_xx{d})
          K_xx{d} = covfuncs{d}(theta{d});
        end
      else
        if isempty(K_pp{d}) || isempty(K_px{d}) || isempty(K_xx{d})
          [K_pp{d}, K_px{d}, K_xx{d}] = covfuncs{d}(theta{d});
        end
      end

      KL = KL + 0; % ???

    end

  end

  % Update q(rho)
  rho = 1;
  logrho = 0;
  KL_rho = 0;

  KL = KL + KL_rho;

  Z = X;
  CovZ = CovX;
  rhoz = rho.*ones(N,1);
  logrhoz = logrho.*ones(N,1);

  %tsplot(X)

  end

  function [c, dc] = rotation_cost(A, U, S, V)
% $$$   if any(options.is_pseudo)
% $$$     error('Rotation cost not implemented for pseudo inputs yet..');
% $$$   end
  % Cost from <log q(X)>
  c = -N * logdet_diag(S);
  dc = -N * inv(A)'; % FIX THIS
  % Cost from -<log p(X)>
  AX = A*X;
  varAX = A.^2*CovX;
  invdiagK = zeros(D,N);
  invK_XA = zeros(N,D);
  for d = 1:D
    if ~options.is_pseudo(d)
      if isempty(linsolve_K{d})
        linsolve_K{d} = get_linsolve_cov(K_xx{d});
      end
      invdiagK(d,:) = 1./diag(K_xx{d});
      invK_XA(:,d) = linsolve_K{d}(AX(d,:)');
    else
      %error('Pseudo cost not maybe correct..');
      if isempty(linsolve_K{d})
        L = chol(K_pp{d}, 'lower');
        Z = linsolve_tril(L, K_px{d});
        h = K_xx{d} - dot(Z,Z,1)';
        h = h + 1e-6*max(h);
        %H = spdiag(h);
        %Z = L(K_xp{d}*diag(h.^(-0.5)));
        H = spdiag(h);
        %plot(h)
        Lambda = speye(size(L)) + Z*spdiag(1./h)*Z';
        L_Lambda = get_linsolve_cov(Lambda);
        linsolve_K{d} = @(x) (H\x - H\(Z'*L_Lambda(Z*(H\x))));

% $$$         L = get_linsolve_cov(K_pp{d});
% $$$         h = K_xx{d} - diagprod(K_xp{d}, L(K_xp{d}')');
        %H = 1e-6*speye(size(H)); % better conditioning...
        %LL = get_linsolve_cov(K_pp{d}+K_xp{d}'*(H\K_xp{d}));
        %linsolve_K{d} = @(x) (H\x - H\(K_xp{d}*LL(K_xp{d}')*(H\x)));
      end
      invK_XA(:,d) = linsolve_K{d}(AX(d,:)');
      invdiagK(d,:) = 1./K_xx{d};
      %invK_XA(:,d) = AX(d,:)' ./ K_xx{d};
    end
    c = c + 0.5 * rho * (AX(d,:)*invK_XA(:,d));
    c = c + 0.5 * varAX(d,:) * invdiagK(d,:)';
  end
  dc = dc + rho * (invK_XA'*X');
  dc = dc + A .* (invdiagK*CovX');
% $$$   dc = dc + 0.5 * (invdiagK*CovX'*A' + A'*invdiagK*CovX');

  end

  function [Z, CovZ] = rotate(A)
  X = A*X;
  CovX = A.^2*CovX;
  Z = X;
  CovZ = CovX;
  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
 function y = ir_phantom_resize(x, nx, ny)
%function y = ir_phantom_resize(x, nx, ny)
%| resize a [mx,my] phantom image to be [nx,ny]
%| by combination of downsampling and cropping

if nargin == 1 && streq(x, 'test'), ir_phantom_resize_test, return, end
if nargin < 3, ir_usage, end

mx = size(x,1);
my = size(x,2);
if nx == mx && ny == my
	y = x;
	return
end

if nx > mx || ny > my
	error 'requested size too large, implement interpolation!'
end

% if power of 2 size, then downsample
if 2^floor(log2(nx)) == nx && nx == ny
	y = downsample2(x, mx/nx);
return
end

nn = 2^max(ceil(log2([nx ny])));
y = downsample2(x, mx/nn);

% trim if non power of 2
if nx < nn
	y = y([1:nx]+round((nn-nx)/2),:);
end
if ny < nn
	y = y(:,[1:ny]+round((nn-ny)/2));
end

function ir_phantom_resize_test
x = zeros(64, 60);
x(20:30,10:15) = 1;
y = ir_phantom_resize(x, 32, 28);
im plc 1 2
im(1, x)
im(2, y)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [img, nSamples] = ellipseInterior(varargin)
% Return an image (matrix) marking the interior of an ellipse
%
% Syntax
%   [img, nSamples] = ellipseInterior(varargin)
%
% Inputs
%   N/A
%
% Optional key/value pairs
%   center  - Center of the ellipse
%   sigma   - Sigma major and minor as a 2-vector
%   theta   - Angle of sigma major from the x-axis (counter clockwise)
%
% Returns
%   img:  1 in the interior and 0 in the exterior
%   nSamples:  Number of interior samples
%
% Description
%   We start with an ellipse with its major axis aligned to the
%   x-axis.  We then rotate (counter clockwise) by an amount theta.
%
%  This is the equation for the contour of the ellipse
%
%     1 = A*(x-h)^2 + 2B*(x-h)*(y-k) + C*(y-k)^2
%
%  The quadratic form is Q = [A B; B C];
%  We test that the parameters are valid (i.e., Q is positive-definite
%  (i.e., Q = M*M') by
%
%    A > 0
%    det(Q) > 0
%
%  We convert from sigma(1), sigma(2), theta to Q and inversely using
%
%        Q = Rotate(theta)'*[diag(a^2,b^2)]*Rotate(theta)
%{
%  [U,S,V] = svd(Q)
%  Figure out theta from U and get a, b from sqrt(S)
%  Since
%      rotMat = [cos(theta) -sin(theta); -sin(theta) -cos(theta)];
%      theta = acos(U(1,1))
%
% Some experimenting
%
% It must be a > b for major/minor axes
a = 3; b = 2; theta = pi/6;
rotMat = [cos(theta) -sin(theta); -sin(theta) -cos(theta)];
Q = rotMat' * diag([a^2, b^2]) * rotMat;

[U,S,V] = svd(Q);
abs(U) - abs(rotMat)  % Not sure why there can be a sign difference
thetaEst = min(acos(U(1,1)),acos(-1*U(1,1)));
theta - thetaEst
%}
%
% Notes:
%
%  https://math.stackexchange.com/questions/264446/the-fastest-way-to-obtain-orientation-%CE%B8-from-this-ellipse-formula
%
% Wandell, January 13 2020
%
% See also
%    ellipsePlot()

% Examples:
%{
   img = ellipseInterior;
   mrvNewGraphWin; imagesc(img); axis equal
%}
%{
   samples = [-3:0.1:3];
   [img, nSamples]  = ellipseInterior('center',[1 2],'spatial samples',samples);
   [img2, nSamples] = ellipseInterior('center',[0.5 2],'spatial samples',samples);

   % mrvNewGraphWin; imagesc(img); axis image
   % mrvNewGraphWin; imagesc(img2); axis image
   overlap = dot(img(:), img2(:))/nSamples;
   fprintf('Overlap is %f\n',overlap);
%}
%{
   samples = [-10:0.1:10];
   sigma = [3,1];
   [img, nSamples]   = ellipseInterior('center',[0 0],'sigma',sigma,'spatial samples',samples);
   [img2, nSamples]  = ellipseInterior('center',[-1 -4],'sigma',sigma,'spatial samples',samples);
   mrvNewGraphWin; imagesc(samples,samples,img); axis image
   mrvNewGraphWin; imagesc(samples,samples,img2); axis image
   overlap = dot(img(:), img2(:))/nSamples;
   fprintf('Overlap is %f\n',overlap);
%}
%{
   img = ellipseInterior('center',[1 1],'quadratic',0.3*eye(2));
   mrvNewGraphWin; imagesc(img);
%}
%% Input parameters

varargin = mrvParamFormat(varargin);

p = inputParser;
p.addParameter('center',[0 0],@isvector);
p.addParameter('spatialsamples',(-3:0.05:3),@isvector);
vFunc = @(x)(length(x) == 2 && x(1) >= x(2) && isnumeric(x));
p.addParameter('sigma',[1 1],vFunc);   % Length in degrees
p.addParameter('theta',0,@isvector);   % Radians

p.parse(varargin{:});
c = p.Results.center;
samples = p.Results.spatialsamples;
sigma   = p.Results.sigma;
theta   = p.Results.theta;

%%  Test that we have a true ellipse formula

rotMat = [cos(theta) -sin(theta); -sin(theta) -cos(theta)];
Q = rotMat' * diag([1/sigma(1)^2, 1/sigma(2)^2]) * rotMat;
if Q(1,1) <= 0 || det(Q) <= 0 || Q(1,2) ~= Q(2,1)
    error('Q is not positive definite');
end

%% Create the spatial samples and compute the values

[X,Y] = meshgrid(samples,samples);

X = X - c(1);
Y = Y - c(2);
V = Q(1,1)*X.^2 + 2*Q(1,2)*X.*Y + Q(2,2)*Y.^2;
% mrvNewGraphWin; mesh(V)

%% Binarize
img = zeros(size(X));
img(V<=1) = 1;

% mrvNewGraphWin; imagesc(img)

if nargout > 1
    nSamples = sum(img(:));
end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function srcPatch = sc_prep_source_patch(img, uvTform, optS)

% SC_PREP_SOURCE_PATCH
%
% Prepare source patches according to uvTform
%
% Input:
%   - img:       input image
%   - uvPixSub:  target patch position
%   - optS:      parameter
% Output:
%   - srcPatch: [pNumPix] x [3] x [numUvPix]

numUvPix =  size(uvTform, 1);

% Prepare source patch sampling position
% srcPatchPos = zeros(optS.pNumPix, 3, numUvPix);

% Get srcPatchPos
c1 = reshape(uvTform(:,1:3)', 1, 3, numUvPix);
c2 = reshape(uvTform(:,4:6)', 1, 3, numUvPix);
c3 = reshape(uvTform(:,7:9)', 1, 3, numUvPix);

% Get the source patch pixel positions
srcPatchPos = bsxfun(@times, optS.refPatchPos(:,1), c1) + ...
              bsxfun(@times, optS.refPatchPos(:,2), c2);
srcPatchPos = bsxfun(@plus, srcPatchPos, c3);

% Convert back to Eucledian coordinate
srcPatchPos = bsxfun(@rdivide, srcPatchPos, srcPatchPos(:, 3, :));

% Grab the color values of source patch using bilinear interpolation
srcPatch = vgg_interp2(img, srcPatchPos(:,1,:), srcPatchPos(:,2,:), 'linear', 0);

% Convert to the format [pNumPix] x [3] x [numUvPix]
srcPatch = permute(srcPatch, [1,3,2]);

% for i = 1 : optS.pNumPix
%     dx = optS.refPatchPos(1,i);
%     dy = optS.refPatchPos(2,i);
%
%     srcPatchPos(i, 1, :) = uvTform(1,:)*dx + uvTform(4,:)*dy + uvTform(7,:);
%     srcPatchPos(i, 2, :) = uvTform(2,:)*dx + uvTform(5,:)*dy + uvTform(8,:);
%     srcPatchPos(i, 3, :) = uvTform(3,:)*dx + uvTform(6,:)*dy + uvTform(9,:);
% end
% srcPatchPos(:, 1:2,:) = bsxfun(@rdivide, srcPatchPos(:, 1:2,:), srcPatchPos(:, 3,:));

% % Avoid sample out of boundary positions
% % srcPatchPos(:,1,:) = sc_clamp(srcPatchPos(:,1,:), 1, size(img,2));
% % srcPatchPos(:,2,:) = sc_clamp(srcPatchPos(:,2,:), 1, size(img,1));
%
% % Sampling source patch
% srcPatch = mirt2D_mexinterp(img, srcPatchPos(:, 1, :), srcPatchPos(:, 2, :));

% srcPatch = permute(srcPatch, [1,3,2]);
end

repo: StructCompletion


text:
function gradWeight = getGradWtFromClusters(viewClusters,i)
    thisCluster = find(cellfun(@(x)(sum(x==i)),viewClusters));
    meanElems = mean(cellfun(@length,viewClusters));
    gradWeight = meanElems/length(viewClusters{thisCluster});
end

repo: CategoryShapes


text:
% ----------------------------------------------------------------------------
% function hfssRotate(fid, ObjList, Axis, Degrees)
%
% Description :
% -------------
% Creates the VB Script necessary to rotate a given set of objects.
%
% Parameters :
% ------------
% fid     - file identifier of the HFSS script file.
% ObjList - a cell array of objects that need to be rotated.
% Axis    - axis of the rotation.
% Degrees - value of the rotation in degrees.
%
% Note :
% ------
%
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ... create some objects here ...
% hfssRotate(fid, {'Dipole1', 'Dipole2'}, 'X', -30);
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% CHANGELOG
%
% 23-Sep-2012: *Initial release (DRP).
% 21-Apr-2021: *Fix rotation of just one element with string name (DRP).
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% Written by Daniel Rodriguez Prado
% danysan@gmail.com
% 23 September 2012
% ----------------------------------------------------------------------------
function hfssRotate(fid, ObjectList, Axis, Degrees)

if (~iscell(ObjectList)) % Fixes rotation when there is only one string.
    ObjectList = {ObjectList};
end

nObjects = length(ObjectList);

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'oEditor.Rotate _\n');
fprintf(fid, 'Array("NAME:Selections", _\n');

% Object Selections.
fprintf(fid, '"Selections:=", "');
for iObj = 1:nObjects
	fprintf(fid, '%s', ObjectList{iObj});
	if (iObj ~= nObjects)
		fprintf(fid, ',');
    end
end
fprintf(fid, '"), _\n');

% Transalation Vector.
fprintf(fid, 'Array("NAME:RotateParameters", _\n');
fprintf(fid, '"RotateAxis:=", "%s", _\n', upper(Axis));
fprintf(fid, '"RotateAngle:=", "%fdeg")\n', Degrees);

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
% lineSegmentFitting extracts LineSegments from each EdgeContour.
% Input     - EdgeContours : The list of Edge Contours.
%           - DeviationThreshold : The threshold chosen to determine how
%                       close the line segments fit on to the edge contour.
% output    - SegmentList : The list of Segments. Each of them corresponds
%                          to each Edge Contour.
function SegmentList = lineSegmentFitting(EdgeContours, DeviationThreshold)

%% Initialize
NumberOfContours = length(EdgeContours);
SegmentList = cell(1,NumberOfContours);

%% Exctrace Line Segemnts from each Edge Contour.
for i =1:NumberOfContours
    EdgePoints = EdgeContours{i};
    LastIndex = length(EdgePoints);
    LastPoint = EdgePoints(LastIndex,:);
    LineSegments = [extractLinearCue(EdgePoints, DeviationThreshold);
        LastPoint];
    SegmentList{i} = LineSegments;
end

end

% extractLinearCue transforms the point-wise Edge Contour into
% the piece-wise linear segments.
% Input     - EdgePoints : The subset points of EdgeContour which will be
%                         transformed into the LineSegment
%           - DeviationThreshold : The threshold chosen to determine how
%                       close the line segments fit on to the edge contour.
% output    - LineSegments : The Line Segments extracted from edge points.
%                           For efficient recursion, It doesn't contain the
%                           last point.
% Example   - Given EdgePoints = [first, p_2, p_3, ..., p_n-1, last],
%             LineSegments could be [first, p_k1, p_k2, ..., p_km], (last)
%             , where k1 < k2 < ... < km
%               and each p_ki belongs to p_j for 2<=j<=n-1
function LineSegments = extractLinearCue(EdgePoints, DeviationThreshold)
%% Initialize
EdgeLength = length(EdgePoints);

% Initialize first and last points.
FirstPoint = EdgePoints(1,:);
LastPoint = EdgePoints(EdgeLength,:);

% If EdgePoints contains only two points, which are first and last,
% return only first point. (Since we remove the last point from the output)
if EdgeLength == 2
    LineSegments = FirstPoint;
else

    max_d = 0;

    % Find maximum deviation
    for i = 2:EdgeLength-1
        CurPoint = EdgePoints(i,:);
        % The deviation di for Point i is denoted :
        % di = | xi(y1-yn)+yi(xn-x1)+ynx1-y1xn |/sqrt((y1-yn)^2+(xn-x1)^2)
        d = abs( CurPoint(2)*(FirstPoint(1)-LastPoint(1)) + ...
            CurPoint(1)*(LastPoint(2)-FirstPoint(2)) + ...
            LastPoint(1)*FirstPoint(2) - FirstPoint(1)*LastPoint(2) ) / ...
            sqrt((FirstPoint(1)-LastPoint(1))^2 + (LastPoint(2)-FirstPoint(2))^2);

        % Compare di with maximum value
        if max_d < d
            max_d = d;
            max_idx = i;

        end
    end

    % If the maximum deviation is smaller than the Threshold,
    % consider FirstPoint-LastPoint as one line segment.
    % Since we remove the last point from the output, return only
    % first point.
    if max_d < DeviationThreshold
        LineSegments = FirstPoint;
    else
        LeftEdgePoints = EdgePoints(1:max_idx,:);
        RightEdgePoints = EdgePoints(max_idx:EdgeLength,:);
        LineSegments = [ extractLinearCue(LeftEdgePoints, DeviationThreshold) ;
            extractLinearCue(RightEdgePoints, DeviationThreshold)];
    end

end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% HLP_ARGSTRUCT2LINEARCELL - Linearize configation output of arg_guipanel
%
% Usage:
%   >> cellval = hlp_argstruct2linearcells( cfg );
%
% Inputs:
%   cfg  - output configuration structure from arg_guipanel
%
% Output:
%   cellval - cell array of output values
%
% Author: Arnaud Delorme, SCCN & CERCO, CNRS, 2013-

% Copyright (C) 2013 Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function cellval = hlp_argstruct2linearcell(cfg);

    cellval = {};
    if isstruct(cfg)
        ff = fieldnames(cfg);
        for iField = 1:length(ff)
            if ~strcmpi(ff{iField}, 'arg_direct')
                if strcmpi(ff{iField}, 'arg_selection')
                    cellval = { cfg.arg_selection cellval{:} };
                else
                    val = cfg.(ff{iField});
                    if isstruct(val)
                        val = hlp_argstruct2linearcell(val);
                        cellval = { cellval{:} ff{iField} val{:} };
                    elseif iscell(val)
                        cellval = { cellval{:} ff{iField} vararg2str(val) };
                    else
                        cellval = { cellval{:} ff{iField} val };
                    end
                end
            end
        end
    else
        cellval = cfg;
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% Visualize faces in a face recognition model.
% The model can be downloaded from
% https://github.com/AlfredXiangWu/face_verification_experiment .
% This model is trained with CASIA-webface, People_ID can be set from 1
% to 10575, to get the memorized face of the corresponding person.

% Based on paper:
% Feng Wang, Haijun Liu, Jian Cheng,
% Visualizing Deep Neural Network by Alternately Image Blurring and Deblurring
caffe.reset_all();
caffe.set_mode_gpu();
gpu_id = 0;  % we will use the first gpu in this demo
caffe.set_device(gpu_id);

% net_model = 'DeepFace_set003_inception.prototxt';% Please, remember to add force_backward:true to this file.
% net_weights = 'D:\project\lfw_face_verification_experiment\model\DeepFace_set003_net_iter.caffemodel';
net_model = 'LightenedCNN_B_deploy.prototxt';% Please, remember to add force_backward:true to this file.
net_weights = 'D:\project\lfw_face_verification_experiment\model\LightenedCNN_B.caffemodel';

mean_file = [];

train_net = caffe.Net(net_model,net_weights,'test');
if isempty(mean_file)
    mean_image = zeros(128,128);
else
    mean_image = caffe.read_mean(mean_file);
end;
for People_ID = 1:10575

% mean_image = mean_image(17:240,17:240,:);
% mean_image = mean_image + randn(size(mean_image));
input_data = zeros(size(mean_image,1), size(mean_image,2), 1, 1, 'single');

mean_file = 'webface_mean.proto';
mean_face = caffe.read_mean(mean_file);
mean_face = mean_face / 256;
input_data(:,:,1,1) = imresize(mean_face,[128,128]);

use_clip = false;
use_cv_norm = false;
use_weight_decay = false;
use_image_blur = false;
use_image_deblur = false;
use_gradient_blur = false;
use_dropout = false;
use_maxGray = false;

H = fspecial('gaussian',[7 7],1.2);
prob = train_net.forward({input_data});
% input_data = input_data - min(input_data(:));
% input_data = input_data / max(input_data(:));
[max_prob,max_idx] = sort(prob{1},'descend');
max_idx = People_ID;%max_idx(1);
this_prob = prob{1}(max_idx);
back_data = ones(size(prob{1}),'single') * -1;
back_data(max_idx) = 1;
back_cell = prob;
back_cell{1} = back_data;
blur_data = zeros(size(input_data));
base_lr = 0.01;
max_lr = 0.01;
lambda1 = 0.00001;
lambda2 = 0.1;
last_prob = -999;
momentum = 0.8;
momentum2 = 0.99;
lastgrad = zeros(size(mean_image));
lastgrad2 = zeros(size(mean_image));
mask = ones(size(mean_image,1), size(mean_image,2));
iter = 1;
dropout = 0.5;

while 1
    lr = base_lr;% * sqrt(this_prob / (1 - this_prob));
    res = train_net.backward(back_cell);

    bak_data = input_data;

    if use_gradient_blur
        res{1} = imfilter(res{1},H,'same');
    end;
    grad = res{1};

    if use_clip
        app_gradient = sum(abs(res{1} .* input_data(:,:,:,1)),3);
        app_gradient = app_gradient < mean(app_gradient(:)) * 0.5;
        clip_grad = reshape(res{1},[size(mean_image,1)*size(mean_image,2) 3]);
        clip_grad(app_gradient==1,:) = 0;
        clip_grad = reshape(clip_grad,size(input_data));
        res{1} = clip_grad;
    end;

    if use_cv_norm
        I = input_data(:,:,:,1);
%         Gx = sign(I(2:end-1,2:end-1,:) - I(1:end-2,2:end-1,:)) - sign(I(3:end,2:end-1,:) - I(2:end-1,2:end-1,:));
%         Gy = sign(I(2:end-1,2:end-1,:) - I(2:end-1,1:end-2,:)) - sign(I(2:end-1,3:end,:) - I(2:end-1,2:end-1,:));
        Gx = smoothL1(I(2:end-1,:,:) - I(1:end-2,:,:)) - smoothL1(I(3:end,:,:) - I(2:end-1,:,:));
        Gx = [smoothL1(I(1,:,:) - I(2,:,:)); Gx; smoothL1(I(end,:,:) - I(end-1,:,:))];
        Gy = smoothL1(I(:,2:end-1,:) - I(:,1:end-2,:)) - smoothL1(I(:,3:end,:) - I(:,2:end-1,:));
        Gy = [smoothL1(I(:,1,:) - I(:,2,:)) Gy smoothL1(I(:,end,:) - I(:,end-1,:))];
%         Gx = sign(I(2:end-1,:,:) - I(1:end-2,:,:)) - sign(I(3:end,:,:) - I(2:end-1,:,:));
%         Gx = [sign(I(1,:,:) - I(2,:,:)); Gx; sign(I(end,:,:) - I(end-1,:,:))];
%         Gy = sign(I(:,2:end-1,:) - I(:,1:end-2,:)) - sign(I(:,3:end,:) - I(:,2:end-1,:));
%         Gy = [sign(I(:,1,:) - I(:,2,:)) Gy sign(I(:,end,:) - I(:,end-1,:))];
            grad = grad - lambda2 * (Gx + Gy);
    end;

    lastgrad = (1 - momentum) * lr * res{1}   + momentum * lastgrad;%/ norm(res{1}(:))
    input_data(:,:,:,1) = input_data(:,:,:,1) + lastgrad;
%     lastgrad = (1 - momentum) * grad   + momentum * lastgrad;%/ norm(res(:))
%     lastgrad2 = (1 - momentum2) * grad.^2 + momentum2 * lastgrad2;%/ norm(res(:))
%     lg_correct = lastgrad ./ (1 - momentum^iter);
%     lg2_correct = lastgrad2 ./ (1 - momentum2^iter);
%     input_data(:,:,:,1) = input_data(:,:,:,1) + lr * lg_correct ./ (sqrt(lg2_correct) + 1e-8);

    if use_weight_decay
        input_data(:,:,:,1) = input_data(:,:,:,1) - lr * lambda1 * I;
    end;
    if use_maxGray
%         if max(input_data(:))>1
            input_data(:,:,:,1) = input_data(:,:,:,1) - min(input_data(:));
            input_data(:,:,:,1) = input_data(:,:,:,1) / max(input_data(:));
%         end;
    end;
%         input_data = (input_data -mean(input_data(:))) / std(input_data(:)) * 30;
%     end;

%     for_forward = reshape(input_data,[size(mean_image,1)*size(mean_image,2) 3]);
%     mask = rand(size(mean_image,1), size(mean_image,2)) < dropout;
%     for_forward(mask==1,:) = 0;
%     for_forward = reshape(for_forward,size(input_data));

    if mod(iter,10) ==0%&&iter<2000
        if mod(iter,20) ~= 0
            H = fspecial('gaussian',[5 5],rand()/2+0.5);
            if use_image_blur
                input_data(:,:,:,1) = imfilter(input_data(:,:,:,1),H,'same');
            end
        else
            if mod(iter,20) == 0
                if use_image_deblur
                    input_data(:,:,:,1) = deconvlucy(input_data(:,:,:,1), H);
                end;
            end;
        end;
    end;
    prob = train_net.forward({input_data});

    this_prob = prob{1}(max_idx);
    fprintf('id=%d,iter=%d,lr=%f,prob1=%f,last_prob=%f\n',People_ID,iter,lr,prob{1}(max_idx),last_prob);
    iter = iter + 1;

    if mod(iter,100)==0
        figure(2);
        % imshow(uint8(mean_image + input_data));
        output = input_data(:,:,:,1);
        output = output';
        imshow(output);
        I = output;
%         Gx = abs(I(2:end-1,2:end-1,:) - I(1:end-2,2:end-1,:)) + abs(I(3:end,2:end-1,:) - I(2:end-1,2:end-1,:));
%         Gy = abs(I(2:end-1,2:end-1,:) - I(2:end-1,1:end-2,:)) + abs(I(2:end-1,3:end,:) - I(2:end-1,2:end-1,:));
%         figure(3);hist(I(:),1000);
%         figure(4);hist(Gy(:),1000);
        if iter == 200
            break;
        end;
    end;
    if this_prob<last_prob
        base_lr = base_lr * 0.99;
%         input_data = bak_data;
    end;
    if this_prob>last_prob&&base_lr<max_lr%&& (this_prob-last_prob) / this_prob < 0.001
        base_lr = base_lr * 1.01;
%         input_data = bak_data;
    end;
%     if this_prob>last_prob
        last_prob = this_prob;
%     end;
    if lr<0.000001
        break;
    end;
end;
output = input_data(:,:,:,1);
output = output';
imwrite(output, ['gallery/' num2str(People_ID) '.png']);
end;

author: happynear

repo: DeepVisualization


text:
function mc1 = moc1 ( a, b, n, f )

%*****************************************************************************80
%
%% MOC1 estimates a function related to the modulus of continuity.
%
%  Discussion;
%
%    The modulus of continuity function MC(T) for a function F(X) over an
%    interval [A,B] is defined as
%
%      MC(T) = max | F(X+DX) - F(X) | for 0 <= DX <= T, and X and X+DX in [A,B].
%
%    If we define the simpler function
%
%      MC1(DX) = max | F(X+DX) - F(X) | for 0 <= DX, and X and X+DX in [A,B],
%
%    then we can evaluate MC(T) as
%
%      MC(T) = max ( MC1(DX) ), for DX <= T.
%
%    This function estimates the MC1 function based on a discrete
%    set of data at N equally spaced points in the interval [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the left and right endpoints of the interval.
%
%    Input, integer N, the number of equally spaced sample points.
%
%    Input, function F(x), a handle to the function.
%
%    Output, real MC1(N), the maximum difference |F(X+DX)-F(X)| estimated at
%    DX = 0 exactly, H exactly, 2*H exactly, ..., (N-1)*H exactly.
%
  if ( nargin < 1 )
    a = 0.0;
  end

  if ( nargin < 2 )
    b = a + 1.0;
  end

  if ( nargin < 3 )
    n = 501;
  end
%
%  Set the evaluation points, and evaluate the function.
%
  x = linspace ( a, b, n );
  fx = f ( x );
%
%  Compute the maximum difference with a separation of 0*H, 1*H, 2*H,
%  ..., (N-1)*H.
%
  mc1 = zeros ( n, 1 );
  for i = 0 : n - 1
    mc1(i+1) = max ( abs ( fx(1+i:n) - fx(1:n-i) ) );
  end

  return
end


text:
function h = ldivide(f,g)
%.\   Pointwise CHEBFUN left divide.
%   F.\G returns a CHEBFUN that represents the function G(x)/F(x).
%
% See also RDIVIDE, MLDIVIDE.

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

% Call RDIVIDE():
h = rdivide(g, f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function DCM = spm_dcm_load(P, save_mem)
% Load a cell array of DCM filenames into a subjects x models cell array
% FORMAT DCM = spm_dcm_load(P)
%
% P        - a DCM's filename or
%          - a GCM's filename (which contains a cell array of DCM files) or
%          - a cell array of DCM filenames or
%          - a character array of DCM filenames
% save_mem - (Optional) if true, only loads priors, posteriors and F for
%            models 2-N
%
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% Peter Zeidman
% $Id: spm_dcm_load.m 7479 2018-11-09 14:17:33Z peter $

if ischar(P)
    if size(P,1) > 1
        % A character array was provided
        P = cellstr(P);
    else
        % A GCM / DCM filename was provided
        x = load(P);
        if isfield(x,'GCM')
            P = x.GCM;
        elseif isfield(x,'DCM')
            P = {x.DCM};
        else
            error('Unknown DCM file type');
        end
    end
end

if nargin < 2
    save_mem = false;
end

DCM = cell(size(P));

fprintf('Loading DCMs...');
for s = 1:size(P,1)
    for m = 1:size(P,2)

        if isstruct(P{s,m})
            % DCM already a structure - leave unchanged
            dcm = P{s,m};
        else
            % Load DCM from filename
            try
                model = load(P{s,m});
                dcm   = model.DCM;
            catch
               fprintf('File: %s\n',P{s,m});
               error('Failed to load model for subject %d model %d', s, m);
            end
        end

        if save_mem && (m>1)
            % Onlt import basic structures
            DCM{s,m} = struct();
            if isfield(dcm,'a'),  DCM{s,m}.a  = dcm.a; end
            if isfield(dcm,'b'),  DCM{s,m}.b  = dcm.b; end
            if isfield(dcm,'c'),  DCM{s,m}.c  = dcm.c; end
            if isfield(dcm,'d'),  DCM{s,m}.d  = dcm.d; end
            if isfield(dcm,'F'),  DCM{s,m}.F  = dcm.F; end
            if isfield(dcm,'Ep'), DCM{s,m}.Ep = dcm.Ep; end
            if isfield(dcm,'Cp'), DCM{s,m}.Cp = dcm.Cp; end
            if isfield(dcm,'options'), DCM{s,m}.options = dcm.options; end
            if isfield(dcm,'M')
                if isfield(dcm.M,'pE'), DCM{s,m}.M.pE = dcm.M.pE; end
                if isfield(dcm.M,'pC'), DCM{s,m}.M.pC = dcm.M.pC; end
            end
        else
            % Import the whole DCM
            DCM{s,m} = dcm;
        end

    end
end
fprintf('Done\n');


text:
%% SOS Testing
%% LP Solve Example SOS1
clc
f = [-1 -1 -3 -2 -2]';
A = [-1 -1 1 1 0;
     1 0 1 -3 0];
b = [30;30];
ub = [40;1;inf;inf;1];

sostype = '12';
sosind = {[1:2]' [3:5]'};
soswt = {[1:2]' [3:5]'};

% [x,fval,ef,stat] = cplexmilp(f,A,b,[],[],sos,sosind,soswt,[],ub)

% SOS 1 OPTI Version
opts = optiset('solver','cbc');

prob = optiprob('f',f,'ineq',A,b,'ub',ub,'sos',sostype,sosind,soswt,'options',opts)

sos = [];
sos.type = sostype;
sos.index = sosind;
sos.weight = soswt;
prob2 = optiprob('f',f,'ineq',A,b,'ub',ub,'sos',sos,'options',opts)

Opt = opti(prob,opts)

[x,fval,ef,stat] = solve(Opt)

%% LP Solve Example SOS2
clc
f = [-1 -1 -3 -2 -2]';
A = sparse([-1 -1 1 1 0;
             1 0 1 -3 0]);
b = [30;30];

lb = zeros(5,1);
ub = [40;1;inf;inf;1];

sos = '1';
sosind = (1:5)';
soswt = (1:5)';

% [x,fval,ef,stat] = cplexmilp(f,A,b,[],[],sos,sosind,soswt,lb,ub,'IIIII');
% x
% stat

% SOS 2 OPTI Version
opts = optiset('solver','cbc');

prob = optiprob('f',f,'ineq',A,b,'bounds',lb,ub,'sos',sos,sosind,soswt,'options',opts)

Opt = opti(prob,opts)

[x,fval,ef,stat] = solve(Opt)

% [x,fval,ef,stat] = solve(Opt);
% x
% fval
% stat

%%
clc
% Objective (f'x)
f = [-1 -1 -3 -2 -2]';

% Row Constraints (rl <= A*x <= ru)
A = sparse([-1 -1  1  1  0;     %sparse A
             1  0  1 -3  0]);
rl = [-Inf;-Inf];
ru = [30;30];

% Bounds (lb <= x <= ub)
lb = [0;0;0;0;0];
ub = [40;1;Inf;Inf;1];

% SOS
sos_type = '1';
sos_index = [1 2 3 4 5]';
sos_weight = [1 2 3 4 5]';

% Setup Options
opts = optiset('solver','cbc'); %CBC is a SOS constraint solver

% Build OPTI Object
Opt = opti('f',f,'lin',A,rl,ru,'bounds',lb,ub,'sos',sos_type,sos_index,...
           sos_weight,'options',opts)

% Solve Problem
[x,fval] = solve(Opt)

%%
% SOS Structure
sos = [];
sos.type = '12'; %each character represents a SOS set
sos.index = {[1 2]' [3:5]'};
sos.weight = {[1 2]' [1:3]'};

% Build OPTI Object
Opt = opti('f',f,'lin',A,rl,ru,'bounds',lb,ub,'sos',sos,'options',opts)

% Solve Problem
[x,fval] = solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function subGraphs=getAllSubGraphs(model)
% getAllSubGraphs
%   Get all metabolic subgraphs in a model. Two metabolites
%   are connected if they share a reaction.
%
%   Input:
%   model         a model structure
%
%   Output:
%   subGraphs     a boolean matrix where the rows correspond to the metabolites
%                 and the columns to which subgraph they are assigned to. The
%                 columns are ordered so that larger subgraphs come first
%
%   Usage: subGraphs=getAllSubGraphs(model)

%Generate the connectivity graph. Metabolites are connected through
%reactions. This is not a bipartite graph with the reactions.
G=model.S*model.S';
G(G~=0)=1;

%Keeps track of which mets have been assigned to a subgraph
isAssigned=false(numel(model.mets),1);

%Allocate space for subgraphs, initially one graph for each met
subGraphs=false(numel(model.mets),numel(model.mets));

%Main loop continues until all mets have been assigned to a subgraph
counter=1;
while ~all(isAssigned)
    currentSG=false(numel(model.mets),1);
    %Find the first non-assigned metabolite and assign it to the current SG
    I=find(isAssigned==false,1);
    currentSG(I)=true;
    isAssigned(I)=true;
    %Then iteratively add all mets that are connected to each other, until
    %no more such mets can be found
    while true
        J=sum(G(:,currentSG),2);

        %Get the new mets for this SG. Also assign them to the current SG
        newOnes=J~=0 & currentSG==false;
        isAssigned(newOnes)=true;
        currentSG(newOnes)=true;

        %If there are no new mets, then abort
        if ~any(newOnes)
            subGraphs(currentSG,counter)=true;
            counter=counter+1;
            break;
        end
    end
end
subGraphs=subGraphs(:,1:counter);

[~, I]=sort(sum(subGraphs),'descend');

subGraphs=subGraphs(:,I);

%Also remove empty subgraphs (can happen when metabolites are never used)
subGraphs(:,sum(subGraphs)==0)=[];
end

author: SysBioChalmers

repo: RAVEN


text:
function nMC=grMinEdgeCover(E)
% Function nMC=grMinEdgeCover(E) solve the minimal edge cover problem.
% Input parameter:
%   E(m,2) or (m,3) - the edges of graph and their weight;
%     1st and 2nd elements of each row is numbers of vertexes;
%     3rd elements of each row is weight of edge;
%     m - number of edges.
%     If we set the array E(m,2), then all weights is 1.
% Output parameter:
%   nMC - the list of the numbers of edges included
%     in the minimal (weighted) edge cover.
% Uses the reduction to integer LP-problem.
% Required the Optimization Toolbox v.3.0.1 or over.
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

% ============= Input data validation ==================
if nargin<1,
  error('There are no input data!')
end
[m,n,E] = grValidation(E); % E data validation

% ============= Parameters of integer LP problem ==========
A=zeros(n,m); % for incidence matrix
A(E(:,1:2)+repmat(([1:m]'-1)*n,1,2))=1; % we fill the incidence matrix
options=optimset('bintprog'); % the default options
options.Display='off'; % we change the output

% ============= We solve the integer LP problem ==========
xmin=bintprog(E(:,3),-A,-ones(n,1),[],[],[],options);
nMC=find(round(xmin)); % the answer - numbers of edges
return

author: HuangCongQing


text:
function [x, featDescription] = getBoundaryFeatures(ind, occ, pb, pb2)
% [x, featDescription] = getBoundaryFeatures(ind, occ, pb, pb2)

if islogical(ind)
    ind = find(ind);
end

featDescription = cell(1, 7);
x = zeros(numel(ind), 7);

[imh, imw, no] = size(occ.po_all);
npix = imh*imw;
%n = zeros(numel(ind), 1);
for f = 1:4
    ind2 = ind + npix*(f-1);
    x(:, f) = occ.po_all(ind2);
    featDescription{f} = ['pocc' num2str(f)];
    %n = n + double(x(:, f)>0);
end
f = f+1;
x(:, f) = max(x(:, 1:f-1), [], 2); % occ.po(ind);
featDescription{f} = 'pocc_max';
%f = f+1;
%x(:, f) = n;
if exist('pb', 'var') && ~isempty(pb)
    f = f+1;
    pb_soft = max(pb.pb_soft, [], 3);
    x(:, f) = pb_soft(ind);
    featDescription{f} = 'pb1';
end
if exist('pb2', 'var') && ~isempty(pb2)
    f = f+1;
    x(:, f) = pb2.pb_soft(ind);
    featDescription{f} = 'pb2';
end

author: Cloud-CV

repo: object-proposals


text:
% Sample Code to test the Tensor Completion program :

clear;

n1 = 50;
n2 = 50;
n3 = 50;
r = 5;
tol   = [];
nitr  = [];
ninit = [];

% generate random low-rank orthogonal tensor
U01 = randn(n1,r);U02 = randn(n2,r);U03 = randn(n3,r);
[U1 temp] = qr(U01);[U2 temp] = qr(U02);[U3 temp] = qr(U03);
U1=U1(:,1:r);U2=U2(:,1:r);U3=U3(:,1:r);
T = zeros(n1,n2,n3);
for i=1:n3
    T(:,:,i) = U1*diag(U3(i,:))*U2';
end

% sample entries
p = 2*(r^0.5*log(n1*n2*n3))/sqrt(n1*n2*n3);
E = ceil(rand(n1,n2,n3)-1+p);
TE = T.*E;

% exact tensor completion
fprintf(1,'exact tensor completion using alternating minimization\n');
[V1 V2 V3 S dist] = TenALS(TE, E, r, ninit, nitr, tol);
rmse=0;
for i3=1:n3
   A1 = U1;
   A2 = U2*diag(U3(i3,:));
   B1 = V1;
   B2 = V2*diag(V3(i3,:).*S(:)');
   rmse = rmse + trace(A1'*A1*A2'*A2) + trace(B1'*B1*B2'*B2) -2*trace(B1'*A1*A2'*B2) ;
end
fprintf(1,'root mean squared error = %e\n',sqrt(rmse/r));

% noisy tensor completion
fprintf(1,'\nnoisy tensor completion using alternating minimization\n');
[V1 V2 V3 S dist] = TenALS(TE+(0.0001/sqrt(n1*n2*n3)*randn(n1,n2,n3).*E), E, r, ninit, nitr, tol);
rmse=0;
for i3=1:n3
   A1 = U1;
   A2 = U2*diag(U3(i3,:));
   B1 = V1;
   B2 = V2*diag(V3(i3,:).*S(:)');
   rmse = rmse + trace(A1'*A1*A2'*A2) + trace(B1'*B1*B2'*B2) -2*trace(B1'*A1*A2'*B2) ;
end
fprintf(1,'root mean squared error = %e\n',sqrt(rmse/r));

author: andrewssobral


text:
Iref=(im2double(imread('test.jpg')));
tau=[0.2,0.3,0.1,0.2,200,100];
% tau=[0,0,0,0,50,100];
It=Iref;
for i=1:3
It(:,:,i)=warpImg(Iref(:,:,i),tau);
end

imshow([Iref,It])

tic
[ImTrans,tau] = align(It,Iref, zeros(6,1),2);
toc;
imshow([It,Iref,ImTrans,10*abs(ImTrans-Iref)])

author: csjcai

repo: RealSR

save_path: github-repos/MATLAB/csjcai-RealSR


text:
classdef RegressionLoss < dagnn.Loss

  properties
    lossType = 1.
  end

  methods
    function outputs = forward(obj, inputs, params)
        X = inputs{1};
        c = inputs{2};
        assert(numel(X) == numel(c));
        c = reshape(c,[1,1,size(c)]);
        switch obj.lossType
            case 'CE'
                X = vl_nnsigmoid(X);
                Y = sum(squeeze(-c.*log(X) - (1-c).*log(1.001-X)));
            case 'MSE'
                Y = sum(squeeze((X - c).^2));
            case 'Huber'
                delta = 1.0;
                a = abs(squeeze(X - c));
                y1 = 0.5*sum(a(a<=delta).^2);
                y2 = sum((a(a>delta)-0.5*delta)*delta);
                Y = y1 + y2;
        end
        outputs{1} = Y;
        if obj.ignoreAverage, return; end;
        n = obj.numAveraged ;
        m = n + size(inputs{2},2);
        obj.average = (n * obj.average + gather(outputs{1})) / m ;
        obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
        X = inputs{1};
        c = inputs{2};
        assert(numel(X) == numel(c));
        c = reshape(c,[1,1,size(c)]);
        switch obj.lossType
            case 'CE'
                X = vl_nnsigmoid(X);
                Y = X - c;
            case 'MSE'
                Y = X - c;
            case 'Huber'
                delta = 1;
                a = X - c;
                Y = X - c;
                Y(a>delta) = delta;
                Y(a<-delta) = - delta;
        end
        derInputs = {Y, []};
        derParams = {};
    end

    function obj = RegressionLoss(varargin)
      obj.load(varargin) ;
    end
  end
end

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
classdef AscendingNodeTermCondition < LatitudeTermCondition
    %AscendingNodeTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties

    end

    methods
        function obj = AscendingNodeTermCondition()
            obj = obj@LatitudeTermCondition(0.0);
            obj.direction = 1.0;
        end

        function name = getName(obj)
            name = sprintf('Ascending Node');
        end

        function params = getTermCondUiStruct(obj)
            params = getTermCondUiStruct@LatitudeTermCondition(obj);

            params.paramName = 'Ascending Node';
            params.paramUnit = '';
            params.useParam = 'off';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = [];
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function optVar = getNewOptVar(obj)
            optVar = AbstractOptimizationVariable.empty(0,1);
        end

        function optVar = getExistingOptVar(obj)
            optVar = AbstractOptimizationVariable.empty(0,1);
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = AscendingNodeTermCondition();
        end
    end
end


text:
classdef PTKAirwayForContext < PTKPlugin
    % PTKAirwayForContext. Plugin for returning the main bronchi serving the lung regions for a given context
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKAirwayForContext returns a set of bronchi which correspond to the
    %     specified context
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Bronchus for<br>context'
        ToolTip = 'Plugin for returning the main bronchi serving the lung regions for a given context'
        Category = 'Analysis'

        Context = PTKContextSet.Any
        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'DoNothing'
        HidePluginInDisplay = true
        FlattenPreviewImage = false
        PTKVersion = '2'
        ButtonWidth = 6
        ButtonHeight = 1
        GeneratePreview = false
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, context, reporting)
            bronchi = [];
            switch context
                case {PTKContext.Lungs, PTKContext.LungROI, PTKContext.OriginalImage}
                    airways = dataset.GetResult('PTKAirwayCentreline');
                    bronchi =  airways.AirwayCentrelineTree;
                case {PTKContext.LeftLung}
                    airways = dataset.GetResult('PTKAirwaysLabelledByLobe');
                    bronchi = airways.StartBranches.Left;
                case {PTKContext.RightLung}
                    airways = dataset.GetResult('PTKAirwaysLabelledByLobe');
                    bronchi = airways.StartBranches.Right;
                case {PTKContext.LeftUpperLobe}
                    airways = dataset.GetResult('PTKAirwaysLabelledByLobe');
                    bronchi = airways.StartBranches.LeftUpper;
                case {PTKContext.LeftLowerLobe}
                    airways = dataset.GetResult('PTKAirwaysLabelledByLobe');
                    bronchi = airways.StartBranches.LeftLower;
                case {PTKContext.RightUpperLobe}
                    airways = dataset.GetResult('PTKAirwaysLabelledByLobe');
                    bronchi = airways.StartBranches.RightUpper;
                case {PTKContext.RightLowerLobe}
                    airways = dataset.GetResult('PTKAirwaysLabelledByLobe');
                    bronchi = airways.StartBranches.RightLower;
                case {PTKContext.RightMiddleLobe}
                    airways = dataset.GetResult('PTKAirwaysLabelledByLobe');
                    bronchi = airways.StartBranches.RightMid;
                case {PTKContext.R_AP, PTKContext.R_P, PTKContext.R_AN, PTKContext.R_L, ...
                        PTKContext.R_M, PTKContext.R_S, PTKContext.R_MB, PTKContext.R_AB, PTKContext.R_LB, ...
                        PTKContext.R_PB, PTKContext.L_APP, PTKContext.L_APP2, PTKContext.L_AN, PTKContext.L_SL, ...
                        PTKContext.L_IL, PTKContext.L_S, PTKContext.L_AMB, PTKContext.L_LB, PTKContext.L_PB}
                    airways = dataset.GetResult('PTKSegmentsByNearestBronchus');
                    bronchi = PTKAirwayForContext.FindSegmentalBronchus(airways.AirwaysBySegment.Trachea, context);
            end
            for branch = bronchi
                branch.GenerateBranchParameters;
            end

            results = [];
            results.AirwayForContext = bronchi;
        end

        function bronchus = FindSegmentalBronchus(airways, context)
            segment_label = uint8(PTKPulmonarySegmentLabels.(char(context)));
            airways_to_do = CoreStack(airways);
            while ~airways_to_do.IsEmpty
                next_airways = airways_to_do.Pop;
                if next_airways.SegmentIndex == segment_label
                    bronchus = next_airways;
                    return;
                end
                airways_to_do.Push(next_airways.Children);
            end

            bronchus = [];
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function SO3VF = cross(SO3VF1, SO3VF2, varargin)
% pointwise cross product
%
% Syntax
%   SO3VF = cross(SO3VF1, SO3VF2)
%   SO3VF = cross(v, SO3VF2)
%   SO3VF = cross(SO3VF1, v)
%
% Input
%   SO3VF1, SO3VF2 - @SO3VectorField
%   v - @vector3d
%
% Output
%   SO3VF - @SO3VectorField
%

if isa(SO3VF2, 'vector3d') && length(SO3VF2)==1
  SO3VF = SO3VectorFieldHandle(@(rot) cross(SO3VF1.eval(rot),SO3VF2),SO3VF1.CS,SO3VF1.SS);
  if SO3VF2.antipodal
    SO3VF = abs(SO3VF);
  end
  return
end

if isa(SO3VF1, 'vector3d')
  SO3VF = -cross(SO3VF2,SO3VF1);
  return
end

ensureCompatibleSymmetries(SO3VF1,SO3VF2)
SO3VF = SO3VectorFieldHandle(@(rot) cross(SO3VF1.eval(rot),SO3VF2.eval(rot)),SO3VF1.CS,SO3VF1.SS);

end

author: mtex-toolbox

repo: mtex

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


text:
function p = genpathGIT(d)
%GENPATH Generate recursive toolbox path.
%   P = GENPATH returns a new path string by adding all the subdirectories
%   of MATLABROOT/toolbox, including empty subdirectories.
%
%   P = GENPATH(D) returns a path string starting in D, plus, recursively,
%   all the subdirectories of D, including empty subdirectories.
%
%   NOTE 1: GENPATH will not exactly recreate the original MATLAB path.
%
%   NOTE 2: GENPATH only includes subdirectories allowed on the MATLAB
%   path.
%
%   See also PATH, ADDPATH, RMPATH, SAVEPATH.

%   Copyright 1984-2006 The MathWorks, Inc.
%   $Revision: 1.13.4.5 $ $Date: 2008/06/24 17:12:06 $
%------------------------------------------------------------------------------

if nargin==0,
  p = genpath(fullfile(matlabroot,'toolbox'));
  if length(p) > 1, p(end) = []; end % Remove trailing pathsep
  return
end

% initialise variables
classsep = '@';  % qualifier for overloaded class directories
packagesep = '+';  % qualifier for overloaded package directories
p = '';           % path to be returned

% Generate path based on given root directory
files = dir(d);
if isempty(files)
  return
end

% Add d to the path even if it is empty.
p = [p d pathsep];

% set logical vector for subdirectory entries in d
isdir = logical(cat(1,files.isdir));
%
% Recursively descend through directories which are neither
% private nor "class" directories.
%
dirs = files(isdir); % select only directory entries from the current listing

for i=1:length(dirs)
   dirname = dirs(i).name;
   if    ~strcmp( dirname,'.')          && ...
         ~strcmp( dirname,'..')         && ...
         ~strcmp( dirname,'.git')       && ...
         ~strncmp( dirname,classsep,1) && ...
         ~strncmp( dirname,packagesep,1) && ...
         ~strcmp( dirname,'private')
      p = [p genpathGIT(fullfile(d,dirname))]; % recursive calling of this function.
   end
end

%------------------------------------------------------------------------------

author: ZJULearning

repo: MatlabFunc


text:
function value = p25_autonomous ( )

%*****************************************************************************80
%
%% P25_AUTONOMOUS reports whether problem P25 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function scanList = selectScans(vw,title,~)
% scanList = selectScans(vw,[title],[scanDir]);
%
%   Gather a list of scans available in Inplane/TSeries
%   and query the user for a sub-selection.
%
%   An alternate function chooseScans uses numScans(view)
%   to determine the number of scans to choose from
%   Use selectScans if you will be analyzing the tSeries.
%   Use chooseScans, if your code does not depend on the
%   presence/absence of the tSeries files.
%
% Input:
%  scanDir (optional): specify the dir that contains
%    'Scan*' subdir. Default is tSeriesDir(view).
% Output:
%  scanList: list of scans selected.
%
% 4/16/99  dbr Initial code
% 3/30/2001, djh, added optional title string
% 10/21/2004, MA, added optional scanDir
% 11/3/04, JL, make scanDir really optional

if ~exist('title','var'), title = 'Choose scans'; end;
%if ~exist('scanDir','var'), scanDir = tSeriesDir(vw,0); end;
%We will no longer use scanDir

scanNum = viewGet(vw,'N Scans');

scanList = 1:scanNum;

scanNames = cell(length(scanList),1);
for i=scanList
    scanNames{i} = num2str(i);
end

% Which scans to analyze?
iSel = buttondlg(title, scanNames);
scanList = scanList(find(iSel));

return;

author: vistalab

repo: vistasoft


text:
%Function to place dose into the running CERR plan.
%Just edit the strings below.
%JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

% [CTUniform3D, CTUniformInfoS] = getUniformizedCTScan;
function planC = addNewDose(dose3D, planC)

register = 'UniformCT';  %Currently only option supported.  Dose has the same shape as the uniformized CT scan.
doseError = [];
fractionGroupID = 'CERR test';
doseEdition = 'CERR test';
description = 'Test PB distribution.'
overWrite = 'no';  %Overwrite the last CERR dose?
dose2CERR(dose3D,doseError,fractionGroupID,doseEdition,description,register,[],overWrite)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout=bss_decomp_tvfilt(varargin)

% decompose an estimated source into target/interference/noise/artefacts components, assuming the admissible distortion is a time-varying filter.
%
% Usage:
%
% [s_target,e_interf[,e_noise],e_artif]=bss_decomp_tvfilt(se,index,S[,N],tvshape,tvstep,L)
%
% Input:
%   - se: row vector of length T containing the estimated source,
%   - index: points which component of S se has to be compared to,
%   - S: n x T matrix containing the original sources,
%   - N: m x T matrix containing the noise on the obseravtions (if any).
%   - tvshape : row vector of length V at most T containing the shape of the elementary
%     allowed time variations of the filter coefficients
%   - tvstep  : hop size (in number of samples) between two consecutive
%     variations of the filter coefficients
%   - L: the number of lags
%
% Output:
%   - s_target: row vector of length T containing the target source(s)
%   contribution,
%   - e_interf: row vector of length T containing the interferences
%   contribution,
%   - e_noise: row vector of length T containing the noise contribution (if
%   any),
%   - e_artif: row vector of length T containing the artifacts
%   contribution.
%
% Developers:  - Cedric Fevotte (fevotte@tsi.enst.fr) - Emmanuel Vincent
% (emmanuel.vincent@irisa.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

se=varargin{1}; index=varargin{2}; S=varargin{3};

switch nargin
    case 6
        N=[]; tvshape = varargin{4}; tvstep = varargin{5}; L = varargin{6};
    case 7
        N=varargin{4}; tvshape = varargin{5}; tvstep = varargin{6}; L = varargin{7};
    otherwise
        disp('Wrong number of arguments.')
end

[ne,Te]=size(se);
[n,T]=size(S);

%%%%%%%%%% WARNINGS %%%%%%%%%%%%%
switch isempty(N)
    case 1
        if n>T || ne>Te, disp('Watch out: signals must be in rows.'), return; end
        if ne~=1, disp('Watch out: se must contain only one row.'), return; end
        if T~=Te, disp('Watch out: se and S have different lengths.'), return; end
    case 0
        [m,Tm]=size(N);
        if n>T || ne>Te || m>Tm, disp('Watch out: signals must be in rows.'), return; end
        if ne~=1, disp('Watch out: se must contain only one row.'), return; end
        if T~=Te, disp('Watch out: S and Se have different lengths.'), return; end
        if T~=Tm, disp('Watch out: N, S and Se have different lengths.'), return; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create the space of target source(s)
target_space = bss_make_lags(S(index,:),L);
% Create the space of sources
sources_space= bss_make_lags(S,L);
% Create the noise space
noise_space  = bss_make_lags(N,L);

%%% Target source(s) contribution %%%
s_target = bss_tvproj(se,target_space,tvshape,tvstep);

%%% Interferences contribution %%%
P_S_se = bss_tvproj(se,sources_space,tvshape,tvstep);
e_interf = P_S_se - s_target;

if isempty(noise_space) % No noise
    %%% Artifacts contribution %%%
    e_artif= se - P_S_se;

    %%% Output %%%
    varargout{1}=s_target;
    varargout{2}=e_interf;
    varargout{3}=e_artif;

else % Noise
    %%% Noise contribution %%%
    P_SN_se= bss_tvproj(se,[sources_space;noise_space],tvshape,tvstep);
    e_noise=P_SN_se-P_S_se;

    %%% Artifacts contribution %%%
    e_artif=se-P_SN_se;

    %%% Output %%%
    varargout{1}=s_target;
    varargout{2}=e_interf;
    varargout{3}=e_noise;
    varargout{4}=e_artif;
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function a = i4vec_permute ( n, a, p )

%*****************************************************************************80
%
%% I4VEC_PERMUTE permutes an I4VEC in place.
%
%  Discussion:
%
%    This routine permutes an array of integer "objects", but the same
%    logic can be used to permute an array of objects of any arithmetic
%    type, or an array of objects of any complexity.  The only temporary
%    storage required is enough to store a single object.  The number
%    of data movements made is N + the number of cycles of order 2 or more,
%    which is never more than N + N/2.
%
%  Example:
%
%    Input:
%
%      N = 5
%      P = (   2,   4,   5,   1,   3 )
%      A = (   1,   2,   3,   4,   5 )
%
%    Output:
%
%      A    = (   2,   4,   5,   1,   3 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of objects.
%
%    Input, integer A(N), the array to be permuted.
%
%    Input, integer P(N), the permutation.  P(I) = J means
%    that the I-th element of the output array should be the J-th
%    element of the input array.  Pmust be a legal permutation
%    of the integers from 1 to N, otherwise the algorithm will
%    fail catastrophically.
%
%    Output, integer A(N), the permuted array.
%
  base = 1;

  missing = perm_check ( n, p, base );

  if ( missing ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4VEC_PERMUTE - Fatal error!\n' );
    fprintf ( 1, '  The input array does not represent\n' );
    fprintf ( 1, '  a proper permutation.  In particular, the\n' );
    fprintf ( 1, '  array is missing the value %d\n', missing );
    error ( 'I4VEC_PERMUTE - Fatal error!' );
  end
%
%  Search for the next element of the permutation that has not been used.
%
  for istart = 1 : n

    if ( p(istart) < 0 )

      continue;

    elseif ( p(istart) == istart )

      p(istart) = -p(istart);
      continue;

    else

      a_temp = a(istart);
      iget = istart;
%
%  Copy the new value into the vacated entry.
%
      while ( 1 )

        iput = iget;
        iget = p(iget);

        p(iput) = -p(iput);

        if ( iget < 1 | n < iget )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'I4VEC_PERMUTE - Fatal error!\n' );
          fprintf ( 1, '  A permutation index is out of range.\n' );
          fprintf ( 1, '  P(%d) = %d\n', iput, iget );
          error ( 'I4VEC_PERMUTE - Fatal error!' );
        end

        if ( iget == istart )
          a(iput) = a_temp;
          break;
        end

        a(iput) = a(iget);

      end

    end

  end
%
%  Restore the signs of the entries.
%
% p(1:n) = -p(1:n);

  return
end