text:
function N = delaunay_neighbors(q,P)
  % DELAUNAY_NEIGHBORS Find the Delaunay neighbors of q in a list of points P.
  %
  % N = delaunay_neighbors(q,P)
  %
  % Experimentally for uniformly randomly distributed points the number of
  % iterations is O(|P|) and specifically seems to do 1.355*|P| "iterations".
  % There's a sort on the angles so at best this implementation is O(|P|log|P|),
  % but it's very simple... so maybe the constants and overheads are much
  % smaller than doing a full Delaunay triangulation.
  %
  % Inputs:
  %   q  2D position of query point
  %   P  #P by 2 list of nearby points
  % Outputs:
  %   N  #N list of indices into P of neighboring points
  %

  % Example:
  %   q = [0 0];
  %   P = rand(60,2)*2-1;
  %   N = delaunay_neighbors(q,P);
  %   Nm = find(adjacency_matrix(delaunay([P;q]))*sparse(size(P,1)+1,1,1,size(P,1)+1,1));
  %   assert(isempty(setxor(N,Nm)));

  % 2----1
  % |θ₂ / \
  % |  /   \
  % | /     \
  % |/θq   θ₀\
  % q---------0
  %% Angle at p0-q-p1
  %thq = th(p1)-th(p0);
  %% Law of sines: θq/|01| = θ₀/|q1|
  %% θ₀ = |q1|*θq/|01|
  %th0 = l(p1)*thq/sqrt(sum((P(p0,:)-P(p1,:)).^2,2));
  % Vector from p0 to p1

  % Inputs:
  %   P  #P by 2 list of nearby points
  %   l  #P by 2 list of distances to origin: l = normrow(P)
  %   p0  first 2d point
  %   p1  second 2d point
  %   p2  third 2d point
  % Outputs:
  %   flag  whether edge from origin to p1 forms a Delaunay edge between points
  %     p1 and p2.
  %
  isdel = @(P,l,p0,p1,p2) ...
    acos((-P(p0,:)/l(p0))*(((P(p1,:)-P(p0,:))/sqrt(sum((P(p1,:)-P(p0,:)).^2,2)))')) + ...
    acos((-P(p2,:)/l(p2))*(((P(p1,:)-P(p2,:))/sqrt(sum((P(p1,:)-P(p2,:)).^2,2)))')) <= pi;

  %% Or perhaps: http://stackoverflow.com/a/8523979
  %function f = isdel(P,l,p0,p1,p2)
  %  f = acos((-P(p0,:)/l(p0))*(((P(p1,:)-P(p0,:))/sqrt(sum((P(p1,:)-P(p0,:)).^2,2)))')) + ...
  %  acos((-P(p2,:)/l(p2))*(((P(p1,:)-P(p2,:))/sqrt(sum((P(p1,:)-P(p2,:)).^2,2)))')) <= pi;
  %  % Shewchuk's incircle test via determinant
  %  h(1,:) = [P(p0,:) 1];
  %  h(2,:) = [P(p1,:) 1];
  %  h(3,:) = [P(p2,:) 1];
  %  d = det([ ...
  %    q(1) q(2) q(1)^2+q(2)^2 1 ; ...
  %    h(1,1) h(1,2) h(1,1)^2+h(1,2)^2 h(1,3) ; ...
  %    h(2,1) h(2,2) h(2,1)^2+h(2,2)^2 h(2,3) ; ...
  %    h(3,1) h(3,2) h(3,1)^2+h(3,2)^2 h(3,3) ; ...
  %    ]);
  %  if f ~= (d<0)
  %    fprintf('%d %d\n',f,d<0);
  %    pause
  %  end
  %  %if any(l([p0 p1 p2]) > 1e3)
  %  %%if all((l([p0 p1 p2])>1e3)==[0;1;0]) && f
  %  %  fprintf('%d %d %d --> %d\n',l([p0 p1 p2]) > 1e5,f);

  %  %  clf;
  %  %  hold on;
  %  %  tsurf([1 2 3;1 3 4],[0 0;P([p0 p1 p2],:)],'FaceColor','none','LineWidth',2,'EdgeColor','b');
  %  %  tsurf([1 2 4;1 4 2],[0 0;P([p0 p1 p2],:)],'FaceColor','none','LineWidth',1,'EdgeColor','r');
  %  %  scatter([0;P([p0 p1 p2],1)],[0;P([p0 p1 p2],2)],'.','SizeData',500);
  %  %  text([0;P([p0 p1 p2],1)],[0;P([p0 p1 p2],2)],['q';num2str((0:2)')],'FontSize',20);
  %  %  hold off;
  %  %  axis equal;
  %  %  axis([-3 3 -3 3]);
  %  %  drawnow;

  %  %  pause
  %  %end
  %end

  % Subtract q from P and q
  P = bsxfun(@minus,P,q);
  q = [0 0];
  % Compute distances to origin
  l = sqrt(sum(P.^2,2));

  % Handle case where q is on the convex hull.
  % TODO: Shouldn't _always_ have to added bounding points.
  max_l = max(l);
  % This is HUGE number is a hack. Should handle boundary "protectors" explicitly then.
  % Something like angle between:
  %   finite-∞-finite --> 0
  %   ∞-finite-∞ --> pi
  %   finite-finite-∞ --> ?
  % Maybe this is a case for homogenous coordinates?
  s = 1e10*max_l;
  extra = [size(P,1)+(1:4)];
  P = [P;s*[1 1;1 -1;-1 -1;-1 1]];
  l = [l;repmat(s*sqrt(2),4,1)];
  assert(size(P,1)==size(l,1));

  % Number of input points
  n = size(P,1);
  % Compute angle with x-axis
  th = atan2(P(:,2),P(:,1));
  [th,I] = sort(th);
  [~,closest] = min(l);

  % Reorder so that closest comes first
  I = I([find(I==closest):end 1:find(I==closest)-1]);
  % Assumption: closest point must be Delaunay neighbor
  N = [I(1)];
  % Consider next two points
  p1i = 2;
  p2i = 3;
  p1 = I(p1i);
  p2 = I(p2i);
  p0 = N(end);
  k = 0;
  while true
    k = k+1;
    % is the edge q-p1 between p0 and p2 is delaunay
    if isdel(P,l,p0,p1,p2)
      % push p1 onto "keepers stack" N
      N = [N;p1];
    else
      % edge q-p1 between p0 and p2 is **not** delaunay
      % first try to work backward until finding delaunay
      while numel(N)>1 && ~isdel(P,l,N(end-1),N(end),p2)
        % MOVE BACKWARD
        k = k+1;
        p0 = N(end-1);
        p1 = N(end);
        N = N(1:end-1);
      end
    end
    if p1i == n
      break;
    end
    % MOVE FORWARD
    p1i = mod(p1i,n)+1;
    p2i = mod(p2i,n)+1;
    p1 = I(p1i);
    p2 = I(p2i);
    p0 = N(end);
  end
  N = setdiff(N,size(P,1)-(0:3));
  P = P(1:end-4,:);
end

author: alecjacobson

repo: gptoolbox


text:
function diff = hist_dist( hist1, hist2, method )
    switch method
        case 'x2'
            % diff = 0.5 * sum((hist1 - hist2).^2) ./ sum(hist1 + hist2 + eps);
            diff = 0.5 * sum( (hist1 - hist2).^2 ./ (hist1 + hist2 + eps) );
        case 'jsd'      % Jensen-Shannon Divergence
            diff = 0.5*(sum(hist1.*log((hist1+eps)./(hist2+eps))) + sum(hist2.*log((hist2+eps)./(hist1+eps))));
        otherwise
            error( 'unknown type for computing histogram distance' );
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [q,qd,q2d] = mixed_traj(t,C,A,B,w,N)
% ---------------------------------------------------------------------
% This function computes "mixed trajectory" meaning the trajectory
% that consistes of finite fourier series and fifth order polynomial
% Inputs:
%   t - time instants at which trajectory should be evaluated
%   C - coefficients of the fifth order polynomail
%   A - coeffincients of the sine in finite fourier series
%   B - coefficinets of the cosine in finitne fourier series
%   w - fundamental frequency
%   N - number of harmonics
% ---------------------------------------------------------------------
% finite fourier series
[qh,qhd,qh2d] = fourier_series_traj(t,zeros(6,1),A,B,w,N);

% fifth order polynomail trajectory
qp = C(:,1) + C(:,2).*t + C(:,3).*t.^2 + C(:,4).*t.^3 + ...
                                    C(:,5).*t.^4 + C(:,6).*t.^5;
qpd = C(:,2) + 2*C(:,3).*t + 3*C(:,4).*t.^2 + ...
                                4*C(:,5).*t.^3 + 5*C(:,6).*t.^4;
qp2d = 2*C(:,3) + 6*C(:,4).*t + 12*C(:,5).*t.^2 + 20*C(:,6).*t.^3;

q = qh + qp;
qd = qhd + qpd;
q2d = qh2d + qp2d;

author: shamilmamedov


text:
% Compute the relative perceptual saliency (PS)
function E = Relative_PS(imgIR, imgTV)

E = PS(imgIR)/PS(imgTV);

end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [ u, l ] = vand2_inverse_ul ( n, x )

%*****************************************************************************80
%
%% VAND2_INVERSE_UL returns the UL factors of the VAND2 inverse matrix.
%
%  Discussion:
%
%    inverse ( A ) = U * L
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Turner,
%    Inverse of the Vandermonde Matrix with Applications,
%    NASA Technical Note TN D-3547, 1966.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real X(N), the values that define A.
%
%    Output, real U(N,N), L(N,N), the UL factors of inverse(A).
%
  u = zeros ( n, n );

  for i = 1 : n
    u(i,i) = 1.0;
    if ( i == 1 )
      for j = i + 1 : n
        u(i,j) = - u(i,j-1) * x(j-1);
      end
    else
      for j = i + 1 : n
        u(i,j) = u(i-1,j-1) - u(i,j-1) * x(j-1);
      end
    end
  end

  l = zeros ( n, n );

  for i = 1 : n
    l(i,1:i) = 1.0;
    for j = 1 : i
      for k = 1 : i
        if ( j ~= k )
          l(i,j) = l(i,j) / ( x(j) - x(k) );
        end
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function g = robOneDynamicsLogLikeGradients(model)

% ROBONEDYNAMICSLOGLIKEGRADIENTS Gradients of the robot one dynamics wrt parameters.

% FGPLVM

g = [];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function matrix2html(matrix, filename, varargin)
% matrix2html Generates a HTML table for a given matrix
%
% The function generates a HTML table based on a given matrix or cell array.
%
% Input
% - matrix (matrix): A 2 dimensional numerical or cell array
% - filename (string): A valid filename or a file handle
% - varargs[RowLabels] (array): Can be used to label the rows of the table
% - varargs[ColumnLabels] (array): Can be used to label the columns of the table
% - varargs[Format] (string): Format of the numeric input data in the printf format.
% - varargs[Class] (string): CSS class of the table.
% - varargs[Title] (string): Title of the table.
% - varargs[EmbedNumers] (boolean): Embed numerical data as data attributes.
%

    rowLabels = [];
    colLabels = [];
    format = '%.4f';
    css_class = [];
    title = [];
    embed_numbers = true;
    if (rem(nargin,2) == 1 || nargin < 2)
        error('Incorrect number of arguments to %s.', mfilename);
    end

    for i = 1:2:length(varargin)
        switch lower(varargin{i})
            case 'rowlabels'
                rowLabels = varargin{i+1};
            case 'columnlabels'
                colLabels = varargin{i+1};
            case 'format'
                format = lower(varargin{i+1});
            case 'class'
                css_class = [css_class, ' ', varargin{i+1}];
            case 'title'
                title = varargin{i+1};
            case 'embednumbers'
                embed_numbers = varargin{i+1};
            otherwise
                error(['Unknown switch ', varargin{i}, '!']) ;
        end
    end

    rowLabels = matrix2cells(rowLabels, 'th', format, false);
    colLabels = matrix2cells(colLabels, 'th', format, false);
    matrix = matrix2cells(matrix, 'td', format, embed_numbers);

    if (ischar(filename))
        fid = fopen(filename, 'w');
        close_file = 1;
    else
        fid = filename;
        close_file = 0;
    end;

    if ~isempty(colLabels) && ~isempty(rowLabels)

        stub = cell(size(colLabels, 1), size(rowLabels, 2));
        stub{1, 1} = sprintf(' 	', ...
                size(rowLabels, 2), size(colLabels, 1));

        head = cat(2, stub, colLabels);

        body = cat(2, rowLabels, matrix);

    elseif ~isempty(colLabels)

        head = colLabels;

        body = matrix;

    elseif ~isempty(rowLabels)

        head = [];
        body = cat(2, rowLabels, matrix);

    end;

    if(~isempty(title))
        fprintf(fid, '
%s
\n', title);
    end

    if(~isempty(css_class))
        fprintf(fid, '
\n', css_class);
    else
        fprintf(fid, '
\n');
    end

    if ~isempty(head)

        skip = cellfun(@(x) isempty(x), head, 'UniformOutput', true);
        fprintf(fid, '
');
        for i = 1 : size(head, 1)
            fprintf(fid, '
');
            cellfun(@(x) fwrite(fid, x), head(i, ~skip(i, :)), 'UniformOutput', true);
            fprintf(fid, '
');
        end;
        fprintf(fid, '
');
    end;

    skip = cellfun(@(x) isempty(x), body, 'UniformOutput', true);

    fprintf(fid, '
');
    for i = 1 : size(body, 1)
        fprintf(fid, '
');
        cellfun(@(x) fwrite(fid, x), body(i, ~skip(i, :)), 'UniformOutput', true);
        fprintf(fid, '
');
    end;
    fprintf(fid, '
');

    fprintf(fid, '
\r\n');

    if (close_file)
        fclose(fid);
    end;

end

function [cells] = matrix2cells(matrix, element, format, embed_numbers)

    if isnumeric(matrix)
        matrix = num2cell(matrix);
    end;
    [cells, rowspan, colspan, attributes] = cellfun(@(x) cell2cell(x, format, embed_numbers), ...
        matrix, 'UniformOutput', false);
    rowspan = cell2mat(rowspan);
    colspan = cell2mat(colspan);

    width = size(matrix, 2);
    height = size(matrix, 1);

    for h=1:height
        for w=1:width
            if rowspan(h, w) == 0 || colspan(h, w) == 0
                cells{h, w} = '';
            elseif rowspan(h, w) == 1 && colspan(h, w) == 1
                cells{h, w} = sprintf('<%s %s>%s</%s>', element, ...
                    attributes{h, w}, cells{h, w}, element);
            elseif rowspan(h, w) == 1
                cells{h, w} = sprintf('<%s colspan="%d" %s>%s</%s>', ...
                    element, colspan(h, w), attributes{h, w}, cells{h, w}, element);
            elseif colspan(h, w) == 1
                cells{h, w} = sprintf('<%s rowspan="%d" %s>%s</%s>', ...
                    element, rowspan(h, w), attributes{h, w}, cells{h, w}, element);
            else
                cells{h, w} = sprintf('<%s colspan="%d" rowspan="%d" %s>%s</%s>', ...
                    element, colspan(h, w), rowspan(h, w), attributes{h, w}, cells{h, w}, element);
            end;
        end
    end

%     max_width = max(sum(colspan));
%     max_height = max(sum(rowspan));
%
%     if max_width < width
%         cells = cells(:, 1:max_width);
%     end
%
%     if max_height < height
%         cells = cells(:, 1:max_height);
%     end

end

function [str, rowspan, colspan, cell_attributes] = cell2cell(s, format, embed_numbers)

    rowspan = 1;
    colspan = 1;

    cell_attributes = '';

    if isnumeric(s)
        str = num2str(s, format);
        if embed_numbers
            cell_attributes = [cell_attributes, sprintf(' data-numeric="%f" ', s)];
        end;
        return;
    elseif ~isstruct(s)
        str = s;
        return;
    end;

    if ~isfield(s, 'text')
        rowspan = 0;
        colspan = 0;
        str = [];
        return;
    end;

    if isfield(s, 'format')
        format = s.format;
    end;

    if isfield(s, 'columns')
        colspan = s.columns;
    end

    if isfield(s, 'rows')
        rowspan = s.rows;
    end

    if isnumeric(s.text)
        if embed_numbers
            cell_attributes = [cell_attributes, sprintf(' data-numeric="%f" ', s.text)];
        end;
        str = num2str(s.text, format);
    else
        str = s.text;
    end;

    element = 'span';
    attributes = '';

    if isfield(s, 'class')
        cell_attributes = [cell_attributes, ' class="', s.class, '" '];
    end

    if isfield(s, 'tooltip')
        attributes = [attributes, ' title="', s.tooltip, '" '];
    end

    if isfield(s, 'url')
        element = 'a';
        attributes = [attributes, ' href="', s.url, '" '];
    end

    if ~isempty(attributes)
       str = sprintf('<%s %s>%s</%s>', element, attributes, str, element);
    end

end

author: votchallenge

repo: toolkit-legacy


text:
function [xc,good,bad,type] = cornerfinder_saddle_point(xt,I,wintx,winty,wx2,wy2);

%[xc] = cornerfinder_saddle_point(xt,I,wintx,winty);
%
%Finds the sub-pixel corners on the image I with initial guess xt
%xt and xc are 2xN matrices. The first component is the x coordinate
%(horizontal) and the second component is the y coordinate (vertical)
%
%Based on Harris corner finder method
%
%Finds corners to a precision below .1 pixel!
%Oct. 14th, 1997 - UPDATED to work with vertical and horizontal edges as well!!!
%Sept 1998 - UPDATED to handle diverged points: we keep the original points
%good is a binary vector indicating wether a feature point has been properly
%found.
%
%Add a zero zone of size wx2,wy2
%July 15th, 1999 - Bug on the mask building... fixed + change to Gaussian mask with higher
%resolution and larger number of iterations.

% California Institute of Technology
% (c) Jean-Yves Bouguet -- Oct. 14th, 1997

xt = xt';
xt = fliplr(xt);

if nargin < 4,
   winty = 5;
   if nargin < 3,
      wintx = 5;
   end;
end;

if nargin < 6,
   wx2 = -1;
   wy2 = -1;
end;

%mask = ones(2*wintx+1,2*winty+1);
mask = exp(-((-wintx:wintx)'/(wintx)).^2) * exp(-((-winty:winty)/(winty)).^2);

if (wx2>0) & (wy2>0),
   if ((wintx - wx2)>=2)&((winty - wy2)>=2),
      mask(wintx+1-wx2:wintx+1+wx2,winty+1-wy2:winty+1+wy2)= zeros(2*wx2+1,2*wy2+1);
   end;
end;

offx = [-wintx:wintx]'*ones(1,2*winty+1);
offy = ones(2*wintx+1,1)*[-winty:winty];

resolution = 0.005;

MaxIter = 10;

[nx,ny] = size(I);
N = size(xt,1);

xc = xt; % first guess... they don't move !!!

type = zeros(1,N);

for i=1:N,

   v_extra = resolution + 1; 		% just larger than resolution

   compt = 0; 				% no iteration yet

   while (norm(v_extra) > resolution) & (compt<MaxIter),

       cIx = xc(i,1); 			%
       cIy = xc(i,2); 			% Coords. of the point
       crIx = round(cIx); 		% on the initial image
       crIy = round(cIy); 		%
       itIx = cIx - crIx; 		% Coefficients
       itIy = cIy - crIy; 		% to compute
       if itIx > 0, 			% the sub pixel
           vIx = [itIx 1-itIx 0]'; 	% accuracy.
       else
           vIx = [0 1+itIx -itIx]';
       end;
       if itIy > 0,
           vIy = [itIy 1-itIy 0];
       else
           vIy = [0 1+itIy -itIy];
       end;

      % What if the sub image is not in?

      if (crIx-wintx-2 < 1), xmin=1; xmax = 2*wintx+5;
      elseif (crIx+wintx+2 > nx), xmax = nx; xmin = nx-2*wintx-4;
      else
          xmin = crIx-wintx-2; xmax = crIx+wintx+2;
      end;

      if (crIy-winty-2 < 1), ymin=1; ymax = 2*winty+5;
      elseif (crIy+winty+2 > ny), ymax = ny; ymin = ny-2*winty-4;
      else
          ymin = crIy-winty-2; ymax = crIy+winty+2;
      end;

      SI = I(xmin:xmax,ymin:ymax); % The necessary neighborhood
      SI = conv2(conv2(SI,vIx,'same'),vIy,'same');
      SI = SI(2:2*wintx+4,2:2*winty+4); % The subpixel interpolated neighborhood

      px = cIx + offx;
      py = cIy + offy;

      if 1, %~saddle,
          [gy,gx] = gradient(SI); 		% The gradient image
          gx = gx(2:2*wintx+2,2:2*winty+2); % extraction of the useful parts only
          gy = gy(2:2*wintx+2,2:2*winty+2); % of the gradients
          gxx = gx .* gx .* mask;
          gyy = gy .* gy .* mask;
          gxy = gx .* gy .* mask;

          bb = [sum(sum(gxx .* px + gxy .* py)); sum(sum(gxy .* px + gyy .* py))];

          a = sum(sum(gxx));
          b = sum(sum(gxy));
          c = sum(sum(gyy));

          dt = a*c - b^2;

          xc2 = [c*bb(1)-b*bb(2) a*bb(2)-b*bb(1)]/dt;
      else

          SI = SI(2:2*wintx+2,2:2*winty+2);
          A =  repmat(mask(:),1,6) .* [px(:).^2 px(:).*py(:) py(:).^2 px(:) py(:) ones((2*wintx+1)*(2*winty+1),1)];
          param = inv(A'*A)*A'*( mask(:).*SI(:));
          xc2 = (-inv([2*param(1) param(2) ; param(2) 2*param(3) ]) * param(4:5))';

      end;

      v_extra = xc(i,:) - xc2;

      xc(i,:) = xc2;

      compt = compt + 1;

  end;

  if 1,

      cIx = xc(i,1); 			%
      cIy = xc(i,2); 			% Coords. of the point
      crIx = round(cIx); 		% on the initial image
      crIy = round(cIy); 		%
      itIx = cIx - crIx; 		% Coefficients
      itIy = cIy - crIy; 		% to compute
      if itIx > 0, 			% the sub pixel
          vIx = [itIx 1-itIx 0]'; 	% accuracy.
      else
          vIx = [0 1+itIx -itIx]';
      end;
      if itIy > 0,
          vIy = [itIy 1-itIy 0];
      else
          vIy = [0 1+itIy -itIy];
      end;

      % What if the sub image is not in?

      if (crIx-wintx-2 < 1), xmin=1; xmax = 2*wintx+5;
      elseif (crIx+wintx+2 > nx), xmax = nx; xmin = nx-2*wintx-4;
      else
          xmin = crIx-wintx-2; xmax = crIx+wintx+2;
      end;

      if (crIy-winty-2 < 1), ymin=1; ymax = 2*winty+5;
      elseif (crIy+winty+2 > ny), ymax = ny; ymin = ny-2*winty-4;
      else
          ymin = crIy-winty-2; ymax = crIy+winty+2;
      end;

      SI = I(xmin:xmax,ymin:ymax); % The necessary neighborhood
      SI = conv2(conv2(SI,vIx,'same'),vIy,'same');
      SI = SI(2:2*wintx+4,2:2*winty+4); % The subpixel interpolated neighborhood
      px = cIx + offx;
      py = cIy + offy;

      SI = SI(2:2*wintx+2,2:2*winty+2);
      A =  repmat(mask(:),1,6) .* [px(:).^2 px(:).*py(:) py(:).^2 px(:) py(:) ones((2*wintx+1)*(2*winty+1),1)];
      param = inv(A'*A)*A'*( mask(:).*SI(:));
      xc2 = (-inv([2*param(1) param(2) ; param(2) 2*param(3) ]) * param(4:5))';

      v_extra = xc(i,:) - xc2;

      xc(i,:) = xc2;
  end;

end;

% check for points that diverge:

delta_x = xc(:,1) - xt(:,1);
delta_y = xc(:,2) - xt(:,2);

%keyboard;

bad = (abs(delta_x) > wintx) | (abs(delta_y) > winty);
good = ~bad;
in_bad = find(bad);

% For the diverged points, keep the original guesses:

xc(in_bad,:) = xt(in_bad,:);

xc = fliplr(xc);
xc = xc';

bad = bad';
good = good';

repo: ekfmonoslam


function [sims_with_endopath,EPSn] = cforecasts(endo_path,endo_index,forecast_data,Phi,Sigma)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'cforecasts' performs the forecast conditional on a path for one or more
% endogenous variables using all the shocks of the VAR
% references: Waggoner and Zha (1999) and Jarocinski (2010)

% Inputs:
% - endo_path, path for the endogenous variables
% - endo_index, order of the endogenous variables with a contional path
% - forecast_data, last data
% - Phi, AR parameters of the VAR
% - Sigma, Covariance matrix of the reduced form VAR shocks

% Output:
% - sims_with_endopath
% 1st dimension:   horizon
% 2nd dimension:   variable
% - EPSn, shocks generating the path

% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018
% Revised, 9/11/2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[fhor,Ncondvar] = size(endo_path);
Nvar            = size(Sigma,2);

if size(Ncondvar,2) ~= length(endo_index)
    error('something went wrong: the number of path and the number of vars are not consistent.')
end
if size(endo_index,1) > size(endo_index,2)
    error('the index needs to be a row vector');
end
Nres = size(endo_index,2)*fhor;

[F,~,~,const,~,lags]     = var2ss(Phi,Sigma);

% endo_path_deviation = endo_path - repmat(const',1,fhor);

% no shocks forecasts
[sims_no_shock,~]   = forecasts(forecast_data,Phi,Sigma,fhor,lags);

% restrictions
% err                 = endo_path_deviation - sims_no_shock(:,endo_path_index);
% err = zeros(Nvar*fhor);
err = endo_path - sims_no_shock(:,endo_index);
err = reshape(err',size(endo_index,1)*fhor,1);

% 1 standard deviation increase (if unit = eye(N))
[C] = chol(Sigma,'lower');

R    = zeros(Nvar * fhor);
tmp0 = [];
for ff = 1 :fhor
    tmp  = F^(ff-1);
    tmp0 = [tmp(1:Nvar,1:Nvar) * C tmp0];
    R((ff-1)*Nvar + 1 : (ff)*Nvar, 1: size(tmp0,2) ) = tmp0;
    index(ff,:) = (ff-1)*Nvar + endo_index;
end
index = reshape(index',Nres,1);
% index= [4 10];
Rtilde = R(index,:);

[U,D,V] = svd(Rtilde);
V1 = V(:,1:Nres);
V2 = V(:,Nres+1:end);

eps = V1*inv(D(1:Nres,1:Nres))*U'*err + V2 * randn(size(R,1)-Nres,1);

% orthogonal schoks
EPSn = reshape(eps,Nvar,fhor);
% reduced form shocks
EPS  = C*EPSn;

[~,sims_with_endopath] = forecasts(forecast_data,Phi,Sigma,fhor,lags,EPS');


text:
function runAllTests(varargin)
%RUNALLTESTS  Explore all subdirectories, and run 'runtests' function
%
%   output = runAllTests(input)
%
%   Example
%   runAllTests
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-07-20,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% Add matGeom to path
addpath(genpath(strrep(mfilename('fullpath'),['tests' filesep mfilename],'matGeom')))

directories = {...
    'geom2d', ...
    'polygons2d', ...
    'graphs', ...
    'geom3d', ...
    'meshes3d', ...
    };

for i = 1:length(directories)
    name = directories{i};
    disp(['Running tests for directory: ' name]);

    cd(name);
    runtests;
    cd('..');
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function mean = rayleigh_mean ( a )

%*****************************************************************************80
%
%% RAYLEIGH_MEAN returns the mean of the Rayleigh PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the PDF.
%    0.0 < A.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = a * sqrt ( 0.5 * pi );

  return
end


text:
function urband(varargin)

	query = reshape(transpose([strvcat(varargin{:}),[repmat('+',size(strvcat(varargin{:}),1)-1,1);' ']]),1,[]);

	url = ['''http://www.urbandictionary.com/define.php?term=',query,''''];

	url(url==' ') = [];

	eval(['web ', url ,' -browser'])

end

author: Sable

repo: mcbench-benchmarks


text:
% ir_mex_build_table
% run matlab's "mex" command to "compile" the table interpolation code
% into mex files.

dir_current = pwd;
dir_nufft = path_find_dir('nufft');
dir_table = [dir_nufft filesep 'table'];
cd(dir_table)

% compile mex files using ir_mex_fun2 ("2" since there are two arguments)
fun = ir_mex_fun();
fun('interp1_table_adj_mex.c',	'interp1_table1_adj.c')
fun('interp1_table_mex.c',	'interp1_table1_for.c')
fun('interp2_table_adj_mex.c',	'interp2_table1_adj.c')
fun('interp2_table_mex.c',	'interp2_table1_for.c')
fun('interp3_table_adj_mex.c',	'interp3_table1_adj.c')
fun('interp3_table_mex.c',	'interp3_table1_for.c')

cd(dir_current)

% "test" mex files by running them (each should display usage)
try
	interp1_table_adj_mex
	interp1_table_mex
    interp2_table_adj_mex
    interp2_table_mex
    interp3_table_adj_mex
    interp3_table_mex
catch
    disp(['WARNING! Error compiling mex files in ', mfilename])
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [tSpan,x,u] = unPackDecVar(z,dim)
% function [tSpan,x,u] = unPackDecVar(z,dim)
%
% This function unpacks the decision variables for
% trajectory optimization into the duration (t),
% state (x), and control (u) matricies
%
% INPUTS:
%   z = [2+ns*ms+nc*nc, 1] = vector of decision variables
%   dim = struct with the size of the state and control matricies
%       .nState = size(x);
%       .nControl = size(u);
%
% OUTPUTS:
%   tSpan = [1, 2] = [t0, tF] = time span
%   x = [ns, ms] = state matrix = [dimension in state space, grid point]
%   u = [nc, mc] = control matrix = [dimension in control space, grid point]
%
% See Also: PACKDECVAR

ns = dim.nState;  % ns = size(x)
nc = dim.nControl; % nc = size(u)

is = 1:(ns(1)*ns(2));
ic = 1:(nc(1)*nc(2));

tSpan = [z(1), z(2)];
x = reshape(z(2+is),ns(1),ns(2));
u = reshape(z(2+is(end)+ic), nc(1),nc(2));

end

author: MatthewPeterKelly

repo: dscTutorials


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

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

% Last Modified by GUIDE v2.5 13-Mar-2004 18:07:02

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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

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

% Update all temperature values
deg_f = str2num( get(hObject,'String') );
deg_f = max( [ 32 deg_f] );
deg_f = min( [212 deg_f] );
deg_c = to_c(deg_f);

% Now update the fields
set (handles.Edit1,'String',sprintf('%.1f',deg_f));
set (handles.Edit2,'String',sprintf('%.1f',deg_c));
set (handles.Slider1,'Value',deg_c);

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

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

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

% Update all temperature values
deg_c = str2num( get(hObject,'String') );
deg_c = max( [  0 deg_c] );
deg_c = min( [100 deg_c] );
deg_f = to_f(deg_c);

% Now update the fields
set (handles.Edit1,'String',sprintf('%.1f',deg_f));
set (handles.Edit2,'String',sprintf('%.1f',deg_c));
set (handles.Slider1,'Value',deg_c);

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

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

% --- Executes on slider movement.
function Slider1_Callback(hObject, eventdata, handles)
% hObject    handle to Slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Update all temperature values
deg_c = get(hObject,'Value');
deg_f = to_f(deg_c);

% Now update the fields
set (handles.Edit1,'String',sprintf('%.1f',deg_f));
set (handles.Edit2,'String',sprintf('%.1f',deg_c));
set (handles.Slider1,'Value',deg_c);

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

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


text:
function [ y, m, d, ierror ] = ymd_check_roman ( y, m, d )

%*****************************************************************************80
%
%% YMD_CHECK_ROMAN checks a Roman YMD date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, the YMD date, which may
%    be corrected if necessary and possible.
%
%    Output, integer IERROR, is 0 if the date is legal.
%
  ierror = 0;
%
%  Check the year.
%
  if ( y <= 0 )
    ierror = 1;
    return
  end
%
%  Check the month.
%
  [ y, m ] = month_borrow_roman ( y, m );

  [ y, m ] = month_carry_roman ( y, m );
%
%  Check the day.
%
  [ y, m, d ] = day_borrow_roman ( y, m, d );

  [ y, m, d ] = day_carry_roman ( y, m, d );

  return
end


text:
function v = lift_catreg_v2(r,t,c,theta,cBound,tBound)
%% Lifting for category registration
%% Heng Yang, July 05, 2021
r           = r(:);
t           = t(:);
c           = c(:);
K           = length(c);
theta       = theta(:);
cBoundSq    = cBound^2;
tBoundSq    = tBound^2;

x           = [r;t;c];
v1          = [1;x;theta;kron(theta,x)];

if tBoundSq < t'*t
    v2      = [1;theta] * 0;
else
    v2      = [1;theta] * sqrt(tBoundSq - t'*t);
end

if cBoundSq < c'*c
    tmp     = [1;theta] * 0;
else
    tmp     = [1;theta] * sqrt(cBoundSq - c'*c);
end

v           = {v1;v2;tmp};

for k = 1:K
    if c(k) < 0
        tmp = [1;theta] * 0;
    else
        tmp = [1;theta] * sqrt(c(k));
    end
    v       = [v;{tmp}];
end

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function [prob,sol,fmin] = nls_prob(varargin)
%NLS_PROB  Return an OPTI NLS
%
%   prob = nls_prob(no) return a pre-built optiprob of a saved NLS.
%
%   [prob,sol,fmin] = NLS_prob(no) returns the optimum solution and function
%   eval at the optimum
%
%   no = nls_prob() returns the number of problems available for testing.

%   (C) 2012 Jonathan Currie (IPL)

% Functions are taken from:
% More, J. J., B. S. Garbow, and K. E. Hillstrom. "Testing Unconstrained
% Optimization Software." ACM Transactions on Mathematical Software 7,
% no. 1 (1981): 17-41.

%Check if just returning no problems
if(nargin < 1)
    prob = 20; sol = []; fmin = [];
    return;
else
    no = varargin{1};
end

%Big switch yard
switch(no)
    case 1
        fun = @(x) [10*(x(2)-x(1)^2); 1 - x(1)];
        ydata = zeros(2,1);
        x0 = [-1.2;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1,1];
        fmin = 0;

    case 2
        fun = @(x) [-13 + x(1) + ((5 - x(2))*x(2) - 2)*x(2);
                    -29 + x(1) + ((x(2) + 1)*x(2) - 14)*x(2)];
        ydata = zeros(2,1);
        x0 = [0.5;-2];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [5,4];
        fmin = 0;

    case 3
        fun = @(x) [1e4*x(1)*x(2) - 1;
                    exp(-x(1)) + exp(-x(2)) - 1.0001];
        ydata = zeros(2,1);
        x0 = [0;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1.098e-5;9.106];
        fmin = 0;

    case 4
        fun = @(x) [x(1) - 1e6; x(2) - 2e-6; x(1)*x(2) - 2];
        ydata = zeros(3,1);
        x0 = [1;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1e6,2e-6];
        fmin = 0;

    case 5
        fun = @(x) [1.5 - x(1)*(1 - x(2)); 2.25 - x(1)*(1 - x(2)^2); 2.625 - x(1)*(1 - x(2)^3) ];
        ydata = zeros(3,1);
        x0 = [1;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [3;0.5];
        fmin = 0;

    case 6
        i = (1:10)';
        fun = @(x) 2 + 2*i - (exp(x(1)*i) + exp(x(2)*i));
        ydata = zeros(10,1);
        x0 = [0.3;0.4];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [0.2578;0.2578];
        fmin = 124.36226865;

    case 7
        fun = @more7;
        ydata = zeros(3,1);
        x0 = [-1;0;0];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1;0;0];
        fmin = 0;

    case 8
        u = (1:15)';
        v = 16 - u;
        w = zeros(15,1);
        for i = 1:15
            w(i) = min(u(i),v(i));
        end
        y = [0.14;0.18;0.22;0.25;0.29;0.32;0.35;0.39;0.37;0.58;0.73;0.96;1.34;2.10;4.39];
        fun = @(x) (x(1) + u ./ (v*x(2) + w*x(3)));
        x0 = [1;1;1];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 8.21487e-3;

    case 9
        i = (1:15)';
        y = [0.0009;0.0044;0.0175;0.0540;0.1295;0.2420;0.3521;0.3989;0.3521;0.2420;0.1295;0.0540;0.0175;0.0044;0.0009];
        t = (8 - i)./2;
        fun = @(x) x(1)*exp((-x(2)*(t - x(3)).^2)./2);
        x0 = [0.4;1;0];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 1.12793e-8;

    case {10,11}
        i = (1:16)';
        y = [34780;28610;23650;19630;16370;13720;11540;9744;8261;7030;6005;5147;4427;3820;3307;2872];
        t = 45 + 5*i;
        fun = @(x) x(1)*exp(x(2)./(t + x(3)));
        x0 = [0.02;4000;250];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 87.9458;

%     case 11 %not working??
%         m = 100;
%         i = [1:m]';
%         t = i./100;
%         y = 25 + (-50 * log(t)).^(2/3);
%         fun = @(x) exp(-(abs(y*100.*i*x(2)).^x(3))./x(1)) - t;
%         x0 = [5;2.5;0.15];
%         ydata = zeros(m,1);
%         prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
%         sol = [50;25;1.5];
%         fmin = 0;

    case 12
        m = 3;
        i = (1:m)';
        t = 0.1*i;
        fun = @(x) exp(-t*x(1)) - exp(-t*x(2)) - x(3)*exp(-t) - exp(-10*t);
        x0 = [0;10;20];
        ydata = zeros(m,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1;10;1];
        fmin = 0;

    case 13
        fun = @(x) [x(1) + 10*x(2);
                    sqrt(5)*(x(3) - x(4));
                    (x(3) - 2*x(3))^2;
                    sqrt(10)*(x(1)-x(4))^2];
        x0 = [3;-1;0;1];
        ydata = zeros(4,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [0;0;0;0];
        fmin = 0;

    case 14
        fun = @(x) [10*(x(2)-x(1)^2);
                    1 - x(1);
                    sqrt(90)*(x(4)-x(3)^2);
                    1 - x(3);
                    sqrt(10)*(x(2) + x(4) - 2);
                    10^(-0.5)*(x(2)-x(4))];
        x0 = [3;-1;-3;-1];
        ydata = zeros(6,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1;1;1;1];
        fmin = 0;

    case 15
        y = [0.1957;0.1947;0.1735;0.16;0.0844;0.0627;0.0456;0.0342;0.0323;0.0235;0.0246];
        u = [4;2;1;0.5;0.25;0.167;0.125;0.1;0.0833;0.0714;0.0625];
        fun = @(x) (x(1)*(u.^2 + u*x(2))) ./ (u.^2 + u*x(3) + x(4));
        x0 = [0.25;0.39;0.415;0.39];
        prob = optiprob('fun',fun,'ydata',-y,'x0',x0);
        sol = [];
        fmin = 3.07505e-4;

    case 16
        m = 20;
        i = (1:m)';
        t = i/5;
        fun = @(x) (x(1) + t*x(2) - exp(t)).^2 + (x(3) + x(4)*sin(t) - cos(t)).^2;
        x0 = [25;5;-5;-1];
        ydata = zeros(m,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [];
        fmin = 85822.2016;

    case 17
        m = 33;
        i = (1:m)';
        t = 10*(i - 1);
        y = [0.844;0.908;0.932;0.936;0.925;0.908;0.881;0.850;0.818;0.784;0.751;0.718;0.685;0.658;0.628;0.603;...
             0.58;0.558;0.538;0.522;0.506;0.490;0.478;0.467;0.457;0.448;0.438;0.431;0.424;0.42;0.414;0.411;0.406];
        fun = @(x) (x(1) + x(2)*exp(-t*x(4)) + x(3)*exp(-t*x(5)));
        x0 = [0.5;1.5;-1;0.01;0.02];
        prob = optiprob('fun',fun,'ydata',-y,'x0',x0);
        sol = [];
        fmin = 5.46489e-5;

    case 18
        m = 13;
        i = (1:m)';
        t = 0.1*i;
        y = exp(-t) - 5*exp(-10*t) + 3*exp(-4*t);
        fun = @(x) x(3)*exp(-t*x(1)) - x(4)*exp(-t*x(2)) + x(6)*exp(-t*x(5));
        x0 = [1;2;1;1;1;1];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 0; %5.65565e-3; seems wrong?

    case 19
        m = 65;
        i = (1:m)';
        t = (i-1)/10;
        y = [1.366;1.191;1.112;1.013;0.991;0.885;0.831;0.847;0.786;0.725;0.746;0.679;0.608;0.655;0.616;0.606;0.602;0.626;0.651;0.724;0.649;0.649;...
            0.694;0.644;0.624;0.661;0.612;0.558;0.533;0.495;0.5;0.423;0.395;0.375;0.372;0.391;0.396;0.405;0.428;0.429;0.523;0.562;0.607;0.653;...
            0.672;0.708;0.633;0.668;0.645;0.632;0.591;0.559;0.597;0.625;0.739;0.710;0.729;0.720;0.636;0.581;0.428;0.292;0.162;0.098;0.054];
        fun = @(x) (x(1)*exp(-t*x(5)) + x(2)*exp(-(t-x(9)).^2*x(6)) + x(3)*exp(-(t-x(10)).^2*x(7)) + x(4)*exp(-(t-x(11)).^2*x(8)));
        x0 = [1.3;0.65;0.65;0.7;0.6;3;5;7;2;4.5;5.5];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 4.01377e-2;

    case 20
        m = 31; n = 12;
        fun = @more20;
        x0 = zeros(n,1);
        ydata = zeros(m,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [];
        fmin = 4.72238e-10;

    otherwise
        error('Problem not available or not implemented yet');
end

function j = more7(x)

if(x(1) > 0)
    theta = 1/(2*pi)*atan(x(2)/x(1));
else
    theta = 1/(2*pi)*atan(x(2)/x(1)) + 0.5;
end

j(1,1) = 10*(x(3) - 10*theta);
j(2,1) = 10*(x(1)^2 + x(2)^2)^0.5 - 1;
j(3,1) = x(3);

function f = more20(x)
n = length(x);
j = (2:n)';
j2 = (1:n)';
t = 0.1*(1:29)';
f = zeros(31,1);
for k = 1:29
    f(k) = sum( (j-1).*x(2:end).*t(k).^(j-2) ) - sum( x.*t(k).^(j2-1) ) - 1;
end
f(30) = x(1);
f(31) = x(2) - x(1)^2 - 1;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [nodes, edges, faces] = gcontour3d(img)
%GCONTOUR3D Create contour graph of a 3D binary image.
%
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 28/06/2004.
%

nodes = zeros([0 3]);   % 3 coordinates vertices
edges = zeros([0 2]);   % first node and second nodes
faces = zeros([0 4]);   % indices of 4 corners of each square face

D1 = size(img, 1);
D2 = size(img, 2);
D3 = size(img, 3);

% first direction for image
for y = 1:D2
    for z = 1:D3
        % find transitions between the two phases
        ind = find(img(1:D1-1, y, z)~=img(2:D1, y, z));

        % process each transition in direction 1
        for i2 = 1:length(ind)

            % coordinates of each node
            n1 = [ind(i2)+.5 y-.5 z-.5];
            n2 = [ind(i2)+.5 y-.5 z+.5];
            n3 = [ind(i2)+.5 y+.5 z+.5];
            n4 = [ind(i2)+.5 y+.5 z-.5];

            % add the face (and edges) with the 4 given nodes
            [nodes, edges, faces] = addFace(nodes, edges, faces, [n1; n2; n3; n4]);
        end
    end
end

% second direction for image
for x=1:D1
    for z=1:D3
        % find transitions between the two phases
        ind = find(img(x, 1:D2-1, z)~=img(x, 2:D2, z));

        % process each transition in direction 1
        for i2 = 1:length(ind)

            % coordinates of each node
            n1 = [x-.5 ind(i2)+.5 z-.5];
            n2 = [x-.5 ind(i2)+.5 z+.5];
            n3 = [x+.5 ind(i2)+.5 z+.5];
            n4 = [x+.5 ind(i2)+.5 z-.5];

            % add the face (and edges) with the 4 given nodes
            [nodes, edges, faces] = addFace(nodes, edges, faces, [n1; n2; n3; n4]);
        end
    end
end

% third direction for image
for x=1:D1
    for y=1:D2
        % find transitions between the two phases
        ind = find(img(x, y, 1:D3-1)~=img(x, y, 2:D3));

        % process each transition in direction 1
        for i2 = 1:length(ind)

            % coordinates of each node
            n1 = [x-.5 y-.5 ind(i2)+.5];
            n2 = [x-.5 y+.5 ind(i2)+.5];
            n3 = [x+.5 y+.5 ind(i2)+.5];
            n4 = [x+.5 y-.5 ind(i2)+.5];

            % add the face (and edges) with the 4 given nodes
            [nodes, edges, faces] = addFace(nodes, edges, faces, [n1; n2; n3; n4]);
        end
    end
end

return;

function [nodes, edges, faces] = addFace(nodes, edges, faces, faceNodes)
% add given nodes and coresponding face to the graph.

n1 = faceNodes(1,:);
n2 = faceNodes(2,:);
n3 = faceNodes(3,:);
n4 = faceNodes(4,:);

% search indices of each nodes
ind1 = find(ismember(nodes, n1, 'rows'));
ind2 = find(ismember(nodes, n2, 'rows'));
ind3 = find(ismember(nodes, n3, 'rows'));
ind4 = find(ismember(nodes, n4, 'rows'));

% if nodes are not in the list, we add them
if isempty(ind1)
    nodes = [nodes; n1];
    ind1 = size(nodes, 1);
end
if isempty(ind2)
    nodes = [nodes; n2];
    ind2 = size(nodes, 1);
end
if isempty(ind3)
    nodes = [nodes; n3];
    ind3 = size(nodes, 1);
end
if isempty(ind4)
    nodes = [nodes; n4];
    ind4 = size(nodes, 1);
end

% add current face to the list
faces(size(faces, 1)+1, 1:4) = [ind1(1) ind2(1) ind3(1) ind4(1)];

% create edges of the face
% (first index is the smallest one, by convention)
e1 = [min(ind1, ind2) max(ind1, ind2)];
e2 = [min(ind2, ind3) max(ind2, ind3)];
e3 = [min(ind3, ind4) max(ind3, ind4)];
e4 = [min(ind4, ind1) max(ind4, ind1)];

 % if nodes are not in the list, we add them
if isempty(find(ismember(edges, e1, 'rows'), 1))
    edges = [edges; e1];
end
if isempty(find(ismember(edges, e2, 'rows'), 1))
    edges = [edges; e2];
end
if isempty(find(ismember(edges, e3, 'rows'), 1))
    edges = [edges; e3];
end
if isempty(find(ismember(edges, e4, 'rows'), 1))
    edges = [edges; e4];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function C = getcolor01(A)
%------------------------------------------------------------------------------
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: June 6, 1999.
% (c) 1999-2002 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
[n, m] = size(A);
C=A(1:2:n, 2:2:m);
%------------------------------------------------------------------------------


text:
function value = ellipse_area_2d ( r1, r2 )

%*****************************************************************************80
%
%% ELLIPSE_AREA_2D returns the area of an ellipse in 2D.
%
%  Discussion:
%
%    The ellipse is defined as points (X,Y) such that
%
%      ( ( X - XC ) / R1 )^2 + ( ( Y - YC ) / R2 )^2 <= 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the major and minor semi-axes.
%
%    Output, real ELLIPSE_AREA_2D, the area of the ellipse.
%
  value = pi * r1 * r2;

  return
end


text:
function [l,L_f,L_lf] = fromFrameIdpLin(F,lf)

% FROMFRAMEIDPLIN  Transforms IDP line from local frame to global frame.
%   I = FROMFRAMEIDPLIN(F,LF) transforms the Inverse Depth line LF from the
%   local frame F to the global frame. The frame F must be specified via a
%   structure containing at least the fields F.t, F.q, F.R and F.Rt
%   (translation, quaternion, rotation matrix and its transpose).
%
%   [I,L_f,L_lf] = FROMFRAMEIDPLIN(...) returns the Jacobians wrt F and IF.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout == 1
    [p1f,p2f] = idpLin2idpPnts(lf);

    p1 = fromFrameIdp(F,p1f);
    p2 = fromFrameIdp(F,p2f);

    l  = [p1;p2(4:6,:)];

else
    % idp parts
    xf  = lf(1:3);
    w1f = lf(4:5);
    r1  = lf(6);
    w2f = lf(7:8);
    r2  = lf(9);

    % dir. vectors
    [m1f, M1F_w1f] = py2vec(w1f);
    [m2f, M2F_w2f] = py2vec(w2f);

    % from frame
    [x, X_f, X_xf] = fromFrame(F,xf);
    [m1, M1_f, M1_m1f]   = fromFrameVec(F,m1f);
    [m2, M2_f, M2_m2f]   = fromFrameVec(F,m2f);

    % angle vectors
    [w1, W1_m1] = vec2py(m1);
    [w2, W2_m2] = vec2py(m2);

    % partial Jacobians
    W1_f   = W1_m1*M1_f;
    W2_f   = W2_m2*M2_f;
    W1_w1f = W1_m1*M1_m1f*M1F_w1f;
    W2_w2f = W2_m2*M2_m2f*M2F_w2f;
    R_f    = zeros(1,7);

    % new idp line
    l = [x;w1;r1;w2;r2];

    % Jacobians
    L_f  = [X_f;W1_f;R_f;W2_f;R_f];
    L_lf = [...
        X_xf        zeros(3)          zeros(3)
        zeros(2,3)  W1_w1f     [0;0]  zeros(2,3)
        0 0 0       0 0         1     0 0        0
        zeros(2,3)  zeros(2,3)        W2_w2f    [0;0]
        0 0 0       0 0         0     0 0        1    ] ;

end

return

%% jac

syms x y z a b c d X Y Z A1 B1 R1 A2 B2 R2 real
F.x = [x;y;z;a;b;c;d];
F   = updateFrame(F);
l_F = [X;Y;Z;A1;B1;R1;A2;B2;R2];

[l,L_f,L_lf] = fromFrameIdpLin(F,l_F);

simplify(L_f  - jacobian(l,F.x))
simplify(L_lf - jacobian(l,l_F))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function dim_num = p03_dim_num ( )

%*****************************************************************************80
%
%% P03_DIM_NUM returns the spatial dimension for problem p03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DIM_NUM, the spatial dimension of the
%    dependent variables.
%
  dim_num = 2;

  return
end


text:
% Author : FUAT COGUN
% Date   : 16.07.2009
%
%
% Gets the dimensions of the screen.
%
%
% Inputs  :
% ======
%
% Outputs : -Width of the screen                                (Integer)
% =======   -Height of the screen                               (Integer)

function [screenWidth screenHeight] = getScreenSize()

    screenSize = get(0,'ScreenSize');
    screenWidth = screenSize(3);
    screenHeight = screenSize(4);


text:
function x_extract = extractsolution(momentdata,options)
%EXTRACTSOLUTIONS Tries to extract solutions from moment matrices
%
%
%  xoptimal = extractsolution(momentstructure)
%
%   xoptimal    : Extracted solutions
%
%   momentdata : Problem data, obtained from SOLVEMOMENT (or SOLVESDP)
%   options    : Options structure from SDPSETTINGS
%
%   See also SOLVEMOMENT, SDPSETTINGS

moment = momentdata.moment;
x = momentdata.x;
monomials = momentdata.monomials;
n = momentdata.n;
d = momentdata.d;

[U,S,V,ranks] = numranks(moment);

if options.moment.extractrank>0
    % We try a extraction from highest order moment no matter what
    flat = length(moment);
    ranks(flat) = options.moment.extractrank;
else
    % Find a flat extension
    flat = d+min(find(ranks(1+d:end)-ranks(1:end-d)==0));
end

if ~isempty(flat)

    % Find a basis
    r = ranks(flat);
    V = U{flat}(:,1:r)*sqrt(diag(diag(S{flat}(1:r,1:r))));
    if options.moment.rceftol >= 0
        [R,pivot] = rref(V',options.moment.rceftol);R = R';
    else
        % Try to find a reasonable tolerance by avoiding severly badly
        % conditioned R. Hack.., but seem to behave rather robustly.
        cV = cond(V);
        tol = 1e-10;
        [R,pivot] = rref(V',tol);R = R';
        while tol<1 & (cond(R)/cond(cV)>1e4)
            tol = tol*5;
            [R,pivot] = rref(V',tol);R = R';
        end
    end

    % Figure out multiplying matrices using YALMIP code
    w = monomials(pivot);
    for i = 1:n
        xw = x(i)*w;
        k = [];
        for j = 1:length(xw)
            k = [k;find(ismember(xw(j),monomials))];
        end
        N{i} = R(k,:);
    end

    % Things missing in the basis...
    if ~all(cellfun('prodofsize',N)==length(w)^2)
       x_extract = {[]};
       return;
    end

    % Create random convex combination
    rands = rand(n,1);rands = rands/sum(rands);
    M = 0;
    for i = 1:n
        M = M + rands(i)*N{i};
    end

    [Q,T] = schur(M);
    % Extract solution
    for i = 1:r
        for j = 1:n
            x_extract{i}(j,1) =  Q(:,i)'*N{j}*Q(:,i);
        end
    end

    % Refine solutions v-Rw = e(x)=0
    if options.moment.refine>0
        xtemp = double(x);
        e = monomials(1:size(R,1))-R*w;
        dedx = jacobian(e,x);
        for j = 1:r
            assign(x,x_extract{j});
            for i = 1:options.moment.refine
                assign(x,double(x)-double(dedx)\double(e));
            end
            x_extract{j} = double(x);
        end
        assign(x,xtemp);
    end

else
    x_extract = {};
end

% Somewhat more stable rank detection
% looking for sharp drops in singular value
function [U,S,V,ranks] = numranks(moment)
for i = 1:length(moment)
    [U{i},S{i},V{i}] = svd(moment{i});
    s = diag(S{i});
    decay = s(2:end)./(eps+s(1:end-1));
    r = min(find(decay<1e-3));
    if isempty(r)
        ranks(i) = rank(moment{i},1e-8);
    else
        ranks(i) = min(r,rank(moment{i},1e-8));
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function t_opf_dc_glpk(quiet)
%T_OPF_DC_GLPK  Tests for DC optimal power flow using GLPK solver.

%   MATPOWER
%   Copyright (c) 2004-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

algs  = [1; 1; 2];
alg_names = {
    'primal simplex',
    'dual simplex',
    'interior',
};
if have_feature('octave')
    if have_feature('octave', 'vnum') < 3.007
        dual = [0; 1; 1];
    else
        dual = [1; 2; 1];
    end
else
    dual = [0; 2; 1];
end
num_tests = 43 * length(algs);

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end
if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
end

mpopt = mpoption('out.all', 0, 'verbose', verbose);
mpopt = mpoption(mpopt, 'opf.dc.solver', 'GLPK');

%% run DC OPF
if have_feature('glpk')
  for k = 1:length(algs)
    mpopt = mpoption(mpopt, 'glpk.opts.lpsolver', algs(k), 'glpk.opts.dual', dual(k));
    if algs(k) == 2     %% interior point sometimes fails on GLPK 4.6x w/o this
        mpopt = mpoption(mpopt, 'glpk.opts.scale', 1);
    end
    t0 = sprintf('DC OPF (GLPK %s): ', alg_names{k});

    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved DC power flow case from MAT-file
    load soln9_dcopf;       %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    t = t0;
    [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with angle difference limits  -----
    t = [t0 'w/angle diff lims : '];
    mpc = loadcase(casefile);
    mpc.branch(4, ANGMAX) = 3;
    mpc.branch(7, ANGMIN) = -4.5;
    r = rundcopf(mpc, mpopt);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(   f, 6456.7213, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,PG        ),    [99.98497;89.35133;125.66371], 4, [t 'gen dispatch']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    e = zeros(size(branch, 1), 1);
    e(4) = 297.83776;
    e(7) = -26.94788;
    t_is(branch(:,MU_ANGMAX )-branch(:,MU_ANGMIN ), e, 4, [t 'branch ang diff mu']);

    t = [t0 'w/ignored angle diff lims : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    r = rundcopf(mpc, mpopt1);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% two new z variables
    %%      0 <= z1, P3 - P1 <= z1
    %%      0 <= z2, P3 - P2 <= z2
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[10;12;13;12;11;14],[-1;1;-1;1;-1;-1],2,14);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [13;14], [1;1], 2, 14);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 1 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    %% with A and N sized for AC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[19;21;25;21;20;26],[-1;1;-1;1;-1;-1],2,26);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [25;26], [1;1], 2, 26);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 2 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    t = [t0 'infeasible : '];
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1], [10;11], [1;1], 1, 14);   %% Pg1 + Pg2
    mpc.u = Inf;
    mpc.l = 600;
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(~success, [t 'no success']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc = loadcase(casefile);
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = rundcopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end
  end
else
    t_skip(num_tests, 'GLPK not available');
end

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function mv = mv_pca(mv,centered);
%
% mv = mv_pca([mv],[centered]);
%
% Perform PCA on multi-voxel data.
%
%
%
% ras 05/05.
if ieNotDefined('mv')
    mv = get(gcf,'UserData');
end

if ieNotDefined('centered')
    % this is currently not used,
    % since I always do it both ways
    centered = 1;
end

% get the voxels x conditions amplitudes matrix
amps = mv_amps(mv);

% apply PCA
[pcs score k Uk Lk] = pca_all(amps); % kalanit's code, uncentered
eigenvals = Lk(sub2ind(size(Lk),1:k,1:k));      % diagonal entries of gamma mat
varExplained = 100 .* eigenvals ./ sum(eigenvals);     % percent variance explained

% redo centered, using MATLAB command
[pcs score latent tsquare] = princomp(amps);

% add to mv struct
mv.pca.pcs = pcs;
mv.pca.score = score;
mv.pca.latent = latent;
mv.pca.tsquare = tsquare;
mv.pca.k = k;
mv.pca.Uk = Uk;
mv.pca.Lk = Lk;

% if a UI exists, visualize, set as user data
if isfield(mv.ui,'fig') & ishandle(mv.ui.fig)
    set(mv.ui.fig,'UserData',mv);
    figure(mv.ui.fig)
    multiVoxelUI; % refresh UI
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function net = fcnInitializeModel16sGeneric(labelCount, net)
%FCNINITIALIZEMODEL16S Initialize the FCN-16S model from FCN-32

%% Remove the last layer
net.removeLayer('deconv32') ;

%% Add the first Deconv layer
filters = single(bilinear_u(4, 1, labelCount)) ;
net.addLayer('deconv2', ...
  dagnn.ConvTranspose('size', size(filters), ...
                      'upsample', 2, ...
                      'crop', [1 1 1 1], ...
                      'hasBias', false), ...
             'x38', 'x39', 'deconv1f') ;
f = net.getParamIndex('deconv1f') ;
net.params(f).value = filters ;
net.params(f).learningRate = 0 ;
net.params(f).weightDecay = 1 ;

%% Add skip layers on top of pool4
net.addLayer('skip4', ...
     dagnn.Conv('size', [1 1 512 labelCount], 'pad', 0), ...
     'x24', 'x40', {'skip4f','skip4b'});

f = net.getParamIndex('skip4f') ;
net.params(f).value = zeros(1, 1, 512, labelCount, 'single') ;
net.params(f).learningRate = 0.1 ;
net.params(f).weightDecay = 1 ;

f = net.getParamIndex('skip4b') ;
net.params(f).value = zeros(1, 1, labelCount, 'single') ;
net.params(f).learningRate = 2 ;
net.params(f).weightDecay = 1 ;

%% Add the summation layer
net.addLayer('sum1', dagnn.Sum(), {'x39', 'x40'}, 'x41') ;

%% Add deconvolutional layer implementing bilinear interpolation
filters = single(bilinear_u(32, labelCount, labelCount)) ;
net.addLayer('deconv16', ...
  dagnn.ConvTranspose('size', size(filters), ...
                      'upsample', 16, ...
                      'crop', 8, ...
                      'numGroups', labelCount, ...
                      'hasBias', false), ...
             'x41', 'prediction', 'deconv16f') ;

f = net.getParamIndex('deconv16f') ;
net.params(f).value = filters ;
net.params(f).learningRate = 0 ;
net.params(f).weightDecay = 1 ;

% Make sure that the output of the bilinear interpolator is not discared for
% visualization purposes
net.vars(net.getVarIndex('prediction')).precious = 1 ;

% empirical test
if 0
  figure(100) ; clf ;
  n = numel(net.vars) ;
  for i=1:n
    vl_tightsubplot(n,i) ;
    showRF(net, 'input', net.vars(i).name) ;
    title(sprintf('%s', net.vars(i).name)) ;
    drawnow ;
  end
end

author: nightrome

repo: matconvnet-calvin


text:
classdef L2RegressionLoss < dagnn.Loss
    methods
        function outputs = forward(obj, inputs, params)
            sz = size(inputs{2});
            mass = sz(1) * sz(2) + 1;

            outputs{1} = vl_nnloss_regression(inputs{1}, inputs{2}, [], ...
                'loss', obj.loss, ...
                'instanceWeights', 1./mass) ;
            n = obj.numAveraged ;
            m = n + size(inputs{1},4) ;
            obj.average = (n * obj.average + double(gather(outputs{1}))) / m ;
            obj.numAveraged = m ;
        end

        function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
            sz = size(inputs{2});
            mass = sz(1) * sz(2) + 1;
            derInputs{1} = vl_nnloss_regression(inputs{1}, inputs{2}, derOutputs{1}, ...
                'loss', obj.loss, ...
                'instanceWeights', 1./mass) ;
            derInputs{2} = [] ;
            derParams = {} ;
        end

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

author: aimerykong


text:
function [net, info] = cnn_train(net, imdb, getBatch, varargin)
%CNN_TRAIN  An example implementation of SGD for training CNNs
%    CNN_TRAIN() is an example learner implementing stochastic
%    gradient descent with momentum to train a CNN. It can be used
%    with different datasets and tasks by providing a suitable
%    getBatch function.
%
%    The function automatically restarts after each training epoch by
%    checkpointing.
%
%    The function supports training on CPU or on one or more GPUs
%    (specify the list of GPU IDs in the `gpus` option). Multi-GPU
%    support is relatively primitive but sufficient to obtain a
%    noticable speedup.

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

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;
opts.profile = false ;

opts.conserveMemory = true ;
opts.backPropDepth = +inf ;
opts.sync = false ;
opts.cudnn = true ;
opts.errorFunction = 'multiclass' ;
opts.errorLabels = {} ;
opts.plotDiagnostics = false ;
opts.plotStatistics = true;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end
if isnan(opts.val), opts.val = [] ; end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

net = vl_simplenn_tidy(net); % fill in some eventually missing values
net.layers{end-1}.precious = 1; % do not remove predictions, used for error
vl_simplenn_display(net, 'batchSize', opts.batchSize) ;

evaluateMode = isempty(opts.train) ;

if ~evaluateMode
  for i=1:numel(net.layers)
    if isfield(net.layers{i}, 'weights')
      J = numel(net.layers{i}.weights) ;
      for j=1:J
        net.layers{i}.momentum{j} = zeros(size(net.layers{i}.weights{j}), 'single') ;
      end
      if ~isfield(net.layers{i}, 'learningRate')
        net.layers{i}.learningRate = ones(1, J, 'single') ;
      end
      if ~isfield(net.layers{i}, 'weightDecay')
        net.layers{i}.weightDecay = ones(1, J, 'single') ;
      end
    end
  end
end

% setup GPUs
numGpus = numel(opts.gpus) ;
if numGpus > 1
  if isempty(gcp('nocreate')),
    parpool('local',numGpus) ;
    spmd, gpuDevice(opts.gpus(labindex)), end
  end
elseif numGpus == 1
  gpuDevice(opts.gpus)
end
if exist(opts.memoryMapFile), delete(opts.memoryMapFile) ; end

% setup error calculation function
hasError = true ;
if isstr(opts.errorFunction)
  switch opts.errorFunction
    case 'none'
      opts.errorFunction = @error_none ;
      hasError = false ;
    case 'multiclass'
      opts.errorFunction = @error_multiclass ;
      if isempty(opts.errorLabels), opts.errorLabels = {'top1err', 'top5err'} ; end
    case 'binary'
      opts.errorFunction = @error_binary ;
      if isempty(opts.errorLabels), opts.errorLabels = {'binerr'} ; end
    otherwise
      error('Unknown error function ''%s''.', opts.errorFunction) ;
  end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  load(modelPath(start), 'net', 'info') ;
  net = vl_simplenn_tidy(net) ; % just in case MatConvNet was updated
end

for epoch=start+1:opts.numEpochs

  % train one epoch and validate
  learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  train = opts.train(randperm(numel(opts.train))) ; % shuffle
  val = opts.val ;

  if numGpus <= 1
    [net,stats.train,prof] = process_epoch(opts, getBatch, epoch, train, learningRate, imdb, net) ;
    [~,stats.val] = process_epoch(opts, getBatch, epoch, val, 0, imdb, net) ;
    if opts.profile
      profile('viewer') ;
      keyboard ;
    end
  else
    fprintf('%s: sending model to %d GPUs\n', mfilename, numGpus) ;
    spmd(numGpus)
      [net_, stats_train_,prof_] = process_epoch(opts, getBatch, epoch, train, learningRate, imdb, net) ;
      [~, stats_val_] = process_epoch(opts, getBatch, epoch, val, 0, imdb, net_) ;
    end
    net = net_{1} ;
    stats.train = sum([stats_train_{:}],2) ;
    stats.val = sum([stats_val_{:}],2) ;
    if opts.profile
      mpiprofile('viewer', [prof_{:,1}]) ;
      keyboard ;
    end
    clear net_ stats_train_ stats_val_ ;
  end

  % save
  if evaluateMode, sets = {'val'} ; else sets = {'train', 'val'} ; end
  for f = sets
    f = char(f) ;
    n = numel(eval(f)) ;
    info.(f).speed(epoch) = n / stats.(f)(1) * max(1, numGpus) ;
    info.(f).objective(epoch) = stats.(f)(2) / n ;
    info.(f).error(:,epoch) = stats.(f)(3:end) / n ;
  end
  if ~evaluateMode
    fprintf('%s: saving model for epoch %d\n', mfilename, epoch) ;
    tic ;
    save(modelPath(epoch), 'net', 'info') ;
    fprintf('%s: model saved in %.2g s\n', mfilename, toc) ;
  end

  if opts.plotStatistics
    switchfigure(1) ; clf ;
    subplot(1,1+hasError,1) ;
    if ~evaluateMode
      semilogy(1:epoch, info.train.objective, '.-', 'linewidth', 2) ;
      hold on ;
    end
    semilogy(1:epoch, info.val.objective, '.--') ;
    xlabel('training epoch') ; ylabel('energy') ;
    grid on ;
    h=legend(sets) ;
    set(h,'color','none');
    title('objective') ;
    if hasError
      subplot(1,2,2) ; leg = {} ;
      if ~evaluateMode
        plot(1:epoch, info.train.error', '.-', 'linewidth', 2) ;
        hold on ;
        leg = horzcat(leg, strcat('train ', opts.errorLabels)) ;
      end
      plot(1:epoch, info.val.error', '.--') ;
      leg = horzcat(leg, strcat('val ', opts.errorLabels)) ;
      set(legend(leg{:}),'color','none') ;
      grid on ;
      xlabel('training epoch') ; ylabel('error') ;
      title('error') ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end
end

% -------------------------------------------------------------------------
function err = error_multiclass(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
[~,predictions] = sort(predictions, 3, 'descend') ;

% be resilient to badly formatted labels
if numel(labels) == size(predictions, 4)
  labels = reshape(labels,1,1,1,[]) ;
end

% skip null labels
mass = single(labels(:,:,1,:) > 0) ;
if size(labels,3) == 2
  % if there is a second channel in labels, used it as weights
  mass = mass .* labels(:,:,2,:) ;
  labels(:,:,2,:) = [] ;
end

m = min(5, size(predictions,3)) ;

error = ~bsxfun(@eq, predictions, labels) ;
err(1,1) = sum(sum(sum(mass .* error(:,:,1,:)))) ;
err(2,1) = sum(sum(sum(mass .* min(error(:,:,1:m,:),[],3)))) ;

% -------------------------------------------------------------------------
function err = error_binary(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
error = bsxfun(@times, predictions, labels) < 0 ;
err = sum(error(:)) ;

% -------------------------------------------------------------------------
function err = error_none(opts, labels, res)
% -------------------------------------------------------------------------
err = zeros(0,1) ;

% -------------------------------------------------------------------------
function  [net_cpu,stats,prof] = process_epoch(opts, getBatch, epoch, subset, learningRate, imdb, net_cpu)
% -------------------------------------------------------------------------

% move the CNN to GPU (if needed)
numGpus = numel(opts.gpus) ;
if numGpus >= 1
  net = vl_simplenn_move(net_cpu, 'gpu') ;
  one = gpuArray(single(1)) ;
else
  net = net_cpu ;
  net_cpu = [] ;
  one = single(1) ;
end

% assume validation mode if the learning rate is zero
training = learningRate > 0 ;
if training
  mode = 'train' ;
  evalMode = 'normal' ;
else
  mode = 'val' ;
  evalMode = 'test' ;
end

% turn on the profiler (if needed)
if opts.profile
  if numGpus <= 1
    prof = profile('info') ;
    profile clear ;
    profile on ;
  else
    prof = mpiprofile('info') ;
    mpiprofile reset ;
    mpiprofile on ;
  end
end

res = [] ;
mmap = [] ;
stats = [] ;
start = tic ;

for t=1:opts.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d: ', mode, epoch, ...
          fix(t/opts.batchSize)+1, ceil(numel(subset)/opts.batchSize)) ;
  batchSize = min(opts.batchSize, numel(subset) - t + 1) ;
  numDone = 0 ;
  error = [] ;
  for s=1:opts.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+opts.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
    [im, labels] = getBatch(imdb, batch) ;

    if opts.prefetch
      if s==opts.numSubBatches
        batchStart = t + (labindex-1) + opts.batchSize ;
        batchEnd = min(t+2*opts.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
      getBatch(imdb, nextBatch) ;
    end

    if numGpus >= 1
      im = gpuArray(im) ;
    end

    % evaluate the CNN
    net.layers{end}.class = labels ;
    if training, dzdy = one; else, dzdy = [] ; end
    res = vl_simplenn(net, im, dzdy, res, ...
                      'accumulate', s ~= 1, ...
                      'mode', evalMode, ...
                      'conserveMemory', opts.conserveMemory, ...
                      'backPropDepth', opts.backPropDepth, ...
                      'sync', opts.sync, ...
                      'cudnn', opts.cudnn) ;

    % accumulate training errors
    error = sum([error, [...
      sum(double(gather(res(end).x))) ;
      reshape(opts.errorFunction(opts, labels, res),[],1) ; ]],2) ;
    numDone = numDone + numel(batch) ;
  end % next sub-batch

  % gather and accumulate gradients across labs
  if training
    if numGpus <= 1
      [net,res] = accumulate_gradients(opts, learningRate, batchSize, net, res) ;
    else
      if isempty(mmap)
        mmap = map_gradients(opts.memoryMapFile, net, res, numGpus) ;
      end
      write_gradients(mmap, net, res) ;
      labBarrier() ;
      [net,res] = accumulate_gradients(opts, learningRate, batchSize, net, res, mmap) ;
    end
  end

  % collect and print learning statistics
  time = toc(start) ;
  stats = sum([stats,[0 ; error]],2); % works even when stats=[]
  stats(1) = time ;
  n = t + batchSize - 1 ; % number of images processed overall
  speed = n/time ;
  fprintf('%.1f Hz%s\n', speed) ;

  m = n / max(1,numlabs) ; % num images processed on this lab only
  fprintf(' obj:%.3g', stats(2)/m) ;
  for i=1:numel(opts.errorLabels)
    fprintf(' %s:%.3g', opts.errorLabels{i}, stats(i+2)/m) ;
  end
  fprintf(' [%d/%d]', numDone, batchSize);
  fprintf('\n') ;

  % collect diagnostic statistics
  if training & opts.plotDiagnostics
    switchfigure(2) ; clf ;
    diag = [res.stats] ;
    barh(horzcat(diag.variation)) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTickLabel',horzcat(diag.label), ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1]) ;
    drawnow ;
  end

end

% switch off the profiler
if opts.profile
  if numGpus <= 1
    prof = profile('info') ;
    profile off ;
  else
    prof = mpiprofile('info');
    mpiprofile off ;
  end
else
  prof = [] ;
end

% bring the network back to CPU
if numGpus >= 1
  net_cpu = vl_simplenn_move(net, 'cpu') ;
else
  net_cpu = net ;
end

% -------------------------------------------------------------------------
function [net,res] = accumulate_gradients(opts, lr, batchSize, net, res, mmap)
% -------------------------------------------------------------------------
if nargin >= 6
  numGpus = numel(mmap.Data) ;
else
  numGpus = 1 ;
end

for l=numel(net.layers):-1:1
  for j=1:numel(res(l).dzdw)

    % accumualte gradients from multiple labs (GPUs) if needed
    if numGpus > 1
      tag = sprintf('l%d_%d',l,j) ;
      tmp = zeros(size(mmap.Data(labindex).(tag)), 'single') ;
      for g = setdiff(1:numGpus, labindex)
        tmp = tmp + mmap.Data(g).(tag) ;
      end
      res(l).dzdw{j} = res(l).dzdw{j} + tmp ;
    end

    if j == 3 && strcmp(net.layers{l}.type, 'bnorm')
      % special case for learning bnorm moments
      thisLR = net.layers{l}.learningRate(j) ;
      net.layers{l}.weights{j} = ...
        (1-thisLR) * net.layers{l}.weights{j} + ...
        (thisLR/batchSize) * res(l).dzdw{j} ;
    else
      % standard gradient training
      thisDecay = opts.weightDecay * net.layers{l}.weightDecay(j) ;
      thisLR = lr * net.layers{l}.learningRate(j) ;
      net.layers{l}.momentum{j} = ...
        opts.momentum * net.layers{l}.momentum{j} ...
        - thisDecay * net.layers{l}.weights{j} ...
        - (1 / batchSize) * res(l).dzdw{j} ;
      net.layers{l}.weights{j} = net.layers{l}.weights{j} + ...
        thisLR * net.layers{l}.momentum{j} ;
    end

    % if requested, collect some useful stats for debugging
    if opts.plotDiagnostics
      variation = [] ;
      label = '' ;
      switch net.layers{l}.type
        case {'conv','convt'}
          variation = thisLR * mean(abs(net.layers{l}.momentum{j}(:))) ;
          if j == 1 % fiters
            base = mean(abs(net.layers{l}.weights{j}(:))) ;
            label = 'filters' ;
          else % biases
            base = mean(abs(res(l+1).x(:))) ;
            label = 'biases' ;
          end
          variation = variation / base ;
          label = sprintf('%s_%s', net.layers{l}.name, label) ;
      end
      res(l).stats.variation(j) = variation ;
      res(l).stats.label{j} = label ;
    end
  end
end

% -------------------------------------------------------------------------
function mmap = map_gradients(fname, net, res, numGpus)
% -------------------------------------------------------------------------
format = {} ;
for i=1:numel(net.layers)
  for j=1:numel(res(i).dzdw)
    format(end+1,1:3) = {'single', size(res(i).dzdw{j}), sprintf('l%d_%d',i,j)} ;
  end
end
format(end+1,1:3) = {'double', [3 1], 'errors'} ;
if ~exist(fname) && (labindex == 1)
  f = fopen(fname,'wb') ;
  for g=1:numGpus
    for i=1:size(format,1)
      fwrite(f,zeros(format{i,2},format{i,1}),format{i,1}) ;
    end
  end
  fclose(f) ;
end
labBarrier() ;
mmap = memmapfile(fname, 'Format', format, 'Repeat', numGpus, 'Writable', true) ;

% -------------------------------------------------------------------------
function write_gradients(mmap, net, res)
% -------------------------------------------------------------------------
for i=1:numel(net.layers)
  for j=1:numel(res(i).dzdw)
    mmap.Data(labindex).(sprintf('l%d_%d',i,j)) = gather(res(i).dzdw{j}) ;
  end
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchfigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
%
%   [NSx, NEV] = combineNSxNEV(filename1, filename2)
%
%   This function loads two NSx and NEV files and it will combine them
%   together into a sinlge NSx and NEV structure in MATLAB. To merge two
%   NSx and NEV files into individual files see mergeNSxNEV. The time
%   difference between the two sets of recordings is removed. To determine
%   the time differnce between the two data files, use
%   NSx.MetaTags.DateTimeRaw or NEV.MetaTags.DateTimeRaw variables.
%
%
%   filename1:  The name of the first NSx file. This input is optional. In
%               its absense, a dialog will open and will prompt the user to
%               select an NSx file.
%               (OPTIONAL)
%
%   filename2:  The name of the second NSx file. This input is also
%               optional. In its absense, a dialog will open and will
%               prompt the user to select an NSx file.
%               (OPTIONAL)
%
%
%   Example:
%
%   [NSx, NEV] = combineNSxNEV('c:\data\saveddata1.ns5', 'c:\data\saveddata2.ns5');
%
%   The above example reads the two files (full path needed)
%   c:\data\saveddata1.ns5 and c:\data\saveddata2.ns5 and their corresponding
%   NEV files (saveddata1.nev and saveddata2.nev) in the same folder and
%   combines them into single variables NSx and NEV into MATLAB workspace.
%
%   Kian Torab
%   ktorab@blackrockmicro.com
%   Blackrock Microsystems
%
%   Version 1.1.2.0

function [NSx1, NEV1] = combineNSxNEV(filename1, filename2)

% Openning NSx files
if exist('filename1', 'var') && exist('filename2', 'var')
    disp('Load the first NSx file.');
    NSx1 = openNSx('read', filename1);
    disp('Load the second NSx file.');
    NSx2 = openNSx('read', filename2);
    if NSx1.MetaTags.SamplingFreq ~= NSx2.MetaTags.SamplingFreq;
        disp('The sampling frequencies are not the same.');
        return;
    end
else
    NSx1 = openNSx('read');
    NSx2 = openNSx('read');
end

% Determining length of the first NSx file
conversionFactor = 30000/NSx1.MetaTags.SamplingFreq;
NSx1DataLength = NSx1.MetaTags.DataPoints * conversionFactor;

% Combining NSx files
NSx1.Data = [NSx1.Data, NSx2.Data];

% Opening NEV files
fileNameNEV1 = [NSx1.MetaTags.FilePath '/' NSx1.MetaTags.Filename(1:end-3), 'nev'];
fileNameNEV2 = [NSx1.MetaTags.FilePath '/' NSx2.MetaTags.Filename(1:end-3), 'nev'];
clear NSx2;

if (exist(fileNameNEV1, 'file') == 2) && (exist(fileNameNEV2, 'file') ==2)
    disp('Openning corresponding NEV files...');
    NEV1 = openNEV('read', fileNameNEV1);
    NEV2 = openNEV('read', fileNameNEV2);
else
    disp('Load the first NEV file.');
    NEV1 = openNEV('read');
    disp('Load the second NEV file.');
    NEV2 = openNEV('read');
end

% Adjusting the timestamp on the second NEV file
NEV2.Data.Comments.TimeStamp = NEV2.Data.Comments.TimeStamp + NSx1DataLength;
NEV2.Data.Comments.TimeStampSec = NEV2.Data.Comments.TimeStampSec + double(NSx1DataLength)/30;
NEV2.Data.SerialDigitalIO.TimeStamp = NEV2.Data.SerialDigitalIO.TimeStamp + NSx1DataLength;
NEV2.Data.SerialDigitalIO.TimeStampSec = NEV2.Data.SerialDigitalIO.TimeStampSec + double(NSx1DataLength)/30;
NEV2.Data.Spikes.TimeStamp = NEV2.Data.Spikes.TimeStamp + NSx1DataLength;
NEV2.Data.VideoSync.TimeStamp = NEV2.Data.VideoSync.TimeStamp + NSx1DataLength;
if ~isempty(NEV2.Data.Tracking)
    trackingFieldNames = fieldnames(NEV2.Data.Tracking);
    for idx = 1:size(trackingFieldNames, 1)
        NEV2.Data.Tracking.(trackingFieldNames{idx}).TimeStamp = NEV2.Data.Tracking.(trackingFieldNames{idx}).TimeStamp + NSx1DataLength;
    end
end
NEV2.Data.PatientTrigger.TimeStamp = NEV2.Data.PatientTrigger.TimeStamp + NSx1DataLength;
NEV2.Data.Reconfig.TimeStamp = NEV2.Data.Reconfig.TimeStamp + NSx1DataLength;

% Combining the two NEV files
NEV1.Data.Spikes.Electrode      = [NEV1.Data.Spikes.Electrode, NEV2.Data.Spikes.Electrode];
NEV1.Data.Spikes.TimeStamp      = [NEV1.Data.Spikes.TimeStamp, NEV2.Data.Spikes.TimeStamp];
NEV1.Data.Spikes.Unit           = [NEV1.Data.Spikes.Unit, NEV2.Data.Spikes.Unit];
NEV1.Data.Spikes.Waveform       = [NEV1.Data.Spikes.Waveform, NEV2.Data.Spikes.Waveform];
NEV1.Data.Comments.TimeStamp    = [NEV1.Data.Comments.TimeStamp, NEV2.Data.Comments.TimeStamp];
NEV1.Data.Comments.TimeStampSec = [NEV1.Data.Comments.TimeStampSec, NEV2.Data.Comments.TimeStampSec];
NEV1.Data.Comments.CharSet      = [NEV1.Data.Comments.CharSet, NEV2.Data.Comments.CharSet];
NEV1.Data.Comments.Color        = [NEV1.Data.Comments.Color, NEV2.Data.Comments.Color];
NEV1.Data.Comments.Text         = [NEV1.Data.Comments.Text; NEV2.Data.Comments.Text];
if ~isempty(NEV2.Data.Tracking)
    for idx = 1:size(trackingFieldNames, 1)
        NEV1.Data.Tracking.(trackingFieldNames{idx}).TimeStamp = [NEV1.Data.Tracking.(trackingFieldNames{idx}).TimeStamp, NEV2.Data.Tracking.(trackingFieldNames{idx}).TimeStamp];
        NEV1.Data.Tracking.(trackingFieldNames{idx}).TimeStampSec = [NEV1.Data.Tracking.(trackingFieldNames{idx}).TimeStampSec, NEV2.Data.Tracking.(trackingFieldNames{idx}).TimeStampSec];
        NEV1.Data.Tracking.(trackingFieldNames{idx}).ParentID = [NEV1.Data.Tracking.(trackingFieldNames{idx}).ParentID, NEV2.Data.Tracking.(trackingFieldNames{idx}).ParentID];
        NEV1.Data.Tracking.(trackingFieldNames{idx}).NodeCount = [NEV1.Data.Tracking.(trackingFieldNames{idx}).NodeCount, NEV2.Data.Tracking.(trackingFieldNames{idx}).NodeCount];
        NEV1.Data.Tracking.(trackingFieldNames{idx}).MarkerCount = [NEV1.Data.Tracking.(trackingFieldNames{idx}).MarkerCount, NEV2.Data.Tracking.(trackingFieldNames{idx}).MarkerCount];
        NEV1.Data.Tracking.(trackingFieldNames{idx}).MarkerCoordinates = [NEV1.Data.Tracking.(trackingFieldNames{idx}).MarkerCoordinates NEV2.Data.Tracking.(trackingFieldNames{idx}).MarkerCoordinates];
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Chapter 7: Statistical estimation
%
%  counting_problem_poisson.m - Section 7.1.1: Counting problems with Poisson distribution
%  ML_covariance_est.m        - Section 7.1.1: Covariance estimation for Gaussian variables
%  probbounds.m               - Section 7.4.3: Probability bounds example with Voronoi diagram
%  expdesign.m                - Section 7.5.2: Experiment design
%  logistics.m                - Figure 7.1: Logistic regression
%  logistics_gp.m             - Figure 7.1: Logistic regression (GP version)
%  maxent.m                   - Example 7.2: Maximum entropy distribution
%  detector2.m                - Example 7.4: Binary hypothesis testing
%  cheb.m                     - Computes Chebyshev lower bounds on probability vectors
%  cher.m                     - Computes Chernoff upper bounds on probability
%  montecarlo.m               - Probability estimate using Monte Carlo methods
help Contents

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX