text:
%        NOTE        NOTE              NOTE

% NOTE: THIS CODE IS the work of 'Pascal Getreuer ' which converts all the
% possible spaces to various different spaces
%Source : Mathworks

function varargout = colorspace(Conversion,varargin)
%%% Input parsing %%%
if nargin < 2, error('Not enough input arguments.'); end
[SrcSpace,DestSpace] = parse(Conversion);

if nargin == 2
   Image = varargin{1};
elseif nargin >= 3
   Image = cat(3,varargin{:});
else
   error('Invalid number of input arguments.');
end

FlipDims = (size(Image,3) == 1);

if FlipDims, Image = permute(Image,[1,3,2]); end
if ~isa(Image,'double'), Image = double(Image)/255; end
if size(Image,3) ~= 3, error('Invalid input size.'); end

SrcT = gettransform(SrcSpace);
DestT = gettransform(DestSpace);

if ~ischar(SrcT) && ~ischar(DestT)
   % Both source and destination transforms are affine, so they
   % can be composed into one affine operation
   T = [DestT(:,1:3)*SrcT(:,1:3),DestT(:,1:3)*SrcT(:,4)+DestT(:,4)];
   Temp = zeros(size(Image));
   Temp(:,:,1) = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
   Temp(:,:,2) = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
   Temp(:,:,3) = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
   Image = Temp;
elseif ~ischar(DestT)
   Image = rgb(Image,SrcSpace);
   Temp = zeros(size(Image));
   Temp(:,:,1) = DestT(1)*Image(:,:,1) + DestT(4)*Image(:,:,2) + DestT(7)*Image(:,:,3) + DestT(10);
   Temp(:,:,2) = DestT(2)*Image(:,:,1) + DestT(5)*Image(:,:,2) + DestT(8)*Image(:,:,3) + DestT(11);
   Temp(:,:,3) = DestT(3)*Image(:,:,1) + DestT(6)*Image(:,:,2) + DestT(9)*Image(:,:,3) + DestT(12);
   Image = Temp;
else
   Image = feval(DestT,Image,SrcSpace);
end

%%% Output format %%%
if nargout > 1
   varargout = {Image(:,:,1),Image(:,:,2),Image(:,:,3)};
else
   if FlipDims, Image = permute(Image,[1,3,2]); end
   varargout = {Image};
end

return;

function [SrcSpace,DestSpace] = parse(Str)
% Parse conversion argument

if ischar(Str)
   Str = lower(strrep(strrep(Str,'-',''),'=',''));
   k = find(Str == '>');

   if length(k) == 1         % Interpret the form 'src->dest'
      SrcSpace = Str(1:k-1);
      DestSpace = Str(k+1:end);
   else
      k = find(Str == '<');

      if length(k) == 1      % Interpret the form 'dest<-src'
         DestSpace = Str(1:k-1);
         SrcSpace = Str(k+1:end);
      else
         error(['Invalid conversion, ''',Str,'''.']);
      end
   end

   SrcSpace = alias(SrcSpace);
   DestSpace = alias(DestSpace);
else
   SrcSpace = 1;             % No source pre-transform
   DestSpace = Conversion;
   if any(size(Conversion) ~= 3), error('Transformation matrix must be 3x3.'); end
end
return;

function Space = alias(Space)
Space = strrep(strrep(Space,'cie',''),' ','');

if isempty(Space)
   Space = 'rgb';
end

switch Space
case {'ycbcr','ycc'}
   Space = 'ycbcr';
case {'hsv','hsb'}
   Space = 'hsv';
case {'hsl','hsi','hls'}
   Space = 'hsl';
case {'rgb','yuv','yiq','ydbdr','ycbcr','jpegycbcr','xyz','lab','luv','lch'}
   return;
end
return;

function T = gettransform(Space)
% Get a colorspace transform: either a matrix describing an affine transform,
% or a string referring to a conversion subroutine
switch Space
case 'ypbpr'
   T = [0.299,0.587,0.114,0;-0.1687367,-0.331264,0.5,0;0.5,-0.418688,-0.081312,0];
case 'yuv'
   % sRGB to NTSC/PAL YUV
   % Wikipedia: http://en.wikipedia.org/wiki/YUV
   T = [0.299,0.587,0.114,0;-0.147,-0.289,0.436,0;0.615,-0.515,-0.100,0];
case 'ydbdr'
   % sRGB to SECAM YDbDr
   % Wikipedia: http://en.wikipedia.org/wiki/YDbDr
   T = [0.299,0.587,0.114,0;-0.450,-0.883,1.333,0;-1.333,1.116,0.217,0];
case 'yiq'
   % sRGB in [0,1] to NTSC YIQ in [0,1];[-0.595716,0.595716];[-0.522591,0.522591];
   % Wikipedia: http://en.wikipedia.org/wiki/YIQ
   T = [0.299,0.587,0.114,0;0.595716,-0.274453,-0.321263,0;0.211456,-0.522591,0.311135,0];
case 'ycbcr'
   % sRGB (range [0,1]) to ITU-R BRT.601 (CCIR 601) Y'CbCr
   % Wikipedia: http://en.wikipedia.org/wiki/YCbCr
   % Poynton, Equation 3, scaling of R'G'B to Y'PbPr conversion
   T = [65.481,128.553,24.966,16;-37.797,-74.203,112.0,128;112.0,-93.786,-18.214,128];
case 'jpegycbcr'
   % Wikipedia: http://en.wikipedia.org/wiki/YCbCr
   T = [0.299,0.587,0.114,0;-0.168736,-0.331264,0.5,0.5;0.5,-0.418688,-0.081312,0.5]*255;
case {'rgb','xyz','hsv','hsl','lab','luv','lch','cat02lms'}
   T = Space;
otherwise
   error(['Unknown color space, ''',Space,'''.']);
end
return;

function Image = rgb(Image,SrcSpace)
% Convert to sRGB from 'SrcSpace'
switch SrcSpace
case 'rgb'
   return;
case 'hsv'
   % Convert HSV to sRGB
   Image = huetorgb((1 - Image(:,:,2)).*Image(:,:,3),Image(:,:,3),Image(:,:,1));
case 'hsl'
   % Convert HSL to sRGB
   L = Image(:,:,3);
   Delta = Image(:,:,2).*min(L,1-L);
   Image = huetorgb(L-Delta,L+Delta,Image(:,:,1));
case {'xyz','lab','luv','lch','cat02lms'}
   % Convert to CIE XYZ
   Image = xyz(Image,SrcSpace);
   % Convert XYZ to RGB
   T = [3.2406, -1.5372, -0.4986; -0.9689, 1.8758, 0.0415; 0.0557, -0.2040, 1.057];
   R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3);  % R
   G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3);  % G
   B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3);  % B
   % Desaturate and rescale to constrain resulting RGB values to [0,1]
   AddWhite = -min(min(min(R,G),B),0);
   R = R + AddWhite;
   G = G + AddWhite;
   B = B + AddWhite;
   % Apply gamma correction to convert linear RGB to sRGB
   Image(:,:,1) = gammacorrection(R);  % R'
   Image(:,:,2) = gammacorrection(G);  % G'
   Image(:,:,3) = gammacorrection(B);  % B'
otherwise  % Conversion is through an affine transform
   T = gettransform(SrcSpace);
   temp = inv(T(:,1:3));
   T = [temp,-temp*T(:,4)];
   R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
   G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
   B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
   Image(:,:,1) = R;
   Image(:,:,2) = G;
   Image(:,:,3) = B;
end

% Clip to [0,1]
Image = min(max(Image,0),1);
return;

function Image = xyz(Image,SrcSpace)
% Convert to CIE XYZ from 'SrcSpace'
WhitePoint = [0.950456,1,1.088754];

switch SrcSpace
case 'xyz'
   return;
case 'luv'
   % Convert CIE L*uv to XYZ
   WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
   WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
   L = Image(:,:,1);
   Y = (L + 16)/116;
   Y = invf(Y)*WhitePoint(2);
   U = Image(:,:,2)./(13*L + 1e-6*(L==0)) + WhitePointU;
   V = Image(:,:,3)./(13*L + 1e-6*(L==0)) + WhitePointV;
   Image(:,:,1) = -(9*Y.*U)./((U-4).*V - U.*V);                  % X
   Image(:,:,2) = Y;                                             % Y
   Image(:,:,3) = (9*Y - (15*V.*Y) - (V.*Image(:,:,1)))./(3*V);  % Z
case {'lab','lch'}
   Image = lab(Image,SrcSpace);
   % Convert CIE L*ab to XYZ
   fY = (Image(:,:,1) + 16)/116;
   fX = fY + Image(:,:,2)/500;
   fZ = fY - Image(:,:,3)/200;
   Image(:,:,1) = WhitePoint(1)*invf(fX);  % X
   Image(:,:,2) = WhitePoint(2)*invf(fY);  % Y
   Image(:,:,3) = WhitePoint(3)*invf(fZ);  % Z
case 'cat02lms'
    % Convert CAT02 LMS to XYZ
   T = inv([0.7328, 0.4296, -0.1624;-0.7036, 1.6975, 0.0061; 0.0030, 0.0136, 0.9834]);
   L = Image(:,:,1);
   M = Image(:,:,2);
   S = Image(:,:,3);
   Image(:,:,1) = T(1)*L + T(4)*M + T(7)*S;  % X
   Image(:,:,2) = T(2)*L + T(5)*M + T(8)*S;  % Y
   Image(:,:,3) = T(3)*L + T(6)*M + T(9)*S;  % Z
otherwise   % Convert from some gamma-corrected space
   % Convert to sRGB
   Image = rgb(Image,SrcSpace);
   % Undo gamma correction
   R = invgammacorrection(Image(:,:,1));
   G = invgammacorrection(Image(:,:,2));
   B = invgammacorrection(Image(:,:,3));
   % Convert RGB to XYZ
   T = inv([3.2406, -1.5372, -0.4986; -0.9689, 1.8758, 0.0415; 0.0557, -0.2040, 1.057]);
   Image(:,:,1) = T(1)*R + T(4)*G + T(7)*B;  % X
   Image(:,:,2) = T(2)*R + T(5)*G + T(8)*B;  % Y
   Image(:,:,3) = T(3)*R + T(6)*G + T(9)*B;  % Z
end
return;

function Image = hsv(Image,SrcSpace)
% Convert to HSV
Image = rgb(Image,SrcSpace);
V = max(Image,[],3);
S = (V - min(Image,[],3))./(V + (V == 0));
Image(:,:,1) = rgbtohue(Image);
Image(:,:,2) = S;
Image(:,:,3) = V;
return;

function Image = hsl(Image,SrcSpace)
% Convert to HSL
switch SrcSpace
case 'hsv'
   % Convert HSV to HSL
   MaxVal = Image(:,:,3);
   MinVal = (1 - Image(:,:,2)).*MaxVal;
   L = 0.5*(MaxVal + MinVal);
   temp = min(L,1-L);
   Image(:,:,2) = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
   Image(:,:,3) = L;
otherwise
   Image = rgb(Image,SrcSpace);  % Convert to sRGB
   % Convert sRGB to HSL
   MinVal = min(Image,[],3);
   MaxVal = max(Image,[],3);
   L = 0.5*(MaxVal + MinVal);
   temp = min(L,1-L);
   S = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
   Image(:,:,1) = rgbtohue(Image);
   Image(:,:,2) = S;
   Image(:,:,3) = L;
end
return;

function Image = lab(Image,SrcSpace)
% Convert to CIE L*a*b* (CIELAB)
WhitePoint = [0.950456,1,1.088754];

switch SrcSpace
case 'lab'
   return;
case 'lch'
   % Convert CIE L*CH to CIE L*ab
   C = Image(:,:,2);
   Image(:,:,2) = cos(Image(:,:,3)*pi/180).*C;  % a*
   Image(:,:,3) = sin(Image(:,:,3)*pi/180).*C;  % b*
otherwise
   Image = xyz(Image,SrcSpace);  % Convert to XYZ
   % Convert XYZ to CIE L*a*b*
   X = Image(:,:,1)/WhitePoint(1);
   Y = Image(:,:,2)/WhitePoint(2);
   Z = Image(:,:,3)/WhitePoint(3);
   fX = f(X);
   fY = f(Y);
   fZ = f(Z);
   Image(:,:,1) = 116*fY - 16;    % L*
   Image(:,:,2) = 500*(fX - fY);  % a*
   Image(:,:,3) = 200*(fY - fZ);  % b*
end
return;

function Image = luv(Image,SrcSpace)
% Convert to CIE L*u*v* (CIELUV)
WhitePoint = [0.950456,1,1.088754];
WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));

Image = xyz(Image,SrcSpace); % Convert to XYZ
Denom = Image(:,:,1) + 15*Image(:,:,2) + 3*Image(:,:,3);
U = (4*Image(:,:,1))./(Denom + (Denom == 0));
V = (9*Image(:,:,2))./(Denom + (Denom == 0));
Y = Image(:,:,2)/WhitePoint(2);
L = 116*f(Y) - 16;
Image(:,:,1) = L;                        % L*
Image(:,:,2) = 13*L.*(U - WhitePointU);  % u*
Image(:,:,3) = 13*L.*(V - WhitePointV);  % v*
return;

function Image = lch(Image,SrcSpace)
% Convert to CIE L*ch
Image = lab(Image,SrcSpace);  % Convert to CIE L*ab
H = atan2(Image(:,:,3),Image(:,:,2));
H = H*180/pi + 360*(H < 0);
Image(:,:,2) = sqrt(Image(:,:,2).^2 + Image(:,:,3).^2);  % C
Image(:,:,3) = H;                                        % H
return;

function Image = cat02lms(Image,SrcSpace)
% Convert to CAT02 LMS
Image = xyz(Image,SrcSpace);
T = [0.7328, 0.4296, -0.1624;-0.7036, 1.6975, 0.0061; 0.0030, 0.0136, 0.9834];
X = Image(:,:,1);
Y = Image(:,:,2);
Z = Image(:,:,3);
Image(:,:,1) = T(1)*X + T(4)*Y + T(7)*Z;  % L
Image(:,:,2) = T(2)*X + T(5)*Y + T(8)*Z;  % M
Image(:,:,3) = T(3)*X + T(6)*Y + T(9)*Z;  % S
return;

function Image = huetorgb(m0,m2,H)
% Convert HSV or HSL hue to RGB
N = size(H);
H = min(max(H(:),0),360)/60;
m0 = m0(:);
m2 = m2(:);
F = H - round(H/2)*2;
M = [m0, m0 + (m2-m0).*abs(F), m2];
Num = length(m0);
j = [2 1 0;1 2 0;0 2 1;0 1 2;1 0 2;2 0 1;2 1 0]*Num;
k = floor(H) + 1;
Image = reshape([M(j(k,1)+(1:Num).'),M(j(k,2)+(1:Num).'),M(j(k,3)+(1:Num).')],[N,3]);
return;

function H = rgbtohue(Image)
% Convert RGB to HSV or HSL hue
[M,i] = sort(Image,3);
i = i(:,:,3);
Delta = M(:,:,3) - M(:,:,1);
Delta = Delta + (Delta == 0);
R = Image(:,:,1);
G = Image(:,:,2);
B = Image(:,:,3);
H = zeros(size(R));
k = (i == 1);
H(k) = (G(k) - B(k))./Delta(k);
k = (i == 2);
H(k) = 2 + (B(k) - R(k))./Delta(k);
k = (i == 3);
H(k) = 4 + (R(k) - G(k))./Delta(k);
H = 60*H + 360*(H < 0);
H(Delta == 0) = nan;
return;

function Rp = gammacorrection(R)
Rp = zeros(size(R));
i = (R <= 0.0031306684425005883);
Rp(i) = 12.92*R(i);
Rp(~i) = real(1.055*R(~i).^0.416666666666666667 - 0.055);
return;

function R = invgammacorrection(Rp)
R = zeros(size(Rp));
i = (Rp <= 0.0404482362771076);
R(i) = Rp(i)/12.92;
R(~i) = real(((Rp(~i) + 0.055)/1.055).^2.4);
return;

function fY = f(Y)
fY = real(Y.^(1/3));
i = (Y < 0.008856);
fY(i) = Y(i)*(841/108) + (4/29);
return;

function Y = invf(fY)
Y = fY.^3;
i = (Y < 0.008856);
Y(i) = (fY(i) - 4/29)*(108/841);
return;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ segmentedImage, binaryImage ] = ColorSegments( regionMatrix )
%COLORSEGMENTS Summary of this function goes here
%   Detailed explanation goes here

uniqueRegionLabels = unique(regionMatrix);
uniqueRegionLabelCount = numel(uniqueRegionLabels);

labelColors = [];

for ii = 1 : uniqueRegionLabelCount

    rValue = randi(255);

    gValue = randi(255);

    bValue = randi(255);

    labelColors = [labelColors; rValue gValue bValue];
end

frequentLabel = mode(mode(regionMatrix));

[rows, cols] = size(regionMatrix);

segmentedImage = uint8(zeros(rows, cols, 3));

binaryImage = uint8(zeros(rows, cols, 1));

for ii = 1 : rows

    for jj = 1 : cols

        regionLabel = regionMatrix(ii, jj);

        colorIndex = find(uniqueRegionLabels == regionLabel);

        segmentedImage(ii, jj, 1) = labelColors(colorIndex, 1);

        segmentedImage(ii, jj, 2) = labelColors(colorIndex, 2);

        segmentedImage(ii, jj, 3) = labelColors(colorIndex, 3);

        if regionLabel == frequentLabel
            binaryImage(ii, jj) = 0;
        else
            binaryImage(ii, jj) = 255;
        end

    end

end

end


text:
function complainif_notposscalar(var,varname,callfun)

if nargin<3
    callfun = mfilename;
end

if  isempty(var) || ~isscalar(var) || ~isnumeric(var) || var<=0 || ...
        ( isinf(var) )
   error('%s: %s should be a positive scalar.',upper(callfun),varname);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% By Sherif A. Tawfik, Faculty of Engineering, Cairo University
% [x,val,status]=IP1(f,A,b,Aeq,beq,lb,ub,M,e)
% this function solves the following mixed-integer linear programming problem
%   min f*x
%  subject to
%        A*x <=b
%        Aeq * x = beq
%        lb <= x <= ub
%        M is a vector of indeces for the variables that are constrained to be integers
%        e is the integarilty tolerance
% the return variables are :
% x : the solution
% val: value of the objective function at the optimal solution
% status =1 if successful
%        =0 if maximum number of iterations reached in he linprog function
%        =-1 if there is no solution
% Example:
%        maximize 17 x1 + 12 x2
%        subject to
% 	             10 x1 + 7 x2 <=40
%                   x1 +   x2 <= 5
%                   x1, x2 >=0 and are integers
% f=[-17, -12]; %take the negative for maximization problems
% A=[ 10  7; 1 1];
% B=[40; 5];
% lb=[0 0];
% ub=[inf inf];
% M=[1,2];
% e=2^-24;
% [x v s]= IP(f,A,B,[],[],lb,ub,M,e)

function [x,val,status]=IP1(f,A,b,Aeq,beq,lb,ub,M,e)
options = optimset('display','off');
bound=inf; % the initial bound is set to +ve infinity
[x0,val0]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);
[x,val,status,b]=rec(f,A,b,Aeq,beq,lb,ub,x0,val0,M,e,bound); % a recursive function that processes the BB tree

function [xx,val,status,bb]=rec(f,A,b,Aeq,beq,lb,ub,x,v,M,e,bound)
options = optimset('display','off');
% x is an initial solution and v is the corressponding objective function value

% solve the corresponding LP model with the integarily constraints removed
[x0,val0,status0]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);

% if the solution is not feasible or the value of the objective function is
% higher than the current bound return with the input intial solution
if status0<=0 | val0 > bound
    xx=x; val=v; status=status0; bb=bound;
    return;
end

% if the integer-constraint variables turned to be integers within the
% input tolerance return
ind=find( abs(x0(M)-round(x0(M)))>e );
if isempty(ind)
    status=1;
    if val0 < bound    % this solution is better than the current solution hence replace
        x0(M)=round(x0(M));
        xx=x0;
        val=val0;
        bb=val0;
    else
        xx=x;  % return the input solution
        val=v;
        bb=bound;
    end
    return
end

% if we come here this means that the solution of the LP relaxation is
% feasible and gives a less value than the current bound but some of the
% integer-constraint variables are not integers.
% Therefore we pick the first one that is not integer and form two LP problems
% and solve them recursively by calling the same function (branching)

% first LP problem with the added constraint that Xi <= floor(Xi) , i=ind(1)
br_var=M(ind(1));
br_value=x(br_var);
if isempty(A)
    [r c]=size(Aeq);
else
    [r c]=size(A);
end
A1=[A ; zeros(1,c)];
A1(end,br_var)=1;
b1=[b;floor(br_value)];

% second LP problem with the added constraint that Xi >= ceil(Xi) , i=ind(1)
A2=[A ;zeros(1,c)];
A2(end,br_var)=-1;
b2=[b; -ceil(br_value)];

% solve the first LP problem
[x1,val1,status1,bound1]=rec(f,A1,b1,Aeq,beq,lb,ub,x0,val0,M,e,bound);
status=status1;
if status1 >0 & bound1<bound % if the solution was successfull and gives a better bound
   xx=x1;
   val=val1;
   bound=bound1;
   bb=bound1;
else
    xx=x0;
    val=val0;
    bb=bound;
end

% solve the second LP problem
[x2,val2,status2,bound2]=rec(f,A2,b2,Aeq,beq,lb,ub,x0,val0,M,e,bound);

if status2 >0 & bound2<bound % if the solution was successfull and gives a better bound
    status=status2;
    xx=x2;
    val=val2;
    bb=bound2;
end

author: Sable

repo: mcbench-benchmarks


text:
function [ c_est, r_est, c_err, r_err ] = sphere_random_centralize ( ...
  m, n, c, r )

%*****************************************************************************80
%
%% SPHERE_RANDOM_CENTRALIZE estimates the center of a sphere.
%
%  Discussion:
%
%    We generate N points on the M dimensional sphere of center C and radius R.
%    We seek to estimate the values of C and R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of random points to generate.
%
%    Input, real C(M,1), the center.
%
%    Input, real R, the radius.
%
%    Output, real C_EST(M,1), R_EST, the estimated center and radius.
%
%    Output, real C_ERR, R_ERR, the errors in the estimates.
%
  c = c(:);
  x = sphere_surface_sample ( m, n, c, r );

  c_est = sum ( x, 2 ) / n;
  c_err = norm ( c - c_est );

  r_mat = x - repmat ( c_est, 1, n );
  r_mat = r_mat.^2;
  r_mat = sum ( r_mat, 1 );
  r_mat = sqrt ( r_mat );
  r_est = sum ( r_mat ) / n;
  r_err = abs ( r - r_est );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
Here you can see part of the example running log for MV-CNN and RGB CNN.

I use K40 to run these codes.
Please run the following code before running fast_classification_UCF101_spatial.m and fast_classification_UCF101_temporal.m:
nvidia-smi -pm 1
nvidia-smi -ac 3004,875

The compuational time for temporal CNN (MV CNN) should around 0.16-0.17 sec.
The compuational time for spatial CNN (RGB CNN) should around 0.10 sec.

Since 250 optical flow and 25 rgb are processed, the speed should be 925 fps.

author: zbwglory

repo: MV-release

save_path: github-repos/MATLAB/zbwglory-MV-release


text:
classdef LineSpecEnum < matlab.mixin.SetGet
    %LineSpecEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        SolidLine('Solid Line', '-')
        DashedLine('Dashed Line', '--')
        DottedLine('Dotted Line', ':')
        DashedDotLine('Dashed-dot Line', '-.')
    end

    properties
        name char = '';
        linespec char = '';
    end

    methods
        function obj = LineSpecEnum(name, linespec)
            obj.name = name;
            obj.linespec = linespec;
        end
    end

    methods(Static)
        function [listBoxStr, m] = getListboxStr()
            m = enumeration('LineSpecEnum');
            listBoxStr = {m.name};
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('LineSpecEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('LineSpecEnum');
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef DMDicomLibraryInterface < CoreBaseClass
    % DMDicomLibraryInterface Interface for a class that parses Dicom
    % files.
    %
    % This interface allows you to use your own Dicom parsing libraries
    % with DicoMat functions, by wrapping them up in a class that
    % implements this interface.
    %
    %
    %     Licence
    %     -------
    %     Part of DicoMat. https://github.com/tomdoel/dicomat
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the BSD 3-Clause license. Please see the file LICENSE for details.
    %

    methods (Abstract)

        % Tests whether a file is in DICOM format
        isDicom = isdicom(obj, fileName)

        % Reads the metaheader data from a Dicom file
        metaheader = dicominfo(obj, varargin)

        % Reads the image data from a Dicom file
        imageData = dicomread(obj, fileName_or_metaHeader)

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio]...
    = cwmr_var(fid, data, varargins, opts)
% This program starts the CWMR-Var algorithm
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ....
%           = cwmr_var(fid, data, varargins, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%          = cwmr_var(fid, data, {2, 0.5, 0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
phi =varargins{1};     % Confidence parameter
epsilon =varargins{2};     %
tc = varargins{3};      % transaction cost fee rate

% Run the CWMR-Var algorithm
[cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = cwmr_var_run(fid, data, phi,  epsilon, tc, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function initializeLengthFilterPanel( horizontalBar, leftBottom, panelSize, hdls )
%INITIALIZELENGTHFILTERPANEL Initialize length filter panel with 2 edit boxes for min and max length of eddy tracks

lengthFilterPanel = uipanel(horizontalBar, ...
    'Units', 'normalize', ...
    'Position', [leftBottom panelSize]);

uicontrol(lengthFilterPanel, ...
    'String', 'Min length', ...
    'Style', 'text',...
    'FontSize', 11, ...
    'Units', 'normalized', ...
    'Position', [0 .5 .5 .5]);

uicontrol(lengthFilterPanel, ...
    'String', 'Max length', ...
    'Style', 'text',...
    'FontSize', 11, ...
    'Units', 'normalized', ...
    'Position', [.5 .5 .5 .5]);

uicontrol(lengthFilterPanel, ...
    'Style', 'edit', ...
    'Units', 'normalized',...
    'Position', [0 0 .5 .5], ...
    'String', 'any', ...
    'Callback', {@minLengthCallback, hdls});

uicontrol(lengthFilterPanel, ...
    'Style', 'edit', ...
    'Units', 'normalized',...
    'Position', [0.5 0 .5 .5], ...
    'String', 'any', ...
    'Callback', {@maxLengthCallback, hdls});

hdls.add_attribute_filter('lifetimes');

end

function minLengthCallback(hObject, ~, hdls)
% Update track min length and generate data again

old_pointer = start_busy_pointer(hdls);

minLength = get(hObject, 'String');
if isempty(minLength) || strcmp(minLength, 'any')
    hdls.set_attribute_filter_min_value('lifetimes', -inf);
else
    minLength = str2double(minLength);
    if(isnan(minLength)) % Not valid string
        errordlg('Min length must be an integer, "any" or empty', 'Min length input error');
        return;
    end
    hdls.set_attribute_filter_min_value('lifetimes', minLength);
end
hdls.generate_data();

end_busy_pointer(hdls, old_pointer);

end

function maxLengthCallback(hObject, ~, hdls)
% Update track max length and generate data again

old_pointer = start_busy_pointer(hdls);

maxLength = get(hObject, 'String');
if(isempty(maxLength) || strcmp(maxLength, 'any'))
    hdls.set_attribute_filter_max_value('lifetimes', inf);
else
    maxLength = str2double(maxLength);
    if(isnan(maxLength)) % Not valid string
        errordlg('Max length must be an integer, "any" or empty', 'Max length input error');
        return;
    end
    hdls.set_attribute_filter_max_value('lifetimes', maxLength);
end
hdls.generate_data();

end_busy_pointer(hdls, old_pointer);

end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
clear all; close all;
I=imread('rice.png');
J=I>120;
[width, height]=size(I);
for i=1:width
    for j=1:height
        if (I(i, j)>130)
            K(i, j)=1;
        else
            K(i, j)=0;
        end
    end
end
figure;
subplot(121);  imshow(J);
subplot(122);  imshow(K);

author: HuangCongQing


text:
function rotError = get_rot_error(R_est,R_gt,d)

nrNodes = size(R_est,1) / d;
rotError = zeros(nrNodes, 1);
for i=1:nrNodes
    Rdiff = R_gt(blkIndices(i,d),:)' * R_est(blkIndices(i,d),:);
    rotError(i, 1) = abs( atan2(Rdiff(2, 1), Rdiff(1, 1)) );
end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function rasterSegs = structUnion(rasterSegs1, rasterSegs2, scanNum, planC)
%"structUnion"
%   Find the union between two structures given their entire set of
%   rasterSegments.  Performs this intersect one slice at a time to save
%   memory. Returns rasterSegments representing the new structure.
%
%   By JRA 10/1/03
%
%   rasterSegs1    : rasterSegments of first structure
%   rasterSegs2    : rasterSegments of second structure
%   scanNum        : single scanNumber that both raster segs are defined on.
%   planC          : CERR planC
%
%   rasterSegs     : rasterSegments of union of structures 1 & 2.
%
%Usage:
%   rasterSegs = structUnion(rasterSegs1, rasterSegs2, scanNum, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
rasterSegs = [];

%sort input rasterSegments by CTSliceValue
rasterSegs1 = sortrows(rasterSegs1, 6);
rasterSegs2 = sortrows(rasterSegs2, 6);

%get list of CTSlices to iterate over.
slices1 = unique(rasterSegs1(:,6));
slices2 = unique(rasterSegs2(:,6));

%for union, need to worry about slices where either structure has segments.
slicesToCalculate = union(slices1, slices2);

%for each slice we are calculating on, create a mask for each structure and
%intersect them. Then convert from that mask to raster segments.
for i=1:length(slicesToCalculate)
    sliceNum = slicesToCalculate(i);
    rasterIndices = find(rasterSegs1(:,6) == sliceNum);
    mask1 = rasterToMask(rasterSegs1(rasterIndices,:), scanNum, planC);
    rasterIndices = find(rasterSegs2(:,6) == sliceNum);
    mask2 = rasterToMask(rasterSegs2(rasterIndices,:), scanNum, planC);
    unionMask = mask1 | mask2;
    rasterSegs = [rasterSegs;maskToRaster(unionMask, sliceNum, scanNum, planC)];
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ output_args ] = gmm_gradient_descent(X,dlog_pi,dlog_mu,dlog_Sigmas,K,D)
%GMM_GRADIENT_DESCENT
%
%   Maximise the log-likelihood of the data via gradient descent
%
%    input ----------------------------------------------------------------
%
%       o X : (D x N), set of Data
%
%       o dlog_pi : function handle, returns derivatives of the priors
%
%       o dlog_mu : function handle, returns derivaties of the Mu's
%
%       o dlog_Sigma : function handle, returns derivatives of the sCovariance.
%
%

end

author: epfl-lasa


text:
function vw=inplaneMotionCompSeries(vw,scan,slice,baseFrame)
%
% vw=inplaneMotionCompSeries(vw,scan,slice,baseFrame)
%
% Loads tSeries, estimates motion between each frame and the baseFrame,
% and warps the image to compensate for the motion.  Holds the warped
% tSeries in vw.tSeries.
%
% Copies the original tSeries.dat file to a origTSeries.dat.
% Saves the new (warped) tSeries in tSeries.dat.
%
% To revert back to the original tSeries, use revertMotionComp.
%
%   scan: default current scan
%   slice: default current slice
%   baseFrame: default 1
%
% djh, 4/16/99

if ~exist('scan','var')
   scan = getCurScan(vw);
end
if ~exist('slice','var')
   slice = viewGet(vw, 'Current Slice');
end
if ~exist('baseFrame','var')
   baseFrame=1;
end

% Load tSeries
[~,nii] = loadtSeries(vw,scan,slice);

dims = niftiGet(nii,'Dim');
data = single(niftiGet(nii,'Data'));

%Now, let us take the tSeries data and transform it into the same
%format as previously saved

nFrames = dims(4);
voxPerSlice = prod(dims(1:2));

tSeries = squeeze(data(:,:,slice,:)); % rows x cols x time
tSeries = reshape(tSeries, [voxPerSlice nFrames])'; % time x voxels

nFrames = numFrames(vw,scan);
dims = sliceDims(vw,scan);
baseIm = reshape(tSeries(baseFrame,:),dims);

% Compute motion estimates and warped tSeries
wtSeries = zeros(size(tSeries));
waitStr = ['Computing motion estimates. scan:' num2str(scan) ' slice:' num2str(slice)];
waitHandle = mrvWaitbar(0,waitStr);
for frame = 1:nFrames
  mrvWaitbar(frame/nFrames)
  im = reshape(tSeries(frame,:),dims);
  M = estMotionIter2(baseIm,im,2,eye(3),1);
  warpedIm = warpAffine2(im,M);
  wtSeries(frame,:) = warpedIm(:)';
end
close(waitHandle)

% Return warped tSeries in vw.tSeries and save it to the tSeries file

size = viewGet(vw,'Functional Slice Dim');
newwtSeries = reshape(wtSeries,size);

%Now let's update the nifti
data(:,:,slice,:) = newwtSeries;
nii = niftiSet(nii,'Data',data);
dim = size(niftiGet(nii,'Data'));
nii = niftiSet(nii,'Dim',dim);

savetSeries(wtSeries,vw,scan,slice,nii);

return

%%% Debug
mrLoadRet
INPLANE{1} = inplaneMotionComp(INPLANE{1});
INPLANE{1} = revertMotionComp(INPLANE{1});

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function imout=padimage(im, sizeimout)

sizeim = size(im);

sizedifhalf = floor((sizeimout-sizeim)/2);

imout = padarray(im, sizedifhalf);

sizedif = sizeimout-size(imout);

if sizedif(1)>0
    imout = [imout; zeros(1,size(imout,2))];
end
if sizedif(2)>0
    imout = [imout, zeros(size(imout,1),1)];
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [spath,mpath,stdpath,fwd,stdfwd,xpath,epath,r1,r2]=schwartzsmithsim(start,param,ftemp,T,rmatrix,sim)
kappa=param(1);
mue=param(2);
sigmax=param(3);
sigmae=param(4);
lambdae=param(5);
lambdax=param(6);
pxe=param(7);
x0=param(8);
e0=param(9);

iter=start;
for k=1:length(T)
    loc=mod(iter,12);

    if loc==0
        f(k)=ftemp(12);
    else
        f(k)=ftemp(loc);
    end
    iter=iter+1;
end

for i=1:length(T)
    explns(i,1)=exp(-kappa*T(i))*x0+e0-((1-exp(-kappa*T(i)))*lambdax/kappa)+(mue-lambdae)*T(i);
    varlns(i,1)=((1-exp(-2*kappa*T(i)))*sigmax^2/(2*kappa))+sigmae^2*T(i)+2*(1-exp(-kappa*T(i)))*pxe*sigmax*sigmae/kappa;
    fwd(i,1)=f(i)*exp(explns(i,1)+.5*varlns(i,1));
    varlnfwd(i,1)= T(i)*(sigmae)^2 + (1-exp(-2*kappa*T(i)))*(sigmax^2)/(2*kappa) + 2*(1-exp(-kappa*T(i)))*pxe*sigmax*sigmae/kappa;
    stdfwd(i,1)=f(i)*sqrt((exp(varlnfwd(i,1))-1)*exp(2*explns(i,1)+varlnfwd(i,1)));
    stdspot(i,1)=f(i)*sqrt((exp(varlnfwd(i,1))-1)*exp(2*explns(i,1)+varlnfwd(i,1)));
    ve(i,1)=sigmae^2*T(i);
    vx(i,1)=(1-exp(-2*kappa*T(i)))*(.5*sigmax^2)/kappa;
end

m=length(T);

r1(:,:)=rmatrix;
r2(:,:)=pxe.*r1(:,:)+sqrt(1-pxe^2).*randn(sim,m);

for i=1:sim

x(1)=x0-lambdax*T(1)-kappa*x0*T(1)+sigmax*sqrt(T(1))*r1(i,1);
e(1)=e0+(mue-lambdae)*T(1)+sigmae*sqrt(T(1))*r2(i,1);
s(1)=f(1)*exp(e(1)+x(1));

for t=1:m-1
    x(t+1)=x(t)-lambdax*(T(t+1)-T(t))-kappa*x(t)*(T(t+1)-T(t))+sigmax*sqrt((T(t+1)-T(t)))*r1(i,t+1);
    e(t+1)=e(t)+(mue-lambdae)*(T(t+1)-T(t))+sigmae*sqrt((T(t+1)-T(t)))*r2(i,t+1);
    s(t+1)=f(t+1)*exp(e(t+1)+x(t+1));
end
epath(i,:)=e(1:end);
xpath(i,:)=x(1:end);
spath(i,:)=s(1:end);
lspath(i,:)=log(s(1:end));
c(i)=corr(r1(i,:)',r2(i,:)');
c2(i)=corr(xpath(i,:)',epath(i,:)');
end
evar = var(epath)';
xvar = var(xpath)';

%E(S(T))=E(F(T,T))
mpath=mean(spath)';
%V(ln(S(T)))=V[ln(F(T,T)]
varlnpath=var(lspath)';
%E(ln(S(T)))=E[ln(F(T,T)]
mlnpath=mean(lspath)';
stdpath=sqrt((exp(varlnpath)-1).*exp(2*mlnpath+varlnpath));


text:
function [c, ceq] = lotkaConstraintFCN(u,x,Ts,N,p)
%% Constraint function of nonlinear MPC for Lotka-Volterra system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon
%
% Output:
%   c:      inequality constraints applied across prediction horizon
%   ceq:    equality constraints (empty)
%

%% Nonlinear MPC design parameters
% Predator population size always positive: >0, min pop size
zMin = [10];

%% Inequality constraints calculation
c = zeros(N,1);
% Apply N population size constraints across prediction horizon, from time
% k+1 to k+N
xk = x;
uk = u(1);
for ct=1:N
    % obtain new cart position at next prediction step
    xk1 = rk4u(@sparseGalerkinControl_Discrete,xk,uk,Ts,1,[],p);

    % -z + zMin < 0 % lower bound
    c(ct) = -xk1(2)+zMin; %c(2*ct-1)
    % z - zMax < 0 % upper bound
    %c(2*ct) = xk1(1)-zMax;
    % update plant state and input for next step
    xk = xk1;
    if ct<N
        uk = u(ct+1);
    end
end
%% No equality constraints
ceq = [];

author: eurika-kaiser

repo: SINDY-MPC


text:
function [predicted,memberships, numhits] = fknn(data, labels, test, ...
	    testlabels, k_values, info, fuzzy)
% FKNN Fuzzy k-nearest neighbor classification algorithm.
% 	Y = FKNN(DATA, LABELS, TEST, TESTLABELS, K, INFO) Runs fuzzy
% 	k-nearest neighbors on the given data. DATA is a N-by-D data matrix
% 	consisting of N patterns each of which is D-dimensional. LABELS is a
% 	N-vector containing the class labels (1,2,...,C) for each pattern.
%	TEST is a M-by-D matrix consisting of M test patterns. TESTLABELS
%	is an optional M-vector for the true class labels of the given test
%	data. If you don't have true labels, just give an empty matrix for this
%	TESTLABELS.
%	K is the number of nearest neighbors to look at.
%	The algorithm will print an information line at every INFO test
%	patterns, if INFO>0. If INFO is zero, nothing will be printed.
%	Y is a M-vector containing the predicted class labels for the given test
%	patterns.
%
% 	[Y,MEMS,HITS] = FKNN(DATA, LABELS, TEST, TESTLABELS, K, INFO) returns
% 	the fuzzy class-memberships values in MEMS, for each test pattern. It is
% 	a M-by-C matrix, C being the number of classes.
%	HITS is the number of correctly predicted test patterns. Note that HITS
%	is meaningless if TESTLABELS is not provided.
%
% 	[Y,MEMS,HITS] = FKNN(DATA, LABELS, TEST, TESTLABELS, K, INFO, FUZZY) If
% 	you don't want to do "fuzzy" k-nn, then give FUZZY as 'false'.
%
%	This m-file is capable of testing several k-values simultaneously. If
%	you pass a vector of k-values, rather than a single scalar, in K, then
%	each output variable is populated accordingly. So, if you give K as
%	[5 10 15], then Y becomes M-by-3, MEMS M-by-C-by-3 and HITS  3-by-1.
%
%	References:
%	J. M. Keller, M. R. Gray, and J. A. Givens, Jr.,
%	"A Fuzzy K-Nearest Neighbor Algorithm",
%	IEEE Transactions on Systems, Man, and Cybernetics,
%	Vol. 15, No. 4, pp. 580-585.
%
% TODO: Generalize this m-file to Lp norm
%
% Emre Akbas [eakbas2 @ uiuc.edu]  Dec 2006.
%

if nargin<7
    fuzzy = true;
end

num_train = size(data,1);
num_test  = size(test,1);

% scaling factor for fuzzy weights. see [1] for details
m = 2;

% convert class labels to unary membership vectors (of 1s and 0s)
max_class = max(labels);
temp = zeros(length(labels),max_class);
for i=1:num_train
    temp(i,:) = [zeros(1, labels(i)-1) 1 zeros(1,max_class - labels(i))];
end
labels = temp;
clear temp;

% allocate space for storing predicted labels
predicted = zeros(num_test, length(k_values));

% allocate space for 'numhits'. This will only be used if 'testlabels' is
% provided
numhits = zeros(length(k_values),1);

% will the memberships be stored? if yes, allocate space
store_memberships = false;
if nargout > 1,
    store_memberships=true;
    memberships = zeros(num_test, max_class, length(k_values));
end

%% BEGIN kNN
% for each test point, do:
t0=clock;
tstart = t0;
for i=1:num_test
    distances = (repmat(test(i,:), num_train,1) - data).^2;
    % for efficiency, no need to take sqrt since it is a non-decreasing function
    distances = sum(distances,2)';

    % sort the distances
    [junk, indeces] = sort(distances);

    for k=1:length(k_values)
	neighbor_index = indeces(1:k_values(k));
	weight = ones(1,length(neighbor_index));
	if fuzzy,
 	    % originally, this weight calculation should be:
 	    % weight = distances(neighbor_index).^(-2/(m-1));
 	    % but since we didn't take sqrt above and the inverse 2th power
 	    % the weights are:
 	    % weight = sqrt(distances(neighbor_index)).^(-2/(m-1));
	    % which is equaliavent to:
 	    weight = distances(neighbor_index).^(-1/(m-1));

 	    % set the Inf (infite) weights, if there are any, to  1.
 	    if max(isinf(weight))
 		warning(['Some of the weights are Inf for sample: ' ...
			num2str(i) '. These weights are set to 1.']);
 		weight(isinf(weight))=1;
 	    end
	end
	test_out = weight*labels(neighbor_index,:)/(sum(weight));

	if store_memberships, memberships(i,:,k) = test_out; end;

	% find predicted class (the one with the max. fuzzy vote)
	[junk, index_of_max] = max(test_out');

	predicted(i,k) = index_of_max;

	% compute current hit rate, if test labels are given
	if ~isempty(testlabels) && predicted(i,k)==testlabels(i)
	    numhits(k) = numhits(k)+1;
	end
    end

    % print info
    if mod(i,info)==0
	elapsed = etime(clock, t0);
	fprintf(1,['%dth sample done.  Elapsed (from previous info): %.2f' ...
	    ' sn.  Estimated left: %.2f sn.\n\tHit rate(s) so far:   '], ...
	    i, elapsed, etime(clock, tstart)*((num_test-i)/i) );
	for k=1:length(k_values)
	    fprintf(1,'%3d: %.3f\t',k_values(k), 100*numhits(k)/i);
	end
	fprintf(1,'\n');

	t0=clock; % start timer again
    end
end


text:
function [c] = power(a,b,varargin)
%C=A.^B
%   [C]=POWER(A,B) Computes A.^B for a QTT-Tucker A and natural B.
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
if isa(a,'qtt_tucker') && numel(b) == 1
    c=a;
    b=b-1;
    while(b>0)
      c=c.*a;
      b=b-1;
    end

else
    error('Incorrect usage for TT-power');
end

author: oseledets

repo: TT-Toolbox


text:
function value = r8_big ( )

%*****************************************************************************80
%
%% R8_BIG returns a "big" real number.
%
%  Discussion:
%
%    The value returned by this function is NOT required to be the
%    maximum representable R8.
%    We simply want a "very large" but non-infinite number.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, a huge number.
%
  value = 1.0E+30;

  return
end


text:
% Build and initialize the computational graph for regression based speech
% enhancement/dereverberation
%
function [layer, para] = Build_EnhanceNet_Gaussian(Data_tr, para, stage)
para.output = 'tmp';

layer = genNetworkDereverb_Gaussian(para.topology, stage);     % generate the network graph
para.preprocessing{1} = {};                     % optional preprocessing for each data stream
para.preprocessing{2} = {};
para.cost_func.layer_idx = length(layer);       % specify which layers are cost function layers
para.cost_func.layer_weight = [1];              % set the weights of each cost function layer
para.IO.inputFeature = [1 1];
para.IO.isTensor = [1 1];
para = ParseOptions2(para);

% generating the scaling factor for the input, as we will need to use a
% small constant in the logarithm. We need to make sure that the power of
% speech are larger than this constant most of the time.
scale = 1e4;        % we hard code the scale to be a constant so that all network will use the same number
scale = scale/2^16; % note that we are using int16 to store waveform samples, so need to scale down
if para.topology.useWav
    layer = InitWavScaleLayer(layer, scale);
end

if stage>1 && isfield(para.topology, 'initModel')
    modelfiles = findFiles(['nnet/' para.topology.initModel], 'mat');
    modelfiles = sort_nnet_by_dev(modelfiles);
    fprintf('Use initial model %s\n', modelfiles{1});
    dnnInit = load(modelfiles{1});

    % find the shared, mean, and var subnets
    [sharedLayers, meanLayers, varLayers] = findSubnets(layer);
    if stage == 3.3     % initialize with regression model, use the initial model only for mean prediction
        [sharedLayersInit, meanLayersInit] = findSubnetsRegression(dnnInit.layer);
    elseif stage == 3.2     % initialize with regression model. share the LSTM between mean and variance prediction
        [sharedLayersInit, meanLayersInit] = findSubnetsRegression(dnnInit.layer);
    else            % initialize with stage 1 or 2 model
        [sharedLayersInit, meanLayersInit, varLayersInit] = findSubnets(dnnInit.layer);
    end
    layer = CopyNetWeightsByLayerIdx(dnnInit.layer, [sharedLayersInit meanLayersInit], layer, [sharedLayers meanLayers]);
    if floor(stage) ==3 && stage~=3.1 && stage~=3.2 && stage~=3.3   % for stage 2, we don't initilize var subnet
        layer = CopyNetWeightsByLayerIdx(dnnInit.layer, varLayersInit, layer, varLayers);
    end
else
    if strcmpi(para.topology.RegressionNetType, 'DNN')    % if use DNN, splice the frames
        idx = ReturnLayerIdxByName(layer, 'splice');
    else
        idx = ReturnLayerIdxByName(layer, 'delta'); % if use LSTM, use dynamic features
    end
    fft_net_length = idx(1);
    fft_net = layer(1:fft_net_length);
    paraTmp = para;
    paraTmp.out_layer_idx = fft_net_length;
    paraTmp.IO = RemoveIOStream(paraTmp.IO, [2]);
    paraTmp.IO.nStream = 1;
    paraTmp.IO.inputFeature = paraTmp.IO.DynamicDistortion.inputFeature(1);
    if isfield(paraTmp.IO, 'fileReader'); paraTmp.IO = rmfield(paraTmp.IO, 'fileReader'); end
    paraTmp.IO.fileReader(1)= paraTmp.IO.DynamicDistortion.fileReader(1);
    fprintf('Generate global MVN weights for mask subnet - %s\n', datestr(now));
    [layer{fft_net_length+1}.W, layer{fft_net_length+1}.b] = computeGlobalCMVN(Data_tr(1), 100, paraTmp, fft_net);
    VerifyPreprocessingTree(layer(1:fft_net_length+1), Data_tr(1), paraTmp, 100);
end
end

function [sharedLayers, meanLayers, varLayers] = findSubnets(layer)

delta_idx = ReturnLayerIdxByName(layer, 'delta');
input_affine_idx = delta_idx(1)+1;
tmp = length(layer) + layer{end}.prev;
mean_idx = tmp(1);
var_idx = tmp(2);
var_start = mean_idx+1;
sharedOut_idx = layer{var_start}.prev + var_start;
sharedLayers = (input_affine_idx:sharedOut_idx);
meanLayers = (sharedOut_idx+1:var_start-1);
varLayers = (var_start:var_idx);

end

function [sharedLayers, meanLayers] = findSubnetsRegression(layer)

delta_idx = ReturnLayerIdxByName(layer, 'delta');
input_affine_idx = delta_idx(1)+1;
tmp = length(layer) + layer{end}.prev;
mean_idx = delta_idx(2);

sharedOut_idx = mean_idx-2;
sharedLayers = (input_affine_idx:sharedOut_idx);
meanLayers = (sharedOut_idx+1:mean_idx);

end

author: singaxiong

repo: SignalGraph


text:
function [dat] = read_plexon_ddt(filename, begsample, endsample)

% READ_PLEXON_DDT reads header or data from a Plexon *.ddt file,
% which is a Plexon continuous data file optimized for continuous
% (streaming) recording where every channel is continuously recorded
% without gaps and the recording includes any dead time between spikes.
%
% Use as
%   [hdr] = read_plexon_ddt(filename)
%   [dat] = read_plexon_ddt(filename, begsample, endsample)
%
% samples start counting at 1
% returned values are in mV

% Copyright (C) 2005-2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: read_plexon_ddt.m 945 2010-04-21 17:41:20Z roboos $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version 100: Samples are assumed to be 12 bits.  All channels have the same NIDAQ gain, and preamp
% gain is not saved.
%     int     Version;         // =100
%     int     DataOffset;      // Offset into the file where the data starts
%     double  Freq;            // Digitization frequency
%     int     NChannels;       // Number of channels
%     int     Year;            // Time/date when the data was acquired
%     int     Month;
%     int     Day;
%     int     Hour;
%     int     Minute;
%     int     Second;
%     int     Gain;            // NIDAQ gain (the same gain for all channels
%     char    Comment[128];    // User-supplied comment
%     unsigned char Padding[256];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dat = [];
fid = fopen(filename, 'rb', 'ieee-le');

Version    = fread(fid, 1, 'int'    );
dat.Version    = Version;

% this is common for all subformats
dat.DataOffset = fread(fid, 1, 'int'    );
dat.Freq       = fread(fid, 1, 'double' );
dat.NChannels  = fread(fid, 1, 'int'    );
dat.Year       = fread(fid, 1, 'int'    );
dat.Month      = fread(fid, 1, 'int'    );
dat.Day        = fread(fid, 1, 'int'    );
dat.Hour       = fread(fid, 1, 'int'    );
dat.Minute     = fread(fid, 1, 'int'    );
dat.Second     = fread(fid, 1, 'int'    );
dat.Gain       = fread(fid, 1, 'int'    );
dat.Comment        = char(fread(fid, 128, 'char')');

if Version==100
  % the remainder of the header contains bytes for padding
  dat.Padding    = fread(fid, 256, 'char' );
elseif Version==101
  % Version 101: A field was added to indicate bits-per-sample (12 or 16).  All channels have the same NIDAQ
  % gain, and preamp gain is not saved.
  dat.BitsPerSample = fread(fid,   1, 'char');
  dat.Padding       = fread(fid, 255, 'char');
elseif Version==102
  % Version 102: A byte array of per-channel NIDAQ gains was added.  The Gain field now indicates the
  % preamp gain,   a value of 1 usually indicating that no preamp gain was specified in the application that
  % recorded the DDT.
  dat.BitsPerSample = fread(fid,   1, 'char');
  dat.ChannelGain   = fread(fid,  64, 'char');
  dat.Padding       = fread(fid, 191, 'char');
elseif Version==103
  % Version 103: A field was added to specify ADC maximal input voltage. The value is integer number in
  % millivolts.
  dat.BitsPerSample  = fread(fid,   1, 'char');
  dat.ChannelGain    = fread(fid,  64, 'char');
  dat.MaxMagnitudeMV = fread(fid,   1, 'short');
  dat.Padding        = fread(fid, 189, 'char');
else
  error('unsupported version of ddt file');
end

% determine the number of samples by looking at the length of the datafile
offset = ftell(fid);
fseek(fid, 0, 'eof');
dat.NSamples = (ftell(fid) - offset) / (2 * dat.NChannels);

if nargin==1
  % do not read any data, only return the header
  fclose(fid);
  return
end

% A DDT file contains a file header followed by an array of A/D samples stored as 16-bit integers regardless
% of whether values are 12 or 16-bits. The samples are multiplexed.

% read the desired samples of the data
offset = dat.DataOffset + (begsample-1)*dat.NChannels;
nsampl = endsample - begsample + 1;
fseek(fid, offset, 'bof');
dat.data = fread(fid, [dat.NChannels nsampl], 'int16');

fclose(fid);

% apply the calibration to obtain the signal in physical units
if Version==100
  PreAmpGain = 1000;
  dat.data = dat.data * 5000 / (2048 * dat.Gain * PreAmpGain);
elseif Version==101
  PreAmpGain = 1000;
  dat.data = dat.data * 5000 / (0.5 * (2^dat.BitsPerSample) * dat.Gain * PreAmpGain);
elseif Version==102
  for i=1:dat.NChannels
    dat.data(i,:) = dat.data(i,:) * 5000 ./ (0.5 * (2^dat.BitsPerSample) * dat.ChannelGain(i));
  end
elseif Version==103
  % I am not sure whether the calibration like this is correct, since the
  % Plexon documentation does not explicitely specify how to do it for the
  % 104 file format. I presume that it is identical to the 103 format.
  for i=1:dat.NChannels
    dat.data(i,:) = dat.data(i,:) * 5000 ./ (0.5 * (2^dat.BitsPerSample) * dat.ChannelGain(i));
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function d = p09_fd ( p )

%*****************************************************************************80
%
%% P20_FD is a signed distance function for problem 9.
%
%  Discussion:
%
%    ICAM is represented by 21 regions: 20 rectangles and one polygon.
%
%    Region 12 is not a room, and should be excluded.
%    Region 13 is a staircase to the second floor, and might be excluded.
%
%    The building is contained in the rectangle
%
%      -7 <= X <= 34
%       0 <= Y <= 98
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P, one or more points.
%
%    Output, real D, the signed distance of each point to the boundary of the region.
%
  v = [ ...
     7.25,  0.25; ... % Lower left corner of R1.
    14.75,  0.25; ...
    14.75,  7.75; ...
    13.00,  7.75; ...
    13.00,  8.25; ...
    25.00,  8.25; ...
    25.00,  7.75; ...
    15.25,  7.75; ...
    15.25,  0.25; ... % Lower left corner of R2.
    28.75,  0.25; ...
    28.75,  7.75; ...
    28.75,  8.25; ...
    29.25,  8.25; ...
    29.25,  7.75; ...
    29.25,  0.25; ... % Lower left corner of R3.
    40.75,  0.25; ...
    40.75,  7.75; ...
    33.00,  7.75; ...
    33.00,  8.25; ...
    33.25,  8.25; ... % Lower left corner of R5.
    40.75,  8.25; ...
    40.75, 17.75; ...
    33.25, 17.75; ...
    33.25, 12.00; ...
    32.75, 12.00; ...
    32.75, 21.75; ...
    20.75, 21.75; ...
    20.75, 22.25; ...
    21.25, 22.25; ... % Lower left corner of R8.
    32.75, 22.25; ...
    32.75, 25.25; ...
    35.75, 25.25; ...
    35.75, 34.75; ...
    32.75, 34.75; ...
    32.75, 37.75; ...
    21.25, 37.75; ...
    21.25, 26.00; ...
    20.75, 26.00; ...
    20.75, 38.25; ...
    25.75, 38.25; ...
    26.25, 38.25; ... % Lower left corner of R11.
    32.75, 38.25; ...
    32.75, 45.75; ...
    26.25, 45.75; ...
    26.25, 41.75; ...
    25.75, 41.75; ...
    25.75, 46.25; ...
    26.25, 46.25; ... % Lower left corner of R15.
    32.75, 46.25; ...
    32.75, 53.75; ...
    26.25, 53.75; ...
    26.25, 49.25; ...
    25.75, 49.25; ...
    25.75, 53.75; ...
    23.75, 53.75; ...
    23.75, 62.75; ...
    23.75, 97.75; ...
    16.25, 97.25; ...
    16.25, 88.00; ...
    15.75, 88.00; ...
    15.75, 97.25; ...
     0.25, 97.25; ...
     0.25, 75.25; ... % Lower left corner of R21
    15.75, 75.25; ...
    15.75, 85.00; ...
    16.25, 85.00; ...
    16.25, 74.75; ...
     0.25, 74.75; ...
     0.25, 63.25; ... % Lower left corner of R19
    16.25, 63.25; ...
    16.25, 62.75; ...
    16.25, 58.25; ... % Lower left corner of R18
    19.25, 58.25; ...
    19.25, 54.25; ...
    20.25, 54.25; ...
    20.25, 53.75; ...
    20.25, 50.75; ... % Lower left corner of R16.
    23.25, 50.75; ...
    23.25, 41.75; ...
    17.25, 41.75; ...
    16.75, 41.75; ...
    16.75, 53.75; ...
     7.25, 53.75; ...
     7.25, 38.25; ...
    17.25, 38.25; ...
    17.25, 37.75; ...
     7.25, 37.75; ...
     7.25, 22.25; ...
    16.75, 22.25; ...
    16.75, 34.25; ...
    17.25, 34.25; ...
    17.25, 21.75; ...
     7.25, 21.75; ...
     7.25,  8.25; ... % Lower left corner of R4.
     9.00,  8.25; ...
     9.00,  7.75; ...
     7.25,  7.75; ...
     7.25,  0.25 ];

  d = dpoly ( p, v );

  return
end


text:
function plots(d,xp,yp,hyper,methods)

% PLOTS  plotting for spider
%
% plots(d,xp,yp,[hyper],[methods])
%
% Inputs:
%
%  d          -- results in group of data objects
%  xp,yp      -- x-axis and y-axis variables, can be hyperparameters
%                of algorithms, or the type of loss, e.g xp='C', y='class_loss'
%  hyper      -- string of hyperparameters
%  methods    -- cell array of strings of names methods of interest
%
% Hyperparameters:
%  scatter=0      -- set to 1 for a scatter plot, otherwise tries to draw lines
%  log_scale=0    -- if set to 1 tries to guess when to use a log scale
%  logs=[0 0]     -- set to 1 for x-axis (1st component) or y-axis (2nd)
%                    to be on a log scale
%  no_error_bar=0 -- set to 1 to not draw error bars when given std. dev.
%  ax=[]          -- if nonempty calls axis(ax) to set axis dimensions
%
%  colors=['bkrgymcrgbdkym']          -- default colors to cycle through
%  dots= {'-*','--o','-.',':'};       -- default line styles
%  scats={'*','o','.','+','x','d','s'};- default scatter plot styles
%
% Examples:
%  r=train(cv(param(svm(kernel('rbf')),'rbf',2.^[-5:5])),gen(toy('l=60')));
%  plots(get_mean(r),'rbf','class_loss','log_scale')
%  plots(r,'fold','class_loss')
%
%  [tr a]=train( chain({ param(toy,'l',[10:10:100])  group({svm knn}) }) )
%  r=test(a); plots(r,'l','class_loss')

if nargin<3 yp='class_loss'; end;
if nargin<4 hyper=[]; end;
if nargin<5 methods=[]; end;

d=group(group2cell(d));
if (~isempty(findstr(xp,'loss')) | ~isempty(findstr(yp,'loss')) )
  if isempty(findstr(d.child{1}.name,'loss'))
    disp('[assuming class_loss]');
    d=loss(d,'class_loss');
  end
end

%% ------------------ basic defaults -------------------------
clf
colors=['bkrgymcrgbdkym'];
dots={'-*','--o','-.',':'};
scats={'*','o','.','+','x','d','s'};
assume=1;                      % make assumptions on methods and losses
scatter=0;                     % scatter plot
no_error_bars=0;
log_scale=0; logs=[];
xname=xp; yname=yp;
leg_pos=0; leg_txt=[];
ax=[];
%% ------------------- intelligent defaults -------------------
if strcmp(xp,'fold') | strcmp(yp,'fold') assume=0; scatter=1; end;

% labeling and titles
if strcmp(xname,'feat') xname='features'; end;
if strcmp(yname,'feat') yname='features'; end;
xname(xname=='_')=' '; yname(yname=='_')=' ';
%% ------------------------------------------------------------
evalobject=0;
eval_hyper;             %% evalute hyper params
%% ------------------------------------------------------------

%if nargin<2 calc_loss=1; loss_type='default'; else calc_loss=1; end;
%if isempty(loss_type) calc_loss=0; end;

d=d.child; dx=d; dy=d;

[m o t u a ind]=get_methods(d);
if isempty(methods) methods=m; end;
m=[]; for i=1:length(methods)
  methods{i}(methods{i}=='_')='-';
  m=[m ' ' methods{i}]; end;
disp(['[assuming methods:' m ']']);
if isempty(leg_txt)
  leg_txt=methods;
  leg_err_txt=[];
  for i=1:length(methods)
    leg_err_txt{(i-1)*2+1}=[methods{i} ' error bars'];
    leg_err_txt{(i-1)*2+2}=[methods{i} ];
  end
end;

%if assume & ~isempty(strmatch('cv',u)) & strmatch('fold',t{1})
% dx=get_mean(dx); dy=get_mean(dy);
%end

%if length(dx)~=length(dy)  %% problem: grouping occured in one, not the other
%  if length(dx)<length(dy) dy=dx; else dx=dy; end;
%end

m=length(methods); % number of methods
r=[]; index=[];
for i=[1:m]
  x=[];y=[]; xe=[]; ye=[];
  for k=[1:length(dy)]
   if a(k,ind(i))      %% if method i is used in data k this is TRUE

     x1=get_data_value(dx{k},xp);
     y1=get_data_value(dy{k},yp);

     if ~isempty(x1) & ~isempty(y1)
       %% error bars
       xe1=0;ye1=0;
       if size(x1,2)>1 xe1=x1(1,2); x1=x1(1,1); end;
       if size(y1,2)>1 ye1=y1(1,2); y1=y1(1,1); end;
       xe=[xe xe1]; ye=[ye ye1];
       x=[x x1]; y=[y y1];
     end
   end
  end
  if i==1
    hold off;
    if log_scale & isempty(logs)   %% log scale stuff
       if max(x)>10 logs(1)=1; end;
       if max(y)>10 logs(2)=1; end;
    end
    if length(logs)<2 logs(2)=0; end;
  end;
  if logs(1) x=log(x); end; if logs(2) y=log(y); end;

  c=colors(mod(i-1,length(colors))+1);
  if ~scatter
   d=dots{mod(i-1,length(dots))+1};
  else
   d=scats{mod(i-1,length(scats))+1};
  end
  if sum(ye)==0 | no_error_bars
    plot(x,y,[c d ]); hold on;
  else
    errorbar2(x,y,ye,-ye,[c d ]); hold on;
    leg_txt=leg_err_txt;
  end
end

if logs(1) xname=['log ' xname]; end;
if logs(2) yname=['log ' yname]; end;
xlabel(xname); ylabel(yname);

legend(char(leg_txt),leg_pos);
if ~isempty(ax) axis(ax); end;

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

function val=get_data_value(d,v)

vv=length(v);
if  vv>4 & ~isempty(strcmp('_loss',v(vv-5:vv)))
  val=d.Y;
else
  v=[v '=']; n=d.name;
  f=max(findstr(n,v))+length(v);
  if isempty(f) val=[]; return; end;
  p=[];
  while 1
    if f>length(n) break; end;
    if n(f)~='=' & n(f)~=' '
      p=[p n(f)];
    else
      break;
    end
    f=f+1;
  end
  val=str2num(p);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function F = in_fread_plexon(sFile, SamplesBounds, iChannels, precision)
% IN_FREAD_PLEXON Read a block of recordings from a Plexon file
%
% USAGE:  F = in_fread_plexon(sFile, SamplesBounds=[], iChannels=[])

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

% ===== INSTALL PLEXON SDK =====
[isInstalled, errMsg] = bst_plugin('Install', 'plexon');
if ~isInstalled
    error(errMsg);
end

% ===== PARSE INPUTS =====
if (nargin < 4) || isempty(precision)
    precision = 'double';
elseif ~ismember(precision, {'single', 'double'})
    error('Unsupported precision.');
end
if (nargin < 3) || isempty(iChannels)
    iChannels = 1:sFile.header.ChannelCount;
end
if (nargin < 2) || isempty(SamplesBounds)
    % Read entire recording
    SamplesBounds = round((sFile.prop.times - sFile.header.FirstTimeStamp) .* sFile.prop.sfreq) + 1;
else
    % Readjust the samples call based on the starting time value
    SamplesBounds = SamplesBounds - round(sFile.header.FirstTimeStamp* sFile.prop.sfreq) + 1;
end

% Read the PLX file and assign it to the Brainstorm format
iSelectedChannels = sFile.header.chan_headers;
nChannels = length(iSelectedChannels);
nSamples  = diff(SamplesBounds) + 1;

% Initialize Brainstorm output
F = zeros(nChannels, nSamples, precision);

for iChannel = 1:nChannels
    % plx_ad_span_v returns values in mV
    [adfreq, n, data] = plx_ad_span_v(sFile.filename, iSelectedChannels(iChannel)-1, SamplesBounds(1), SamplesBounds(2));
    F(iChannel,:) = data./1000; % Convert to V
end
end

author: brainstorm-tools


text:
function [y] = temporal_agg_p(z,op1,sc)
% PURPOSE: Temporal aggregation of a time series preserving its dimension
% ------------------------------------------------------------
% SYNTAX: y = temporal_agg_p(z,op1,sc);
% ------------------------------------------------------------
% OUTPUT: y: nx1 temporally aggregated series, missing=0
% ------------------------------------------------------------
% INPUT:  z: nx1 ---> vector of high frequency data
%         op1: type of temporal aggregation
%         op1=1 ---> sum (flow)
%         op1=2 ---> average (index)
%         op1=3 ---> last element (stock) ---> interpolation
%         op1=4 ---> first element (stock) ---> interpolation
%         sc: number of high frequency data points
%            for each low frequency data points
%         Note: n = sc x N
% ------------------------------------------------------------
% LIBRARY: aggreg
% ------------------------------------------------------------

% written by:
%  Enrique M. Quilis
%  Macroeconomic Research Department
%  Ministry of Economy and Competitiveness
%  <enrique.quilis@mineco.es>
% Version 1.0 [October 2010]

[n,m] = size(z);

% ------------------------------------------------------------
% Computes the number of low frequency points.
% Low frequency periods should be complete

N = fix(n/sc);
C = aggreg_p(op1,n,sc);

% -----------------------------------------------------------
% Expanding the aggregation matrix to perform
% extrapolation if needed.

if (n > sc * N)
   pred = n - sc*N;           % Number of required extrapolations
   C=[C zeros(N,pred)];
else
   pred = 0;
end

y = C*z;

author: Sable


text:
function [sampled, label] = seg_label_provider(options,  segpath, sample_size, flip)

% first get the mask
try,

	cnn_input_size = options.cnn_input_size;
	bmask = imread(segpath);

	if(flip)
		bmask = flipdim(bmask,2);
	end

	bmask = imresize(bmask,...
	 [cnn_input_size, cnn_input_size], 'nearest');

	bmask = reshape(bmask, 1, cnn_input_size*cnn_input_size);
	tmask = find(bmask < 255);

	if((length(tmask)>0) && (length(tmask) < sample_size))
		tmask = repmat(tmask, 1, ceil(sample_size/length(tmask)));
	end

	% then sample from each
	sampled = zeros(3, sample_size, 'single');
	label = zeros(1, sample_size, 'single');
	tmask = randsample(tmask,sample_size);
	[y,x] = ind2sub([cnn_input_size, cnn_input_size],tmask);
	label(:) = bmask(tmask);

	%
	sampled(2,:) = y - 1;
	sampled(3,:) = x - 1;
catch,
	keyboard;
end

end

author: aayushbansal

repo: PixelNet


text:
function gaps = findRootNPmets(model, findNCmets)
% Finds the root no production (and no consumption)
% metabolites in a model, used by `gapFind`
%
% USAGE:
%
%    gaps = findRootNPmets(model, findNCmets)
%
% INPUT:
%    model          a COBRA model
%
% OPTIONAL INPUT:
%    findNCmets:    find no consumption mets as well as no production (default = false)
%
% OUTPUT:
%    gaps:          all root no production metabolites
%
% .. Author: - Jeff Orth 7/15/09

if nargin < 2
    findNCmets = false;
end

isRootNPmet = zeros(length(model.mets),1);

for i = 1:length(model.mets)
    row = find(model.S(i,:)); %which rxns this met participates in
    rowR = ismember(row,find(model.lb < 0)); %reversible rxns
    if any(model.S(i,row) > 0) %if met is produced by any reaction
        %don't do anything
    elseif any(rowR) %if met is in any reverible rxns
        %don't do anything
    else
        isRootNPmet(i) = 1;
    end
end

if findNCmets

    isRootNCmet = zeros(length(model.mets),1);

    for i = 1:length(model.mets)
        row = find(model.S(i,:)); %which rxns this met participates in
        rowR = ismember(row,find(model.lb < 0)); %reversible rxns
        if any(model.S(i,row) < 0) %if met is consumed by any reaction
            %don't do anything
        elseif any(rowR) %if met is in any reverible rxns
            %don't do anything
        else
            isRootNCmet(i) = 1;
        end
    end
end

if findNCmets
    gaps = model.mets((isRootNPmet+isRootNCmet)>=1);
else
    gaps = model.mets(isRootNPmet==1);
end

author: opencobra

repo: cobratoolbox


text:
function gX = ggXgaussianKernGradX(ggKern, gaussianKern, X, X2, covGrad)

% GGXGAUSSIANKERNGRADX Compute gradient between the GG and GAUSSIAN
% kernels wrt the input locations
% FORMAT
% DESC computes the gradient between the GG and GAUSSIAN kernels with
% respect to the input positions where both the row positions and column
% positions are provided separately.
% RETURN g : the returned gradients.
% ARG kern : kernel structure for which gradients are being computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
%
% SEEALSO : gaussianKernParamInit, kernGradX, gaussianKernDiagGradX
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009, 2010

% KERN

if nargin < 5,
    covGrad = X2;
    X2 = X;
else
    U = X;
    X = X2;
    X2 = U;
end

[K, Kbase, Pqrinv, Prinv, P] = ggXgaussianKernCompute(ggKern, ....
    gaussianKern, X2, X);

if ggKern.isArd
    PX = X*diag(P);
    PX2 = X2*diag(P);
else
    PX = P*X;
    PX2 = P*X2;
end

temp = covGrad'.*K;

gX2 = zeros(size(X));
for i=1:size(X,2),
    mPX2 = PX2(:,i);
    MPX2 = mPX2(:,ones(1,size(X,1)));
    mPX = PX(:,i)';
    MPX = mPX(ones(size(X2,1),1),:);
    gX2(:,i) = sum(temp.*(MPX2-MPX),1)';
end

gX = gX2;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function logp = logNormalizeVec(logp)
%function logp = logNormalizeRows(logp)
logp = log(normalize_vec(exp(logp)));
end

author: jacobeisenstein

repo: SAGE


text:
%  CORRESPONDENCE INFORMATION
%  This code is written by Gaofeng MENG
%
%  Gaofeng MENG:
%  National Laboratory of Pattern Recognition,
%  Institute of Automation, Academy of Sciences, Beijing 100190
%  Comments and bug reports are welcome.  Email to gfmeng@nlpr.ia.ac.cn
%
%  WORK SETTING:
%  This code has been compiled and tested by using MATLAB R2009a
%
%  For more detials, please see our paper:
%  Gaofeng MENG, Ying WANG, Jiangyong DUAN, Shiming XIANG, Chunhong PAN.
%  Efficient Image Dehazing with Boundary Constraint and Contextual Regularization,
%  ICCV, Sydney, Australia, pp.617-624, 3-6 Dec., 2013.
%
%  Last Modified: Feb. 14, 2014, By Gaofeng MENG
%

function A = DH_bccr_Airlight(HazeImg, method, wsz)
% estimating the global airlight
%
if strcmp(method, 'manual')
    h = figure, imshow(HazeImg, []);
    title('manual airlight estimation: left click to pick a most hazy pixel. ')
    [x, y] = ginput(1);
    A = HazeImg(round(y), round(x), :);
    A = double(A) -1;
    A = min(A, 255);
    close(h);
elseif strcmp(method, 'he')
    A = airlight_he(HazeImg, wsz);
elseif strcmp(method, 'our')
    A = airlight_our(HazeImg, wsz);
else
    error('parameter error.');
end

%%
function A = airlight_our(HazeImg, wsz)
% estimating A channel by channel separately
for k = 1 : 3
    minImg = ordfilt2(double(HazeImg(:, :, k)), 1, ones(wsz), 'symmetric');
    A(k) = max(minImg(:));
end

%%
function A = airlight_he(HazeImg, wsz)
% estimating A using He's method
hsv = rgb2hsv(HazeImg);
GrayImg = hsv(:, :, 3);
[nRows, nCols, bt] = size(HazeImg);

% computing dark channel
DarkImg = min(double(HazeImg), [], 3);
DarkImg = ordfilt2(DarkImg, 1, ones(wsz), 'symmetric');
%
topDark = sort(DarkImg(:), 'descend');
idx = round(0.001 * length(topDark));
val = topDark(idx);
id_set = find(DarkImg >= val);  % the top 0.1% brightest pixels in the dark channel
BrightPxls = GrayImg(id_set);
iBright = find(BrightPxls >= max(BrightPxls));
id = id_set(iBright); id = id(1);
row = mod(id, nRows);
col = floor(id / nRows) + 1;

% A is a vector
A = HazeImg(row, col, :);
A = double(A(:));

author: AomanHao

repo: Matlab-Image-Dehaze-Enhance


text:
function [vectorTo] = niftiCreateStringInplane(vectorFrom, sliceDimension)
%
%Given a vector string and the slice dimension, return the correct string
%orientation in the Vistalab Inplane format, which does not touch the slice
%dimension

%We are looking to return something that will keep the slice dimension the
%same, but make it be the 3rd dimension, then if, L-R are there, make it
%the second dimension, otherwise, leave the other two dimensions as they are

vectorTo = vectorFrom; %We will operate purely on vectorTo

if (sliceDimension ~= 3)
    %We need to move the slice dimension to 3
    tmpChar = vectorTo(sliceDimension);
    vectorTo(sliceDimension) = vectorTo(3);
    vectorTo(3) = tmpChar;
end %if

%We now have a vector that has the 3rd dimension be the slice dimension
%Now, let's move LR dimension to 1st position

LRdim = niftiFindDimOfString(vectorTo,'R'); %This will find either R or L

if (LRdim ~= 3 && LRdim ~= 2) %i.e. it is not the slice dimension nor is it already the second dimension
    tmpChar = vectorTo(2);
    vectorTo(2) = vectorTo(LRdim);
    vectorTo(LRdim) = tmpChar;
end %if

tmpChar = vectorTo(3); %Get the slice dimension and reapply it after, to preserve it
vectorTo = niftiCreateComplementString(vectorTo);
vectorTo(3) = tmpChar;

return


text:
function varargout = mtimes(varargin)
%*   Pointwise multiplication for DISKFUN objects.
%   c*F or F*c multiplies a DISKFUN F by a scalar c.
%
%   F*G computes the integral of F(s,t)G(l,s) over s, and this is the
%   continuous analogue of matrix-matrix multiplication.
%
% See also DISKFUN/TIMES.

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

[varargout{1:nargout}] = mtimes@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fx2 = p16_fx2 ( x )

%*****************************************************************************80
%
%% P16_FX2 evaluates the second derivative of the function for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX2, the second derivative of the function at X.
%
  e = 0.8;

  fx2 = e * pi * pi * sin ( pi * x / 180.0  ) / 180.0 / 180.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = minus(f, g)
%-	Subtraction of two CLASSICFUN objects.
%   F - G subtracts G from F, where F and G may be CLASSICFUN objects or scalars.
%
%   If F and G are both CLASSICFUN objects, they are assumed to have the same domain.
%   The method gives no warning if their domains don't agree, but the output of
%   the method will be meaningless.
%
% See also PLUS, UMINUS.

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

% To avoid duplication, we note that (f - g) = f + (-g):
f = plus(f, uminus(g));

end


text:
for i = 1 : length(missed)
    if ~isempty(missed{i})
        I = imread(fullfile(param.dataset_root, val2db.list{i}));
        hold on
        r = missed{i};
        r=[r(:,[1,2]) r(:,3)-r(:,1) r(:,4)-r(:,2)];
        I = insertShape(I, 'rectangle', r);
        imshow(I)
        waitforbuttonpress;
    end
end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function [f,G,H] = self_collision_barrier(V,E,tol)
  % SELF_COLLISION_BARRIER Compute the a barrier function (and its derivatives)
  % for point-edge collisions of a given *strictly feasible** line-complex in
  % 2D. See "Bijective Parameterization with Free Boundaries" [Smith and
  % Schaefer 2015] or for a similar barrier "Incremental Potential Contact:
  % Intersection- and Inversion-free, Large-Deformation Dynamics" [Li et al.
  % 2020].
  %
  % [f,G,H] = self_collision_barrier(V,E,tol)
  %
  % Inputs:
  %   V  #V by 2 list of input vertex positions
  %   E  #E by 2 list of edge indices into rows of V
  %   tol  distance at which barrier term becomes positive {1e-3}
  % Outputs:
  %   f  scalar total objective value
  %   G  #V*2 gradient
  %   H  #V*2 by #V*2 sparse Hessian matrix
  %
  % Note: this function uses a special symbolic library trick which generates
  % files self_collision_barrier_cap_sym.m and self_collision_barrier_line_sym.m
  % if they don't already exist. If your change the symbolic-math part of this
  % file, you must delete these automatically generated files.
  %

  function [sqrD,T] = point_segment_squared_distance(P,A,B)
    PA = P-A;
    BA = B-A;
    T = min(max(sum(PA.*BA,2)./sum(BA.^2,2),0),1);
    % vector to closest point
    PC = PA-BA.*T;
    sqrD = sum(PC.^2,2);
  end

  function [f,G,H] = self_collision_barrier_cap(V,E,IJ)
    % vertex i colliding with edge jk

    n = size(V,1);
    % [x x x y y y]
    IJ = [IJ  n+IJ];
    X = V(IJ);

    %% Sanity check
    %p = X(:,1:2:end);
    %c = X(:,2:2:end);
    %pc = p-c;
    %d = sqrt(sum(pc.^2,2));
    %if ~isempty(d) && max(d) > tol && nargout > 1
    %  error
    %end

    % Writing the symbol
    path = mfilename('fullpath');
    aux = [path '_cap_sym.m'];
    % should also check date...
    if ~exist(aux,'file')
      % From here, only touch X
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      sX = sym('X',[1 4]);
      stol = sym('tol',[1 1]);
      sp = sX(1:2:end);
      sc = sX(2:2:end);
      spc = sp-sc;
      sd = sqrt(sum(spc.^2,2));
      sf = barrier(sd,stol);

      hess = @(sf,sX) cell2sym(arrayfun(@(g) gradient(g,sX),gradient(sf,sX),'UniformOutput',false));

      aux_handle = ...
        matlabFunction(sf,gradient(sf,sX),hess(sf,sX),'vars',{sX,stol},'File',aux);
    else
      aux_name = [mfilename('func') '_cap_sym'];
      aux_handle = str2func(aux_name);
    end

    faux=[];gaux = [];Haux = [];
    switch nargout
    case {0,1}
      [faux] = aux_handle(X,tol);
    case 2
      [faux,gaux] = aux_handle(X,tol);
    case 3
      [faux,gaux,Haux] = aux_handle(X,tol);
    end

    % unnecessary indirection?
    f_fun = @(X) faux;
    dfdX_fun = @(X) gaux;
    d2fdX2_fun = @(X) Haux;
    f = sum(f_fun(X));
    if nargout<=1
      return;
    end
    dfdX = dfdX_fun(X);
    G = full_sparse(IJ,ones(size(IJ)),reshape(dfdX,size(IJ)),2*n,1);
    if nargout<=2
      return;
    end
    d2fdX2 = double(d2fdX2_fun(X));

    d2fdX2 = reshape(d2fdX2,[],4*4);
    if psd_project
      d2fdX2 = psd_project_rows(d2fdX2);
    end

    HI = repmat(IJ,[1 1 size(IJ,2)]);
    HJ = permute(repmat(IJ,[1 1 size(IJ,2)]),[1 3 2]);
    H = fast_sparse(HI(:),HJ(:),d2fdX2(:),2*n,2*n);
  end

  function [f,G,H] = self_collision_barrier_line(V,E,IJK)

    n = size(V,1);
    % [x x x y y y]
    IJK = [IJK n+[IJK]];
    X = V(IJK);

    %% Sanity check
    %p = X(:,1:3:end);
    %a = X(:,2:3:end);
    %b = X(:,3:3:end);
    %pa = p-a;
    %ba = b-a;
    %t = sum(pa.*ba,2)./sum(ba.^2,2);
    %pc = pa - ba.*t;
    %d = sqrt(sum(pc.^2,2));
    %if ~isempty(d) && max(d) > tol && nargout > 1
    %  error
    %end

    % Writing the symbol
    path = mfilename('fullpath');
    aux = [path '_line_sym.m'];
    % should also check date...
    if ~exist(aux,'file')
      % From here, only touch X
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      sX = sym('X',[1 6]);
      stol = sym('tol',[1 1]);
      sp = sX(1:3:end);
      sa = sX(2:3:end);
      sb = sX(3:3:end);
      spa = sp-sa;
      sba = sb-sa;
      st = sum(spa.*sba,2)./sum(sba.^2,2);
      spc = spa - sba.*st;
      sd = sqrt(sum(spc.^2,2));
      sf = barrier(sd,stol);

      hess = @(sf,sX) cell2sym(arrayfun(@(g) gradient(g,sX),gradient(sf,sX),'UniformOutput',false));
      aux_handle = ...
        matlabFunction(sf,gradient(sf,sX),hess(sf,sX),'vars',{sX,stol},'File',aux);
    else
      aux_name = [mfilename('func') '_line_sym'];
      aux_handle = str2func(aux_name);
    end

    faux=[];gaux = [];Haux = [];
    switch nargout
    case {0,1}
      [faux] = aux_handle(X,tol);
    case 2
      [faux,gaux] = aux_handle(X,tol);
    case 3
      [faux,gaux,Haux] = aux_handle(X,tol);
    end

    % unnecessary indirection?
    f_fun = @(X) faux;
    dfdX_fun = @(X) gaux;
    d2fdX2_fun = @(X) Haux;
    f = sum(f_fun(X));
    if nargout<=1
      return;
    end
    dfdX = dfdX_fun(X);
    G = full_sparse(IJK,ones(size(IJK)),reshape(dfdX,size(IJK)),2*n,1);
    if nargout<=2
      return;
    end
    d2fdX2 = double(d2fdX2_fun(X));

    d2fdX2 = reshape(d2fdX2,[],6*6);
    if psd_project
      d2fdX2 = psd_project_rows(d2fdX2);
    end

    HI = repmat(IJK,[1 1 size(IJK,2)]);
    HJ = permute(repmat(IJK,[1 1 size(IJK,2)]),[1 3 2]);
    H = fast_sparse(HI(:),HJ(:),d2fdX2(:),2*n,2*n);

  end

  psd_project = true;
  % The max is not needed because we're handling that explicitly
  %
  % Smith and Schaefer
  %barrier = @(d,tol) (tol./d - 1).^2;
  % IPC
  barrier = @(D,tol) -(D-tol).^2.*log(D./tol);

  b = unique(E);
  [b1,b2] = box_each_element(V,b);
  [E1,E2] = box_each_element(V,E);
  I = box_intersect(b1-tol,b2+tol,E1-tol,E2+tol);
  I = I(~any(b(I(:,1))==E(I(:,2),:),2),:);
  BC = barycenter(V,E);

  [sqrD,T] = point_segment_squared_distance(V(b(I(:,1)),:),V(E(I(:,2),1),:),V(E(I(:,2),2),:));

  keep = find(sqrD<tol^2);

  I = I(keep,:);
  T = T(keep);

  %BC = barycenter(V,E);
  %plot_edges(V,E,'-ok');
  %hold on;
  %sct(V(b(I(:,1)),:),'or','LineWidth',2);
  %I(:,2)
  %size(E)
  %sct(BC(I(:,2),:),'og','LineWidth',2);
  %hold off;

  cap = T==0 | T==1;
  IJKline = [b(I(~cap,1)) E(I(~cap,2),:)];

  sqrD = sqrD(keep);

  Icap = I(cap,:);
  Tcap = T(cap);
  Jcap = E(Icap(:,2),1);
  Jcap(Tcap==1) = E(Icap(Tcap==1,2),2);
  IJcap = unique(sort([b(Icap(:,1)) Jcap],2),'rows');

  switch nargout
  case {0,1}
    [fline] = self_collision_barrier_line(V,E,IJKline);
    [fcap] = self_collision_barrier_cap(V,E,IJcap);
    f = fline + fcap;
  case 2
    [fline,Gline] = self_collision_barrier_line(V,E,IJKline);
    [fcap,Gcap] = self_collision_barrier_cap(V,E,IJcap);
    f = fline + fcap;
    G = Gline + Gcap;
  case 3
    [fline,Gline,Hline] = self_collision_barrier_line(V,E,IJKline);
    [fcap,Gcap,Hcap] = self_collision_barrier_cap(V,E,IJcap);
    f = fline + fcap;
    G = Gline + Gcap;
    H = Hline + Hcap;
  end

end

author: alecjacobson

repo: gptoolbox


text:
function [C]=structComplete(A,B,emptyFixOpt)

% function [C]=structComplete(A,B,emptyFixOpt)
% ------------------------------------------------------------------------
% This function fills in the missing data in the structure A with the
% content from B. The structure B can be seen as a default input structure
% for instance and A might be an incomplete set of inputs. The optional
% parameter emptyFixOpt (0=no, 1=yes) determines whether empty entries in A
% are overwritten by the corresponding default values in B.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2017/11/20
%------------------------------------------------------------------------

%%

fieldNameSet=fieldnames(B);

C=A; %Initialize C as the same as A
for q=1:1:numel(fieldNameSet) %Loop over field names
    fieldNameNow=fieldNameSet{q}; %Current field name
    if isfield(A,fieldNameNow) %If A contains the field in B

        if isempty(A.(fieldNameNow)) %if the field in A is empty
            if emptyFixOpt==1 %if empty field fixing is on
                C.(fieldNameNow)=B.(fieldNameNow); %Replace empty in C by default
            end
        end

        if isstruct(A.(fieldNameNow)) && isstruct(B.(fieldNameNow)) %If the field in A is a structure, check structure recursively
            [C.(fieldNameNow)]=structComplete(A.(fieldNameNow),B.(fieldNameNow),emptyFixOpt);
        end

    else %If the field is missing, add it
        C.(fieldNameNow)=B.(fieldNameNow);
    end
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 Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

% for i=1: size(X,1);
%   x = X(i, :)';
%   for k=1: K;
%     Z(i, k) = x' * U(:, k);
%   end
% end

U_reduce = U(:, 1:K);
Z = X * U_reduce;

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

end

author: benoitvallon

repo: coursera-machine-learning


text:
function h = imdisp(im, opt)

% imdisp -- Display an image
%
% Usage:
%       h = imdisp(im, opt)
%
% Input:
%       im          Image to display
%       opt         Display options
%
% Output:
%       h           Handle to image object
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2014-05-29
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 2,
  opt = [];
end

% Check for box plot region options
brgn = [];
if isfield(opt,'bregion'),
  brgn = opt.bregion;
end

% Check for zoom region options
zrgn = [];
if isfield(opt,'zregion'),
  zrgn = opt.zregion;
elseif isfield(opt,'zgap') && isfield(opt,'bregion'),
  zrgn = opt.bregion + opt.zgap*[-1 -1 1 1];
end
if ~isempty(zrgn),
  im = im(zrgn(1):zrgn(3),zrgn(2):zrgn(4),:);
  if ~isempty(brgn),
    brgn = brgn - [zrgn(1:2) zrgn(1:2)];
  end
end

% Handle image offset and scaling
if 0,
if isfloat(im) && min(im(:)) < 0,
  im = im - min(im(:));
end
if isfloat(im) && max(im(:)) > 1,
    im = im / max(im(:));
end
else
  if isfloat(im),
    im = im - min(im(:));
    im = im / max(im(:));
  end
end

% Different behaviour for greyscale and RGB images
if size(im,3) == 1,
  hh = imagesc(im);
  colormap(gray);
else
  hh = image(im);
end
axis image; axis off;

% Plot box region if option supplied
if ~isempty(brgn),
  bx = brgn([2 2 4 4 2]) + [-0.5 -0.5 0.5  0.5 -0.5] + 1;
  by = brgn([1 3 3 1 1]) + [-0.5  0.5 0.5 -0.5 -0.5] + 1;
  bl = 'r-';
  if isfield(opt,'bline'),
    bl = opt.bline;
  end
  hold on;
  plot(bx, by, bl);
  hold off;
end

% Check for title option
if isfield(opt, 'title'),
  title(opt.title);
end

if nargout > 0
    h = hh;
end

return

author: thfylsty


text:
function value = cc_abscissa ( order, i )

%*****************************************************************************80
%
%% CC_ABSCISSA returns the I-th abscissa of the Clenshaw Curtis rule.
%
%  Discussion:
%
%    Our convention is that the abscissas are numbered from left to
%    right.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, integer I, the index of the desired abscissa.  1 <= I <= ORDER.
%
%    Output, real VALUE, the value of the I-th abscissa in the
%    rule of order ORDER.
%
  if ( order < 1 )
    value = - Inf;
  elseif ( i < 1 | order < i )
    value = - Inf;
  elseif ( order == 1 )
    value = 0.0;
  elseif ( 2 * ( order - i ) == order - 1 )
    value = 0.0;
  else
    value = cos ( ( order - i ) * pi / ( order - 1 ) );
  end

  return
end


text:
function normals = vertexNormal(vertices, faces)
%VERTEXNORMAL Compute normals to a mesh vertices.
%
%   N = vertexNormal(V, F)
%   Computes vertex normals of the mesh given by vertices V and F.
%   V is a vertex array with 3 columns, F is either a NF-by-3 or NF-by-4
%   index array, or a cell array with NF elements.
%
%   Example
%     % Draw the vertex normals of a sphere
%     s = [10 20 30 40];
%     [v f] = sphereMesh(s);
%     drawMesh(v, f);
%     view(3);axis equal; light; lighting gouraud;
%     normals = vertexNormal(v, f);
%     drawVector3d(v, normals);
%
%   See also
%     meshes3d, meshFaceNormals, triangulateFaces

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-12-19, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

nv = size(vertices, 1);
nf = size(faces, 1);

% unit normals to the faces
faceNormals = normalizeVector3d(meshFaceNormals(vertices, faces));

% compute normal of each vertex: sum of normals to each face
normals = zeros(nv, 3);
if isnumeric(faces)
    for i = 1:nf
        face = faces(i, :);
        for j = 1:length(face)
            v = face(j);
            normals(v, :) = normals(v,:) + faceNormals(i,:);
        end
    end
else
    for i = 1:nf
        face = faces{i};
        for j = 1:length(face)
            v = face(j);
            normals(v, :) = normals(v,:) + faceNormals(i,:);
        end
    end
end

% normalize vertex normals to unit vectors
normals = normalizeVector3d(normals);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%Cytokinesis process test case
%
% Author: Jayodita Sanghvi, jayodita@stanford.edu
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/8/2010
classdef Cytokinesis_Test < edu.stanford.covert.cell.sim.ProcessTestCase
    methods
        function this = Cytokinesis_Test(name)
            this = this@edu.stanford.covert.cell.sim.ProcessTestCase(name);
        end

        function loadTestFixture(this)
            this.loadTestFixture@edu.stanford.covert.cell.sim.ProcessTestCase();
            m = this.process;

            c = m.chromosome;
            c.segregated = true;

            g = m.geometry;
            g.width = 200 * 1e-9;
            g.pinchedDiameter = 200 * 1e-9;
        end
    end

    %tests
    methods
        function testNumEdges(this)
            m = this.process;
            g = m.geometry;
            r = m.ftsZRing;
            assertTrue(r.numEdges > 0);
            assertTrue(r.numEdges * r.filamentLengthInNm*1e-9 < g.pinchedCircumference);
        end

        function testBindingUnlimited(this)
             m = this.process;
             r = m.ftsZRing;

             m.substrates(:) = 1e6;
             m.enzymes(m.enzymeIndexs_ftsZ_GTP_polymer) = 1e3;
             m.rateFilamentBindingMembrane = 1;  %force binding
             m.rateFtsZGtpHydrolysis = 0;

             m.evolveState();
             assertEqual(0, r.numEdgesOneStraight);
             assertEqual(r.numEdges, r.numEdgesTwoStraight);
             assertEqual(1e3-2*r.numEdges, m.enzymes(m.enzymeIndexs_ftsZ_GTP_polymer));
             assertEqual(2*r.numEdges, m.boundEnzymes(m.enzymeIndexs_ftsZ_GTP_polymer));
        end

        function testBindingPolymerLimited(this)
             m = this.process;
             r = m.ftsZRing;

             m.substrates(:) = 1e6;
             m.enzymes(m.enzymeIndexs_ftsZ_GTP_polymer) = 1;
             m.rateFilamentBindingMembrane = 1;  %force binding

             m.evolveState();
             assertEqual(1, r.numEdgesOneStraight);
             assertEqual(0, m.enzymes(m.enzymeIndexs_ftsZ_GTP_polymer));
             assertEqual(1, m.boundEnzymes(m.enzymeIndexs_ftsZ_GTP_polymer));
        end

        function testBendingUnlimited(this)
            m = this.process;
            g = m.geometry;
            r = m.ftsZRing;

            r.numEdgesTwoStraight = r.numEdges;
            m.substrates(:) = 1e6;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GTP_polymer) = 2*r.numEdges;
            m.rateFtsZGtpHydrolysis = 1;  %force bending
            m.rateFilamentDissociation = 0;

            initialDiameter = g.pinchedDiameter;

            m.evolveState();
            assertEqual(0, r.numEdgesTwoStraight);
            assertEqual(r.numEdges, r.numEdgesTwoBent);
            assertEqual(0, r.numResidualBent);
            assertEqual(2*r.numEdges, m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer));
            assertEqual(0, m.boundEnzymes(m.enzymeIndexs_ftsZ_GTP_polymer));
            assertEqual(1e6-(2*r.numEdges*r.numFtsZSubunitsPerFilament), m.substrates(m.substrateIndexs_water));
            assertEqual(1e6+(2*r.numEdges*r.numFtsZSubunitsPerFilament), m.substrates(m.substrateIndexs_phosphate));
            assertEqual(1e6+(2*r.numEdges*r.numFtsZSubunitsPerFilament), m.substrates(m.substrateIndexs_hydrogen));
            assertTrue(g.pinchedDiameter < initialDiameter);
        end

        function testBendingRateLimited(this)
            m = this.process;
            g = m.geometry;
            r = m.ftsZRing;

            r.numEdgesTwoStraight = r.numEdges;
            m.substrates(:) = 1e6;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GTP_polymer) = 30;

            initialDiameter = g.pinchedDiameter;

            m.evolveState();
            assertTrue(r.numEdgesTwoBent > 0);
            assertTrue(r.numEdgesTwoBent < r.numEdges);
            assertEqual(initialDiameter, g.pinchedDiameter);

            m.rateFtsZGtpHydrolysis = 1;  %force bending
            m.rateFilamentDissociation = 0;

            m.evolveState();
            assertEqual(0, r.numEdgesTwoStraight);
            assertEqual(r.numEdges, r.numEdgesTwoBent);
            assertTrue(g.pinchedDiameter < initialDiameter);
        end

        function testFirstRingFallingOffUnlimited(this)
            m = this.process;
            r = m.ftsZRing;

            m.substrates(:) = 0;
            m.enzymes(:) = 0;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer) = 30;
            r.numEdgesTwoBent = r.numEdges;
            m.rateFilamentDissociation = 1;  %force falling off

            m.evolveState();
            assertEqual(0, r.numEdgesTwoBent);
            assertEqual(r.numEdges, r.numResidualBent);
            assertEqual(r.numEdges, m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer));
            assertEqual(...
                r.numEdges*r.numFtsZSubunitsPerFilament,...
                m.enzymes(m.enzymeIndexs_ftsZ_GDP));
        end

        function testFirstRingFallingOffRateLimited(this)
            m = this.process;
            r = m.ftsZRing;

            m.substrates(:) = 0;
            m.enzymes(:) = 0;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer) = 30;
            r.numEdgesTwoBent = r.numEdges;

            m.evolveState();
            assertTrue(r.numEdgesTwoBent < r.numEdges);
            assertTrue(r.numEdgesTwoBent > 0);
            assertTrue(r.numResidualBent > 0);
            assertTrue(r.numResidualBent < r.numEdges);

            m.rateFilamentDissociation = 1;  %force falling off

            m.evolveState();
            assertEqual(0, r.numEdgesTwoBent);
            assertEqual(r.numEdges, r.numResidualBent);
        end

        function testSecondRingFallingOffUnlimited(this)
            m = this.process;
            r = m.ftsZRing;

            r.numEdgesOneStraight = r.numEdges;
            r.numResidualBent = 5;
            m.enzymes(:) = 0;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GTP_polymer) = r.numEdgesOneStraight;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer) = r.numResidualBent;

            m.rateFilamentDissociation = 1;

            m.evolveState();
            assertEqual(0, m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer));
            assertEqual(0, r.numResidualBent);
            assertEqual(...
                5 * r.numFtsZSubunitsPerFilament,...
                m.enzymes(m.enzymeIndexs_ftsZ_GDP));
        end

        function testSecondRingFallingOffRateLimited(this)
            m = this.process;
            r = m.ftsZRing;

            r.numEdgesOneStraight = r.numEdges;
            r.numResidualBent = 12;
            m.enzymes(:) = 0;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GTP_polymer) = r.numEdgesOneStraight;
            m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer) = r.numResidualBent;

            m.evolveState();
            assertTrue(r.numResidualBent < 15);
            assertTrue(r.numResidualBent > 0);

            m.rateFilamentDissociation = 1;

            m.evolveState();
            assertEqual(0, m.boundEnzymes(m.enzymeIndexs_ftsZ_GDP_polymer));
            assertEqual(0, r.numResidualBent);
            assertEqual(...
                12 * r.numFtsZSubunitsPerFilament,...
                m.enzymes(m.enzymeIndexs_ftsZ_GDP));
        end

        function testPinchCell(this)
            m = this.process;
            g = m.geometry;
            r = m.ftsZRing;

            m.enzymes(:) = 1e6;
            m.substrates(:) = 1e6;
            m.rateFtsZGtpHydrolysis = 1;
            m.rateFilamentBindingMembrane = 1;
            m.rateFilamentDissociation = 1;

            for j = 1:75
                m.evolveState();
            end

            assertEqual(0, g.pinchedDiameter);
            assertTrue(g.pinched);
            assertEqual(0, r.numEdges);
            assertEqual(0, r.numEdgesOneStraight);
            assertEqual(0, r.numEdgesTwoStraight);
            assertEqual(0, r.numEdgesTwoBent);
            assertEqual(0, r.numResidualBent);
        end

        function testNoEnzymes(this)
            m = this.process;
            m.enzymes(:) = 0;
            m.substrates(:) = 0;
            m.rateFtsZGtpHydrolysis = 1;
            m.rateFilamentBindingMembrane = 1;
            m.rateFilamentDissociation = 1;
            for j = 1:10
                m.evolveState();
            end

            g = m.geometry;
            assertEqual(g.width, g.pinchedDiameter);
        end

        function testNoWater(this)
            m = this.process;
            m.enzymes(:) = 1e6;
            m.substrates(:) = 1e6;
            m.substrates(m.substrateIndexs_water) = 0;
            m.rateFtsZGtpHydrolysis = 1;
            m.rateFilamentBindingMembrane = 1;
            m.rateFilamentDissociation = 1;
            for j = 1:10
                m.evolveState();
            end

            g = m.geometry;
            assertEqual(g.width, g.pinchedDiameter);
        end

        function testGeneEssentiality(this)
            m = this.process;
            g = m.geometry;
            g.pinchedDiameter = g.width;
            m.enzymes(:) = 1e6;
            m.substrates(:) = 1e6;
            m.rateFtsZGtpHydrolysis = 1;
            m.rateFilamentBindingMembrane = 1;
            m.rateFilamentDissociation = 1;

            initialPinchedDiameter = g.pinchedDiameter;

            this.helpTestGeneEssentiality(...
                {'MG_224'},... %cell division protein ftsZ
                @(m, i) m.geometry.pinchedDiameter < initialPinchedDiameter,...
                struct('lengthSec',10));
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell