text:
function [Out1,Out2]=qsynth(Action,In1,In2)
% Dick Benson
% Design single terminated elliptic passive LC low pass filters.
% Copyright 2001-2013 The MathWorks, Inc.
% [L,C] = qsynth('getLC');

    if nargin ==0
       Action = 'init';
    end;

    if ~strcmp(Action,'init')
      Hqsynth_ = get(gcf,'userdata');
    end

    % handle index definitions
    i_f1      = 1;
    i_topology   = 2;
    i_order   = 3;
    i_synthpb = 4;
    i_ripl    = 5;
    i_ripv    = 6;
    i_stopl   = 7;
    i_stopv   = 8;
    i_ax1     = 9;
    i_pl1     = 10;
    i_pl2     = 11;
    i_fcl     = 12;
    i_fcv     = 13;
    i_lv      = 14;
    i_cv      = 15;
    i_lcl     = 16;
    i_testpb  = 17;

    if strcmp(Action,'init')

        Hqsynth_(i_f1)= figure('Color',[0 0 0] ,'Name','Qsynth: Single Terminated Ladder Elliptic Filter',...
                          'Position',[20 20 640 440],...
                          'Resize','off',...
                          'Color',[0 0 1],...
                          'NumberTitle','off','visible','on');
                          % ###### userdata contains vector of handles if desired

        Hqsynth_(i_topology)=uicontrol('Style','popup','visible','on',...
                                    'String','RIN=1,RL=inf|RIN=0,RL=1',...
                                    'Position',[10 400 105 20],...
                                    'BackgroundColor',[0 0 0],...
                                    'ForeGroundColor',[0 1 1],...
                                    'Value',1,...
                                    'callback','qsynth(''clear'')',...
                                    'HorizontalAlignment','Left');

        Hqsynth_(i_order)=uicontrol('Style','popup','visible','on',...
                                'String','Order:3|Order:5|Order:7|Order:9|Order:11',...
                                'Position',[10 380 105 20],...
                                'BackgroundColor',[0 0 0],...
                                'ForeGroundColor',[0 1 1],...
                                'Value',4,...
                                'HorizontalAlignment','Left');

        Hqsynth_(i_ripl) =  uicontrol('Style','text','visible','on',...
                                'String','Ripple (dB):',...
                                'Position',[10,355,60,18],...
                                'BackGroundColor',[1 1 1]*0.8,...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','left');

         Hqsynth_(i_ripv) = uicontrol('Style','edit','visible','on',...
                                'Position',[70,355,45,18],...
                                'String','1',...
                                'BackgroundColor',[1 1 0],...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','center',...
                                'Max',1);

         Hqsynth_(i_stopl) =  uicontrol('Style','text','visible','on',...
                                'String','Stop (dB):',...
                                'Position',[10,330,60,18],...
                                'BackGroundColor',[1 1 1]*0.8,...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','left');

         Hqsynth_(i_stopv) = uicontrol('Style','edit','visible','on',...
                                'Position',[70,330,45,18],...
                                'String','60',...
                                'BackgroundColor',[1 1 0],...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','center',...
                                'Max',1);

         Hqsynth_(i_fcl) =  uicontrol('Style','text','visible','on',...
                                'String','Cutoff (Hz):',...
                                'Position',[10,305,60,18],...
                                'BackGroundColor',[1 1 1]*0.8,...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','left');

         Hqsynth_(i_fcv) = uicontrol('Style','edit','visible','on',...
                                'Position',[70,305,45,18],...
                                'String','1.0e6',...
                                'BackgroundColor',[1 1 0],...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','center',...
                                'Max',1);

          Hqsynth_(i_lcl) =  uicontrol('Style','text','visible','on',...
                                'String','              L                         C',...
                                'Position',[10,280,165,18],...
                                'BackGroundColor',[1 1 1]*0.8,...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','left');

        Hqsynth_(i_lv)  =  uicontrol('Style','edit','visible','on',...
                                'Position',[10,35,85,240],...   % was 80
                                'String','',...
                                'BackgroundColor',[1 1 0],...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','left',...
                                'Max',14);

        Hqsynth_(i_cv)  =  uicontrol('Style','edit','visible','on',...
                                'Position',[100,35,85,240],...  % was 95,35,80
                                'String','',...
                                'BackgroundColor',[1 1 0],...
                                'ForeGroundColor',[0 0 0],...
                                'HorizontalAlignment','left',...
                                'Max',14);

        Hqsynth_(i_ax1)=axes('Units','pixels','Position',[240,80,370,300],...
                                'Box','on',...
                                'visible','on',...
                                'NextPlot','add',...
                                'DrawMode','fast',...
                                'Color',[0 0 0],...
                                'TickDir','out',...
                                'YlimMode','auto',...
                                'XlimMode','auto',...
                                'Xcolor',[1 1 1],...
                                'Ycolor',[1 1 1],...
                                'FontSize',12);

        title('Singly Terminated Elliptic Filter ','color',[1 1 1],'fontsize',9);
        xlabel('Frequency in Hz','fontsize',9);
        ylabel('Magnitude in dB','fontsize',9)

        Hqsynth_(i_pl1) = plot(Hqsynth_(i_ax1),0,0,'clipping','on',...
                                  'Color',[0 1 0],...
                                  'erasemode','xor','visible','on');

        Hqsynth_(i_pl2) = plot(Hqsynth_(i_ax1),0,0,'clipping','on',...
                                  'Color',[1 0 0],...
                                  'erasemode','xor','visible','on');

        Hqsynth_(i_synthpb) = uicontrol('Style','Pushbutton',...
                              'Position',[10,10,105,20],...
                              'String','Synthesize',...
                              'Callback','qsynth(''synth'');');

        Hqsynth_(i_testpb) = uicontrol('Style','Pushbutton',...
                              'Position',[120,10,105,20],...
                              'String','Test',...
                              'Callback','qsynth(''test'');');

        set(Hqsynth_(i_f1),'userdata',Hqsynth_);
        zoom on

    elseif strcmp(Action,'synth')
       fmtstr = '%8.5g';

       n  = get(Hqsynth_(i_order),'value')*2+1;          % only odd orders supported
       rp = str2num(get(Hqsynth_(i_ripv),'string'));
       rs = str2num(get(Hqsynth_(i_stopv),'string'));
       wn = 2*pi*str2num(get(Hqsynth_(i_fcv),'string'));

       [b,a]=ellip(n,rp,rs,wn,'s');
       save lpf_coefs b a
       w = 0:(wn/25):10*wn;
       h = freqs(b,a,w);

       set(Hqsynth_(i_ax1),'Ylim',[-rs-20,5],'Xlim',[0,10*wn/(2*pi)]);
       set(Hqsynth_(i_pl1),'xdata',w/(2*pi),'ydata',20*log10(h));

       topology = get(Hqsynth_(i_topology),'value');
       if topology ==1
          %  Rin = 1, Rl = inf , Vsource with 1 ohm Z0, no load
          % must be odd order
          m     = b(2:length(b));
          zeros = roots(m);
          m2 = [];
          n2 = [];
          for i=1:length(a)
              if rem(i,2)
                 n2=[n2,a(i),0];
              else
                 m2=[m2,a(i),0];
              end;
          end;
          m2=m2(1:(length(m2)-1));

          % remove admittance pole @ s=inf
          % z=m2/n2  y = n2/m2 yp = n2/m2-c1s

          C(1) = n2(1)/m2(1);  % shunt capacitor
          n2 = n2 - conv(m2,[C(1),0]);
          % now it gets tricky ....
          for i = 1:(n-1)/2
              % need an xmission zero @ w=wx
              wx = j*abs(zeros(2*i)); % pick an xmission zero....
              % do a partial removal of the impedance pole @ inf
              % first evaluate impedance @ w=wx
              zmag = polyval(m2,wx)/polyval(n2,wx);
              L(2*i) = zmag/wx;

              if i==1
                 n2=n2(3:length(n2));
              end;
              % this partially removes the pole
              m2 = m2-conv(n2,[L(2*i),0]);

              % now remaining admittance has a pole at w=wx
              % remove it by shunting a series LC

              % first must determine L & C values
              [q,r1]= deconv(m2,[1 0 abs(wx)^2]);  % r should be zero
              k    = n2(1);
              n2p  = n2/k;
              n2p  = n2p(1:length(n2p)-1);
              L(2*i+1)=1/(k*(polyval(n2p,wx)/polyval(q,wx)));
              C(2*i+1)= 1/(L(2*i+1)*(abs(wx)^2));

              % then remove the addmittance.... not obvious!
              [n2pp,r2] = deconv(( n2p-q/(k*L(2*i+1)) )*k,[1 0 1/(C(2*i+1)*L(2*i+1))]);
              n2 = conv(n2pp,[1 0]);
              m2 = q;
          end;
          L=L';
          C=C';

       elseif topology==2
          % Rin=0, Rl=1 (Vsource input, 1 ohm load)
          %
          m     = b(2:length(b)); % drop leading 0.0 term
          zeros = roots(m);       %
          m2 = [];
          n2 = [];
          for i=1:length(a)
              if rem(i,2)
                 n2=[n2,a(i),0];
              else
                 m2=[m2,a(i),0];
              end;
          end;

          m2=m2(1:(length(m2)-1));
          for i = 1:(n-1)/2
              % need an xmission zero @ w=wx
              wx = j*abs(zeros(2*i)); % pick an xmission zero....

              % do a partial removal of the impedance pole @ inf
              % first evaluate impedance @ w=wx
              zmag = polyval(n2,wx)/polyval(m2,wx);
              L(2*i-1) = zmag/wx;
              % this partially removes the pole

              n2 = n2-conv(m2,[L(2*i-1),0]);
              % now remaining admittance has a pole at w=wx

              % remove it by shunting a series LC
              % first must determine L & C values
              [q,r1] = deconv(n2,[1 0 abs(wx)^2]);  % r should be zero
              k = polyval(m2,wx)/polyval(conv(q,[1 0]),wx);
              L(2*i) = 1/k;
              C(2*i) = 1/(L(2*i)*(abs(wx)^2));
              [p,r2] = deconv( m2-k*conv(q,[1 0])   ,  [1 0 abs(wx)^2]);
              m2=p;
              n2=q;
          end;
          L(n)= n2(1)/m2;     % final series L .... amazing .... absolutly amazing!
          L = fliplr(L)';     % need to reverse this network
          C = [0,fliplr(C)]'; %

       end;
       set(Hqsynth_(i_ripl),'userdata',L);
       set(Hqsynth_(i_ripv),'userdata',C);

       s='';
       for i=1:(length(L)-1)
           s=[s,sprintf([fmtstr,'\n'],L(i))];
       end;
       s=[s,sprintf(fmtstr,L(length(L)))];
       set(Hqsynth_(i_lv),'string',s);
       s='';
       for i=1:(length(C)-1)
           s=[s,sprintf([fmtstr,'\n'],C(i))];
       end;
       s=[s,sprintf(fmtstr,C(length(C)))];
       set(Hqsynth_(i_cv),'string',s);

    elseif strcmp(Action,'getLC');
      % Out1 = get(Hqsynth_(i_ripl),'userdata'); % inductors
      % Out2 = get(Hqsynth_(i_ripv),'userdata'); % capacitors
      % huh? ,  should get it from strings in edit fields
       caps = get(Hqsynth_(i_cv),'string');
       inds = get(Hqsynth_(i_lv),'string');

       for i=1:length(caps(:,1))
          Out2(i) = str2num(caps(i,:));
       end;

       for i=1:length(inds(:,1))
          Out1(i) = str2num(inds(i,:));
       end;

    elseif strcmp(Action,'test')
       topology = get(Hqsynth_(i_topology),'value');   % network topology
       wn = 2*pi*str2num(get(Hqsynth_(i_fcv),'string'));
       w  = 0:(wn/25):10*wn;
       [L,C]=qsynth('getLC');
       if topology ==1
          % test the synthesis
          % begin with source resistor
          [a,b,c,d,q]=pabcd([],[],[],[],[],'rs',1);
          % then the shunt capacitor
          [a,b,c,d,q]=pabcd(a,b,c,d,q,'cp',C(1));
          for i = 1:(length(L)-1)/2
              % series L
              [a,b,c,d,q]=pabcd(a,b,c,d,q,'ls',L((i*2)));
              % shunt series LC
              [a,b,c,d,q]=pabcd(a,b,c,d,q,'lcs',L((i*2)+1),C((i*2)+1));
          end;

       elseif topology ==2
          a = []; b = []; c = []; d = []; q = [];
          for i = 1:(length(L)-1)/2
              % series L
              [a,b,c,d,q]=pabcd(a,b,c,d,q,'ls',L((i*2-1)));
              % shunt series LC
              [a,b,c,d,q]=pabcd(a,b,c,d,q,'lcs',L((i*2)),C((i*2)));
          end;
          % finish with a series L
          [a,b,c,d,q]=pabcd(a,b,c,d,q,'ls',L(length(L)));
          % and a 1 ohm shunt R load resistor
          [a,b,c,d,q]=pabcd(a,b,c,d,q,'rp',1);
       end;

       % [mag,phase]   =   bode(q,a,w);
       h = freqs(q,a,w);
       set(Hqsynth_(i_pl2),'xdata',w/(2*pi),'ydata',20*log10(h));

    elseif strcmp(Action,'clear')
       set(Hqsynth_(i_ripl),'userdata',[]);
       set(Hqsynth_(i_ripv),'userdata',[]);
       set(Hqsynth_(i_cv),'string','');
       set(Hqsynth_(i_lv),'string','');
       set(Hqsynth_(i_pl2),'xdata',[],'ydata',[]);
       set(Hqsynth_(i_pl1),'xdata',[],'ydata',[]);

    else
       disp([Action,' not reconized in qsynth'])
    end;

% end qsynth

function [aout,bout,cout,dout,qout]=pabcd(ain,bin,cin,din,qin,element,val1,val2,val3)
% [aout,bout,cout,dout,qout]=pabcd(ain,bin,cin,din,qin,element,val1,val2)
% Uses chain matrix parameters to cascades abcd with denominator q (in)
% with an element that has val1 (val2 for LC).
% Useful for computing functions & parameters of ladder networks.
% element code   description
%     rs       = series resistor
%     lcs      = shunt series lc val1=l val2=c
%     ls       = series l
%     cp       = shunt c
%     rp       = shunt r
%     lcp      = series chunt lc val1=l val2=c val3 = r (in series with l for loss)
%
%     a*e2-b*i2 = e1
%     c*e2-d*i2 = i1
%
% Dick Benson
%

      if strcmp(element,'rs')
           % rs  = series resistor
           a = 1;
           b = val1;
           c = 0;
           d = 1;
           q = 1;
      elseif strcmp(element,'lcs')
           % lcs = shunt series lc
           % l=val1 c=val2
           q = [val1*val2, 0 , 1];
           a = q;
           b = 0;
           c = [val2 0];
           d = q;

      elseif strcmp(element,'lcp')
          % l=val1 c=val2
           q = [val1*val2 val2*val3 1];
           a = q;
           b = [val1 val3];
           c = 0;
           d = q;
      elseif strcmp(element,'ls')
           % ls = series l
           a = 1;
           b = [val1 0];
           c = 0;
           d = 1;
           q = 1;
      elseif strcmp(element,'cp')
           % cp = shunt c
           q = 1;
           a = 1;
           b = 0;
           c = [val1 0];
           d = 1;
      elseif strcmp(element,'rp')
           % rp = shunt r
           q = 1;
           a = 1;
           b = 0;
           c = 1/val1;
           d = 1;
      else
         disp([element,' not supported in pabcd '])
      end;

      if isempty(ain)
         aout=a;
         bout=b;
         cout=c;
         dout=d;
         qout=q;
      else
         aout = padd(conv(ain,a),conv(bin,c));
         bout = padd(conv(ain,b),conv(bin,d));
         cout = padd(conv(cin,a),conv(din,c));
         dout = padd(conv(cin,b),conv(din,d));
         qout = conv(qin,q);
      end;

% end function

function pout = padd(p1,p2)
    lp1 = length(p1);
    lp2 = length(p2);
    if lp1==lp2
       pout = p1+p2;
    elseif lp2>lp1
       pout = p2;
       pout((lp2-lp1+1):lp2)=pout((lp2-lp1+1):lp2)+p1;
    else
       pout = p1;
       pout((lp1-lp2+1):lp1)=pout((lp1-lp2+1):lp1)+p2;
    end;
% end padd


text:
function S = size(v)
% SIZE  Size of the expansion coefficients of a BALLFUNV
%   S = SIZE(V) returns a 3 by 3 matrix whose vertical entries represent
%   the sizes of Vx, Vy and Vz
%
%   See also BALLFUN/SIZE.

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

% Return the list [m,n,p] corresponding to the size of the ballfun
% functions in f

if isempty( v )
    S = [];
    return
end

V = v.comp;

S = [ size( V{1} );
      size( V{2} );
      size( V{3} ); ];
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function r = corr ( n, x, m )

%*****************************************************************************80
%
%% CORR computes the sample correlation of a signal sample.
%
%  Discussion:
%
%    The sample correlation is defined, for 0 <= i < N, as
%
%      R(i) = 1/N * sum ( 0 <= j <= N - 1 - i ) X(i+j) * X(j)
%
%    The sample correlation is an estimate of the correlation function.
%
%    It is usually the case that the signal X is assumed to
%    have zero mean.  Here, we compute the mean and adjust the
%    calculation accordingly:
%
%      R(i) = 1/N * sum ( 0 <= j <= N - 1 - i )
%        ( X(i+j) - Xbar ) * ( X(j) - Xbar )
%
%    Experience suggests that only the first 5 or 10 percent of
%    the lags are statistically reliable, so that one might choose
%    M = N / 20 or M = N / 10, for instance.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Sophocles Orfanidis,
%    Introduction to Signal Processing,
%    Prentice-Hall, 1995,
%    ISBN: 0-13-209172-0,
%    LC: TK5102.5.O246.
%
%  Parameters:
%
%    Input, integer N, the number of equally spaced signal
%    samples.
%
%    Input, real X(N), the signal samples.
%
%    Input, integer M, the maximum lag to consider.
%    0 <= M < N.
%
%    Output, real R(1:M+1), the sample correlations.
%
  r(1:m+1) = 0.0;

  xbar = sum ( x(1:n) ) / n;

  for i = 0 : m
    for j = 0 : n - i - 1
      r(i+1) = r(i+1) + ( x(i+j+1) - xbar ) * ( x(j+1) - xbar );
    end
  end

  r(1:m+1) = r(1:m+1) / n;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
port = int16(2000);
client = py.carla.Client('localhost', port);
client.set_timeout(10.0);
world = client.get_world();
blueprint_library = world.get_blueprint_library();
actor_list = py.list(world.get_actors());

id = 0;

for i=1:length(actor_list)
    k = actor_list{i}.type_id;

    if k(1:7) == "vehicle"
        id = i - 1;
    end
end

if (id == 0)
    fprintf('ERROR: Could not find any vehicle!\n');
    fprintf('The program will now exit!\n');
    return;
end

carla_is_running = true;
yolo_detection = true;

% Sensor 1
blueprint = world.get_blueprint_library().find('sensor.camera.rgb');
blueprint.set_attribute('image_size_x', '1920')
blueprint.set_attribute('image_size_y', '1080')
% blueprint.set_attribute('sensor_tick', '.1');
transform = py.carla.Transform(py.carla.Location(pyargs('x',0.8, 'z',1.7)));
sensor = world.spawn_actor(blueprint, transform, pyargs('attach_to',actor_list{id}));

pyModule = py.importlib.import_module('senFunc');
py.senFunc.getData(sensor);

pause(.5);
% data = uint8(py.array.array('d',py.numpy.nditer(ml)));

% Load YOLO model
if exist('yoloml','var') ~= 1
    disp('loading modified network')
    load('yoloml.mat')
end

% Display Initial Frame
image = py.getattr(pyModule, 'array');
% image = py.numpy.ascontiguousarray(image);
imMat = uint8(image);

% Set Yolo parameters
probThresh = 0.175;
iouThresh = 0.5;

myImageSize = size(imMat);

classLabels = ["aeroplane",	"bicycle",	"bird"	,"boat",	"bottle"	,"bus"	,"car",...
"cat",	"chair"	,"cow"	,"diningtable"	,"dog"	,"horse",	"motorbike",	"person",	"pottedplant",...
"sheep",	"sofa",	"train",	"tvmonitor"];

imageHandle = imshow(imMat);

% Detect After how many frames
totalFramesElapsed = 0;
predictEvery = 2;

while carla_is_running
%     pause(1);
%     clc;

     image = py.getattr(pyModule, 'array');
     imMat = uint8(image);

     % Detect Cars
%      [bboxes,scores] = detect(detector,imMat);
%      imMat = insertObjectAnnotation(imMat,'rectangle',bboxes,1);

     if yolo_detection
         redImage = single(imresize(imMat,[448 448]))/255;

         if ~rem(totalFramesElapsed, predictEvery)
            out = predict(yoloml,redImage,'ExecutionEnvironment','gpu');
            % fprintf('%d\n', totalFramesElapsed);
         end
         class = out(1:980);
         boxProbs = out(981:1078);
         boxDims = out(1079:1470);

         outArray = zeros(7,7,30);

         for j = 0:6
             for i = 0:6
                 outArray(i+1,j+1,1:20) = class(i*20*7+j*20+1:i*20*7+j*20+20);
                 outArray(i+1,j+1,21:22) = boxProbs(i*2*7+j*2+1:i*2*7+j*2+2);
                 outArray(i+1,j+1,23:30) = boxDims(i*8*7+j*8+1:i*8*7+j*8+8);
             end
         end

         [cellProb, cellIndex] = max(outArray(:,:,21:22),[],3);
         contain = max(outArray(:,:,21:22),[],3)>probThresh;

         [classMax, classMaxIndex] = max(outArray(:,:,1:20),[],3);

         boxes = [];

         counter = 0;
         for i = 1:7
             for j = 1:7
                 if contain(i,j) == 1
                     counter = counter+1;

                     % Bounding box center relative to cell
                     x = outArray(i,j,22+1+(cellIndex(i,j)-1)*4);
                     y = outArray(i,j,22+2+(cellIndex(i,j)-1)*4);

                     % Yolo outputs the square root of the width and height of the
                     % bounding boxes (subtle detail in paper that took me forver to realize).
                     % Relative to image size.
                     w = (outArray(i,j,22+3+(cellIndex(i,j)-1)*4))^2;
                     h = (outArray(i,j,22+4+(cellIndex(i,j)-1)*4))^2;

                     %absolute values scaled to image size
                     wS = w*448;
                     hS = h*448;
                     xS = (j-1)*448/7+x*448/7-wS/2;
                     yS = (i-1)*448/7+y*448/7-hS/2;

                     % this array will be used for drawing bounding boxes in Matlab
                     boxes(counter).coords = [xS yS wS hS];

                    %save cell indices in the structure
                     boxes(counter).cellIndex = [i,j];

                    %save classIndex to structure
                     boxes(counter).classIndex = classMaxIndex(i,j);

                    % save cell proability to structure
                     boxes(counter).cellProb = cellProb(i,j);

                    % put in a switch for non max which we will use later
                     boxes(counter).nonMax = 1;
                 end
             end
         end

         for i = 1:length(boxes)
             for j = i+1:length(boxes)
                 %calculate intersection over union (can also use bboxOverlapRatio
                 %with proper toolbox
                 intersect = rectint(boxes(i).coords,boxes(j).coords);
                 union = boxes(i).coords(3)*boxes(i).coords(4)+boxes(j).coords(3)*boxes(j).coords(4)-intersect;
                 iou(i,j) = intersect/union;
                 if boxes(i).classIndex == boxes(j).classIndex && iou(i,j) > iouThresh
                     [value(i), dropIndex(i)] = min([boxes(i).cellProb boxes(j).cellProb]);
                     if dropIndex(i) == 1
                         boxes(i).nonMax=0;
                     elseif dropIndex(i) == 2
                         boxes(j).nonMax=0;
                     end
                 end
             end
         end

         %imagesc(imMat);

         for i = 1:length(boxes)
             if boxes(i).nonMax == 1
                 textStr = convertStringsToChars(classLabels(boxes(i).classIndex));

                 if sum(textStr == [ "bicycle", "car", "motorbike",	"person"])
                     position = [(boxes(i).cellIndex(2)-1)*myImageSize(2)/7 (boxes(i).cellIndex(1)-1)*myImageSize(1)/7];
                     %text(position(1),position(2),textStr,'Color',[0 1 0],'fontWeight','bold','fontSize',12);
                     transformedCoor = boxes(i).coords .* [myImageSize(2) myImageSize(1) myImageSize(2) myImageSize(1)] / 448;

                     imMat = insertText(imMat, transformedCoor([1, 2]), textStr, 'AnchorPoint','LeftBottom', 'FontSize', 36);
                     imMat = insertShape(imMat, 'Rectangle', transformedCoor, 'LineWidth',5);

                     % rectangle('Position',transformedCoor, 'EdgeColor','green','LineWidth',2);
                 end
             end
         end
     end

     % Write directly to Cdata. Consideribly faster than imshow
     set(imageHandle,'Cdata',imMat);
     drawnow;

     totalFramesElapsed = totalFramesElapsed + 1;

     pause(0.001);

      % Directly write to Cdata, Noticibly faster than just imshow
     % rectangle('Position',[1 1 30 30], 'EdgeColor','green','LineWidth',2);
     %imagesc(imMat);
     %pause(.05);
     %py.cv2.imshow('Carla', image);
     % pause(0.05);
%      data = uint8(py.array.array('d',py.numpy.nditer(image)));
%
%      data = reshape(data, 4, []);
%
%      imMat(:,:,1) = reshape(data(3,:), 1920, 1080)';
%      imMat(:,:,2) = reshape(data(2,:), 1920, 1080)';
%      imMat(:,:,3) = reshape(data(1,:), 1920, 1080)';
%
%      imshow(imMat);

%      for n=0:1080*1920-1
%         i = fix(n/1920) + 1;
%         j = rem(n,1920) + 1;
%
%         imMat(i, j, 3) = data(4*n + 1);
%         imMat(i, j, 2) = data(4*n + 2);
%         imMat(i, j, 1) = data(4*n + 3);
%      end

% %     currentLocation = actor_list{id}.get_location();
% %     fprintf('Location: [m]\n')
% %     fprintf('x: %.2f \t y: %.2f \t z: %.2f \n\n',currentLocation.x, currentLocation.y, currentLocation.z);
% %
% %     fprintf('Velocity: [m/s]\n')
% %     currentVelocity = actor_list{id}.get_velocity();
% %     fprintf('x: %.2f \t y: %.2f \t z: %.2f \n\n',currentVelocity.x, currentVelocity.y, currentVelocity.z);
% %
% %     fprintf('Acceleration: [m/s^2]\n')
% %     currentAcceleration = actor_list{id}.get_acceleration();
% %     fprintf('x: %.2f \t y: %.2f \t z: %.2f \n\n',currentAcceleration.x, currentAcceleration.y, currentAcceleration.z);
% %
% %     omega = char (hex2dec ( '03C9' ));
% %     fprintf('Angular Velocity: [%s]\n', omega);
% %     currentAngVelocity = actor_list{id}.get_angular_velocity();
% %     fprintf('x: %.2f \t y: %.2f \t z: %.2f \n\n',currentAngVelocity.x, currentAngVelocity.y, currentAngVelocity.z);
%
end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
function yearly_tmax_correction(filenameout)

% this program is design to adjust the daily Tmax generated by WG using
% new yearly Tmax series produced by FFT.

% load the yearly averaged Tmax after FFT
load('Pnew_tmax');
tmax_FFT=Pnew_tmax';

% load WG generated daily Tmax
load(filenameout);
[n,m]=size(gTmax);
gTmax=gTmax';
tmax_WG=reshape(gTmax,[],1);

n=length(tmax_FFT); % years
m=length(tmax_WG);  % days

% calculate yearly Tmax generated by WG, namely (Y)
j=1;
Z=zeros(n,1);
for i=1:365:m
    Z(j,1)=mean(tmax_WG(i:i+365-1,1));
    j=j+1;
end

% calculate the ratio of yearly Tmax between those after FFT and
% generated by weather generator
for i=1:n
    tmax_ratio(i,1)=tmax_FFT(i,1)-Z(i,1);
end

% extend the yearly Tmax ratio to daily scale,the data in each
% year are the same
tmax_extent=zeros(m,1);
j=1;
for i=1:365:m
    tmax_extent(i:i+365-1,1)=tmax_ratio(j,1);
    j=j+1;
end

% adjust the daily Tmax generated by WG using above ratios
tmax_adjust=zeros(size(tmax_WG));
for i=1:m
    tmax_adjust(i,1)=tmax_WG(i,1)+tmax_extent(i,1);
end
yearly_corrected_tmax=tmax_adjust;
save('yearly_corrected_tmax','yearly_corrected_tmax')

author: Sable

repo: mcbench-benchmarks


text:
function F = in_fread_eeg(sFile, sfid, iEpoch, SamplesBounds)
% IN_FREAD_EEG:  Read an epoch from a Neuroscan .eeg file (list of epochs).
%
% USAGE:  F = in_fread_eeg(sFile, sfid, iEpoch, SamplesBounds)
%         F = in_fread_eeg(sFile, sfid, iEpoch)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, 2009-2011

fileSamples = round(sFile.prop.times .* sFile.prop.sfreq);
% Parse inputs
if (nargin < 4) || isempty(SamplesBounds)
    SamplesBounds = fileSamples;
% Check start and stop samples
elseif (SamplesBounds(1) < fileSamples(1)) || (SamplesBounds(1) > SamplesBounds(2)) || (SamplesBounds(2) > fileSamples(2))
    error('Invalid samples range.');
end

% Get some information on the file
nChannels = sFile.header.epochs(iEpoch).datasize(1);
% Position cursor in file to read this data block
startSample = SamplesBounds(1) - fileSamples(1);
pos = double(sFile.header.epochs(iEpoch).datapos) + startSample * nChannels * double(sFile.header.data.bytes_per_samp);
fseek(sfid, double(pos), 'bof');
% Read [nChannels, nSamples]
F = fread(sfid, [nChannels, SamplesBounds(2) - SamplesBounds(1) + 1], sFile.header.data.dataformat);
% Calibrate data
F = neuroscan_apply_calibration(F, sFile.header);

author: brainstorm-tools


text:
function alpha = meshDihedralAngles(vertices, edges, faces)
%MESHDIHEDRALANGLES Dihedral at edges of a polyhedal mesh.
%
%   ALPHA = meshDihedralAngles(V, E, F)
%   where V, E and F represent vertices, edges and faces of a mesh,
%   computes the dihedral angle between the two adjacent faces of each edge
%   in the mesh. ALPHA is a column array with as many rows as the number of
%   edges. The i-th element of ALPHA corresponds to the i-th edge.
%
%   Note: the function assumes that the faces are correctly oriented. The
%   face vertices should be indexed counter-clockwise when considering the
%   supporting plane of the face, with the outer normal oriented outwards
%   of the mesh.
%
%   Example
%   [v, e, f] = createCube;
%   rad2deg(meshDihedralAngles(v, e, f))
%   ans =
%       90
%       90
%       90
%       90
%       90
%       90
%       90
%       90
%       90
%       90
%       90
%       90
%
%   See also
%   meshes3d, polyhedronMeanBreadth, trimeshMeanBreadth, dihedralAngle, meshEdgeFaces
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-10-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% compute normal of each face
normals = meshFaceNormals(vertices, faces);

% indices of faces adjacent to each edge
edgeFaces = meshEdgeFaces(vertices, edges, faces);

% allocate memory for resulting angles
Ne = size(edges, 1);
alpha = zeros(Ne, 1);

% iterate over edges
for i = 1:Ne
    % indices of adjacent faces
    indFace1 = edgeFaces(i, 1);
    indFace2 = edgeFaces(i, 2);

    % normal vector of adjacent faces
    normal1 = normals(indFace1, :);
    normal2 = normals(indFace2, :);

    % compute dihedral angle of two vectors
    alpha(i) = vectorAngle3d(normal1, normal2);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [error_x,error_y] = ...
         navierpost_bc(viscosity,aez,fezx,fezy,elerrorx,elerrory,xy,ev,ebound);
%navierpost_q1p0_bc   postprocesses Poisson error estimator
%   [error_x,error_y] = ...
%   navierpost_q1p0_bc(viscosity,aez,fezx,fezy,elerrorx,elerrory,xy,ev,ebound);;
%   input
%          viscosity            viscosity parameter
%          aez                  elementwise Poisson problem matrices
%          fezx,fezy            elementwise rhs vectors
%          elerrorx, elerrory   elementwise error estimate (without BC imposition)
%          xy                   vertex coordinate vector
%          ev                   element mapping matrix
%          ebound               element edge boundary matrix
%   output
%          error_x, error_y     component elementwise error estimate
%
%   calls function localbc_xycd
%   IFISS function: DJS; 11 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      x=xy(:,1); y=xy(:,2);
      nel=length(ev(:,1));
      lev=[ev,ev(:,1)];
	  error_x=elerrorx;  error_y=elerrory;
%
% recompute contributions from elements with Dirichlet boundaries
      nbde=length(ebound(:,1));
      ebdy = zeros(nel,1);
      edge = zeros(nel,1);
% isolate boundary elements
      for el = 1:nbde
      ee = ebound(el,1);
      ebdy(ee) = ebdy(ee)+1; edge(ee)=ebound(el,2);
      end
%
% two edge elements
      k2=find(ebdy==2);
      nel2b=length(k2);
% loop over two edge elements
      for el = 1:nel2b
      el2e=k2(el);
      kk=find(ebound(:,1) == el2e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
	  ae=squeeze(aez(el2e,1:5,1:5));
      fex=fezx(el2e,:)'; fey=fezy(el2e,:)';
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el2e,:)); yl=y(lev(el2e,:));
	  [bae,fex,fey] = localbc_xycd(viscosity,ae,fex,fey,edges,xl,yl);
% solve local problems
      errx=viscosity*bae\fex;  erry=viscosity*bae\fey;
	  error_x(el2e,1) = errx'*ae*errx; error_y(el2e,1) = erry'*ae*erry;
	  end
% end of element loop
%
% one edge elements
      k1=find(ebdy==1);
      nel1b=length(k1);
% loop over one edge elements
      for el = 1:nel1b
      el1e=k1(el);
      kk=find(ebound(:,1) == el1e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
      fex=fezx(el1e,:)'; fey=fezy(el1e,:)';
	  ae=squeeze(aez(el1e,1:5,1:5));
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el1e,:)); yl=y(lev(el1e,:));
	  [bae,fex,fey] = localbc_xycd(viscosity,ae,fex,fey,edges,xl,yl);
% solve local problems
      errx=viscosity*bae\fex;  erry=viscosity*bae\fey;
	  error_x(el1e,1) = errx'*ae*errx; error_y(el1e,1) = erry'*ae*erry;
      end
% end of element loop
%
      err_x = sqrt(sum(error_x)); error_x = sqrt(error_x);
	  err_y = sqrt(sum(error_y)); error_y = sqrt(error_y);
      fprintf('estimated velocity error (in energy):  (%10.6e,%10.6e) \n',err_x,err_y)
 return


text:
function currentVal=getImageRotate(view)
%
% currentVal=getImageRotate(view)
%
% Returns the value of the image rotate slider in radians
if isequal(view.name, 'hidden')
	% instead of a handle to the slider, the ImageRotate field will
	% directly store the rotations for hidden views:
	currentVal = view.ui.ImageRotate;
else
	% non-hidden view: view.ui.ImageRotate should be a slider struct
	currentVal = get(view.ui.ImageRotate.sliderHandle, 'Value');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [X p ei ej] = chrobak_payne_straight_line_drawing(A,varargin)
% CHROBAK_PAYNE_STRAIGHT_LINE_DRAWING Draw planar graphs with straight lines
%
% X = chrobak_payne_straight_line_drawing(A) generates coordinates for each
% vertex such that a planar graph A can be drawn without any edge
% crossings.  This function reports an error if A is not planar.
%
% [X,p,ei,ej] = ... returns additional information.  p is a
% canonical planar ordering of the vertices, and [ei ej] are additional
% edges required to make A a maximal_planar graph.
%
% ... = chrobak_payne_straight_line_drawing(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   options.is_maximal: is A already a maximal planar graph [{0} | 1]
%
% Note: Be careful with is_maximal=1 and nocheck=1.  If the graph is not
% maximal, then the call will crash Matlab.
%
% Example:
%   [A,xy] = grid_graph(6,5);
%   X = chrobak_payne_straight_line_drawing(A);
%   gplot(A,X,'.-'); hold on; gplot(A,xy*20,'r.-'); hold off
%   % it's still planar, but not obviously a grid!

% David Gleich
% Copyright, Stanford University, 2008

%% History
%  2007-10-06: Initial coding
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

options = struct('is_maximal',0);
options = merge_options(options,varargin{:});
if check
    check_matlab_bgl(A,struct('sym',1));
    if options.is_maximal,
        [i j] = make_maximal_planar(A);
        if ~isempty(i), error('matlab_bgl:checkFailed',...
            'The graph was not a maximal planar but is_maximal was set.'); end
    end
end

[ei ej p X] = planar_drawing_mex(A,options.is_maximal,0);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function R = Ehn_GF(img)

mi = min(img(:));
ma = max(img(:));
img1 = (img-mi)/(ma-mi)*255+0;

log_img = log(img1+1);

[m,n] = size(img1);
r = floor(0.04*max(m,n));
eps = 0.1;
% ***using fast guided filter
base = 255*fastguidedfilter_md(img1/255, img1/255, r, eps^2, max(1,r/4));
% base = 255*guidedfilter(img1/255, img1/255, r, eps^2);
log_base = log(base+1);

log_detail = log_img - log_base;

target_contrast = log(4);
compressfactor = target_contrast/(max(log_base(:))-min(log_base(:)));
log_absolute_scale = (1-compressfactor)*max(log_base(:));
log_output = compressfactor*log_base + log_detail + log_absolute_scale;
R = exp(log_output);
R = min(R, 255);
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [vert,econ,tria] = cfmtri2(vert,econ)
%CFMTRI2 compute a conforming 2-simplex Delaunay triangulat-
%ion in the two-dimensional plane.
%   [VERT,CONN,TRIA]=CFMTRI2(VERT,CONN) computes the confor-
%   ming Delaunay trianguation, given the points VERT, and
%   edge constraints CONN. New points are inserted to bisect
%   edge constraints until all are recovered. VERT is a
%   V-by-2 array of XY coordinates to be triangulated, TRIA
%   is a T-by-3 array of vertex indexing, with each row
%   defining a triangle, such that VERT(TRIA(II,1),:),
%   VERT(TRIA(II,2),:) and VERT(TRIA(II,3),:) are the coord-
%   inates of the II-TH triangle. CONN is a C-by-2 array of
%   constraining edges, where each row defines an edge, as
%   per the TRIA array.
%
%   See also DELTRI2, DELAUNAYN

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 07/07/2017

%---------------------------------------------- basic checks
    if ( ~isnumeric(vert) || ...
         ~isnumeric(econ) )
        error('cfmtri2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(vert) ~= +2 || ndims(econ) ~= +2)
        error('cfmtri2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    if (size(vert,2)~= +2 || size(econ,2)~= +2)
        error('cfmtri2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

%-- the DELAUNAYN routine is *not* well-behaved numerically,
%-- so explicitly re-scale the problem about [-1,-1; +1,+1].
    vmax = max(vert,[],1) ;
    vmin = min(vert,[],1) ;

    vdel = vmax - vmin;
    vdel = mean(vdel) ;
    vdel = vdel * +.5 ;

    vmid = vmax + vmin;
    vmid = vmid * +.5 ;

    vert = vert - vmid;
    vert = vert / vdel;

%-- keep bisecting edge constraints until they are all reco-
%-- vered!
    while (true)

    %----------------- un-constrained delaunay triangulation
        tria = delaunay2(vert) ;

        nv = size(vert,+1);
        nt = size(tria,+1);

    %----------------------------- build non-unique edge-set
        ee = zeros(nt*3,2);
        ee((1:nt)+nt*0,:) = tria(:,[1,2]);
        ee((1:nt)+nt*1,:) = tria(:,[2,3]);
        ee((1:nt)+nt*2,:) = tria(:,[3,1]);

    %----------------- find constraints within tria-edge set
       [in] = setset2(econ,ee) ;

    %----------------------------- done when have contraints
        if (all(in)), break; end

    %----------------------------- un-recovered edge centres
        vm = vert(econ(~in,1),:) ...
           + vert(econ(~in,2),:) ;
        vm = vm * +.5 ;

    %----------------------------- un-recovered edge indexes
        ev = nv+(1:size(vm,1))';
        en = [econ(~in,+1), ev;
              econ(~in,+2), ev];

    %----------------------------- push new vert/edge arrays
        vert = [vert( :,:); vm];
        econ = [econ(in,:); en];

    end

%--------------------------------- undo geomertic re-scaling
    vert = vert * vdel ;
    vert = vert + vmid ;

end

function [tria] = delaunay2(vert)
%DELAUNAY2 thin wrapper for DELAUNAYN, so that we can have a
%   more efficient version in OCTAVE...

    isoctave = exist( ...
        'OCTAVE_VERSION','builtin')>+0;

    if (isoctave)

    %-- call QHULL and then filter zero-volume simplexes via
    %-- vectorised area comparisons.

    %-- note silliness re. EVAL, so that MATLAB doesn't com-
    %-- plain re. OCTAVE '__' names.

        tria = eval( ...
          '__delaunayn__(vert)') ;

        ab = vert(tria(:,2),:) ...
           - vert(tria(:,1),:) ;
        ac = vert(tria(:,3),:) ...
           - vert(tria(:,1),:) ;

        aa = ab(:,1).* ac(:,2) ...
           - ab(:,2).* ac(:,1) ;

        lb = sumsq(ab,2) ;
        lc = sumsq(ac,2) ;

        ll = max (lb,lc) ;

        keep = abs(aa) >= ll * eps^.8 ;

        tria = tria(keep,:);

    else

    %-- the default call in MATLAB seems to be fast enough!!

        tria = ...
        delaunay (vert(:,1),vert(:,2));

    end

end

author: dengwirda


text:
function [ears,ear_opp,flops,flop_opp] = find_ears(F)
  % FIND_EARS  Find all ears (faces with two boundary edges) in a given mesh
  %
  % [ears,ear_opp,flops,flop_opp] = find_ears(F)
  %
  % Inputs:
  %   F  #F by 3 list of triangle mesh indices
  % Outputs:
  %   ears  #ears list of indices into F of ears
  %   ear_opp  #ears list of indices indicating which edge is non-boundary
  %     (connecting to flops)
  %   flops  #ears list of indices into F of faces incident on ears
  %   flop_opp  #ears list of indices indicating which edge is non-boundary
  %     (connecting to ears)
  %
  % Known limitation: (V,F) must be manifold near ears for flops and flop_opp
  % outputs.
  %

  if nargout<=2
    % Works on non-manifold meshes...
    [~,B] = on_boundary(F);
    ears = find(sum(B,2) == 2);
    [~,ear_opp] = min(B(ears,:),[],2);
  else
    % Must be manifold (actually only ears need to be)
    [Fp, Fi] = triangle_triangle_adjacency(F);
    % Definition: an ear has 2 boundary **edges**
    ears = find(sum(Fp==-1,2)==2);
    [flops,ear_opp] = max(Fp(ears,:),[],2);
    flop_opp = Fi(sub2ind(size(Fp),ears,ear_opp));
  end
end

author: alecjacobson

repo: gptoolbox


text:
function [rollback, ev, nefv] = pne_detect_events(reg_ev, nefv, cefv, step)
%PNE_DETECT_EVENTS  Detect events from event function values
%   [ROLLBACK, CRITICAL_EVENTS, NEFV] = ...
%       PNE_DETECT_EVENTS(REG_EV, NEFV, CEFV, STEP)
%
%   Determines whether any events have been detected based on the
%   corresponding "current" and "next" event function values. If any events
%   are detected, the relevant information is returned in CRITICAL_EVENTS,
%   along with a ROLLBACK flag indicating whether to roll back the proposed
%   next step.
%
%   Inputs:
%       REG_EV : struct containing info about registered event fcns
%       NEFV : cell array of Next Event Function Values
%       CEFV : cell array of Current Event Function Values
%       STEP : current step size
%
%   Outputs:
%       ROLLBACK : flag indicating whether any event has requested a
%           rollback step, if ROLLBACK is true, then CRITICAL_EVENTS will
%           have length 1
%       CRITICAL_EVENTS : struct array containing information about any
%           detected events, with fields:
%           eidx        : event index, index in REG_EV, 0 if no event detected
%           name        : name of event function, empty if none detected
%           zero        : 1 if zero has been detected, 0 otherwise
%                           (i.e. interval detected or no event detected)
%           idx         : index(es) of critical elements in event function
%           step_scale  : linearly interpolated estimate of scaling factor
%                         for current step size required to reach event zero
%           log         : 1 log the event in the results, 0 don't log the event
%                         (set to 1 for zero events, 0 otherwise, can be
%                           modified by callbacks)
%           msg         : event message, set to something generic like
%                           'ZERO detected for TARGET_LAM event' or
%                           'INTERVAL detected for QLIM(3) event', but intended
%                         to be changed/updated by callbacks
%       NEFV : cell array of Next Event Function Values
%           (possibly updated to be exactly zero to avoid re-detection of
%           an interval following a zero detection)

%   MP-Opt-Model
%   Copyright (c) 2016-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Shrirang Abhyankar, Argonne National Laboratory
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% initialize result variables
rollback = 0;
ev = struct( ...            %% initialize struct for storing critical events
    'eidx', 0, ...
    'zero', 0, ...
    'step_scale', 1, ...
    'log', 0, ...
    'name', '', ...
    'idx', 0, ...
    'msg', '' ...
);

%% other initialization
i = 1;              %% index into ev struct
nef = length(nefv); %% number of event functions

%% detect events, first look for event intervals for events requesting rollback
for eidx = 1:nef
    if ~reg_ev(eidx).locate     %% if event doesn't request rollback to locate zero
        continue;               %%   skip to next event
    end

    %% current and previous event function signs
    c_sign = sign(nefv{eidx});
    p_sign = sign(cefv{eidx});

    %% if there's been a sign change and we aren't within event tolerance ...
    idx = find( abs(c_sign) == 1 & c_sign == -p_sign & ...
                abs(nefv{eidx}) > reg_ev(eidx).tol  );
    if ~isempty(idx)
        if step == 0    %% if it's a "repeat" step
            %% (e.g. after warmstart with possible fcn change)
            %% ... make this one the critical one and call it a ZERO event
            ev.eidx = eidx;
            ev.zero = 1;
            ev.step_scale = 1;
            ev.log = 1;
            ev.name = reg_ev(eidx).name;
            ev.idx = idx;
            ev.msg = 'ZERO (BIFURCATION)';
            i = i + 1;
            break;
        else
            %% ... compute step size scaling factors and find index of smallest one
            [step_scale, j] = ...
                min(cefv{eidx}(idx) ./ (cefv{eidx}(idx) - nefv{eidx}(idx)) );

            %% if it's smaller than the current critical one ...
            if step_scale < ev.step_scale
                %% ... make this one the critical one
                ev.eidx = eidx;
                ev.zero = 0;
                ev.step_scale = step_scale;
                ev.log = 0;
                ev.name = reg_ev(eidx).name;
                ev.idx = idx(j);
                ev.msg = 'INTERVAL';
                rollback = 1;   %% signal that a rollback event has been detected
            end
        end
    end
end

%% if no rollback events were detected
if rollback == 0
    %% search for event zeros
    for eidx = 1:nef
        %% if there's an event zero ...
        idx = find( abs(nefv{eidx}) <= reg_ev(eidx).tol );
        if ~isempty(idx)
            %% set event function to exactly zero
            %% (to prevent possible INTERVAL detection again on next step)
            nefv{eidx}(idx) = 0;

            %% ... make this one the critical one
            ev(i).eidx = eidx;
            ev(i).zero = 1;
            ev(i).step_scale = 1;
            ev(i).log = 1;
            ev(i).name = reg_ev(eidx).name;
            ev(i).idx = idx;
            ev(i).msg = 'ZERO';
            i = i + 1;
        end
    end

    %% and if no zeros were detected
    if i == 1
        %% search for intervals for non-rollback events
        for eidx = 1:nef
            %% current and previous event function signs
            c_sign = sign(nefv{eidx});
            p_sign = sign(cefv{eidx});

            %% if there's been a sign change ...
            idx = find( abs(c_sign) == 1 & c_sign == -p_sign );
            if ~isempty(idx)
                %% ... compute step size scaling factors ...
                step_scale = cefv{eidx}(idx) ./ (cefv{eidx}(idx) - nefv{eidx}(idx));

                %% ... and save the info as an interval detection
                ev(i).eidx = eidx;
                ev(i).zero = 0;
                ev(i).step_scale = step_scale;
                ev(i).log = 0;
                ev(i).name = reg_ev(eidx).name;
                ev(i).idx = idx;
                ev(i).msg = 'INTERVAL';
                i = i + 1;
            end
        end
    end
end

%% update msgs
if ev(1).eidx
    for i = 1:length(ev)
        if length(nefv{ev(i).eidx}) > 1
            s1 = sprintf('(%d)', ev(i).idx);
        else
            s1 = '';
        end
        if rollback
            s2 = sprintf(' : ROLLBACK by %g', ev(i).step_scale);
        else
            s2 = '';
        end
        ev(i).msg = sprintf('%s detected for %s%s event%s', ...
            ev(i).msg, ev(i).name, s1, s2);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [L1, S1, statsPCP] = fastpcp(V, lambda, loops, rank0, rankThreshold, lambdaFactor )

if( nargin < 6 )
  lambdaFactor = 1.0;
  if( nargin < 5 )
    rankThreshold = 0.01;
    if( nargin < 4 )
      rank0 = 1;                % initial rank
      if( nargin < 3 )
        loops = 2;              % total number of outer loops
      end
    end
  end
end

if isempty(rankThreshold)
	rankThreshold = 0.01;
end

if isempty(lambdaFactor)
	lambdaFactor = 1.0;
end

% Data size
%[Nrows,Ncols] = size(V);

% Set flag (increments the rank plus one at each iteration)
inc_rank = 1;

% ---------------------------------
% >>>  measure time performance <<<
% ---------------------------------
%t = tic;

% ------------------------
% --- First outer loop ---
  rank = rank0;                     % current rank
  statsPCP.rank(1) = rank;          % save current rank

  % Partial SVD
  %[Ulan Slan Vlan] = lansvd(V, rank, 'L');
  %[Ulan,Slan,Vlan] = svds(V, rank);
  [Ulan,Slan,Vlan] = svdsecon(V, rank); % fastest

  % Current low-rank approximation
  L1 = Ulan*Slan*Vlan';

  % Shrinkage
  S1 = shrink(V-L1, lambda);

% ------------------------
% ---    Outer loops   ---

for k = 2:loops,

  if(inc_rank == 1)
     lambda = lambda * lambdaFactor;         % modify Lambda at each iteration
     rank = rank + 1;                        % increase rank
  end

  % low rank (partial SVD)
  %[Ulan Slan Vlan] = lansvd(V-S1, rank, 'L');
  %[Ulan,Slan,Vlan] = svds(V-S1, rank);
  [Ulan,Slan,Vlan] = svdsecon(V-S1, rank); % fastest

  currentEvals = diag(Slan);                                    % extract current evals
  statsPCP.rank(k) = length( currentEvals );                    % save current rank
  statsPCP.rho(k) = currentEvals(end) / sum( currentEvals(1:end-1) );       % relative contribution of the last evec

  % simple rule to keep or increase the current rank's value
  if(statsPCP.rho(k) < rankThreshold )
     inc_rank = 0;
  else
     inc_rank = 1;
  end

  % Current low-rank approximation
  L1 = Ulan*Slan*Vlan';

  % Shrinkage
  S1 = shrink(V-L1, lambda);
end

% ---------------------------------
% >>>  measure time performance <<<
% ---------------------------------
%statsPCP.time = toc(t);

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

function u = shrink(v, lambda)
  u = sign(v).*max(0, abs(v) - lambda);
return

author: andrewssobral

repo: lrslibrary


text:
%nav formulas on a (non rotating flat earth)
%vtyp==0->n frame vel, vtyp==1->b frame vel
function [Cbn_new, Vx_new, pos_new]=strapdown_pln_dcm(Cbn, Vx, pos, a, w, g, dt, vtyp)

%%%Update attitude
%Part I:Body frame update
rot=w*dt;
rot_norm=norm(rot);
sr_a=1-(rot_norm^2/6)+(rot_norm^4/120);
sr_b=(1/2)-(rot_norm^2/24)+(rot_norm^4/720);
mx_a=eye(3)+sr_a*skew(rot)+sr_b*skew(rot)*skew(rot);

Cbn_new=Cbn*mx_a;

%%Update Velocity and position
if (vtyp==0) %vel in n frame
    vel_inc=(Cbn*(a*dt))+[0;0;g]*dt;
    Vx_new=Vx+vel_inc;

    pos_new=pos+Vx*dt;
elseif (vtyp==1) %vel in b frame
    vel_inc1=(a+(Cbn'*[0;0;g]))*dt;
    vel_inc2=(cross(Vx,w))*dt;
    Vx_new=Vx+vel_inc1+vel_inc2;

    pos_new=pos+Cbn*Vx*dt;
end

repo: ekfmonoslam


text:
function [P,N,check]=plane_intersect(N1,A1,N2,A2)
%plane_intersect computes the intersection of two planes(if any)
% Inputs:
%       N1: normal vector to Plane 1
%       A1: any point that belongs to Plane 1
%       N2: normal vector to Plane 2
%       A2: any point that belongs to Plane 2
%
%Outputs:
%   P    is a point that lies on the interection straight line.
%   N    is the direction vector of the straight line
% check is an integer (0:Plane 1 and Plane 2 are parallel'
%                              1:Plane 1 and Plane 2 coincide
%                              2:Plane 1 and Plane 2 intersect)
%
% Example:
% Determine the intersection of these two planes:
% 2x - 5y + 3z = 12 and 3x + 4y - 3z = 6
% The first plane is represented by the normal vector N1=[2 -5 3]
% and any arbitrary point that lies on the plane, ex: A1=[0 0 4]
% The second plane is represented by the normal vector N2=[3 4 -3]
% and any arbitrary point that lies on the plane, ex: A2=[0 0 -2]
%[P,N,check]=plane_intersect([2 -5 3],[0 0 4],[3 4 -3],[0 0 -2]);

%This function is written by :
%                             Nassim Khaled
%                             Wayne State University
%                             Research Assistant and Phd candidate
%If you have any comments or face any problems, please feel free to leave
%your comments and i will try to reply to you as fast as possible.
P=[0 0 0];
N=cross(N1,N2);

%  test if the two planes are parallel
if norm(N) < 10^-7                % Plane 1 and Plane 2 are near parallel
    V=A1-A2;
        if (dot(N1,V) == 0)
            check=1;                    % Plane 1 and Plane 2 coincide
           return
        else
            check=0;                   %Plane 1 and Plane 2 are disjoint
            return
        end
end

 check=2;

% Plane 1 and Plane 2 intersect in a line
%first determine max abs coordinate of cross product
maxc=find(abs(N)==max(abs(N)));

%next, to get a point on the intersection line and
%zero the max coord, and solve for the other two

d1 = -dot(N1,A1);   %the constants in the Plane 1 equations
d2 = -dot(N2, A2);  %the constants in the Plane 2 equations

switch maxc
    case 1                   % intersect with x=0
        P(1)= 0;
        P(2) = (d2*N1(3) - d1*N2(3))/ N(1);
        P(3) = (d1*N2(2) - d2*N1(2))/ N(1);
    case 2                    %intersect with y=0
        P(1) = (d1*N2(3) - d2*N1(3))/ N(2);
        P(2) = 0;
        P(3) = (d2*N1(1) - d1*N2(1))/ N(2);
    case 3                    %intersect with z=0
        P(1) = (d2*N1(2) - d1*N2(2))/ N(3);
        P(2) = (d1*N2(1) - d2*N1(1))/ N(3);
        P(3) = 0;
end


text:
%% Face and Eyes Detection
%
% In this demo, we will learn the basics of face detection using Haar
% Feature-based Cascade Classifiers, and how the same extends for eye
% detection, etc.
%
% This program demonstrates the use of |cv.CascadeClassifier| class to detect
% objects (face + eyes). You can use Haar or LBP features. This classifier can
% detect many kinds of rigid objects, once the appropriate classifier is
% trained. It's most known use is for faces.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/facedetect.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/tapi/ufacedetect.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/gpu/cascadeclassifier.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/facedetect.py>
% * <https://docs.opencv.org/3.2.0/db/d28/tutorial_cascade_classifier.html>
% * <https://docs.opencv.org/3.2.0/d7/d8b/tutorial_py_face_detection.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/objectDetection/objectDetection.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/objectDetection/objectDetection2.cpp>
%

%% Theory
%
% Object Detection using Haar feature-based cascade classifiers is an
% effective object detection method proposed by Paul Viola and Michael Jones
% in their paper, "Rapid Object Detection using a Boosted Cascade of Simple
% Features" in 2001. It is a machine learning based approach where a cascade
% function is trained from a lot of positive and negative images. It is then
% used to detect objects in other images.
%
% Here we will work with face detection. Initially, the algorithm needs a lot
% of positive images (images of faces) and negative images (images without
% faces) to train the classifier. Then we need to extract features from it.
% For this, haar features shown in the below image are used. They are just like
% our convolutional kernel. Each feature is a single value obtained by
% subtracting sum of pixels under the white rectangle from sum of pixels under
% the black rectangle.
%
% <<https://docs.opencv.org/3.4.0/haar_features.jpg>>
%
% Now all possible sizes and locations of each kernel are used to calculate
% lots of features. (Just imagine how much computation it needs? Even a
% 24x24 window results over 160000 features). For each feature calculation, we
% need to find the sum of the pixels under white and black rectangles. To
% solve this, they introduced the integral image. However large your image, it
% reduces the calculations for a given pixels, to an operation involving just
% four pixels. Nice, isn't it? It makes things super-fast.
%
% But among all these features we calculated, most of them are irrelevant. For
% example, consider the image below. The top row shows two good features. The
% first feature selected seems to focus on the property that the region of the
% eyes is often darker than the region of the nose and cheeks. The second
% feature selected relies on the property that the eyes are darker than the
% bridge of the nose. But the same windows applied to cheeks or any other
% place is irrelevant. So how do we select the best features out of 160000+
% features? It is achieved by *Adaboost*.
%
% <<https://docs.opencv.org/3.4.0/haar.png>>
%
% For this, we apply each and every feature on all the training images. For
% each feature, it finds the best threshold which will classify the faces to
% positive and negative. Obviously, there will be errors or misclassifications.
% We select the features with minimum error rate, which means they are the
% features that most accurately classify the face and non-face images. (The
% process is not as simple as this. Each image is given an equal weight in the
% beginning. After each classification, weights of misclassified images are
% increased. Then the same process is done. New error rates are calculated.
% Also new weights. The process is continued until the required accuracy or
% error rate is achieved or the required number of features are found).
%
% The final classifier is a weighted sum of these weak classifiers. It is
% called weak because it alone can't classify the image, but together with
% others forms a strong classifier. The paper says even 200 features provide
% detection with 95% accuracy. Their final setup had around 6000 features.
% (Imagine a reduction from 160000+ features to 6000 features. That is a big
% gain).
%
% So now you take an image. Take each 24x24 window. Apply 6000 features to it.
% Check if it is face or not. Wow.. Isn't it a little inefficient and time
% consuming? Yes, it is. The authors have a good solution for that.
%
% In an image, most of the image is non-face region. So it is a better
% idea to have a simple method to check if a window is not a face region. If
% it is not, discard it in a single shot, and don't process it again. Instead,
% focus on regions where there can be a face. This way, we spend more time
% checking a possible face region.
%
% For this they introduced the concept of *Cascade of Classifiers*. Instead of
% applying all 6000 features on a window, the features are grouped into
% different stages of classifiers and applied one-by-one. (Normally the first
% few stages will contain very many fewer features). If a window fails the
% first stage, discard it. We don't consider the remaining features on it. If
% it passes, apply the second stage of features and continue the process. The
% window which passes all stages is a face region. How is that plan!
%
% The authors' detector had 6000+ features with 38 stages with 1, 10, 25, 25
% and 50 features in the first five stages. (The two features in the above
% image are actually obtained as the best two features from Adaboost).
% According to the authors, on average, 10 features out of 6000+ are evaluated
% per sub-window.
%
% So this is a simple intuitive explanation of how Viola-Jones face detection
% works. Read the paper for more details or check out the following references:
%
% * Video Lecture on
%   <http://www.youtube.com/watch?v=WfdYYNamHZ8 Face Detection and Tracking>
% * An interesting interview regarding Face Detection by
%   <http://www.makematics.com/research/viola-jones/ Adam Harvey>
%
% OpenCV comes with a trainer as well as detector. If you want to train your
% own classifier for any object like car, planes etc. you can use OpenCV to
% create one. See the OpenCV docs for full details on Cascade Classifier
% Training.
%
% Here we will deal with detection. OpenCV already contains many pre-trained
% classifiers for face, eyes, smiles, etc. Those XML files are stored in
% the |opencv/data/haarcascades/| folder.
%

%% Code
%
% In this example, we will create a face and eyes detector with OpenCV:
%
% * First we need to load the required XML classifiers.
% * Then load our input image (or video) in grayscale mode.
% * Now we find the faces in the image. If faces are found, it returns the
%   positions of each detected faces as a rectangle |[x,y,w,h]|. Once we get
%   these locations, we can create a ROI for the face and apply eye detection
%   on this ROI (since eyes are always on the face!).
%

%% Options

% this is the primary trained classifier such as frontal face
cascadeName = fullfile(mexopencv.root(),'test','haarcascade_frontalface_alt.xml');
% this an optional secondary classifier such as eyes
nestedCascadeName = fullfile(mexopencv.root(),'test','haarcascade_eye_tree_eyeglasses.xml');
% image scale greater or equal to 1, try 1.3 for example
scale = 1.0;
% attempts detection of flipped image as well
tryflip = false;

%% Initialization

% download XML files if missing
download_classifier_xml(cascadeName);
download_classifier_xml(nestedCascadeName);

% load cacade classifiers
cascade = cv.CascadeClassifier(cascadeName);
assert(~cascade.empty(), 'Could not load classifier cascade');
nestedCascade = cv.CascadeClassifier();
if ~nestedCascade.load(nestedCascadeName)
    disp('Could not load classifier cascade for nested objects');
end
scale = max(scale, 1.0);

%% Main loop
% (either video feed or a still image)
if false
    % read an image
    frame = cv.imread(fullfile(mexopencv.root(),'test','lena.jpg'), 'Color',true);
    % detect faces/eyes and draw detections
    frame = detectAndDraw(frame, cascade, nestedCascade, scale, tryflip);
    imshow(frame);
else
    % prepare video input
    cap = cv.VideoCapture();
    pause(1);
    assert(cap.isOpened());

    % prepare figure
    frame = cap.read();
    assert(~isempty(frame));
    hImg = imshow(frame);

    % video feed
    while ishghandle(hImg)
        % read frame
        frame = cap.read();
        if isempty(frame), break; end

        % detect faces/eyes and draw detections
        frame = detectAndDraw(frame, cascade, nestedCascade, scale, tryflip);

        % update
        set(hImg, 'CData',frame);
        drawnow;
    end
    cap.release();
end

%% Processing function

function img = detectAndDraw(img, cascadeF, cascadeE, scale, tryflip)
    % downscale image and preprocess it
    fx = 1/scale;
    gray = cv.cvtColor(img, 'RGB2GRAY');
    gray = cv.resize(gray, fx, fx);
    gray = cv.equalizeHist(gray);
    [h,w] = size(gray);

    % detection options
    detectOpts = {
        'ScaleFactor',1.1, ...
        'MinNeighbors',2, ...
        ... 'FindBiggestObject',true, ...
        ... 'DoRoughSearch',true, ...
        'ScaleImage',true, ...
        'MinSize',[30 30]
    };

    % detect faces
    tic
    faces = cascadeF.detect(gray, detectOpts{:});
    if tryflip
        faces2 = cascadeF.detect(cv.flip(gray, 1), detectOpts{:});
        faces2 = cellfun(@(r) [w-r(1)-r(3) r(2:4)], faces2, 'UniformOutput',false);
        faces = [faces(:); faces2(:)];
    end
    toc

    % draw
    clrs = uint8(255 * lines(7));
    for i=1:numel(faces)
        r = faces{i};
        ii = mod(i-1, size(clrs,1)) + 1;
        drawOpts = {'Color',clrs(ii,:), 'Thickness',3};

        % draw faces
        aspect_ratio = r(3)/r(4);
        if 0.75 < aspect_ratio && aspect_ratio < 1.3
            center = round((r(1:2) + r(3:4)*0.5) * scale);
            radius = round((r(3) + r(4)) * 0.25*scale);
            img = cv.circle(img, center, radius, drawOpts{:});
        else
            pt1 = round(r(1:2) * scale);
            pt2 = round((r(1:2) + r(3:4) - 1) * scale);
            img = cv.rectangle(img, pt1, pt2, drawOpts{:});
        end

        if ~cascadeE.empty()
            % detect nested objects (eyes)
            if false && mexopencv.require('images')
                grayROI = imcrop(gray, [r(1:2)+1 r(3:4)]);
            else
                grayROI = cv.Rect.crop(gray, r);
            end
            nestedObjs = cascadeE.detect(grayROI, detectOpts{:});

            % draw eyes
            for j=1:numel(nestedObjs)
                nr = nestedObjs{j};
                center = round((r(1:2) + nr(1:2) + nr(3:4)*0.5) * scale);
                radius = round((nr(3) + nr(4)) * 0.25*scale);
                img = cv.circle(img, center, radius, drawOpts{:});
            end
        end
    end
end

%% Helper function

function download_classifier_xml(fname)
    if exist(fname, 'file') ~= 2
        % attempt to download trained Haar/LBP/HOG classifier from Github
        url = 'https://cdn.rawgit.com/opencv/opencv/3.4.0/data/';
        [~, f, ext] = fileparts(fname);
        if strncmpi(f, 'haarcascade_', length('haarcascade_'))
            url = [url, 'haarcascades/'];
        elseif strncmpi(f, 'lbpcascade_', length('lbpcascade_'))
            url = [url, 'lbpcascades/'];
        elseif strncmpi(f, 'hogcascade_', length('hogcascade_'))
            url = [url, 'hogcascades/'];
        else
            error('File not found');
        end
        fprintf('Downloading cascade classifier "%s"...\n', [f ext]);
        url = [url f ext];
        urlwrite(url, fname);
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function Wop = normalize_wavelet_factory_1d(N, filter_options, scat_options, epsilon)

    if nargin <4
        epsilon = 2^(-20);
    end

    filters = filter_bank(N, filter_options);

    for m = 0:scat_options.M
        filt_ind = min(numel(filters), m+1);
        Wop{m+1} = @(X)(wavelet_renorm(X, m));
    end

    function [S, Utilde] = wavelet_renorm(U, m)
        filt_ind = min(numel(filters), m+1);
        [S, U] = wavelet_layer_1d(U, filters{filt_ind}, scat_options);
        Utilde = renorm_low_pass_layer_1d(S, U, epsilon);
    end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function p = JDDLDR_UP3(X,pi,ZA,ZB,D,C,gamma1,gamma2,MaxIteration_num)
% X: original training data

beta = 0.001;
A  = [];
for ci = 1:size(C,2)
    A = [A D(ci).M*C(ci).M];
end

S = gamma1*ZA*ZA'+gamma2*ZB*ZB';
% S = gamma1*ZA*ZA';
p = pi;
numcomps = size(p,2);

iter_num_sub= 1;
phi1 = X*X';
phi2 = 2*X*A';

while iter_num_sub<MaxIteration_num % subject iteration for updating p, iteration ends when the contiguous function get close enough or reach teh max_iteration
    phi = (X-p*A)*(X-p*A)';
%     phi = eye(size(phi));
    [U,V] = eig(phi-S); p1 = U(:,1:numcomps);
      sum(sum(abs(p+beta*(p1-p) - p)));
      p = p+beta*(p1-p);
%     p = p1;
%     clear U V W phi
    iter_num_sub=iter_num_sub+1;
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
% downloading all models

addpath('../functions/');

mkdir_if_missing('../weight/');
mkdir_if_missing('./temp');

%% for detection files
fprintf('Downloading faster_rcnn_final_model...\n');
urlwrite('https://onedrive.live.com/download?resid=D7AF52BADBA8A4BC!114&authkey=!AERHoxZ-iAx_j34&ithint=file%2czip', ...
    './temp/faster_rcnn_final_model.zip');

fprintf('Unzipping...\n');
unzip('./temp/faster_rcnn_final_model.zip', './temp/');

movefile ./temp/output/faster_rcnn_final/faster_rcnn_VOC0712_vgg_16layers/detection_final  ../weight/detection_final_VGG
movefile ./temp/output/faster_rcnn_final/faster_rcnn_VOC0712_ZF/detection_final  ../weight/detection_final_ZF
fprintf('Detection Models Done.\n');

rmdir('./temp/', 's');

%% for segmentation files
fprintf('Downloading fcn models...\n')
urlwrite('http://dl.caffe.berkeleyvision.org/fcn8s-heavy-pascal.caffemodel', ...
    '../weight/fcn8s-heavy-pascal.caffemodel');

urlwrite('http://dl.caffe.berkeleyvision.org/fcn-alexnet-pascal.caffemodel', ...
    '../weight/fcn-alexnet-pascal.caffemodel');

fprintf('Segmentation Models Done.\n');

author: cihangxie

repo: DAG

save_path: github-repos/MATLAB/cihangxie-DAG


text:
function test_suite = test_sample_unique
% tests for cosmo_sample_unique
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_sample_unique_basics
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_sample_unique(varargin{:}),'');

    for k=from_range_randomly(10,5)
        for n=from_range_randomly(10,5);
            for count=[0,from_range_randomly(100,10)];
                if count==0
                    args={k,n};
                    ncol=1;
                else
                    args={k,n,count};
                    ncol=count;
                end

                if k>n
                    % n too large; should throw error
                    aet(args{:});
                    continue;
                else
                    i=cosmo_sample_unique(args{:});
                    assertEqual(sort(i,1),i);

                    % check size
                    assert(isequal(size(i),[k,ncol]))

                    % check contents
                    i_vec=i(:);
                    assert(all(i_vec>=1 | i_vec<=n | round(i_vec)==i_vec));

                    % each element must be select about equally often,
                    % and each column must have unique elements
                    i_s=sort(i,1);

                    % each element unique
                    assert(all(all(diff(i_s,1,1)>0)));

                    % check counts
                    counts=zeros(n,1);
                    for col=1:ncol
                        counts(i(:,col))=counts(i(:,col))+1;
                    end

                    % counts differ at most by one
                    assert(min(counts)+1>=max(counts));
                end
            end
        end
    end

function test_sample_unique_full_randperm()
    k=2+from_range_randomly(10);
    n=k;
    c=5+from_range_randomly(10);

    i=cosmo_sample_unique(k,n,c);
    assertEqual(sort(i,1),repmat((1:k)',1,c));

function test_sample_unique_difference_sequences()
    k=2+from_range_randomly(10);
    n=k+from_range_randomly(10);
    c=5+from_range_randomly(10);

    seed=1+from_range_randomly(1e5);
    last_args={{'seed',seed},{'seed',[]'},{}};
    for k=1:numel(last_args)
        all_same=true;

        last_arg=last_args{k};
        use_seed=numel(last_arg)==2 && ~isempty(last_arg{2});

        for tries=1:5
            i=cosmo_sample_unique(k,n,c,last_arg{:});
            if tries==1
                i_first=i;
            elseif ~isequal(i_first,i)
                all_same=false;
                break;
            end
        end

        if xor(all_same,use_seed)
            assertFalse(true,'randomness mismatch with respect to seed');
        end
    end

function test_sample_unique_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_sample_unique(varargin{:}),'');

    bad_args={-1, 1.5, [1 1], struct(), false};
    for k=1:numel(bad_args)
        bad_arg=bad_args{k};

        for j=1:3
            all_args={1,1,1};
            all_args{j}=bad_arg;
            aet(all_args{:});
        end
    end

function vs=from_range_randomly(n, count)
    if nargin<=2
        count=1;
    end

    assert(n>=count);
    vals=1:n;

    rp=randperm(n);
    vs=vals(rp(1:count));

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [x,fVal,exitCode]=divRectOpt(f,lowerBounds,upperBounds,options)
%%DIVRECTOPT The dividing rectangles optimization algorithm. This is a
%           derivative-free global optimization algorithm where the
%           multidimensional search space can be bounded. This algorithm
%           tries to minimize a given function of a vector parameter.
%
%INPUTS: f A handle to the function over which the minimization is to be
%          performed. The function [fVal,violatesConstraints]=f(x) takes
%          the NX1 xvector and returns the real scalar function value
%          fVal. The parameter violatesConstraints is used to mark
%          forbidden areas, for example, where the function might not be
%          defined. violatesConstraints=false if there are not constraint
%          violations and it is true if there is a constraint violation.
% lowerBounds,upperBounds NX1 vectors of lower bounds defining lower and
%          upper bounds on the components of x at the optimal points. This
%          defines the search space. The elements must be finite and
%          cannot be NaNs.
%  options A structure where all elements are optional. The elements set
%          options that affect the algorithm and its convergence. Default
%          values are used if a field is missing. Possible fields in the
%          structure are:
%          algorithm (default 0) This selects which of the two direct
%                    optimization algorithms should be posed. Possible
%                    values are
%                    0: Use the algorithm of [1], which is probably
%                       better for functions with many local minima.
%                    1: Use the algorithm of [2], which favors local
%                       search and thus might be better for functions
%                       with few local minima.
%             maxDiv (default 5000) The maximum number of divisions of
%                    the hyperrectangles allowed. This must be >0.
%            maxIter (default 100) The maximum number of iterations.
%           maxFEval (default 500) The maximum number of function
%                     evaluations.
%            epsilon (default 1e-4) The Jones' factor. It is the epsilon
%                    term in Definitions 3.1 and 4.1 in [1]. It affects
%                    the convergence rate of the algorithm by avoiding
%                    oversampling near points with low function values
%                    and biasing the sampling towards a global search.
%                    If a positive value of epsilon is specified, then
%                    the same value is used for all iterations. In [1],
%                    it is suggested that epsilon be set to the desired
%                    solution accuracy (in f) for good convergence,
%                    though other convergence conditions are discussed in
%                    [2]. If a negative value is specified, then epsilon
%                    is adaptively changed using the formula
%                    epsilon = max(1.D-4*abs(fMin),abs(epsilon)).
%                    The magnitude of epislon must be between 0 and 1,
%                    not inclusive of 1, though 0 tends to work.
%         epsilonAbs (default 0) Definitions 3.1 and 4.1 in [1] involve
%                    comparisons to fMin-epsilon*abs(fMin), where fMin is
%                    the current minimum function value found. If
%                    epsilonAbs is nonzero, then the comparison is
%                    actually taken with respect to
%                    min(fMin-epsilon*abs(fMin),epsilonAbs)
%                    epsilonAbs must be positive.
%       volumeRelTol (default 1e-10) Convergence is declared if the
%                    volume of the search region is this fraction of the
%                    original volume (must be between 0 and 1)
%        sigmaRelTol (default 0) Convergence is declared if the hypercube
%                    measure is smaller than this value. The notion of
%                    the hypercube measure sigma(S) is discussed in [2].
%            fGlobal If the actual minimum value of f is known, even
%                    though the value of x providing that minimum is
%                    unknown, then it can be specified. Otherwise, this
%                    parameter is not used.
%      fGlobalRelTol (default 1e-12 if fGlobal is provided) This is only
%                    used if fGlobal is provided. Convergence is declared
%                    if
%                    (fMin - fGlobal)/max(1,abs(fGlobal)) < fGlobalRelTol
%
%OUTPUTS: x The NX1 optimal point found.
%      fVal The value of the function at the optimal point found.
%  exitCode A return value indicating the status of the algorithm on
%          termination. Possible values are:
%          1 Function terminated, because  the number of function
%            evaluations done is larger than maxFEval.
%          2 Function terminated, because the number of iterations equals
%            maxIter.
%          3 The function value found is within fGlobalRelTol of the
%            global optimum given in fGlobal.
%          4 The volume of the hyperrectangle with fVal at its center is
%            less than volumeRelTol times the volume of the original
%            hyperrectangle.
%          5 The measure of the hyperrectangle with fVal at its center is
%            less than sigmaRelTol.
%
%Global optimization algorithms such as this are very good at handling
%functions with many local minima. Often, one might want to use a global
%optimization algorithm to get an initial estimate close to the global
%solution and then do a few iterations of, for example, Newton's method,
%to find the global optimum.
%
%This function is an interface to the DIRECT library that is part of the
%NLOpt set of nonlinear optimization routines available at [4], which
%implements [1] and [2].  Though the aforementioned code online provides
%a Matlab interface, it was decided to include the dividing rectangles
%algorithm without the rest of the library, in part because some of the
%other functions include weak copyleft provisions in the form of the LGPL.
%Additionally, the implementation above did not allow one to specify the
%maximum number of hyperrectangle divisions (it was statically fixed) nor
%did it use the memory allocation routines in Matlab that mex files are
%supposed to (but do not have to) use. Thus, the code was slightly
%modified (but is still included with the rest of the third party code).
%
%Though [1] and [2] use a number of exmaple, functions, they formulae are
%not explicitely given. Rather, one must consult a few references, which
%can be tedious to obtain. On the other hand, numrous exmaples functions
%used in [1] are given explicitely in [3]. Examples of the algorithm are
%thus:
% %Branin's RCOS function
% %Note that the deal function is used to make an anonymous function have
% %two outputs.
% f=@(x)deal((x(2)-(5/(4*pi^2))*x(1)^2+(5/pi)*x(1)-6)^2+10*(1-(1/(8*pi)))*cos(x(1))+10,0);
% lowerBounds=[-5;10];
% upperBounds=[0;15];
% [x,fVal,exitCode]=divRectOpt(f,lowerBounds,upperBounds)
%The globally optimal value is given in [3] as 0.397887357729739. The
%returned fVal is close.
%
% %Six-Hump Camel
% f=@(x)deal((x(1)^2*(4-2.1*x(1)^2+x(1)^4/3)+x(1)*x(2)+x(2)^2*(-4+4*x(2)^2)),0);
% lowerBounds=[-3;-2];
% upperBounds=[3;2];
% [x,fVal,exitCode]=divRectOpt(f,lowerBounds,upperBounds)
%The globally optimal value is given in [3] as -1.0316284535. The returned
%fVal is close.
%
% %Two-Dimensional Schubert
% f=@(x)deal(sum((1:5).*cos(((1:5)+1)*x(1)+(1:5)))*sum((1:5).*cos(((1:5)+1)*x(2)+(1:5))),0)
% lowerBounds=[-10;-10];
% upperBounds=[10;10];
% options.maxFEval=5000;
% options.maxIter=200;
% [x,fVal,exitCode]=divRectOpt(f,lowerBounds,upperBounds,options)
%Here, increasing the number of function evaluations and iterations is
%important to getting a solution that is close to the optimal value as
%this function has 760 local minima (and 18 global minima). The global
%minimum is -186.730908831024 as per [3].
%
%REFERENCES:
%[1] D. R. Jones, C. D. Peritunen, and B. E. Stuckman, "Lipschitzian
%    optimization without the Lipschitz constant," Journal of Optimization
%    Theory and Application, vol. 79, no. 1, pp. 157-181, Oct. 1993.
%[2] J. M. Gablonsky and C. T. Kelley, "A locally-biased form of the
%    DIRECT algorithm," Journal of Global Optimization, vol. 21, no. 1,
%    pp. 27-37, Sep. 2001.
%[3] M. Bj�kman and K. Holmstr�m, "Global optimization using the DIRECT
%    algorithm in Matlab," The Electronic International Journal Advanced
%    Modeling and Optimization, vol. 1, no. 2, pp. 17-37, 1999.
%    [Online]. Available: http://camo.ici.ro/journal/v1n2.htm
%[4] S. G. Johnson. (2014, 20 May) The NLopt nonlinear-optimization
%    package. [Online]. Available: http://ab-initio.mit.edu/nlopt
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function update_editfield_date(src,~)
    %update_editfield_date general function to update x.Value with a datenum when x.String is updated, where x is an editfield
    % use as the callback
    %
    %  initialdate = '2015-01-23 14:31:55.01'
    %  uicontrol('Style','edit','string',initialdate,'Value',datenum(initialdate),'Callback',@update_editfield_date);
    % can interpret as decimal year or full date

    isTextdate = contains(src.String,{':',' ','/','-'});

    if isTextdate
        src.Value = datenum(datetime(src.String)); %provides extra parsing
    else
        d = str2double(get(src,'String'));
        if d >= 1800 && d < 3000
            %treat as year
            set(src,'String',[get(src,'String'), '-01-01 00:00:00']);
            update_editfield_date(src);
            return
        else
            try
                src.Value=str2double(src.String);
            catch ME
                warning(ME)
            end
        end
    end
end

author: CelsoReyes


text:
function [geph,stat]=decode_geph(ver,sat,toc,data)
global glc gls

geph=gls.geph;
stat=1;

[sys,~]=satsys(sat);
if sys~=glc.SYS_GLO
    stat=0;return;
end

geph.sat=sat;

% toc rounded by 15 min in utc
[week,tow]=time2gpst(toc);
toc=gpst2time(week,floor((tow+450)/900)*900);
dow=fix(floor(tow/86400));

% time of frame in utc
if ver<=2.99
    tod=data(3);
else
    tod=rem(data(3),86400);
end
tof=gpst2time(week,tod+dow*86400);
tof=adjday(tof,toc);

geph.toe=utc2gpst(toc);
geph.tof=utc2gpst(tof);

geph.iode=fix(rem(tow+10800,86400)/900+0.5);

geph.taun=-data(1);
geph.gamn=data(2);

geph.pos=[data(4);data(8);data(12)]*1e3;
geph.vel=[data(5);data(9);data(13)]*1e3;
geph.acc=[data(6);data(10);data(14)]*1e3;

geph.svh=fix(data(7));
geph.frq=fix(data(11));
geph.age=fix(data(15));

% some receiver output >128 for minus frequency number
if geph.frq>128
    geph.frq=geph.frq-256;
end

return

repo: GINav


text:
function [fb] = fbCreate(numOrient,startSigma,numScales,scaling,elong)
% function [fb] = fbCreate(numOrient,startSigma,numScales,scaling,elong)
%
% Create a filterbank containing numOrient even and odd-symmetric
% filters and one center-surround filter at numScales scales.
%
% The even-symmetric filter is a Gaussian second derivative.
% The odd-symmetric filter is its Hilbert transform.
%
% See also oeFilter, csFilter, fbRun.
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

if nargin<3, numScales = 1; end
if nargin<4, scaling = sqrt(2); end
if nargin<5, elong = 3; end
support = 3;

fb = cell(2*numOrient,numScales);
for scale = 1:numScales,
  sigma = startSigma * scaling^(scale-1);
  for orient = 1:numOrient,
    theta = (orient-1)/numOrient * pi;
    fb{2*orient-1,scale} = oeFilter(sigma*[elong 1],support,theta, 2,0);
    fb{2*orient,scale} = oeFilter(sigma*[elong 1],support,theta,2,1);
  end
  %fb{2*numOrient+1,scale} = csFilter(sigma*[3 1],support);
end

author: zouchuhang

repo: LayoutNet


text:
function cleaned_matches = removeOutliers(I1_l, I1_r, I2_l, I2_r, matches, match_params)
%REMOVEOUTLIERS Performs outliers removal by applying normalized cross-correlation
% score between patches around the matched points.
%
% INPUT:
%
%   - I1_l: left image for time t-1
%   - I1_r: right images for time t-1
%   - I2_l: left image for time t
%   - I2_r: right images for time t
%   - matches (1,N): strcutre of matched points
%       - pt1_l: keypoint in left image at time t-1
%       - pt1_r: keypoint in right image at time t-1
%       - pt2_l: keypoint in left image at time t
%       - pt2_r: keypoint in right image at time t
%   - match_params: structure comprising of following
%       - match_ncc_window: window size of the patch for normalized cross-correlation
%       - match_ncc_tolerance: threshold for normalized cross-correlation
%
% OUTPUT:
%   - cleaned_matches(1, N2): processed structure of matched points
%       - pt1_l: location of keypoint in left image at time t-1
%       - pt1_r: location of keypoint in right image at time t-1
%       - pt2_l: location of keypoint in left image at time t
%       - pt2_r: location of keypoint in right image at time t

% initialize parameters
ncc_window = match_params.match_ncc_window;
ncc_tolerance = match_params.match_ncc_tolerance;

% store indices of matches that are valid
valid_match_indices = ones(length(matches), 1, 'logical');

% for all matches do
parfor i = 1:length(matches)
    match = matches(i);

    % verify match by computing NCC score
    % a) left frame at t-1 with right frame at t-1
    validity = verifyMatchNCC(I1_l, I1_r, match.pt1_l.location, match.pt1_r.location, ...
                                ncc_window, ncc_tolerance);
    if not(validity)
        valid_match_indices(i) = 0;
        continue;
    end

    % b) right frame at t-1 with right frame at t
    validity = verifyMatchNCC(I1_r, I2_r, match.pt1_r.location, match.pt2_r.location, ...
                                ncc_window, ncc_tolerance);
    if not(validity)
        valid_match_indices(i) = 0;
        continue;
    end

    % c) right frame at t with left frame at t
    validity = verifyMatchNCC(I2_r, I2_l, match.pt1_r.location, match.pt2_l.location, ...
                                ncc_window, ncc_tolerance);
    if not(validity)
        valid_match_indices(i) = 0;
        continue;
    end

    % d) left frame at t with left frame at t-1
    validity = verifyMatchNCC(I2_l, I1_l, match.pt2_l.location, match.pt1_l.location, ...
                                ncc_window, ncc_tolerance);
    if not(validity)
        valid_match_indices(i) = 0;
        continue;
    end
end

% cleaned matched points after NCC
cleaned_matches = matches(valid_match_indices);

end

repo: Stereo-Odometry-SOFT


text:
function result = fit_rayleigh_pdf( x,y,W,hAx )
% fit_rayleigh_pdf - Non Linear Least Squares fit of the Rayleigh distribution.
%                    given the samples of the histogram of the samples, finds the
%                    distribution parameter that fits the histogram samples.
%
%    fits data to the probability of the form:
%        p(r)=r*exp(-r^2/(2*s))/s
%    with parameter: s
%
% format:   result = fit_rayleigh_pdf( x,y,W,hAx )
%
% input:    y   - vector, samples of the histogram to be fitted
%           x   - vector, position of the samples of the histogram (i.e. y = f(x,a))
%           W   - matrix or scalar, a square weighting matrix of the size NxN where
%                 N = length(y), or 0 to indicate no weighting is needed.
%           hAx - handle of an axis, on which the fitted distribution is plotted
%                 if h is given empty, a figure is created.
%
% output:   result  - structure with the fields
%                      s   - fitted parameter
%                      VAR - variance of the estimation
%                      type- weighted LS or not weighted LS
%                      iter- number of iteration for the solution
%

%
% Algorithm
% ===========
%
% We use the WLS algorithm to estimate the PDF from the samples.%
% The rayleigh distribution is given by:
%
%    p(x,s) = x*exp(-x^2/(2*s))/s
%
%    note that X is known and therefore, considered a constant vector
%
% The non liner WLS estimator is given by:
%
%    s(n+1) = s(n) + inv(H'*W*H)*(H') * (y-h) = s(n) + G * err
%
%    where:   h = p(x,s)
%             H = diff( p(x,a) ) with respect to "s"
%             W = weighting matrix of size NxN  (N = length(y))
%             s = a single parameter to be estimated
%
% The error estimation is given by:
%
%    VAR( s ) = G * VAR( err ) * (G')
%
%       or when W=I and the noise is a gaussian noise
%
%    VAR( s ) = inv( H' * H )
%

if (nargin<3)
    error( 'fit_rayleigh_pdf - insufficient input arguments' );
end

s       = x(find(y==max(y))).^2;        % initial guess
y       = y(:);                         % both should be column vectors !
x       = x(:);
x2      = x.^2;                         % save computation time
thresh  = 0.995;                        % convergence threshold for the loop
last_cnt= inf;
iter    = 0;

% check weight matrix input
if (size(W,1)==length(y)) & (size(W,2)==length(y))
    weights_flag    = 1;
    type            = 'WLS';
else
    weights_flag    = 0;
    type            = 'LS';
end

% Estimation
% =============
if (weights_flag)
    % loop for convergence (with weighting matrix)
    % =============================================
    while (1)
        iter    = iter + 1;
        h       = x.*exp( -x2/(2*s) )/s;
        H       = h.*(x2/(2*s^2) - 1/s);
        HTW     = H'*W;
        e       = inv( HTW * H ) * HTW * (y-h);
        s       = s + e;
        control = e*e;
        if ( control > (last_cnt * thresh) )
            break;
        else
            last_cnt = control;
        end
    end

    % summarize results
    h       = x.*exp( -x2/(2*s) )/s;
    H       = h.*(x2/(2*s^2) - 1/s);
    HTW     = H'*W;
    G       = inv( HTW * H ) * HTW;
    err     = (y-h);
    VAR     = G * var(err) * (G');
else

    % loop for convergence (without a weighting matrix)
    % ==================================================
    while (1)
        iter    = iter + 1;
        h       = x.*exp( -x2/(2*s) )/s;
        H       = h.*(x2/(2*s^2) - 1/s);
        HT      = H';
        control = inv( HT * H );
        s       = s + control * HT * (y-h);
        if ( control>(last_cnt * thresh) )
            break;
        else
            last_cnt = control;
        end
    end

	% summarize results
    h       = x.*exp( -x2/(2*s) )/s;
    H       = h.*(x2/(2*s^2) - 1/s);
    err     = (y-h);
    VAR     = inv( (H') * H );
end

% finish summarizing results
% ============================
result.s    = s;
result.VAR  = VAR;
result.RMS  = sqrt( (err')*err/ (x(2)-x(1))^2 / (length(err)-1) );
result.iter = iter;
result.type = type;

% plot distribution if asked for
% ===============================
if (nargin>3)
    if ishandle( hAx )
        plot_rayleigh( x,result,hAx,4 );
    else
        figure;
        plot_rayleigh( x,result1,gca,4 );
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function p = matrix_T_pdf(A, M, V, K, n)
% MATRIX_T_PDF Evaluate the density of a matrix under a Matrix-T distribution
% p = matrix_T_pdf(A, M, V, K, n)

% See "Bayesian Linear Regression", T. Minka, MIT Tech Report, 2001

[d m] = size(K);
is = 1:d;
c1 = prod(gamma((n+1-is)/2)) / prod(gamma((n-m+1-is)/2));
c2 = det(K)^(d/2) / det(pi*V)^(m/2); %% pi or 2pi?
p = c1 * c2 * det((A-M)'*inv(V)*(A-M)*K + eye(m))^(-n/2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function r = isnan(p)
%ISNAN        Logical result:  polynomial p contains NaN components, y/n
%
%   r = isnan(p)
%
%Result 1 iff at least one coefficient of p contains some NaN
%

% written  09/16/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  r = any(isnan(p.c));

repo: OpenMAS