text:
function a = tri_upper_inverse ( alpha, n )

%*****************************************************************************80
%
%% TRI_UPPER_INVERSE returns the inverse of the TRI_UPPER matrix.
%
%  Formula:
%
%    if ( I = J )
%      A(I,J) = 1
%    elseif ( I = J - 1 )
%      A(I,J) = -ALPHA
%    elseif ( I < J )
%      A(I,J) = - ALPHA * ( 1-ALPHA)^(J-I-1)
%    else
%      A(I,J) = 0
%
%  Example:
%
%    ALPHA = 3, N = 5
%
%    1 -3  6 -12  24
%    0  1 -3   6 -12
%    0  0  1  -3   6
%    0  0  0   1  -3
%    0  0  0   0   1
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is nonsingular.
%
%    A is upper triangular.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    det ( A ) = 1.
%
%    A is unimodular.
%
%    LAMBDA(1:N) = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, value used on the superdiagonals.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( i == j )
        a(i,j) = 1.0;
      elseif ( i == j - 1 )
        a(i,j) = - alpha;
      elseif ( i < j )
        a(i,j) = - alpha * ( 1.0 - alpha )^(j-i-1);
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
function scale = p37_scale ( neqn )

%*****************************************************************************80
%
%% P37_SCALE returns scale factors for problem p37.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = [ 1.0; 1.0 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function g = kernGradient(kern, x, varargin)

% KERNGRADIENT Compute the gradient wrt the kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO kernDiagGradient, kernGradX

% KERN

fhandle = str2func([kern.type 'KernGradient']);

% varargin contains (optionally) x2 and covGrad.
g = fhandle(kern, x, varargin{:});

% Check if parameters are being optimised in a transformed space.
factors = kernFactors(kern, 'gradfact');
g(factors.index) = g(factors.index).*factors.val;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function By = mfBy(yc,para,flag);

if ~exist('flag','var'), flag = 'By';  end;

Omega = para.Omega;
dim   = length(Omega);
m     = para.m;
h     = Omega./m;
a     = sqrt(para.mu);
b     = sqrt(para.mu+para.lambda);

% Bu for elastic-grad
% yc -> [u1,u2,u3]
%
%     dim == 3,                     dim == 2,
%     | d11         |
%     | d21         |
%     | d31         |   | u1 |          | d11     |
%     |     d12     |   |    |          | d21     |   | u1 |
%  y= |     d22     | * | u2 |      y = |     d12 | * |    |
%     |     d32     |   |    |          |     d22 |   | u2 |
%     |         d13 |   | u3 |          | d11 d22 |
%     |         d23 |
%     |         d33 |
%     | d11 d22 d33 |

flag = ['elastic-stg-',flag,'-',int2str(dim)];

switch flag
  case 'elastic-stg-By-2'
    nc = m(1)*m(2);
    nn = (m(1)+1)*(m(2)+1);
    ns1 = (m(1)+1)*m(2);
    ns2 = m(1)*(m(2)+1);

    I11 = 1:nc;                        % d11 y1
    I21 = I11(end)+(1:nn-2*(m(1)+1));  % d21 y1
    I12 = I21(end)+(1:nn-2*(m(2)+1));  % d12 y2
    I22 = I12(end)+(1:nc);             % d22 y2
    I00 = I22(end)+(1:nc);             % div Y

    By = zeros(I00(end),1);

    y1 = reshape(yc(1:ns1),m(1)+1,m(2));
    y2 = reshape(yc(ns1+(1:ns2)),m(1),m(2)+1);

    %% d11 y1
    By(I11) = reshape(y1(2:end,:)-y1(1:end-1,:),length(I11),1)/h(1);

    %% d21 y1
    By(I21) = reshape(y1(:,2:end)-y1(:,1:end-1),length(I21),1)/h(2);

    %% d12 y2
    By(I12) = reshape(y2(2:end,:)-y2(1:end-1,:),length(I12),1)/h(1);

    %% d22 y2
    By(I22) = reshape(y2(:,2:end)-y2(:,1:end-1),length(I22),1)/h(2);

    %% d22 y2
    By(I00) = b*By(I11) + b*By(I22);
    By(1:I22(end)) = a*By(1:I22(end));

  case 'elastic-stg-By-3'

    n0 = m(1)*m(2)*m(3);
    n1 = (m(1)+1)*m(2)*m(3);
    n2 = m(1)*(m(2)+1)*m(3);
    n3 = m(1)*m(2)*(m(3)+1);

    j1 = m(1)*(m(2)+1)*(m(3)+1);
    j2 = (m(1)+1)*m(2)*(m(3)+1);
    j3 = (m(1)+1)*(m(2)+1)*m(3);

    p1 = n0 + j2 + j3;
    p2 = p1 + n0 + j1 + j3;
    p3 = p2 + n0 + j1 + j2;

    ly = 4*n0+2*(j1+j2+j3);

%     mfilename, keyboard
    [y1,y2,y3] = vec2array(yc,m,'stg');

    clear yc

    By = zeros(ly,1);

    % d11 y1 : n1 -> n0
    dummy = y1(2:end,:,:) - y1(1:end-1,:,:);
    By(1:n0) = reshape(dummy,n0,1)/h(1);

    % d21 y1 : n1 -> j3
    dummy = zeros(m(1)+1,m(2)+1,m(3));
    dummy(:,2:end-1,:) = y1(:,2:end,:) - y1(:,1:end-1,:);
    By(n0+(1:j3)) = reshape(dummy,j3,1)/h(2);

    % d31 y1 : n1 -> j2
    dummy = zeros(m(1)+1,m(2),m(3)+1);
    dummy(:,:,2:end-1) = y1(:,:,2:end) - y1(:,:,1:end-1);
    By(n0+j3+(1:j2)) = reshape(dummy,j2,1)/h(3);

    % d12 y2 : n2 -> j3
    dummy = zeros(m(1)+1,m(2)+1,m(3));
    dummy(2:end-1,:,:) = y2(2:end,:,:) - y2(1:end-1,:,:);
    By(p1+(1:j3)) = reshape(dummy,j3,1)/h(1);

    % d22 y2 : n2 -> n0
    dummy = y2(:,2:end,:) - y2(:,1:end-1,:);
    By(p1+j3+(1:n0)) = reshape(dummy,n0,1)/h(2);

    % d32 y2 : n2 -> j1
    dummy = zeros(m(1),m(2)+1,m(3)+1);
    dummy(:,:,2:end-1) = y2(:,:,2:end) - y2(:,:,1:end-1);
    By(p1+j3+n0+(1:j1)) = reshape(dummy,j1,1)/h(3);

    % d13 y3 : n3 -> j2
    dummy = zeros(m(1)+1,m(2),m(3)+1);
    dummy(2:end-1,:,:) = y3(2:end,:,:) - y3(1:end-1,:,:);
    By(p2+(1:j2)) = reshape(dummy,j2,1)/h(1);

    % d23 y3 : n3 -> j1
    dummy = zeros(m(1),m(2)+1,m(3)+1);
    dummy(:,2:end-1,:) = y3(:,2:end,:) - y3(:,1:end-1,:);
    By(p2+j2+(1:j1)) = reshape(dummy,j1,1)/h(2);

    % d33 y3 : n3 -> n0
    dummy = y3(:,:,2:end) - y3(:,:,1:end-1);
    By(p2+j2+j1+(1:n0)) = reshape(dummy,n0,1)/h(3);

    % div u
    By(end-n0+1:end) = By(1:n0)+By(p1+j3+(1:n0)) + By(p2+j2+j1+(1:n0));

    By(1:end-n0)     = a*By(1:end-n0);
    By(end-n0+1:end) = b*By(end-n0+1:end);

  case 'elastic-stg-BTy-2'

    % the input yc is the result of B*y, therefore it can be splitted into 5
    % parts: d11u1, d21u1,d12u2,d22u2,Divy

    nc  = m(1)*m(2);
    nn  = (m(1)+1)*(m(2)+1);
    ns1 = (m(1)+1)*m(2);
    ns2 = m(1)*(m(2)+1);

    By  = zeros(ns1+ns2,1);
    I11 = 1:nc;                        % d11 y1
    I21 = I11(end)+(1:nn-2*(m(1)+1));  % d21 y1
    I12 = I21(end)+(1:nn-2*(m(2)+1));  % d12 y2
    I22 = I12(end)+(1:nc);             % d22 y2
    I00 = I22(end)+(1:nc);             % div Y

    yd = reshape(yc(I11),m);
    zd = reshape([-yd(1,:);yd(1:end-1,:)-yd(2:end,:);yd(end,:)],ns1,1)/h(1);
    By(1:ns1) = a*zd;

    yd = reshape(yc(I21),m(1)+1,m(2)-1);
    zd = reshape([-yd(:,1),yd(:,1:end-1)-yd(:,2:end),yd(:,end)],ns1,1)/h(2);
    By(1:ns1) = By(1:ns1) + a*zd;

    yd = reshape(yc(I00),m);
    zd = reshape([-yd(1,:);yd(1:end-1,:)-yd(2:end,:);yd(end,:)],ns1,1)/h(1);
    By(1:ns1) = By(1:ns1) + b*zd;

    yd = reshape(yc(I12),m(1)-1,m(2)+1);
    zd = reshape([-yd(1,:);yd(1:end-1,:)-yd(2:end,:);yd(end,:)],ns2,1)/h(1);
    By(1+ns1:end) = a*zd;

    yd = reshape(yc(I22),m);
    zd = reshape([-yd(:,1),yd(:,1:end-1)-yd(:,2:end),yd(:,end)],ns2,1)/h(2);
    By(1+ns1:end) = By(1+ns1:end) + a*zd;

    yd = reshape(yc(I00),m);
    zd = reshape([-yd(:,1),yd(:,1:end-1)-yd(:,2:end),yd(:,end)],ns2,1)/h(2);
    By(1+ns1:end) = By(1+ns1:end) + b*zd;

  case 'elastic-stg-BTy-3'

    % the input yc is the result of B*y, therefore it can be splitted into 10
    % parts: d11y1, d21y1,d31y1,d12y2,d22y2,d32y2,d13y3,d23y3,u33y3,Divy

    n0 = m(1)*m(2)*m(3);
    n1 = (m(1)+1)*m(2)*m(3);
    n2 = m(1)*(m(2)+1)*m(3);
    n3 = m(1)*m(2)*(m(3)+1);

    j1 = m(1)*(m(2)+1)*(m(3)+1);
    j2 = (m(1)+1)*m(2)*(m(3)+1);
    j3 = (m(1)+1)*(m(2)+1)*m(3);

    p1 = n0 + j2 + j3;
    p2 = p1 + n0 + j1 + j3;
    p3 = p2 + n0 + j1 + j2;

    By = zeros(n1+n2+n3,1);

    % extract Divy
    Divy  = reshape(yc(end-n0+1:end),m(1),m(2),m(3));

    % ---------------------------------- part 1 of u
    % extract D11y1, boundary condition for tangential derivative
    uu = reshape(yc(1:n0),m(1),m(2),m(3));
    dummy = zeros(m(1)+1,m(2),m(3));
    dummy(1,:,:)   = -uu(1,:,:);
    dummy(end,:,:) =  uu(end,:,:);
    dummy(2:end-1,:,:) = uu(1:end-1,:,:) - uu(2:end,:,:);
    By(1:n1) = a*reshape(dummy,n1,1)/h(1);

    % extract D21y1, boundary condition for normal derivative
    uu = reshape(yc(n0+(1:j3)),m(1)+1,m(2)+1,m(3));
    uu(:,[1,end],:) = 0;
    dummy = uu(:,1:end-1,:) - uu(:,2:end,:);
    By(1:n1) = By(1:n1) + a*reshape(dummy,n1,1)/h(2);

    % extract D31u1, boundary condition for normal derivative
    uu = reshape(yc(n0+j3+(1:j2)),m(1)+1,m(2),m(3)+1);
    uu(:,:,[1,end]) = 0;
    dummy = uu(:,:,1:end-1) - uu(:,:,2:end);
    By(1:n1) = By(1:n1) + a*reshape(dummy,n1,1)/h(3);

    % Divy, boundary condition for tangential derivative
    dummy = zeros(m(1)+1,m(2),m(3));
    dummy(1,:,:)   = -Divy(1,:,:);
    dummy(end,:,:) =  Divy(end,:,:);
    dummy(2:end-1,:,:) = Divy(1:end-1,:,:) - Divy(2:end,:,:);
    By(1:n1) = By(1:n1) + b*reshape(dummy,n1,1)/h(1);

    % ---------------------------------- part 2 of y
    % extract D12y2, boundary condition for normal derivative
    uu = reshape(yc(p1+(1:j3)),m(1)+1,m(2)+1,m(3));
    uu([1,end],:,:) = 0;
    dummy = uu(1:end-1,:,:) - uu(2:end,:,:);
    By(n1+(1:n2)) = a*reshape(dummy,n2,1)/h(1);

    % extract D22y2, boundary condition for tangential derivative
    uu = reshape(yc(p1+j3+(1:n0)),m(1),m(2),m(3));
    dummy = zeros(m(1),m(2)+1,m(3));
    dummy(:,1,:)   = -uu(:,1,:);
    dummy(:,end,:) =  uu(:,end,:);
    dummy(:,2:end-1,:) = uu(:,1:end-1,:) - uu(:,2:end,:);
    By(n1+(1:n2)) = By(n1+(1:n2)) + a*reshape(dummy,n2,1)/h(2);

    % extract D32y2, boundary condition for normal derivative
    uu = reshape(yc(p1+j3+n0+(1:j1)),m(1),m(2)+1,m(3)+1);
    uu(:,:,[1,end]) = 0;
    dummy = uu(:,:,1:end-1) - uu(:,:,2:end);
    By(n1+(1:n2)) = By(n1+(1:n2)) + a*reshape(dummy,n2,1)/h(3);

    % Divy, boundary condition for tangential derivative
    dummy = zeros(m(1),m(2)+1,m(3));
    dummy(:,1,:)   = -Divy(:,1,:);
    dummy(:,end,:) =  Divy(:,end,:);
    dummy(:,2:end-1,:) = Divy(:,1:end-1,:) - Divy(:,2:end,:);
    By(n1+(1:n2)) = By(n1+(1:n2)) + b*reshape(dummy,n2,1)/h(2);

    % ---------------------------------- part 3 of y
    % extract D13y3, boundary condition for normal derivative
    uu = reshape(yc(p2+(1:j2)),m(1)+1,m(2),m(3)+1);
    uu([1,end],:,:) = 0;
    dummy = uu(1:end-1,:,:) - uu(2:end,:,:);
    By(n1+n2+(1:n3)) = a*reshape(dummy,n3,1)/h(1);

    % extract D23y3, boundary condition for normal derivative
    uu = reshape(yc(p2+j2+(1:j1)),m(1),m(2)+1,m(3)+1);
    uu(:,[1,end],:) = 0;
    dummy = uu(:,1:end-1,:) - uu(:,2:end,:);
    By(n1+n2+(1:n3)) = By(n1+n2+(1:n3)) + a*reshape(dummy,n3,1)/h(2);

    % extract D32y2, boundary condition for tangential derivative
    uu = reshape(yc(p2+j2+j1+(1:n0)),m(1),m(2),m(3));
    dummy = zeros(m(1),m(2),m(3)+1);
    dummy(:,:,1)   = -uu(:,:,1);
    dummy(:,:,end) =  uu(:,:,end);
    dummy(:,:,2:end-1) = uu(:,:,1:end-1) - uu(:,:,2:end);
    By(n1+n2+(1:n3)) = By(n1+n2+(1:n3)) + a*reshape(dummy,n3,1)/h(3);

    % Divy, boundary condition for tangential derivative
    dummy = zeros(m(1),m(2),m(3)+1);
    dummy(:,:,1)   = -Divy(:,:,1);
    dummy(:,:,end) =  Divy(:,:,end);
    dummy(:,:,2:end-1) = Divy(:,:,1:end-1) - Divy(:,:,2:end);
    By(n1+n2+(1:n3)) = By(n1+n2+(1:n3)) + b*reshape(dummy,n3,1)/h(3);

  otherwise, jmerror(flag)
end;

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

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% get image block - 4
function [img1, img2, img3, img4] = getImgBlock(img)

[h, w] = size(img);
h_cen = floor(h/2);
w_cen = floor(w/2);

img1 = img(1:h_cen+2, 1:w_cen+2);
img2 = img(1:h_cen+2, w_cen-1:w);
img3 = img(h_cen-1:h, 1:w_cen+2);
img4 = img(h_cen-1:h, w_cen-1:w);

% figure;imshow(img1);
% figure;imshow(img2);
% figure;imshow(img3);
% figure;imshow(img4);
%
% F1 = cat(2,img1, img2);
% figure;imshow(F1);
% F2 = cat(2,img3, img4);
% figure;imshow(F2);
% F = cat(1,F1, F2);
% figure;imshow(F);

end

author: thfylsty


text:
function H = hessian_2sided_nrows(f,x,k,varargin)
% Computes the last K rows of a 2-sided finite difference Hessian
%
% USAGE:
%   H=hessian_2sided_nrows(FUNC,X,K,VARARGIN)
%
% INPUTS:
%   FUNC         - Function name, fval = func(x,varargin)
%   X            - Vector of parameters (N x 1)
%   K            - Number of rows to compute
%   VARARGIN     - Optional arguments passed to the function
%
% OUTPUTS:
%   H            - Finite difference, 2-sided hessian: N x K
%
% COMMENTS:
%   Modification of hessian_2sided
%   See also HESSIAN_2SIDED

% Code originally from COMPECON toolbox [www4.ncsu.edu/~pfackler]
% documentation modified to fit the format of the Econometrics Toolbox
% by James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jlesage@spatial-econometrics.com
%
% Further modified (to do 2-sided numerical derivs, rather than 1) by:
% Modifications Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 9/1/2005

if size(x,2)>size(x,1)
    x=x';
end
if size(x,2)~=1
    error('X must be a column vector.')
end
n = size(x,1);

try
    feval(f,x,varargin{:});
catch FE
    errtxt=['There was an error evaluating the function.  Please check the arguements.  The specific error was:' FE.message];
    error(errtxt);
end

fx = feval(f,x,varargin{:});

% Compute the stepsize (h)
h = eps.^(1/3)*max(abs(x),1e-2);
xh = x+h;
h = xh-x;

ee = sparse(1:n,1:n,h,n,n);

% Compute forward and backward steps
gp = zeros(n,1);
gm = zeros(n,1);
for i=1:n
    gp(i) = feval(f,x+ee(:,i),varargin{:});
    gm(i) = feval(f,x-ee(:,i),varargin{:});
end

hh=h*h';
Hm=NaN*ones(n);
Hp=NaN*ones(n);
% Compute "double" forward and backward steps
for i=n-k+1:n
    for j=1:n
        Hp(i,j) = feval(f,x+ee(:,i)+ee(:,j),varargin{:});
        Hp(j,i)=Hp(i,j);
        Hm(i,j) = feval(f,x-ee(:,i)-ee(:,j),varargin{:});
        Hm(j,i)=Hm(i,j);
    end
end
%Compute the hessian
H = zeros(n);
for i=(n-k+1):n
    for j=1:n
        H(i,j) = (Hp(i,j)-gp(i)-gp(j)+fx+fx-gm(i)-gm(j)+Hm(i,j))/hh(i,j)/2;
        H(j,i) = H(i,j);
    end
end

%Only the final K rows
H=H((n-k+1):n,:);

author: bashtage

repo: mfe-toolbox


text:
function [ld, UC] = logdet(A, UC)

% LOGDET The log of the determinant when argument is positive definite.
% FORMAT
% DESC returns the log determinant of a positive definite matrix. If
% the matrix isn't quite positive definite the function adds 'jitter'
% to make it positive definite and gives out a warning message (this
% is done through JITCHOL).
% ARG A : the input positive definite matrix for which the log
% determinant is required.
% RETURN d : the log determinant of A computed using Cholesky
% decomposition.
% RETURN U : the Cholesky decomposition of A.
%
% FORMAT
% DESC returns the log determinant of a positive definite matrix given
% the Cholesky decomposition of A. If jitter is used then the
% amount of jitter used is returned.
% ARG A : the input positive definite matrix for which the log
% determinant is required.
% ARG U : the Cholesky decomposition of A.
% RETURN d : the log determinant of A computed using Cholesky
% decomposition.
% RETURN U : the Cholesky decomposition of A.
% RETURN jitter : the amount of jitter added.
%
% SEEALSO : jitChol, pdinv, chol
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2005, 2006

% NDLUTIL

if nargin < 2
  UC=[];
end

% Obtain a Cholesky decomposition.
if isempty(UC)
  if nargout > 2
    [UC, jitter] = jitChol(A);
  else
    UC = jitChol(A);
  end
end

ld = 2*sum(log(diag(UC)));


text:
function [muX, SXX] = marginalize_gaussian(mu, Sigma, X, Y, ns)
% MARGINALIZE_GAUSSIAN Compute Pr(X) from Pr(X,Y) where X and Y are jointly Gaussian.
% [muX, SXX] = marginalize_gaussian(mu, Sigma, X, Y, ns)

[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function stop = optimplotfvalKsptot(~,optimValues,state,varargin)
% OPTIMPLOTFVAL Plot value of the objective function at each iteration.
%
%   STOP = OPTIMPLOTFVAL(X,OPTIMVALUES,STATE) plots OPTIMVALUES.fval.  If
%   the function value is not scalar, a bar plot of the elements at the
%   current iteration is displayed.  If the OPTIMVALUES.fval field does not
%   exist, the OPTIMVALUES.residual field is used.
%
%   Example:
%   Create an options structure that will use OPTIMPLOTFVAL as the plot
%   function
%     options = optimset('PlotFcns',@optimplotfval);
%
%   Pass the options into an optimization problem to view the plot
%     fminbnd(@sin,3,10,options)

%   Copyright 2006-2010 The MathWorks, Inc.

stop = false;
switch state
    case 'iter'
        if isfield(optimValues,'fval')
            if isscalar(optimValues.fval)
                plotscalar(optimValues.iteration,optimValues.fval);
            else
                plotvector(optimValues.iteration,optimValues.fval);
            end
        else
            plotvector(optimValues.iteration,optimValues.residual);
        end
end

function plotscalar(iteration,fval)
% PLOTSCALAR initializes or updates a line plot of the function value
% at each iteration.
persistent plotfval plotfvalParent

if iteration == 0
    plotfval = plot(iteration,fval,'kd','MarkerFaceColor',[1 0 1]);
    plotfvalParent = plotfval.Parent;
    title(getString(message('MATLAB:optimfun:funfun:optimplots:TitleCurrentFunctionValue',sprintf('%g',fval))),'interp','none');
    xlabel(getString(message('MATLAB:optimfun:funfun:optimplots:LabelIteration')),'interp','none');
    set(plotfval,'Tag','optimplotfval');
    ylabel(getString(message('MATLAB:optimfun:funfun:optimplots:LabelFunctionValue')),'interp','none')
else
%     plotfval = findobj(get(gca,'Children'),'Tag','optimplotfval');
%     plotfval.Parent = plotfvalParent;
    newX = [get(plotfval,'Xdata') iteration];
    newY = [get(plotfval,'Ydata') fval];
    set(plotfval,'Xdata',newX, 'Ydata',newY);
    set(get(plotfvalParent,'Title'),'String',getString(message('MATLAB:optimfun:funfun:optimplots:TitleCurrentFunctionValue',sprintf('%g',fval))));
end

function plotvector(iteration,fval)
% PLOTVECTOR creates or updates a bar plot of the function values or
% residuals at the current iteration.
persistent plotfval

if iteration == 0
    xlabelText = getString(message('MATLAB:optimfun:funfun:optimplots:LabelNumberOfFunctionValues0',sprintf('%g',length(fval))));
    % display up to the first 100 values
    if numel(fval) > 100
        xlabelText = {xlabelText,getString(message('MATLAB:optimfun:funfun:optimplots:LabelShowingOnlyFirst100Values'))};
        fval = fval(1:100);
    end
    plotfval = bar(fval);
    title(getString(message('MATLAB:optimfun:funfun:optimplots:TitleCurrentFunctionValues')),'interp','none');
    set(plotfval,'edgecolor','none')
    set(plotfvalParent,'xlim',[0,1 + length(fval)])
    xlabel(xlabelText,'interp','none')
    set(plotfval,'Tag','optimplotfval');
    ylabel(getString(message('MATLAB:optimfun:funfun:optimplots:LabelFunctionValue')),'interp','none')
else
%     plotfval = findobj(get(gca,'Children'),'Tag','optimplotfval');
    % display up to the first 100 values
    if numel(fval) > 100
        fval = fval(1:100);
    end
    set(plotfval,'Ydata',fval);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function T=Entropy_Yen(X,Trange);
hg=hist(X,Trange);
pg=hg./sum(hg);
Ti=0;
TrangeSpan=Trange(pg>0);
JT=zeros(length(TrangeSpan),1);
for Tk=TrangeSpan
    Ti=Ti+1;
    PT=sum(pg(Trange<=Tk));
    Cf=-log10(sum((pg(Trange<=Tk)./PT).^2));
    Cb=-log10(sum((pg(Trange>Tk)./(1-PT)).^2));
    JT(Ti)=Cf+Cb;
end
JT(isinf(JT))=0;
JT(isnan(JT))=0;
[~,Tind]=max(JT);
T=TrangeSpan(Tind);
end

author: radishgiant

repo: ThresholdAndSegment


text:
function [distancemap]=surfacemap(Vertices,faces,Index)

%this function defines how many vertices each vertex is away from the
%specifc vertex with index "Index3".

%the first collumn gives the vertex indices, the second the number of
%vertices each vertex is distant from the Index vertix
%the third collumn gives the "weigth" for each, meaning the number of
%vertices of the same generation devided by the total number of vertices

faces(:,4)=1:length(faces);
distancemap(:,1)=1:length(Vertices(:,1));
distancemap(Index,2)=1;
distancemap(Index,3)=1-length(double(find(distancemap(:,2))))/length(Vertices(:,1));

a=2;
facestemp1(1,1)=Index;

while isempty(faces)==0

    indicesconnectedfaces=vertcat(find(double(ismember(faces(:,1),facestemp1))),find(double(ismember(faces(:,2),facestemp1))),find(double(ismember(faces(:,3),facestemp1))));
    connectedfaces=faces(indicesconnectedfaces,:);
    indicesvertices=unique(reshape(connectedfaces(:,1:3),3*length(connectedfaces(:,1)),1));
    distancemap(indicesvertices,2)=a;
    faces(indicesconnectedfaces,:)=[];

    indicesoldvertces=find(double(ismember(indicesvertices,facestemp1)));
    distancemap(indicesvertices(indicesoldvertces,:),2)=a-1;
    indicesvertices(indicesoldvertces,:)=[];
    facestemp1=indicesvertices;
    distancemap(indicesvertices,3)=1-length(double(find(distancemap(:,2))))/length(Vertices(:,1));
    a=a+1;
end

author: Sable

repo: mcbench-benchmarks


text:
function constant_test ( )

%*****************************************************************************80
%
%% CONSTANT_TEST examines the constant correlation.
%
%  Discussion:
%
%    This code is based substantially on a document by Toby Driscoll.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Toby Driscoll,
%    Mercer's theorem and the Karhunen-Loeve expansion,
%    Oxford University Mathematical Institute,
%    http://www2.maths.ox.ac.uk/chebfun/examples/stats/pdf/MercerKarhunenLoeve.pdf
%
  rmpath ( '/usr/local/matlab/toolbox/datafeed/datafeed' )
  addpath ( '../chebfun' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONSTANT_TEST:\n' );
  fprintf ( 1, '  Demonstrate Mercer''s theorem and the KL expansion\n' );
  fprintf ( 1, '  for the constant kernel.\n' );
%
%  Set the interval.
%
  a = 0.0;
  b = +10.0;
  fprintf ( 1, '  Using interval [%g,%g]\n', a, b );
  s_num = 21;
  s_vec = linspace ( a, b, s_num );
%
%  FRED is a function from the CHEBFUN library.
%
%  It constructs a "chebop" representing the Fredholm integral operator
%  with kernel K for functions in domain [A,B].
%
  F = fred ( @constant_correlation, domain ( [ a, b ] ) );
%
%  EIGS has been extended to be able to compute the eigenvalues Lambda
%  and eigenfunctions Psi of the Fredholm integral operator represented by F.
%
%  The "LM" switch requests that we return the eigenvalues of largest magnitude.
%
%  Each Psi is a CHEBFUN, that is, it takes a real number argument.
%
  eigen_num = 20;
  [ Psi, Lambda ] = eigs ( F, eigen_num, 'lm' );

  eigen_found = length ( Lambda );
  fprintf ( 1, '  Requested %d eigenmodes, computed %d\n', eigen_num, eigen_found );
  eigen_num = min ( eigen_num, eigen_found );
%
%  Print the eigenvalues.
%
  lambda_vec = diag ( Lambda );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I       Lambda(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %2d  %10g\n', i, lambda_vec(i) );
  end
%
%  Plot the eigenvalues.
%
  figure ( 1 )
  clf
  plot ( lambda_vec, 'Linewidth', 2 );
  hold on
  plot ( lambda_vec, 'b.', 'Markersize', 20 );
  title ( 'constant: Mercer eigenvalues' );
  xlabel ( '<--- N --->')
  grid on
  print -dpng 'constant_figure1.png'
%
%  Plot selected eigenfunctions.
%
  figure ( 2 )
  subplot ( 4, 1, 1 )
  plot ( Psi(:,1), 'Linewidth', 2 );
  title ( 'constant: Mercer eigenfunction PSI(1)' )
  grid on
  if ( 2 <= eigen_found )
    subplot ( 4, 1, 2 )
    plot ( Psi(:,2), 'Linewidth', 2 );
    title ( 'constant: Mercer eigenfunction PSI(2)' )
    grid on
  end
  if ( 3 <= eigen_found )
    subplot ( 4, 1, 3 )
    plot ( Psi(:,5), 'Linewidth', 2 );
    title ( 'constant: Mercer eigenfunction PSI(5)' )
    grid on
  end
  if ( 4 <= eigen_found )
    subplot ( 4, 1, 4 )
    plot ( Psi(:,10), 'Linewidth', 2 );
    title ( 'constant: Mercer eigenfunction PSI(10)' )
    grid on
  end
  print -dpng 'constant_figure2.png'
%
%  Orthonormality check.
%
  ptp = Psi' * Psi;
  error_frobenius = r8mat_is_identity ( eigen_num, ptp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of I - Psi'' * Psi = %g\n', error_frobenius );
%
%  K(S,S) should be exactly 1.
%  Because we are using a truncated representation of K, our estimate of K(S,S)
%  will be smaller than 1.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Truncated estimate of K(s,s) = 1 for S in the interval.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        S        K(s,s) estimate\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 21
    s = s_vec(i);
    ptlp = Psi(s,:) * Lambda * Psi(s,:)';
    fprintf ( 1, '  %10g  %14g\n', s, ptlp );
  end
%
%  Look at eigenvalue decay.
%
  x = diff ( log ( ( 1:eigen_num ) ) );
  y = diff ( log ( lambda_vec ) )';
  c = y ./ x;
  figure ( 3 );
  clf
  plot ( c, 'Linewidth', 2 );
  grid on
  hold on
  plot ( c, 'b.', 'Markersize', 20 );
  xlabel ( '<-- N -->' )
  title ( 'constant: Eigenvalue decay rate' );
  print -dpng 'constant_figure3.png'
%
%  Look at eigenvalue sum.
%
%  The trace of K(s,s) over [a,b] should be the integral of 1 over [a,b],
%  that is, b - a.
%
%  We compare the trace to the partial sums of the lambda's to see how much
%  of the variance of the process we have captured.
%
  trace_K = b - a;
  lambda_cum = cumsum ( lambda_vec ) / trace_K;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index    Cumulative Eigenvalue sum\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %5d  %10g\n', i, lambda_cum(i) );
  end
%
%  We decide to use the first 10 eigenfunctions.
%
  eigen_use = min ( 10, eigen_found );
%
%  Find 400 realizations of the process by selecting, for each realization,
%  10 random parameters Z in the truncated KL expansion.
%
  Z = randn ( eigen_use, 400 );
  X = Psi(:,1:eigen_use) * ( sqrt ( Lambda(1:eigen_use,1:eigen_use) ) * Z );
%
%  Plot 40 of the realizations;
%  Plot their mean, computed from all 400.
%
  figure ( 4 )
  clf
  plot ( X(:,1:40) )
  mu = sum ( X, 2 ) / 400;
  hold on;
  plot ( mu, 'k', 'linewidth', 3 );
  title ( 'constant: 40 Random Realizations X(t,omega), and their Mean.' )
  print -dpng 'constant_figure4.png'
%
%  Estimate the covariance from the data.
%
  figure ( 5 )
  clf
  [ S, T ] = meshgrid ( s_vec, s_vec );
  C = cov ( X(s_vec,:)' );
  mesh ( S, T, C );
  hold on;
  D = constant_correlation ( S, T );
  plot3 ( S, T, D, 'k.', 'markersize', 10 )
  title ( 'constant: Covariance K(S,T) (dots), and Estimate from Realizations (mesh)' )
  print -dpng 'constant_figure5.png'
%
%  Using just 10 functions in the expansion,
%  reduce the correlation length,
%  and examine the sum of the lambda's.
%
  if ( 0 )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use a fixed number of eigenfunctions = %d\n', eigen_use );
  fprintf ( 1, '  but vary the correlation length RHOBAR.\n' );
  fprintf ( 1, '  (We used RHOBAR = 1 above.)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The sum of the eigenvalues, divided by (B-A),\n' );
  fprintf ( 1, '  discloses the relative amount of the total variation\n' );
  fprintf ( 1, '  that is captured by the truncated expansion.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      RHOBAR      VARSUM\n' );
  fprintf ( 1, '\n' );
  rhobar = 4.0;
  for i = 1 : 10
    F = fred ( @constant_correlation, domain ( [ a, b ] ) );
    lambda_vec = eigs ( F, 20, 'lm' );
    varsum = sum ( lambda_vec(1:eigen_use) ) / ( b - a );
    fprintf ( 1, '  %10g  %10g\n', rhobar, varsum );
    rhobar = rhobar / 2.0;
  end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONSTANT_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../chebfun' )

  return
end


text:
% shortread() - Read matrix from short file.
%
% Usage:
%  >> A = shortread(filename,size,'format',offset)
%
% Inputs:
%  filename - Read matrix a from specified file while assuming four byte
%             short integers.
%  size     - The vector SIZE determine the number of short elements to be
%             read and the dimensions of the resulting matrix. If the last
%             element of SIZE is INF the size of the last dimension is determined
%             by the file length.
%
% Optional inputs:
% 'format'  - The option FORMAT argument specifies the storage format as
%             defined by fopen(). Default format ([]) is 'native'.
% offset    - The option OFFSET is offset in shorts from the beginning of file (=0)
%             to start reading (2-byte shorts assumed).
%             It uses fseek to read a portion of a large data file.
%
% Author: Sigurd Enghoff, CNL / Salk Institute, La Jolla, 7/1998
%
% See also: floatread(), floatwrite(), fopen()

% Copyright (C) Sigurd Enghoff, CNL / Salk Institute, La Jolla, 7/1998
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% xx-07-98  written by Sigurd Enghoff, Salk Institute
% 04-26-99  modified by Sigurd Enghoff to handle variable-sized and
%           multi-dimensional data.
% 07-08-99  modified by Sigurd Enghoff, FORMAT argument added.
% 02-08-00  help updated for toolbox inclusion -sm
% 02-14-00  added segment arg -sm

function A = shortread(fname,Asize,fform,offset)

if ~exist('fform') | isempty(fform)|fform==0
	fform = 'native';
end

fid = fopen(fname,'rb',fform);
if fid>0
 if exist('offset')
   stts = fseek(fid,2*offset,'bof');
   if stts ~= 0
     error('shortread(): fseek() error.');
     return
   end
 end
 A = fread(fid,prod(Asize),'short');
else
 error('shortread(): fopen() error.');
 return
end
fprintf('   %d shorts read\n',prod(size(A)));

if Asize(end) == Inf
	Asize = Asize(1:end-1);
	A = reshape(A,[Asize length(A)/prod(Asize)]);
else
	A = reshape(A,Asize);
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function t_most_mpopf(quiet)
%T_MOST_MPOPF  Tests of deterministic multiperiod DC OPF problems
%
%   Cases taken from most_ex5_mpopf.m.

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

if nargin < 1
    quiet = 0;      %% verbose by default
end

ntests = 42;
t_begin(ntests, quiet);

if quiet
    verbose = 0;
else
    verbose = 0;
end
% verbose = 2;

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);
    near_sing_matrix_warn_id = 'Octave:nearly-singular-matrix';
else
    near_sing_matrix_warn_id = 'MATLAB:nearlySingularMatrix';
end
s = warning('query', near_sing_matrix_warn_id);
warning('off', near_sing_matrix_warn_id);

casefile = 'ex_case3b';
solnfile =  't_most_mpopf_soln';
soln = load(solnfile);
mpopt = mpoption('verbose', verbose);
mpopt = mpoption(mpopt, 'out.gen', 1);
mpopt = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'most.solver', 'MIPS');
mpopt = mpoption(mpopt, 'most.dc_model', 1);
% mpopt = mpoption(mpopt, 'opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
%         'mips.comptol', 1e-8, 'mips.costtol', 1e-8);
if have_feature('octave') && have_feature('octave', 'vnum') > 7
    mpopt = mpoption(mpopt, 'mips.linsolver', 'LU');
end
if ~verbose
    mpopt = mpoption(mpopt, 'out.all', 0);
end
% mpopt = mpoption(mpopt, 'out.all', -1);

%% load base case file
mpc = loadcase(casefile);

nb = size(mpc.bus, 1);
nl = size(mpc.branch, 1);
ng = size(mpc.gen, 1);

xgd = loadxgendata('ex_xgd_ramp', mpc);
[iwind, mpc, xgd] = addwind('ex_wind', mpc, xgd);
profiles = getprofiles('ex_wind_profile_d', iwind);
profiles = getprofiles('ex_load_profile', profiles);
nt = size(profiles(1).values, 1);       % number of periods

mpc_full = mpc;
xgd_full = xgd;

%%-----  Multiperiod DC OPF (w/ramp)  -----
t = 'Multiperiod DC OPF w/ramp : ';
mdi = loadmd(mpc, nt, xgd, [], [], profiles);
mdo = most(mdi, mpopt);
% EPg1 = mdo.results.ExpectedDispatch;    % active generation
% Elam1 = mdo.results.GenPrices;          % nodal energy price
t_is(mdo.results.success, 1, 12, [t 'success']);
ms = most_summary(mdo);
ex = soln.wramp;
t_is(ms.f, ex.f, 1.8, [t 'f']);
t_is(ms.Pg, ex.Pg, 3.5, [t 'Pg']);
t_is(ms.Rup, ex.Rup, 2.5, [t 'Rup']);
t_is(ms.Rdn, ex.Rdn, 2.5, [t 'Rdn']);
t_is(ms.Pf, ex.Pf, 4, [t 'Pf']);
t_is(ms.u, ex.u, 8, [t 'u']);
t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
t_is(ms.muF, ex.muF, 5, [t 'muF']);
ex = soln.rramp;
t_is(mdo.results.GenPrices, ex.GenPrices, 3, [t 'GenPrices']);
t_is(mdo.results.GenTLMP, ex.GenTLMP, 2, [t 'TLMP']);
% rramp = mdo.results;
% wramp = most_summary(mdo);
% keyboard;

%%-----  Multiperiod DC OPF (w/ramp+wear/tear)  -----
t = 'Multiperiod DC OPF w/ramp+wear/tear : ';
xgd.RampWearCostCoeff(1:3) = 1;
mdi = loadmd(mpc, nt, xgd, [], [], profiles);
mdo = most(mdi, mpopt);
% EPg1 = mdo.results.ExpectedDispatch;    % active generation
% Elam1 = mdo.results.GenPrices;          % nodal energy price
t_is(mdo.results.success, 1, 12, [t 'success']);
ms = most_summary(mdo);
ex = soln.wwear;
t_is(ms.f, ex.f, 0.5, [t 'f']);
t_is(ms.Pg, ex.Pg, 2.8, [t 'Pg']);
t_is(ms.Rup, ex.Rup, 2.5, [t 'Rup']);
t_is(ms.Rdn, ex.Rdn, 2.5, [t 'Rdn']);
t_is(ms.Pf, ex.Pf, 3, [t 'Pf']);
t_is(ms.u, ex.u, 8, [t 'u']);
t_is(ms.lamP, ex.lamP, 2.5, [t 'lamP']);
t_is(ms.muF, ex.muF, 2, [t 'muF']);
ex = soln.rwear;
t_is(mdo.results.GenPrices, ex.GenPrices, 2.8, [t 'GenPrices']);
t_is(mdo.results.GenTLMP, ex.GenTLMP, 1, [t 'TLMP']);
% rwear = mdo.results;
% wwear = most_summary(mdo);
% keyboard;

t = 'build model without solving : ';
mpopt = mpoption(mpopt, 'most.build_model', 1, ...
                        'most.solve_model', 0, ...
                        'most.resolve_new_cost', 0);
mdo = most(mdi, mpopt);
t_ok(isfield(mdo, 'om'), [t '''om'' field']);
t_ok(isfield(mdo, 'QP'), [t '''QP'' field']);
t_ok(~isfield(mdo.QP, 'x'), [t 'no ''QP.x'' field']);
t_ok(~isfield(mdo.QP, 'f'), [t 'no ''QP.f'' field']);
t_ok(~isfield(mdo.QP, 'exitflag'), [t 'no ''QP.exitflag'' field']);
t_ok(~isfield(mdo.QP, 'output'), [t 'no ''QP.output'' field']);
t_ok(isfield(mdo.results, 'SetupTime'), [t '''results.SetupTime'' field']);
t_ok(~isfield(mdo.results, 'success'), [t 'no ''results.success'' field']);

t = 'solve previously built model : ';
mpopt = mpoption(mpopt, 'most.build_model', 0, ...
                        'most.solve_model', 1, ...
                        'most.resolve_new_cost', 1);
mdi1 = mdo;
mdo = most(mdi1, mpopt);
t_is(mdo.results.success, 1, 12, [t 'success']);
ms = most_summary(mdo);
ex = soln.wwear;
t_is(ms.f, ex.f, 0.5, [t 'f']);
t_is(ms.Pg, ex.Pg, 2.8, [t 'Pg']);
t_is(ms.Rup, ex.Rup, 2.5, [t 'Rup']);
t_is(ms.Rdn, ex.Rdn, 2.5, [t 'Rdn']);
t_is(ms.Pf, ex.Pf, 3, [t 'Pf']);
t_is(ms.u, ex.u, 8, [t 'u']);
t_is(ms.lamP, ex.lamP, 2.5, [t 'lamP']);
t_is(ms.muF, ex.muF, 2, [t 'muF']);
ex = soln.rwear;
t_is(mdo.results.GenPrices, ex.GenPrices, 2.8, [t 'GenPrices']);
t_is(mdo.results.GenTLMP, ex.GenTLMP, 1, [t 'TLMP']);
% wwear = most_summary(mdo);
% keyboard;

t = 'output model is copy of input model';
mdo.om.add_var('test', 10);
t_is(mdo.om.var.N, mdi1.om.var.N+10, 12, t);

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

t_end;

% save t_most_mpopf_soln wramp wwear rramp rwear

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   ABB IRB140.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   09/01/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function robot = parameters()

robot.name= 'Epson_Prosix_VT6';
robot.path = 'robots/Epson/Prosix_VT6';

robot.DH.theta= '[q(1)+pi q(2)-pi/2 q(3) q(4) q(5) q(6)]';
robot.DH.d='[0.412 0 0 0.400 0 0.080]';
robot.DH.a='[0.1 0.420 0 0 0 0]';
robot.DH.alpha= '[-pi/2 0 -pi/2 pi/2 -pi/2 0]';
robot.J=[];

robot.inversekinematic_fn = 'inversekinematic_Prosix_VT6(robot, T)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%number of degrees of freedom
robot.DOF = 6;

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-170) deg2rad(170); %Axis 1, minimum, maximum
                deg2rad(-160) deg2rad(65); %Axis 2, minimum, maximum
                deg2rad(-51) deg2rad(190); %Axis 3
                deg2rad(-200) deg2rad(200); %Axis 4: Unlimited (400� default)
                deg2rad(-125) deg2rad(125); %Axis 5
                deg2rad(-360) deg2rad(360)]; %Axis 6: Really Unlimited to (800� default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(183.7); %Axis 1, rad/s
                deg2rad(122.5); %Axis 2, rad/s
                deg2rad(118.8); %Axis 3, rad/s
                deg2rad(271.4); %Axis 4, rad/s
                deg2rad(296.8); %Axis 5, rad/s
                deg2rad(296.8)];%Axis 6, rad/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

% end effectors maximum velocity
robot.linear_velmax = 2.5; %m/s

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [94 104 245]./255;
%for transparency
robot.graphical.draw_transparent=30;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=0.5;
%adjust for a default view of the robot
robot.axis=[-2 2 -2 5 0 2];
%read graphics files
robot = read_graphics(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DYNAMIC PARAMETERS
%   WARNING! These parameters do not correspond to the actual IRB 140
%   robot. They have been introduced to demonstrate the necessity of
%   simulating the robot and should be used only for educational purposes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[16.1 17.16 9.18 6.12 1.53 0.9];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[0       0          0; %(rx, ry, rz) link 1
                     -0.05	 0.006	 0.1; %(rx, ry, rz) link 2
                    -0.0203	-0.0141	 0.070;  %(rx, ry, rz) link 3
                     0       0.019       0;%(rx, ry, rz) link 4
                     0       0           0;%(rx, ry, rz) link 5
                     0       0         0.032];%(rx, ry, rz) link 6

%Inertia matrices of each link with respect to its D-H reference system.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, for each row
robot.dynamics.Inertia=[0      0.35	0   	0	0	0;
    .13     .524	.539	0	0	0;
    .066	.086	.0125	0	0	0;
    1.8e-3	1.3e-3	1.8e-3	0	0	0;
    .3e-3	.4e-3	.3e-3	0	0	0;
    .15e-3	.15e-3	.04e-3	0	0	0];

robot.motors=load_motors([5 5 5 4 4 4]);
%Speed reductor at each joint
% my values are
% G=[230 450 20 1 5 1]
robot.motors.G=[300 300 300 300 300 300];

repo: ARTE


text:
%%Ex. 10 Continuation to next line
summation1 = 1 + 3 + 5 + 7 ...
 + 9 + 11

%Note: The three periods (...) allow continuation to the next line of commands. The two
%  lines in the above example are essentially one line of "summation1 = 1+3+5+7+9+11".

author: TheAlgorithms

repo: MATLAB-Octave


text:
%compute vx_cm
function [data,units]=compute_vx_cm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
vx_cm=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    vx_cm{1,i}=(trx(larva).x_mm(2:end)-trx(larva).x_mm(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=vx_cm;

author: kristinbranson

repo: JAABA


text:
function [varargout]=icolorbar(varargin)

% function [h,hc]=icolorbar(cLim,h)
%-------------------------------------------------------------------------
%
%
% Change log:
% 2018/05/10 Fixed bug in relation to older matlab versions which do not
% support h.Colormap=c for axis specific color maps. Older versions use
% colormap(h,c) instead.
% 2019/08/14 Added colorbar handle output
%-------------------------------------------------------------------------

%% parse input

switch nargin
    case 0
        cLim=caxis;
        h=[];
    case 1
        cLim=varargin{1};
        h=[];
    case 2
        cLim=varargin{1};
        h=varargin{2};
end

if isempty(cLim)
    cLim=caxis;
end

if isempty(h)
    %Get current figure or open a new one
    if isempty(findobj('type','figure'))
        hf=cFigure; %Open a new figure
    else
        hf=gcf; %Get the current figure
    end
   figure(hf);
   if isempty(hf.CurrentAxes) %isempty(get(hf,'CurrentAxes'))
       h=axes;
   else
       h=gca;
   end
else
    switch class(h)
        case {'matlab.ui.Figure'} %A figure handle
            hf=h;
            if isempty(hf.CurrentAxes) %isempty(get(hf,'CurrentAxes'))
                h=axes;
            else
                h=gca;
            end
        case {'matlab.graphics.axis.Axes'} %An axes handle
            hf=h.Parent; %Get figure handle
        otherwise
            error('Handle does not represent an axis or figure handle');
    end
end

%% Get colormap
if isprop(h,'Colormap') %Available in MATLAB 2018
    c=h.Colormap;
    matlabNew=1;
else %Older MATLAB versions
    c=colormap(h);
    matlabNew=0;
end

%% Set limits and tick labels
caxis([cLim(1)-0.5 cLim(2)+0.5]);
hc=colorbar;
hc.Ticks=cLim(1):1:cLim(2);
cn=resampleColormap(c,numel(hc.Ticks));
if matlabNew==1 %axis colormap
    h.Colormap=cn;
else %Figure colormap
    colormap(h,cn);
end

%% Collect output
switch nargout
    case 1
        varargout{1}=h;
    case 2
        varargout{1}=h;
        varargout{2}=hc;
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function sF = mrdivide(sF,s)
% overload the / operator, i.e. one can now write @S2Fun / 2  in order
% to scale an S2Fun.
%
% Syntax
%   sF = sF / 2
%
% Input
%  sF - @S2Fun
%  s  - constant, vector
%
% Output
%  sF - @S2Fun
%

if ~isnumeric(s)
  error('Second argument has to be numeric. Use ./ instead.')
end

sF = times(sF,inv(s));

end

author: mtex-toolbox

repo: mtex

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


text:
function square_symq_rule_test04 ( degree, n )

%*****************************************************************************80
%
%% SQUARE_SYMQ_RULE_TEST04 gets a rule and tests its accuracy.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    01 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired total polynomial degree exactness
%    of the quadrature rule.  0 <= DEGREE <= 50.
%
%    Input, integer N, the number of nodes to be used by the rule.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_SYMQ_RULE_TEST04\n' );
  fprintf ( 1, '  Get a quadrature rule for the symmetric square.\n' );
  fprintf ( 1, '  Test its accuracy.\n' );
  fprintf ( 1, '  Polynomial exactness degree DEGREE = %d\n', degree );
%
%  Retrieve a symmetric quadrature rule.
%
  [ x, w ] = square_symq ( degree, n );

  npols = ( ( degree + 1 ) * ( degree + 2 ) ) / 2;

  rints = zeros(npols,1);

  for i = 1 : n

    z(1) = x(1,i);
    z(2) = x(2,i);

    pols = lege2eva ( degree, z );

    rints(1:npols) = rints(1:npols) + w(i) * pols(1:npols);

  end

  area = 4.0;

  d = 0.0;
  d = ( rints(1) - sqrt ( area ) )^2;
  for i = 2 : npols
    d = d + rints(i)^2;
  end
  d = sqrt ( d ) / npols;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  RMS error = %g\n', d );

  return
end


text:
function b = r8sr_vxm ( n, nz, row, col, diag, off, x )

%*****************************************************************************80
%
%% R8SR_VXM multiplies a vector times a R8SR matrix.
%
%  Discussion:
%
%    The R8SR storage format stores the diagonal of a sparse matrix in DIAG.
%    The off-diagonal entries of row I are stored in entries ROW(I)
%    through ROW(I+1)-1 of OFF.  COL(J) records the column index
%    of the entry in A(J).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer NZ, the number of offdiagonal nonzero elements in A.
%
%    Input, integer ROW(N+1).  The nonzero offdiagonal elements of row I of A
%    are contained in A(ROW(I)) through A(ROW(I+1)-1).
%
%    Input, integer COL(NZ), contains the column index of the element
%    in the corresponding position in A.
%
%    Input, real DIAG(N), the diagonal elements of A.
%
%    Input, real OFF(NZ), the off-diagonal elements of A.
%
%    Input, real X(N), the vector to be multiplies by A.
%
%    Output, real B(N), the product A' * X.
%
  b(1:n) = diag(1:n) .* x(1:n);

  for i = 1 : n
    for k = row(i) : row(i+1)-1
      j = col(k);
      b(j) = b(j) + off(k) * x(i);
    end
  end

  return
end


text:
function [x,y] = intline(x1, x2, y1, y2)
%INTLINE Integer-coordinate line drawing algorithm.
%   [X, Y] = INTLINE(X1, X2, Y1, Y2) computes an
%   approximation to the line segment joining (X1, Y1) and
%   (X2, Y2) with integer coordinates.  X1, X2, Y1, and Y2
%   should be integers.  INTLINE is reversible; that is,
%   INTLINE(X1, X2, Y1, Y2) produces the same results as
%   FLIPUD(INTLINE(X2, X1, Y2, Y1)).
%
%   Function adapted from the 'strel' function of matlab.
%

dx = abs(x2 - x1);
dy = abs(y2 - y1);

% Check for degenerate case
if dx == 0 && dy == 0
    x = x1;
    y = y1;
    return;
end

flip = 0;
if dx >= dy
    if x1 > x2
        % swap coordinates to draw from left to right.
        t = x1; x1 = x2; x2 = t;
        t = y1; y1 = y2; y2 = t;
        flip = 1;
    end

    % compute line slope, and y from x
    m = (y2 - y1) / (x2 - x1);
    x = (x1:x2).';
    y = round(y1 + m * (x - x1));

else
    if y1 > y2
        % swap coordinates to draw from bottom to top.
        t = x1; x1 = x2; x2 = t;
        t = y1; y1 = y2; y2 = t;
        flip = 1;
    end

    % compute line slope, and x from y
    m = (x2 - x1) / (y2 - y1);
    y = (y1:y2).';
    x = round(x1 + m * (y - y1));
end

if flip
    x = flipud(x);
    y = flipud(y);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% CANFIS Grid Feed-Forward operation.
function y = canfis_grid_forward(x,mean,sigma,b,ThetaL4)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        								 							     				                                       %
%   			                 	NETWORK FUNCTIONALITY SECTION					   %
%      																					                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[NumInVars NumInTerms] = size(mean);

NumRules = NumInTerms^NumInVars;

% LAYER 1 - INPUT TERM NODES
In2 = x*ones(1,NumInTerms);
Out1 = 1./(1 + (abs((In2-mean)./sigma)).^(2*b));

% LAYER 2 - PRODUCT NODES
 precond = comb(Out1);
 Out2 = prod(precond,2);
 S_2 = sum(Out2);

% LAYER 3 - NORMALIZATION NODES
 Out3 = Out2./S_2;

% LAYERS 4 - 5: CONSEQUENT NODES - SUMMING NODE
Aux1 = [x; 1]*Out3';

a = reshape(Aux1,(NumInVars+1)*NumRules,1);  % New Input Training Data shaped as a column vector.

y = ThetaL4'*a;


text:
function test_approx_test08 ( )

%*****************************************************************************80
%
%% TEST_APPROX_TEST08 plots a cubic spline interpolant for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
  approx_filename = 'test08_approx.txt';
  data_filename = 'test08_data.txt';
  jmax = 7;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_APPROX_TEST08\n' );
  fprintf ( 1, '  Plot a cubic spline interpolant for problem 7.\n' );

  prob = 7;
%
%  Get the data.
%
  data_num = p00_data_num ( prob );

  [ xdata, ydata ] = p00_dat ( prob, data_num );

  r8vec2_write ( data_filename, data_num, xdata, ydata );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data values stored in "%s"\n', data_filename );
%
%  Set up the interpolation function.
%
  ibcbeg = 0;
  ibcend = 0;
  ybcbeg = 0.0;
  ybcend = 0.0;

  ypp = spline_cubic_set ( data_num, xdata, ydata, ibcbeg, ybcbeg, ibcend, ...
    ybcend );
%
%  Evaluate the interpolation function.
%
  plot = 0;
  nplot = ( jmax - 1 ) * ( data_num - 1 ) + 1;
  xplot = zeros ( nplot, 1 );
  yplot = zeros ( nplot, 1 );

  for i = 1 : data_num - 1

    if ( i == data_num - 1 )
      jhi = jmax;
    else
      jhi = jmax - 1;
    end

   for j = 1 : jhi

      xval = ( ( jmax - j     ) * xdata(i)     ...
             + (        j - 1 ) * xdata(i+1) ) ...
             / ( jmax     - 1 );

      [ yval, ypval, yppval ] = spline_cubic_val ( data_num, xdata, ydata, ...
        ypp, xval );

      plot = plot + 1;
      xplot(plot) = xval;
      yplot(plot) = yval;

    end

  end

  r8vec2_write ( approx_filename, nplot, xplot, yplot );

  fprintf ( 1, '  Approximant values stored in "%s"\n', approx_filename );

  return
end


text:
function thiem_test
% determination of transmissivity  T                    Holzbecher October 2005
% using Thiem's solution
% and fzero for automatic estimation
global rfit sfit reach Q

% specify fitting data
rfit = [0.8 30 90 215];                     % observation well distance
sfit = [2.236 1.088 0.716 0.25];            % observated steady state drawdowns
Q = 788;                                    % pumping rate
rwell = 1;                                  % well radius
reach = 500;                                % reach
T = 500;                                    % guess for transmissivity
test = 0;                                   % test run flag

% -------------------------------------------------------------------------
if test == 1               % input pumptest 'Oude Korendijk'
    rfit = [0.8 30 90 215]; sfit = [2.236 1.088 0.716 0.25];
    Q = 788; rwell = 0.5; T = 400; reach = 500;
end
options = optimset('Display','iter','TolFun',1e-8);
T = fzero(@myfun,T,options);
res = sum(Q*log(rfit/reach)/T/2/pi+sfit);
display (['Residual = ' num2str(res)]);
r = linspace (rwell,reach,100);
h = Q*log(r/reach)/T/2/pi;
plot (rfit,-sfit,'o',r,h,'-');
legend ('measured','ideally',4);
xlabel ('r'); ylabel ('h');
title ('Thiem Solution');
text (10*rwell,-0.1*max(sfit),['T = ' num2str(T)]);

function f = myfun(T);
global rfit sfit reach Q

% calculate Thiem solution
h = Q*log(rfit/reach)/T/2/pi;

% specify function f to vanish
f = (h+sfit)*h'/T;


text:
function [S,m,psi] = ComputeStokesQuad(HH,HV,VH,VV,SmoothSize)
%COMPUTESTOKESQUAD: Alpha/Entropy approximation for quad pol data (used
%Stokes for dual so we're keeping the name for consistency)

%Pauli basis
P_11 = fastrunmean(abs(HH+VV).*abs(HH+VV),[SmoothSize SmoothSize],'mean');
P_22 = fastrunmean(abs(HH-VV).*abs(HH-VV),[SmoothSize SmoothSize],'mean');
P_33 = fastrunmean(abs(HV-VH).*abs(HV-VH),[SmoothSize SmoothSize],'mean');

%span is the trace of the Pauli basis
S = P_11 + P_22 + P_33;

N_11 = P_11./S;
N_22 = P_22./S;
N_33 = P_33./S;

Frob = sqrt(N_11.*N_11 + N_22.*N_22 + N_33.*N_33);
m = 1.5*(1-Frob.*Frob);
psi = acosd(sqrt(N_11));
%psi = (pi/2)*(1-N_11);

end

author: ngageoint


text:
function im = bmpLoader(im_path)
im = imread(im_path, 'bmp');
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function geometry_test200 ( )

%*****************************************************************************80
%
%% TEST200 tests SPHERE_TRIANGLE_SIDES_TO_ANGLES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  r = 10.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST200\n' );
  fprintf ( 1, '  SPHERE_TRIANGLE_SIDES_TO_ANGLES takes the sides of a\n' );
  fprintf ( 1, '  spherical triangle and determines the angles.\n' );

  as = 121.0 + ( 15.4 / 60.0 );
  bs = 104.0 + ( 54.7 / 60.0 );
  cs =  65.0 + ( 42.5 / 60.0 );

  as = degrees_to_radians ( as );
  bs = degrees_to_radians ( bs );
  cs = degrees_to_radians ( cs );

  as = r * as;
  bs = r * bs;
  cs = r * cs;
%
%  Get the spherical angles.
%
  [ a, b, c ] = sphere_triangle_sides_to_angles ( r, as, bs, cs );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A       = %8f (radians)\n', a );
  a = radians_to_degrees ( a );
  fprintf ( 1, '  A       = %8f (degrees)\n', a );
  a = 117.0 + ( 58.0 / 60.0 );
  fprintf ( 1, '  Correct = %8f (radians)\n', a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  B       = %8f (radians)\n', b );
  b = radians_to_degrees ( b );
  fprintf ( 1, '  B       = %8f (degrees)\n', b );
  b = 93.0 + ( 13.8 / 60.0 );
  fprintf ( 1, '  Correct = %8f (radians)\n', b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  C       = %8f (radians)\n', c );
  c = radians_to_degrees ( c );
  fprintf ( 1, '  C       = %8f (degrees)\n', c );
  c = 70.0 + ( 20.6 / 60.0 );
  fprintf ( 1, '  Correct = %8f (radians)\n', c );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [xdx, xdy, ydx, ydy] = jacobian_matrix2d(vectorField)
% JACOBIAND_MATRIX2D Estimates the elements of the jacobian matrix
%
% function [xdx xdy ydx ydy] = jacobian_matrix2d(vectorField)
%
% INPUT ARGUMENTS
% vectorField       - Vectorfield to estimate the Jacobian for
%
% OPTIONAL INPUT ARGUMENTS
% N/A
%
% OUTPUT ARGUMENTS
% xdx ...           - Matrix elements of the Jacobian

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

[xdx xdy] = gradient(vectorField{1});
[ydx ydy] = gradient(vectorField{2});

author: fordanic

repo: image-registration


text:
function z = sandwich3x3(x, y)

% SANDWICH3X3 compute x*y*x' provided y is Hermitian and dimensionality is 3x3xN

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

% FIXME build in check for hermitianity
z     = complex(zeros(size(x)));
xconj = conj(x);
xabs2 = abs(x).^2;

z(1,1,:,:) = xabs2(1,1,:,:).*y(1,1,:,:)  + ...
             xabs2(1,2,:,:).*y(2,2,:,:)  + ...
             xabs2(1,3,:,:).*y(3,3,:,:)  + ...
             2.*real( x(1,2,:,:).*y(2,1,:,:).*xconj(1,1,:,:) ) + ...
             2.*real( x(1,3,:,:).*y(3,1,:,:).*xconj(1,1,:,:) ) + ...
             2.*real( x(1,3,:,:).*y(3,2,:,:).*xconj(1,2,:,:) );

             %(x(1,2,:,:).*y(2,1,:,:) + x(1,3,:,:).*y(3,1,:,:)).*xconj(1,1,:,:) + ...
             %(x(1,1,:,:).*y(1,2,:,:) + x(1,3,:,:).*y(3,2,:,:)).*xconj(1,2,:,:) + ...
             %(x(1,1,:,:).*y(1,3,:,:) + x(1,2,:,:).*y(2,3,:,:)).*xconj(1,3,:,:);

z(2,1,:,:) = (x(2,1,:,:).*y(1,1,:,:) + x(2,2,:,:).*y(2,1,:,:) + x(2,3,:,:).*y(3,1,:,:)).*xconj(1,1,:,:) + ...
             (x(2,1,:,:).*y(1,2,:,:) + x(2,2,:,:).*y(2,2,:,:) + x(2,3,:,:).*y(3,2,:,:)).*xconj(1,2,:,:) + ...
             (x(2,1,:,:).*y(1,3,:,:) + x(2,2,:,:).*y(2,3,:,:) + x(2,3,:,:).*y(3,3,:,:)).*xconj(1,3,:,:);

z(3,1,:,:) = (x(3,1,:,:).*y(1,1,:,:) + x(3,2,:,:).*y(2,1,:,:) + x(3,3,:,:).*y(3,1,:,:)).*xconj(1,1,:,:) + ...
             (x(3,1,:,:).*y(1,2,:,:) + x(3,2,:,:).*y(2,2,:,:) + x(3,3,:,:).*y(3,2,:,:)).*xconj(1,2,:,:) + ...
             (x(3,1,:,:).*y(1,3,:,:) + x(3,2,:,:).*y(2,3,:,:) + x(3,3,:,:).*y(3,3,:,:)).*xconj(1,3,:,:);

z(1,2,:,:) = conj(z(2,1,:,:));

z(2,2,:,:) = xabs2(2,1,:,:).*y(1,1,:,:)  + ...
             xabs2(2,2,:,:).*y(2,2,:,:)  + ...
             xabs2(2,3,:,:).*y(3,3,:,:)  + ...
             2.*real( x(2,2,:,:).*y(2,1,:,:).*xconj(2,1,:,:) ) + ...
             2.*real( x(2,3,:,:).*y(3,1,:,:).*xconj(2,1,:,:) ) + ...
             2.*real( x(2,3,:,:).*y(3,2,:,:).*xconj(2,2,:,:) );

             %(x(2,2,:,:).*y(2,1,:,:) + x(2,3,:,:).*y(3,1,:,:)).*xconj(2,1,:,:) + ...
             %(x(2,1,:,:).*y(1,2,:,:) + x(2,3,:,:).*y(3,2,:,:)).*xconj(2,2,:,:) + ...
             %(x(2,1,:,:).*y(1,3,:,:) + x(2,2,:,:).*y(2,3,:,:)).*xconj(2,3,:,:);

z(3,2,:,:) = (x(3,1,:,:).*y(1,1,:,:) + x(3,2,:,:).*y(2,1,:,:) + x(3,3,:,:).*y(3,1,:,:)).*xconj(2,1,:,:) + ...
             (x(3,1,:,:).*y(1,2,:,:) + x(3,2,:,:).*y(2,2,:,:) + x(3,3,:,:).*y(3,2,:,:)).*xconj(2,2,:,:) + ...
             (x(3,1,:,:).*y(1,3,:,:) + x(3,2,:,:).*y(2,3,:,:) + x(3,3,:,:).*y(3,3,:,:)).*xconj(2,3,:,:);

z(1,3,:,:) = conj(z(3,1,:,:));

z(2,3,:,:) = conj(z(3,2,:,:));

z(3,3,:,:) = xabs2(3,1,:,:).*y(1,1,:,:)  + ...
             xabs2(3,2,:,:).*y(2,2,:,:)  + ...
             xabs2(3,3,:,:).*y(3,3,:,:)  + ...
             2.*real( x(3,2,:,:).*y(2,1,:,:).*xconj(3,1,:,:) ) + ...
             2.*real( x(3,3,:,:).*y(3,1,:,:).*xconj(3,1,:,:) ) + ...
             2.*real( x(3,3,:,:).*y(3,2,:,:).*xconj(3,2,:,:) );

             %(x(3,2,:,:).*y(2,1,:,:) + x(3,3,:,:).*y(3,1,:,:)).*xconj(3,1,:,:) + ...
             %(x(3,1,:,:).*y(1,2,:,:) + x(3,3,:,:).*y(3,2,:,:)).*xconj(3,2,:,:) + ...
             %(x(3,1,:,:).*y(1,3,:,:) + x(3,2,:,:).*y(2,3,:,:)).*xconj(3,3,:,:);
%
%b1 b2 b7    a1 a2' a4'   b1' b3' b5'
%b3 b4 b8    a2 a3  a5'   b2' b4' b6'
%b5 b6 b9    a4 a5  a6    b7' b8' b9'
%
%b1*a1+b2*a2+b7*a4  b1*a2'+b2*a3+b7*a5  b1*a4'+b2*a5'+b7*a6   b1' b3' b5'
%b3*a1+b4*a2+b8*a4  b3*a2'+b4*a3+b8*a5  b3*a4'+b4*a5'+b8*a6   b2' b4' b6'
%b5*a1+b6*a2+b9*a4  b5*a2'+b6*a3+b9*a5  b5*a4'+b6*a5'+b9*a6   b7' b8' b9'
%
%

%b1*a1*b1'+b2*a2*b1'+b1*a2'*b2'+b2*a3*b2' b1*a1*b3'+b2*a2*b3'+b1*a2'*b4'+b2*a3*b4'
%b3*a1*b1'+b4*a2*b1'+b3*a2'*b2'+b4*a3*b2' b3*a1*b3'+b4*a2*b3'+b3*a2'*b4'+b4*a3*b4'
%
%a1*abs(b1)^2 + a2*(b1'*b2) + a2'*(b1*b2') + a3*abs(b2)^2    a1*b1*b3'    + a2*b2*b3'   + a2'*b1*b4'   + a3*b2*b4'
%a1*b1'*b3    + a2*b1'*b4   + a2'*b2'*b3   + a3*b2'*b4       a1*abs(b3)^2 + a2*(b3'*b4) + a2'*(b3*b4') + a3*abs(b4)^2


text:
function appendTextToReadme(comments, dlgFlag);
% function appendTextToReadme(<comments>, <dlgFlag=1>);
%
% Provides an interface for adding user com.ments to Readme file.
% If a string is entered for the 'comments' argument, uses that
% string as the default comments in the dialog. If the dialog flag
% is set to 0 <default 1>, doesn't put up the dialog, and just writes the
% comments to the readme.
%
% djh, 9/4/01
% ras, 6/17/03 (from createReadmeAppendNotes)
% ras, 03/10/06: added input args
if notDefined('comments'),  comments = '';      end
if notDefined('dlgFlag'),   dlgFlag = 1;        end

if dlgFlag==1
    % Dialog box to input comments
    pos = [200 400 450 200];
    topH = figure(...
        'MenuBar', 'none', ...
        'Name', 'Any other comments about this session?', ...
        'NumberTitle','off', ...
        'UserData', '', ...
        'Position', pos ...
        );

    % Create the comments field:
    labelPos = [0.05 0.2 0.9 0.8];
    h = CreateCommentsField(labelPos, topH, 12, comments);

    % install the file-control buttons
    bpos = [0.1 0.05 0.2 0.1];
    uicontrol( ...
        'Style', 'pushbutton', ...
        'String', 'Accept', ...
        'Units','Normalized',...
        'HorizontalAlignment', 'center', ...
        'Callback', 'uiresume', ...
        'FontSize', 14, ...
        'Position', bpos ...
        );

    % Wait until we get a uiresume, then perform an update. Repeat
    % this cycle until the update reports no errors.
    ok = 0;
    while ~ok
      uiwait(topH);
        ok = 1;
    end

    % Get comments
    comments = get(h,'String');
    close(topH);
end

% Append to Readme.txt
[fid, message] = fopen('Readme.txt','a');
if fid == -1
    warndlg(messsage);
    return
end

for line = 1:size(comments, 1)
    wcount = fprintf(fid,'%s\n', comments(line,:));
end

status = fclose(fid);
if status == -1
    warndlg(messsage);
    return
end

return

author: vistalab

repo: vistasoft


text:
function [z,output] = minf_lbfgs(f,g,z0,options)
%MINF_LBFGS Minimize a function by L-BFGS with line search.
%   [z,output] = minf_lbfgs(f,g,z0) starts at z0 and attempts to find a
%   local minimizer of the real-valued function f(z). The input variables z
%   may be a scalar, vector, matrix, tensor or even a (nested) cell array
%   of tensors and its contents may be real or complex.
%
%   If f(x) is a function of real variables x, the function g(x) should
%   compute the partial derivatives of f with respect to the real variables
%   x, i.e. g(xk) := df(xk)/dx. If f(z) is a function of complex variables
%   z, the function g(z) should compute two times the partial derivative
%   of f with respect to conj(z) (treating z as constant), i.e. g(zk) :=
%   2*df(zk)/d(conj(z)) = 2*conj(df(zk)/dz). If g is the empty matrix [],
%   the real gradient or scaled conjugate cogradient is approximated with
%   finite differences. The output of the function g(z) may have the same
%   structure as z (although this is not necessary). The structure output
%   returns additional information:
%
%      output.alpha      - The line search step length in every iteration.
%      output.fevals     - The total number of function/gradient calls.
%      output.fval       - The value of the objective function f in every
%                          iteration.
%      output.info       - The circumstances under which the procedure
%                          terminated:
%                             1: Objective function tolerance reached.
%                             2: Step size tolerance reached.
%                             3: Maximum number of iterations reached.
%      output.infols     - The circumstances under which the line search
%                          terminated in every iteration.
%      output.iterations - The number of iterations.
%      output.relfval    - The difference in objective function value
%                          between every two successive iterates, relative
%                          to its initial value.
%      output.relstep    - The step size relative to the norm of the
%                          current iterate in every iteration.
%
%   minf_lbfgs(f,g,z0,options) may be used to set the following options:
%
%      options.Display = 10      - Displays output information each
%                                  options.Display iterations. Set to 0 to
%                                  disable.
%      options.LineSearch        - The line search used to minimize the
%      = @ls_mt                    objective function in the quasi-Newton
%                                  descent direction.
%      options.LineSearchOptions - The options structure passed to the line
%                                  search routine.
%      options.M                 - The number of L-BFGS updates to store.
%      = min(30,length(z0))
%      options.MaxIter = 500     - The maximum number of iterations.
%      options.TolFun = 1e-6     - The tolerance for output.relfval.
%      options.TolX = 1e-8       - The tolerance for output.relstep.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Unconstrained
%       optimization of real functions in complex variables", SIAM J. Opt.,
%       Vol. 22, No. 3, 2012, pp. 879-898.

% Evaluate the objective function and gradient.
dim = structure(z0);
fval = f(z0);
if ~isa(g,'function_handle') && isempty(g)
    grad = serialize(deriv(f,z0,fval));
else
    grad = serialize(g(z0));
end
z = z0;

% Check the options structure.
isfunc = @(f)isa(f,'function_handle');
if nargin < 4, options = struct; end
if ~isfield(options,'Display'), options.Display = 10; end
if ~isfield(options,'LineSearch') || ~isfunc(options.LineSearch)
    options.LineSearch = @ls_mt;
end
if ~isfield(options,'LineSearchOptions')
    options.LineSearchOptions = struct;
end
if ~isfield(options.LineSearchOptions,'alpha')
    options.LineSearchOptions.alpha = 1;
end
if ~isfield(options.LineSearchOptions,'c2')
    options.LineSearchOptions.c2 = 0.9;
end
if ~isfield(options,'M'), options.M = min(30,length(grad)); end
if ~isfield(options,'MaxIter'), options.MaxIter = 500; end
if ~isfield(options,'TolFun'), options.TolFun = 1e-6; end
if ~isfield(options,'TolX'), options.TolX = 1e-8; end

% Initialize the algorithm.
S = zeros(numel(grad),options.M);
Y = zeros(numel(grad),options.M);
a = zeros(1,options.M);
r = zeros(1,options.M);
m = 0;
midx = 0;

% L-BFGS with line search.
output.alpha = [];
output.fevals = 1;
output.fval = fval;
output.info = false;
output.infols = [];
output.iterations = 0;
output.relfval = [];
output.relstep = [];
while ~output.info

    % Compute the quasi-Newton step pqn = -H*grad.
    pqn = -grad;
    for i = 1:m
        a(i) = r(midx(i))*real(S(:,midx(i))'*pqn);
        pqn = pqn-a(i)*Y(:,midx(i));
    end
    if m > 0
        y1 = Y(:,midx(1));
        y1y1 = y1'*y1;
        gamma = y1y1*r(midx(1));
        pqn = 1/gamma*pqn;
    end
    for i = m:-1:1
        b = r(midx(i))*real(Y(:,midx(i))'*pqn);
        pqn = pqn+(a(i)-b)*S(:,midx(i));
    end

    % Minimize f along z+alpha*pqn.
    state = output; state.grad = grad;
    [alpha,outputls] = options.LineSearch( ...
        f,g,z,deserialize(pqn,dim),state,options.LineSearchOptions);
    output.alpha(:,end+1) = alpha;
    if length(alpha) < 2, alpha(2) = 1; end

    % Update iterate.
    z1 = serialize(z);
    z = alpha(2)*(z1+alpha(1)*pqn);
    if alpha(2) == 1
        s = alpha(1)*pqn;
    else
        s = z-z1;
    end
    z = deserialize(z,dim);

    % Update gradient and Hessian approximation.
    grad1 = grad;
    if isfield(outputls,'grad')
        grad = outputls.grad;
    else
        if ~isa(g,'function_handle') && isempty(g)
            grad = serialize(deriv(f,z,output.fval(end)));
        else
            grad = serialize(g(z));
        end
    end
    y = grad-grad1;
    sy = real(y'*s);
    if sy > 0
        m = min(m+1,options.M);
        midx = [midx(1)+1:-1:1,m:-1:midx(1)];
        S(:,midx(1)) = s;
        Y(:,midx(1)) = y;
        r(:,midx(1)) = 1/sy;
    end

    % Update the output structure.
    if isfield(outputls,'fevals')
        output.fevals = output.fevals+outputls.fevals;
    end
    if isfield(outputls,'fval')
        output.fval(end+1) = outputls.fval;
    else
        output.fval(end+1) = f(z);
    end
    if isfield(outputls,'info')
        output.infols(end+1) = outputls.info;
    end
    output.iterations = output.iterations+1;
    output.relfval(end+1) = ...
        abs(diff(output.fval(end:-1:end-1)))/abs(output.fval(1));
    output.relstep(end+1) = norm(s)/norm(z1);
    if isnan(output.relstep(end)), output.relstep(end) = 0; end
    if output.relfval(end) <= options.TolFun, output.info = 1; end
    if output.relstep(end) <= options.TolX, output.info = 2; end
    if output.iterations >= options.MaxIter, output.info = 3; end

    % Display progress.
    if options.Display > 0 && (output.iterations == 1 || output.info || ...
       mod(output.iterations,options.Display) == 0)
        if output.iterations == 1
            bold = '%s';
            [~,~,~,~,v] = regexp(version('-release'),'([0-9]+)([ab])');
            if usejava('Desktop') && str2double(v{1}{1}) > 2011 || ...
               (str2double(v{1}{1}) == 2011 && strcmpi(v{1}{2},'b'))
                bold = '%s';
            end
        end
        if output.iterations == 1 || ...
           mod(output.iterations,15*options.Display) == 0
            fprintf('\n%7s%s','',sprintf(bold,'fval'));
            fprintf('%13s%s','',sprintf(bold,'relfval'));
            fprintf('%10s%s','',sprintf(bold,'relstep'));
            fprintf('%10s%s','',sprintf(bold,'alpha'));
            fprintf('\n%30s = %4.e %6s = %4.e\n\n', ...
                    'TolFun',options.TolFun,'TolX',options.TolX);
        end
        if output.iterations == 1
            fprintf('%4i: % 14.8e |\n',0,output.fval(1));
        end
        stralpha = [repmat('%10.4e ',1,size(output.alpha,1)) '\n'];
        fprintf(['%4i: % 14.8e | %14.8e | %14.8e | ' stralpha], ...
                output.iterations,output.fval(end), ...
                output.relfval(end),output.relstep(end), ...
                abs(output.alpha(:,end)));
    end

end

% Display termination message.
if options.Display > 0
    ahref = '\n%s\n\n';
    x = round(linspace(0,output.iterations,min(500,output.iterations)));
    if length(bold) > 2
        ahref = sprintf(['\n%%s\n\n'],mat2str(x'), ...
            mat2str([output.fval(x+1)' [nan output.relfval(x(2:end))]' ...
                    [nan output.relstep(x(2:end))]'],3));
    end
    switch output.info
        case 1, fprintf(ahref,'Objective function tolerance reached.');
        case 2, fprintf(ahref,'Step size tolerance reached.');
        case 3, fprintf(ahref,'Maximum number of iterations reached.');
    end
end

end

function [z,offset] = deserialize(z,dim,offset)
    if iscell(dim)
        v = z;
        z = cell(size(dim));
        if nargin < 3, offset = 0; end
        for i = 1:numel(z)
            if iscell(dim{i})
                [z{i},offset] = deserialize(v,dim{i},offset);
            else
                n = prod(dim{i}(:));
                z{i} = reshape(v(offset+(1:n)),dim{i});
                offset = offset+n;
            end
        end
    elseif ~isempty(dim)
        z = reshape(z,dim);
    end
end

function z = serialize(z)
    if iscell(z)
        for i = find(cellfun(@iscell,z(:).'))
            z{i} = serialize(z{i});
        end
        s = cellfun(@numel,z(:)); o = [0; cumsum(s)];
        c = z; z = zeros(o(end),1);
        for i = 1:length(s), z(o(i)+(1:s(i))) = c{i}(:); end
    else
        z = z(:);
    end
end

function dim = structure(z)
    if iscell(z)
        dim = cellfun(@size,z,'UniformOutput',false);
        for i = find(cellfun(@iscell,z(:).'))
            dim{i} = structure(z{i});
        end
    else
        dim = size(z);
        if numel(z) == dim(1), dim = []; end
    end
end

author: andrewssobral

repo: lrslibrary


text:
function varargout = corr_filter(x, y, der_w, varargin)

opts.lambda = nan;
opts = vl_argparse(opts, varargin);

% x is [m1, m2, p, b]
% y is [m1, m2]
% der_w is same size as x

sz = size_min_ndims(x, 4);
n = prod(sz(1:2));

y_f = fft2(y);
x_f = fft2(x);
% k = 1/n sum_i x_i corr x_i + lambda delta
assert(~isnan(opts.lambda), 'lambda must be specified');
k_f = 1/n*sum(conj(x_f).*x_f, 3) + opts.lambda;
% a must satisfy n (k conv a) = y
% The signal a contains a weight per example (shift)
a_f = 1/n*bsxfun(@times, y_f, 1./k_f);

if isempty(der_w)
    % Use same weight a for all channels i.
    % w[i] = a corr x[i]
    w_f = bsxfun(@times, conj(a_f), x_f);
    % w = ifft2(w_f, 'symmetric');
    w = real(ifft2(w_f));
    varargout = {w};
else
    der_w_f = fft2(der_w);
    % a, x -> w
    % w[i] = a corr x[i]
    % dw[i] = da corr x[i] + a corr dx[i]
    % F dw[i] = conj(F da) .* F x[i] + conj(F a) .* F dx[i]
    % <der_w, dw> = sum_i <der_w[i], dw[i]> = sum_i <F der_w[i], F dw[i]>
    %   = sum_i <F der_w[i], conj(F da) .* F x[i] + conj(F a) .* F dx[i]>
    %   = <F da, sum_i conj(F der_w[i]) .* F x[i]> + sum_i <F der_w[i] .* F a, F dx[i]>
    der_a_f = sum(x_f .* conj(der_w_f), 3);
    der_x_f = bsxfun(@times, a_f, der_w_f);
    % k, y -> a
    % k conv a = 1/n y
    % dk conv a + k conv da = 1/n dy
    % dk_f .* a_f + k_f .* da_f = 1/n dy_f
    % <der_a, da> = <der_a_f, da_f>
    %   = <der_a_f, k_f^-1 .* (1/n dy_f - dk_f .* a_f)>
    %   = <1/n der_a_f .* conj(k_f^-1), dy_f> + <-der_a_f .* conj(k_f^-1 .* a_f), dk_f>
    %   = <der_y_f, dy_f> + <der_k_f, dk_f>
    der_y_f = 1/n*sum(der_a_f .* conj(1 ./ k_f), 4); % accumulate gradients over batch
    der_y = real(ifft2(der_y_f));
    der_k_f = -der_a_f .* conj(a_f ./ k_f);
    % x -> k
    % k = 1/n sum_i x_i corr x_i + lambda delta
    % dk = 1/n sum_i {dx[i] corr x[i] + x[i] corr dx[i]}
    % F dk = 1/n sum_i {conj(F dx[i]) .* F x[i] + conj(F x[i]) .* F dx[i]}
    % <der_k, dk> = <der_k, 1/n sum_i {dx[i] corr x[i] + x[i] corr dx[i]}>
    %   = sum_i <F der_k, 1/n conj(F dx[i]) .* F x[i] + conj(F x[i]) .* F dx[i]>
    %   = sum_i <F dx[i], 1/n conj(F der_k) .* F x[i]> + <1/n F der_k .* F x[i], F dx[i]>
    %   = sum_i <F dx[i], 1/n [F der_k + conj(F der_k)] .* F x[i]>
    %   = sum_i <F dx[i], 2/n real(F der_k) .* F x[i]>
    %   = sum_i <F der_x[i], F dx[i]>
    der_x_f = der_x_f + 2/n*bsxfun(@times, real(der_k_f), x_f);
    % der_x = ifft2(der_x_f, 'symmetric');
    der_x = real(ifft2(der_x_f));
    varargout = {der_x, der_y};
end

end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
% -------------------------------------------------------------------------
% Plot individual slopes of regressions, using conf. interval for x range
% -------------------------------------------------------------------------
function igls_plot_slopes(out, X, varargin)

    % data
    % -------------------------------------------------------------------
    N = out.sub; % subjects

    betas = out.beta_indiv([2 1], :);  % format for compatibility with mediation.m
    % slope is first row, intercept is 2nd

    groupbeta = out.beta([2 1]);

    w = ones(N, 1);
    if ~isempty(varargin)
        % we have stats structure with (maybe) weights
        if isfield(out, 'w')
            w = out.w;
            w = w .* N;
        end
    end

    ncols = 1;

    % plot
    % -------------------------------------------------------------------
    fh = create_figure('Slope_Plot');

    subplot(1, ncols, 1);
    linehan = plot_individual_slopes(betas, X, w(:,1));
    set(linehan, 'Color', [.7 .7 .7]);

    linehan_group = plot_individual_slopes(groupbeta, X(:), 3);
    set(linehan_group, 'LineWidth', 4);

    %xlabel(vnames{1}), ylabel(vnames{3});
    title('Slopes');

% %     subplot(1, ncols, 2);
% %     plot_individual_slopes(bbetas, M, w(:,2));
% %     xlabel(vnames{3}), ylabel(vnames{2});
% %     title('b: M->Y controlling X');
% %
% %     subplot(1, ncols, 3);
% %     plot_individual_slopes(cpbetas, X, w(:,3));
% %     xlabel(vnames{1}), ylabel(vnames{2});
% %     title('c'': X->Y controlling M');
% %
% %     subplot(1, ncols, 4);
% %     plot_individual_slopes(cbetas, X, w(:,4));
% %     xlabel(vnames{1}), ylabel(vnames{2});
% %     title('c: X->Y');
end

function linehan = plot_individual_slopes(betas, X, w)

    % sort so that we plot from lowest to highest weights
    [w, sorti] = sort(w);
    betas = betas(:,sorti);
    if iscell(X), X = X(sorti); else X = X(:,sorti); end

    % get colors based on weights
    N = size(betas, 2);
    minwt = .4;     % make sure all lines are visible.
    w = rescale_range(w, [minwt 1]);

    % line widths: median split, top half gets 1, bottom gets 1/2.
    linew = (w(:,1) >= median(w(:,1))) +.5;

    w = 1 - w;  % for colors, 0 is black
    colors = repmat([1 1 1], N, 1) .* repmat(w, 1, 3);
    colors(colors < 0) = 0; % include to remove rounding error

    for i = 1:N
        if iscell(X), x = X{i}; else x = X(:,i); end

        % 95% conf. interval for x
        [nanvec, x] = nanremove(x);
        mx = mean(x);
        s = std(x) * tinv(.975, length(x)-1);
        x = [mx - s mx + s];
        y = betas(2, i) + betas(1, i) * x;
        linehan(i) = plot(x, y, '-', 'Color', colors(i,:), 'LineWidth', linew(i));
    end

    drawnow
end

function rx = rescale_range(x, y)
    % re-scale x to range of y
    m = range(y)./range(x);

    if isinf(m)
        % no range/do not rescale
        rx = x;
    else
        b = y(1) - m * x(1);
        rx = m*x + b;
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [growthRate, minProd, maxProd] = testOptKnockSol(model, targetRxn, deletions)
% Tests an `OptKnock` knockout strain
%
% USAGE:
%
%    [growthRate, minProd, maxProd] = testOptKnockSol(model, targetRxn, deletions)
%
% INPUTS:
%    model:         COBRA model structure
%    targetRxn:     Target reaction (e.g. 'EX_etoh(e)')
%
% OPTIONAL INPUT:
%    deletions:     Set of reaction deletions (e.g. {'PGI', 'TPI'})
%                   (Default = [])
%
% OUTPUTS:
%    growthRate:    Maximim growth rate of the strain
%    minProd:       Minimum production rate at max growth rate
%    maxProd:       Maximum production rate at max growth rate
%
% .. Author - Markus Herrgard 5/23/07

if (nargin < 3)
    deletions = [];
end

tol = 1e-7;

% Number of deletions
nDel = length(deletions);

modelKO = model;
for i = 1:nDel
    modelKO = changeRxnBounds(modelKO,deletions{i},0,'b');
end
% Calculate optimal growth rate
solKO = optimizeCbModel(modelKO);
growthRate = solKO.f;
if (solKO.stat == 1)
    % Max & min production of the metabolite at the optimal growth rate
    grRounded = floor(solKO.f/tol)*tol;
    modelKO = changeRxnBounds(modelKO,modelKO.rxns(modelKO.c==1),grRounded,'l');
    modelKO = changeObjective(modelKO,targetRxn);
    solMax = optimizeCbModel(modelKO,'max');
    solMin = optimizeCbModel(modelKO,'min');
    maxProd = solMax.f;
    minProd = solMin.f;
else
    maxProd = 0;
    minProd = 0;
end

author: opencobra

repo: cobratoolbox


text:
base64 = org.apache.commons.codec.binary.Base64;
tsv_file = 'MsCelebV1-Faces-Cropped.tsv';

root_folder = 'data';
mkdir(root_folder);

tsv_fid = fopen(tsv_file, 'r');
total_lines = 8456240;
image_id = 0;

pdollar_toolbox_path='D:/face project/pdollar-toolbox';
addpath(genpath(pdollar_toolbox_path));

MTCNN_path = 'D:\face project\MTCNN_face_detection_alignment\code\codes\MTCNNv2';
caffe_model_path=[MTCNN_path , '/model'];
addpath(genpath(MTCNN_path));

coord5points = [30.2946, 65.5318, 48.0252, 33.5493, 62.7299; ...
                51.6963, 51.5014, 71.7366, 92.3655, 92.2041];
imgSize = [112, 96];

%caffe.set_mode_cpu();
gpu_id=0;
caffe.set_mode_gpu();
caffe.set_device(gpu_id);
caffe.reset_all();

%three steps's threshold
threshold=[0.6 0.7 0.7]

minsize = 100;

%scale factor
factor=0.709;

%load caffe models
prototxt_dir =strcat(caffe_model_path,'/det1.prototxt');
model_dir = strcat(caffe_model_path,'/det1.caffemodel');
PNet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir = strcat(caffe_model_path,'/det2.prototxt');
model_dir = strcat(caffe_model_path,'/det2.caffemodel');
RNet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir = strcat(caffe_model_path,'/det3.prototxt');
model_dir = strcat(caffe_model_path,'/det3.caffemodel');
ONet=caffe.Net(prototxt_dir,model_dir,'test');
prototxt_dir =  strcat(caffe_model_path,'/det4.prototxt');
model_dir =  strcat(caffe_model_path,'/det4.caffemodel');
LNet=caffe.Net(prototxt_dir,model_dir,'test');
faces=cell(0);

while ~feof(tsv_fid)
    line = fgetl(tsv_fid);%MID,EntityNameString,ImageURL,FaceID,FaceRectangle_Base64Encoded,FaceData_Base64Encoded
    image_id = image_id + 1;
    if image_id < 5000000
        if mod(image_id, 100000)==0
            disp(image_id);
        end;
        continue;
    end;
    C = strsplit(line,'\t');
    folder = fullfile(root_folder,C{1});
    if exist(folder,'dir')==0
        mkdir(folder);
    end;
    filename = fullfile(folder,[C{2} '-FaceId-0.jpg']);
    if exist(filename, 'file')>0
        continue;
    end;
    raw  = base64decode(C{end});
%     img_fid = fopen(filename, 'wb');
%     fwrite(img_fid, raw, 'uint8');
%     fclose(img_fid);
    jImg = javax.imageio.ImageIO.read(java.io.ByteArrayInputStream(raw));
    h = jImg.getHeight;
    w = jImg.getWidth;
    p = typecast(jImg.getData.getDataStorage, 'uint8');
    img = permute(reshape(p, [3 w h]), [3 2 1]);
    img = img(:,:,[3 2 1]);
    if size(img, 3) < 3
       img(:,:,2) = img(:,:,1);
       img(:,:,3) = img(:,:,1);
    end
        disp([num2str(image_id) '/' num2str(total_lines) ' ' filename]);
    [boundingboxes points]=detect_face(img,minsize,PNet,RNet,ONet,LNet,threshold,false,factor);
    if isempty(boundingboxes)
        continue;
    end;
    default_face = 1;
    if size(boundingboxes,1) > 1
        for bb=2:size(boundingboxes,1)
            if abs((boundingboxes(bb,1) + boundingboxes(bb,3))/2 - size(img,2) / 2) + abs((boundingboxes(bb,2) + boundingboxes(bb,4))/2 - size(img,1) / 2) < ...
                    abs((boundingboxes(default_face,1) + boundingboxes(default_face,3))/2 - size(img,2) / 2) + abs((boundingboxes(default_face,2) + boundingboxes(default_face,4))/2 - size(img,1) / 2)
                default_face = bb;
            end;
        end;
    end;
    facial5points = double(reshape(points(:,default_face),[5 2])');
    Tfm =  cp2tform(facial5points', coord5points', 'similarity');
    cropImg = imtransform(img, Tfm, 'XData', [1 imgSize(2)],...
                                  'YData', [1 imgSize(1)], 'Size', imgSize);
    imwrite(cropImg, filename);
	%show detection result
% 	numbox=size(boundingboxes,1);
%     figure(1);
% 	imshow(img)
% 	hold on;
% 	for j=1:numbox
% 		plot(points(1:5,j),points(6:10,j),'g.','MarkerSize',10);
% 		r=rectangle('Position',[boundingboxes(j,1:2) boundingboxes(j,3:4)-boundingboxes(j,1:2)],'Edgecolor','g','LineWidth',3);
%     end
%     hold off;
%     figure(2);
%     imshow(cropImg);
% 	pause
end;
fclose(tsv_fid);

author: happynear

repo: FaceVerification


text:
function [y,deriv] = boost_obj(w,T,weights,logit_prior)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% Weighted binary classifier cross-entropy objective, based on 'boosting'
% proper scoring rule. This rule places more emphasis on extreme scores,
% than the logariothmic scoring rule.
%
%  Differentiable inputs:
%   w: is vector of N detection scores (in log-likelihood-ratio format)
%
%  Fixed parameters:
%   T: is vector of N labels: 1 for target and -1 for non-target.
%   weights: is N-vector of objective function weights, one per trial.
%   logit_prior: is logit(prior), this controls the region of interest

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w)boost_obj(w,T,weights,logit_prior);
    return;
end

if isa(w,'function_handle')
    outer = boost_obj([],T,weights,logit_prior);
    y = compose_mv(outer,w,[]);
    return;
end

w = w(:);
scores = w.';
arg = bsxfun(@plus,scores,logit_prior).*T;
wobj = exp(-arg/2).*weights; % 1*N
y = sum(wobj);

if nargout>1
    deriv = @(dy) deriv_this(dy,wobj(:),T);
end

function [g,hess,linear] = deriv_this(dy,wobj,T)
g0 = -0.5*wobj.*T(:);
g = dy*g0;
linear = false;
hess = @(d) hessianprod(d,dy,g0,wobj);

function [h,Jv] = hessianprod(d,dy,g0,wobj)

h = dy*(0.25*wobj(:).*d(:));

if nargout>1
    Jv = d.'*g0;
end

function test_this()
N = 30;
T = [ones(1,N/3),-ones(1,N/3),zeros(1,N/3)];
scores = randn(1,N);
weights = [rand(1,2*N/3),zeros(1,N/3)];
f = @(w) brier_obj(w,T,weights,-2.23);
f = @(w) boost_obj(w,T,weights,-2.23);
test_MV2DF(f,scores(:));

author: nesl


text:
function plotSpektra(psi,varargin)
% plot the Legendre coefficients of the kernel function

bw = get_option(varargin,'bandwidth',32);
bw = min(bw,length(psi.A)-1);
A = reshape(psi.A(1:bw+1),1,[]);

if check_option(varargin,'logarithmic')
  optiondraw(loglog(0:bw,abs(A)./(2*(0:bw)+1),...
    'marker','o','MarkerSize',5),varargin{:});
else
  optiondraw(semilogx(0:bw,A./(2*(0:bw)+1),...
    'marker','o','MarkerSize',5),varargin{:});
end
set(gcf,'Name',['Legendre coefficients of the kernel ',inputname(1)]);

end


text:
function efv = pne_event_nose(cx, opt)
%PNE_EVENT_NOSE  Event function to detect the limit or nose point
%   EFV = PNE_EVENT_NOSE(CX, OPT)
%
%   PNE_MASTER event function to detect the limit or nose point of the
%   continuation curve, based on the sign of the lambda component of the
%   tangent vector.
%
%   Set OPT.stop_at to 'NOSE' to use this function to trigger termination
%   of the continuation at the limit or nose point.
%
%   Inputs:
%       CX : struct containing info about current point (continuation soln)
%       OPT - PNES_MASTER options struct
%
%   Outputs:
%       EFV : event function value
%
%   See also PNES_MASTER, PNE_REGISTER_EVENTS, PNE_EVENT_TARGET_LAM.

%   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.

%% event function value is dlam, the last element of the
%% normalized tangent vector at the current soln
efv = cx.z(end);

author: MATPOWER


text:
% change in angle to closest point on wall in fly's coordinate system
function [data,units] = compute_dangle2closestroi2(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  % set sign so that negative means going toward 0, positive means going
  % away from 0
  if trx(fly).nframes <= 1,
    data{i} = [];
  else
    data{i} = sign(trx(fly).angle2closestroi2(1:end-1)).*...
      modrange(diff(trx(fly).angle2closestroi2,1,2),-pi,pi)./trx(fly).dt;
  end
end
units = parseunits('rad/s');


text:
function modelRet = modelTest(modelType, numIn, varargin)

% MODELTEST Run some tests on the specified model.
% FORMAT
% DESC runs some tests on the specified model to ensure it is
% correctly implemented.
% ARG modelType : type of model to test. For example, 'linear' or
% 'mlp'.
% RETURN model : the model that was generated for the tests.
%
% SEEALSO : modelCreate
%
% COPYRIGHT : Neil D. Lawrence,  2006

% MLTOOLS

%if exist([modelType 'Test']) == 2
%  feval([modelType 'Test'])
%end
if ~isstruct(modelType)
  if nargin < 2
    numIn = 4;
  end
  numData = 20;
  numOut = 3;

  % Generate some x positions.
  x = randn(numData, numIn);
  x2 = randn(numData, numIn);
  options = modelOptions(modelType);
  model = modelCreate(modelType, numIn, numOut, options);
  model = modelParamInit(model);
  model.X = x;
  model.y = modelOut(model, x2);

  % Set the parameters randomly.
  params = modelExtractParam(model);
  params = randn(size(params))./sqrt(randn(size(params)).^2);
  model = modelExpandParam(model, params);
else
  model = modelType;
end

% Check model parameter gradient.
modelGradientCheck(model, varargin{:});

if exist([model.type 'OutputGrad'])==2;
  epsilon = 1e-6;
  [params, names] = modelExtractParam(model);
  if length(names) == 0
    for i = 1:model.numParams
      names{i} = ['Param ' num2str(i)];
    end
  end
  origParams = params;
  Lplus = zeros(size(model.X, 1), model.numParams, model.outputDim);
  Lminus = zeros(size(model.X, 1), model.numParams, model.outputDim);
  for i = 1:length(params);
    params = origParams;
    params(i) = origParams(i) + epsilon;
    model = modelExpandParam(model, params);
    Lplus(:, i, :) = reshape(modelOut(model, model.X), ...
                             size(model.X, 1), 1, model.outputDim);
    params(i) = origParams(i) - epsilon;
    model = modelExpandParam(model, params);
    Lminus(:, i, :) = reshape(modelOut(model, model.X), ...
                              size(model.X, 1), 1, model.outputDim);
  end
  params = origParams;
  model = modelExpandParam(model, params);
  gLDiff = .5*(Lplus - Lminus)/epsilon;
  g = modelOutputGrad(model, model.X);

  outParamMaxDiff = max(max(max(abs(gLDiff-g))));
  if outParamMaxDiff > 2*epsilon
    l = 0;
    for i = 1:length(names)
      if l < length(names{i})
        l = length(names{i});
      end
    end

    fprintf([char(repmat(32, 1, l)) '\tanalytic   diffs     delta\n']);
    for i = 1:length(names)
      spaceLen = l - length(names{i});
      space = char(repmat(32, 1, spaceLen));
      fprintf([space names{i} ':\t%4.6f\t%4.6f\t%4.6f\n'], ...
              g(i), gLDiff(i), gLDiff(i) - g(i));
    end
  end
  fprintf('Output param max diff: %2.6f.\n', outParamMaxDiff);
else
  fprintf('No grad of output with respect to param implemented.\n');
end

if exist([model.type 'OutputGradX'])==2;
  epsilon = 1e-6;
  X = model.X;
  origX = X;
  Lplus = zeros(size(model.X, 1), size(model.X, 2), model.outputDim);
  Lminus = zeros(size(model.X, 1), size(model.X, 2), model.outputDim);
  for i = 1:size(model.X, 2)
    X = origX;
    X(:, i) = origX(:, i) + epsilon;
    Lplus(:, i, :) = reshape(modelOut(model, X), ...
                             size(model.X, 1), 1, model.outputDim);
    X(:, i) = origX(:, i) - epsilon;
    Lminus(:, i, :) = reshape(modelOut(model, X), ...
                              size(model.X, 1), 1, model.outputDim);
  end
  X = origX;
  gLDiff = .5*(Lplus - Lminus)/epsilon;
  g = modelOutputGradX(model, X);

  outputXMaxDiff = max(max(max(abs(gLDiff-g))));
  if outputXMaxDiff > 2*epsilon
    fprintf('gX\n')
    disp(g)
    fprintf('gXDiff\n')
    disp(gLDiff)
  end
  fprintf('X max diff: %2.6f.\n', outputXMaxDiff);
else
  fprintf('No grad of output with respect to X implemented.\n');
end
if nargout > 0
  modelRet = model;
else
  modelDisplay(model);
end

author: SheffieldML

repo: GPmat


text:
function x = power_mod ( a, n, m )

%*****************************************************************************80
%
%% POWER_MOD computes mod ( A^N, M ).
%
%  Discussion:
%
%    Some programming tricks are used to speed up the computation, and to
%    allow computations in which A**N is much too large to store in a
%    real word.
%
%    First, for efficiency, the power A**N is computed by determining
%    the binary expansion of N, then computing A, A**2, A**4, and so on
%    by repeated squaring, and multiplying only those factors that
%    contribute to A**N.
%
%    Secondly, the intermediate products are immediately "mod'ed", which
%    keeps them small.
%
%    For instance, to compute mod ( A^13, 11 ), we essentially compute
%
%       13 = 1 + 4 + 8
%
%       A^13 = A * A^4 * A^8
%
%       mod ( A^13, 11 ) = mod ( A, 11 ) * mod ( A^4, 11 ) * mod ( A^8, 11 ).
%
%    Fermat's little theorem says that if P is prime, and A is not divisible
%    by P, then ( A^(P-1) - 1 ) is divisible by P.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, the base of the expression to be tested.
%    A should be nonnegative.
%
%    Input, integer N, the power to which the base is raised.
%    N should be nonnegative.
%
%    Input, integer M, the divisor against which the expression is tested.
%    M should be positive.
%
%    Output, integer X, the remainder when A^N is divided by M.
%
  if ( a < 0 )
    x = -1;
    return
  end

  if ( floor ( a ) ~= a )
    x = -1;
    return
  end

  if ( n < 0 )
    x = -1;
    return
  end

  if ( floor ( n ) ~= n )
    x = -1;
    return
  end

  if ( m <= 0 )
    x = -1;
    return
  end

  if ( floor ( m ) ~= m )
    x = -1;
    return
  end
%
%  A contains the successive squares of A.
%
  x = 1;

  while ( 0 < n )

    d = mod ( n, 2 );

    if ( d == 1 )
      x = mod ( x * a, m );
    end

    a = mod ( a * a, m );
    n = floor ( ( n - d ) / 2 );

  end
%
%  Ensure that 0 <= X.
%
  while ( x < 0 )
    x = x + m;
  end

  return
end