text:
function [imageSeg,Rot]= getPlanSeg(XYZworldframeTest,SpaceTest,imgSize)
% input XYZworldframeTest is alinged in Z
%load([ '/n/fs/modelnet/NYUdataSet/NYUdatafeatureNew/' 'feature_' num2str(imageNum) '.mat'],'XYZworldframeTest','rgbTest','SpaceTest','imgDepth');
onPlaneThreshold =0.055;
sizethr =500;
normalAgreeThreshold =0.8;
removeNaN = find(~sum(isnan(XYZworldframeTest),2));
XYZworldframeTest = XYZworldframeTest(removeNaN,:);
[zhist,zc]=hist(XYZworldframeTest(:,3),round(20*(SpaceTest.Rz(2)-SpaceTest.Rz(1))));
[~,zpks] =findpeaks([0,zhist,0],'minpeakdistance',3);
zpks = zpks-1;
% figure,
% bar(bincenter,count);
% hold on;
% plot(bincenter(pks),count(pks),'rx');
%remove
range =[max(min(XYZworldframeTest(:,3)),zc(max(1,zpks-1))'),...
        min(zc(min(size(zc,2),zpks+1)),max(XYZworldframeTest(:,3)))';
        zc(end-2),max(XYZworldframeTest(:,3));
        min(XYZworldframeTest(:,3)),zc(3)];
removepts = inrange(XYZworldframeTest(:,3),range);
XYZworldframeTestRemove=XYZworldframeTest;
XYZworldframeTestRemove(removepts>0,:) =[];

% hough lines
[H,T,R] = votelines(XYZworldframeTestRemove(:,1),XYZworldframeTestRemove(:,2));
% blur it
allangleVote = max(H);
allangleVote = [allangleVote,allangleVote];
shift= 0;
clear rightAngle
for i =1:length(shift)
    rightAngle(i,:) = allangleVote(91+shift(i):180+shift(i))+allangleVote(1:90);
end
[~,AnglepickedLin]= max(rightAngle(:));
[s,Anglepicked]=ind2sub(size(rightAngle),AnglepickedLin);
[~,Rpicked1] = max(H(:,Anglepicked));
[~,Rpicked2] = max(H(:,Anglepicked+90+shift(s)));
P = [R(Rpicked1),T(Anglepicked);R(Rpicked2),T(Anglepicked)+90+shift(s)];

rotation = min(P(:,2));
Rot = getRotationMatrix('z',-1*rotation/180*pi);
XYZworldframeTestNew = [[Rot(1:2,1:2)*XYZworldframeTest(:,[1,2])']', XYZworldframeTest(:,3)];
normals = points2normals(XYZworldframeTestNew);

% find  corner
[xhist,xc]=hist(XYZworldframeTestNew(:,1),round((SpaceTest.Rx(2)-SpaceTest.Rx(1))/0.05));
[~,xpks] =findpeaks([0,xhist,0],'MINPEAKHEIGHT',50);
xpks = xpks-1;
[yhist,yc]=hist(XYZworldframeTestNew(:,2),round((SpaceTest.Ry(2)-SpaceTest.Ry(1))/0.05));
[~,ypks] =findpeaks([0,yhist,0],'MINPEAKHEIGHT',50);
ypks = ypks-1;
[zhist,zc]=hist(XYZworldframeTestNew(:,3),round((SpaceTest.Rz(2)-SpaceTest.Rz(1))/0.05));
[~,zpks] =findpeaks([0,zhist,0],'MINPEAKHEIGHT',50);
zpks = zpks-1;

%project points onto this plan and caculate conected component

gid =1;
imageSeg = zeros(imgSize);
for p =1:size(xpks,2)
    inLine =find(abs(XYZworldframeTestNew(:,1)-xc(xpks(p)))<onPlaneThreshold&abs(normals(1,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            imageSeg(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end

for p =1:size(ypks,2)
    inLine =find(abs(XYZworldframeTestNew(:,2)-yc(ypks(p)))<onPlaneThreshold&abs(normals(2,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            imageSeg(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end

for p =1:size(zpks,2)
    inLine =find(abs(XYZworldframeTestNew(:,3)-zc(zpks(p)))<onPlaneThreshold&abs(normals(3,:)')>normalAgreeThreshold);
    conthisSegMask = zeros(imgSize);
    conthisSegMask(removeNaN(inLine)) =1;
    label = bwlabel(conthisSegMask,8);
    unique_label = unique(label);
    unique_label(unique_label==0)=[];
    for i =1:length(unique_label),
         if sum(label(:)==unique_label(i))>sizethr,
            imageSeg(label==unique_label(i))=gid;
            gid = gid+1;
         end
    end
end
%%

% rotate back
%{
corner_r = [xc(xpks(1)) yc(ypks(1));xc(xpks(end)) yc(ypks(end))];
corner_r = get4corner(corner_r);
corner = [[Rot(1:2,1:2)'*corner_r(:,[1:2])']'];
D1 = cos(P(1,2)*pi/180)*corner(1,1)+sin(P(1,2)*pi/180)*corner(1,2);
D2 = cos(P(2,2)*pi/180)*corner(1,1)+sin(P(2,2)*pi/180)*corner(1,2);
D3 = cos(P(1,2)*pi/180)*corner(3,1)+sin(P(1,2)*pi/180)*corner(3,2);
D4 = cos(P(2,2)*pi/180)*corner(3,1)+sin(P(2,2)*pi/180)*corner(3,2);
P_new = [D1,P(1,2);D2,P(2,2);D3,P(1,2);D4,P(2,2)];
f = figure,
vis_point_cloud(XYZworldframeTest,rgbTest,10,5000);hold on;
for j =1:3
    plot3(corner([j,j+1],1),corner([j,j+1],2),[max(XYZworldframeTest(:,3));max(XYZworldframeTest(:,3))],'-xr','LineWidth',10)
end
plot3(corner([1,4],1),corner([1,4],2),[max(XYZworldframeTest(:,3));max(XYZworldframeTest(:,3))],'-xr','LineWidth',10)

for i =1:gid
    hold on;
    plot3(XYZworldframeTestOrg(imageSeg(:)==i,1),XYZworldframeTestOrg(imageSeg(:)==i,2),XYZworldframeTestOrg(imageSeg(:)==i,3),'+','Color',rand([1,3]));
end

axis equal;
axis tight;
view(30,50)
saveas(f,['./result/' num2str(imageNum) '.fig']);
saveas(f,['./result/' num2str(imageNum) '.jpg']);

for gid =1:length(gropuind)
    imageSeg(removeNaN(gropuind{gid})) = gid;
end
%}

%figure(1),imagesc(imageSeg)
%Boudary = [D1 D2 D3 D4];
%{
if imageNum> 0,
    im = getImagesc(imageSeg);
    mkdir(segpath);
    imwrite(im,sprintf('%s/%04d.jpg',segpath,imageNum));
end
%}
end
function [hough_transform,T,R] = votelines(X,Y)
         thetaResolution = 1;
         rhoResolution = 0.1;
         T = [1:180];
         theta = T/180*pi;
         rho = X(:)*cos(theta)+ Y(:)*sin(theta);
         % quantize the rho
         rhoNorm = max(1,round((rho-min(rho(:)))/rhoResolution));
         R =[1:max(rhoNorm(:))]*rhoResolution+min(rho(:));
         hough_transform = zeros(max(rhoNorm(:)),size(T,2));
         TT = repmat(T,[size(rhoNorm,1),1]);
         hough_transform = accumarray([rhoNorm(:),TT(:)],1,size(hough_transform));
end

author: thusiyuan


text:
% Command list:
% turtleBot('openSerial') - open turtlebot serial port
% turtleBot('setSpeed', trans, rot) - set translational and rotational speed
% pos = turtleBot('getPosition') - get turtleBot position, return value is [x; y; angle]
% turtleBot('goTo', pos) - go to target position
% turtleBot('bye') - close serial and kill all threads

% compile
% mex turtleBot.cpp

% To get connected to the turtleBot, we need to open the serial port
turtleBot('openSerial');

% To simply set speed
% unit: translational speed - m/s, rotational speed - rad/s
turtleBot('setSpeed', trans, rot);

% To obtain the current position
pos = turtleBot('getPosition');
disp(pos);

% Go to a target position, pos = [x y angle heading].
% unit: x - mm y - mm angle - deg heading - TOWARD 0; BACKWARD 1
turtleBot('goTo', pos);

% IMPORTANT: remember to shutdown after finishing everything!
turtleBot('bye');

author: jianxiongxiao

repo: ProfXkit


text:
function h = quiver(sAF,varargin)
% quiver spherical axis field
%
% Syntax
%   quiver3(sAF)
%
% See also
%   S2VectorField/plot
%

% plot the function values
h = plot(sAF,varargin{:});

% remove output if not required
if nargout == 0, clear h; end

end

author: mtex-toolbox

repo: mtex

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


text:
%  INTERNAL FUNCTION: dirichlet distribution
%
%  Note:
%     One would rarely need to call distribution functions directly.
%     Refer to <densities.html> file for list of supported distributions.
%
%

author: jmaih


text:
function [x_g,w_g,phi,p_x,p_y,p_z] = threed_shapeiso(x_local,r,s,t,w)
%-----------------------------------------------------------------------
%  threed_shapeiso.m - computes test functions and derivatives on an
%                      element given element coordinates and Gauss points.
%                      Isoparametric coordinates are used (curved elements)
%
%  Copyright (c) 2002, Jeff Borggaard, Virginia Tech
%  Version: 1.0a
%
%  Usage:    [x_g,w_g,phi,p_x,p_y,p_z] = threed_shapeiso(x_local,r,s,t,w)
%
%  Variables:     x_local
%                        Coordinates of the element nodes
%                 (r,s,t)
%                        Coordinates of Gauss points in unit tetrahedron
%                 w
%                        Gauss weights associated with (r,s,t)
%
%                 x_g
%                        Coordinates of Gauss points in the element
%                 w_g
%                        Gauss weights scaled by the element Jacobian
%                 phi
%                        Value of element shape functions at x_g
%                 p_x
%                 p_y
%                 p_z
%                        First spatial derivatives of phi
%                 p_xx, etc.
%                        Second spatial derivatives of phi
%                        (only applicable with 2nd order or higher)
%                        (currently not implemented)
%-----------------------------------------------------------------------
  x      = x_local;
  [n,t1] = size(x);  % t1 had better be 3
  rule   = length(r);

  if (n == 10)
    % The following assumes isoparametric elements
    c0 =   x( 1,:);                                      % 1
    c1 =-3*x( 1,:) -   x( 2,:) + 4*x( 5,:)            ;  % r
    c2 =-3*x( 1,:) -   x( 3,:) + 4*x( 7,:)            ;  % s
    c3 =-3*x( 1,:) -   x( 4,:) + 4*x(10,:)            ;  % t
    c4 = 2*x( 1,:) + 2*x( 2,:) - 4*x( 5,:)            ;  % r^2
    c5 = 2*x( 1,:) + 2*x( 3,:) - 4*x( 7,:)            ;  % s^2
    c6 = 2*x( 1,:) + 2*x( 4,:) - 4*x(10,:)            ;  % t^2
    c7 = 4*x( 1,:) - 4*x( 5,:) + 4*x( 6,:) - 4*x( 7,:); % rs
    c8 = 4*x( 1,:) - 4*x( 5,:) + 4*x( 8,:) - 4*x(10,:); % rt
    c9 = 4*x( 1,:) - 4*x( 7,:) + 4*x( 9,:) - 4*x(10,:); % st

    x_g(:,1) = c0(1) + c1(1)*r + c2(1)*s + c3(1)*t + c4(1)*r.^2 ...
             + c5(1)*s.^2 + c6(1)*t.^2 + c7(1)*r.*s + c8(1)*r.*t ...
             + c9(1)*s.*t;
    xr = c1(1) + 2*c4(1)*r + c7(1)*s + c8(1)*t;
    xs = c2(1) + 2*c5(1)*s + c7(1)*r + c9(1)*t;
    xt = c3(1) + 2*c6(1)*t + c8(1)*r + c9(1)*s;

    x_g(:,2) = c0(2) + c1(2)*r + c2(2)*s + c3(2)*t + c4(2)*r.^2 ...
             + c5(2)*s.^2 + c6(2)*t.^2 + c7(2)*r.*s + c8(2)*r.*t ...
             + c9(2)*s.*t;
    yr = c1(2) + 2*c4(2)*r + c7(2)*s + c8(2)*t;
    ys = c2(2) + 2*c5(2)*s + c7(2)*r + c9(2)*t;
    yt = c3(2) + 2*c6(2)*t + c8(2)*r + c9(2)*s;

    x_g(:,3) = c0(3) + c1(3)*r + c2(3)*s + c3(3)*t + c4(3)*r.^2 ...
             + c5(3)*s.^2 + c6(3)*t.^2 + c7(3)*r.*s + c8(3)*r.*t ...
             + c9(3)*s.*t;
    zr = c1(3) + 2*c4(3)*r + c7(3)*s + c8(3)*t;
    zs = c2(3) + 2*c5(3)*s + c7(3)*r + c9(3)*t;
    zt = c3(3) + 2*c6(3)*t + c8(3)*r + c9(3)*s;

    % Compute the Jacobian of the (r,s,t) -> (x,y,z) transformation
    jac = ( xr.*ys.*zt + xs.*yt.*zr + xt.*yr.*zs ...
           -xr.*yt.*zs - xs.*yr.*zt - xt.*ys.*zr );

    w_g = jac.*w;

    % Invert derivatives of the mapping
    rx = ( ys.*zt-yt.*zs )./jac;
    ry = ( xt.*zs-xs.*zt )./jac;
    rz = ( xs.*yt-xt.*ys )./jac;

    sx = ( yt.*zr-yr.*zt )./jac;
    sy = ( xr.*zt-xt.*zr )./jac;
    sz = ( xt.*yr-xr.*yt )./jac;

    tx = ( yr.*zs-ys.*zr )./jac;
    ty = ( xs.*zr-xr.*zs )./jac;
    tz = ( xr.*ys-xs.*yr )./jac;

    % Compute shape function and derivatives at Gauss points
    u  = 1 - r  - s  - t ;
    ux =   - rx - sx - tx;
    uy =   - ry - sy - ty;
    uz =   - rz - sz - tz;

    phi = zeros(rule,n);
    phi(:,1)  = u - 2*r.*u - 2*s.*u - 2*t.*u;
    phi(:,2)  = r - 2*r.*u - 2*r.*s - 2*r.*t;
    phi(:,3)  = s - 2*r.*s - 2*s.*u - 2*s.*t;
    phi(:,4)  = t - 2*r.*t - 2*s.*t - 2*t.*u;
    phi(:,5)  = 4*r.*u;
    phi(:,6)  = 4*r.*s;
    phi(:,7)  = 4*s.*u;
    phi(:,8)  = 4*r.*t;
    phi(:,9)  = 4*s.*t;
    phi(:,10) = 4*t.*u;

    p_x = zeros(rule,n);
    p_x(:,1)  = ux - 2*rx.*u - 2*r.*ux - 2*sx.*u ...
                   - 2*s.*ux - 2*tx.*u - 2*t.*ux    ;
    p_x(:,2)  = rx - 2*rx.*u - 2*r.*ux - 2*rx.*s ...
                   - 2*r.*sx - 2*rx.*t - 2*r.*tx    ;
    p_x(:,3)  = sx - 2*rx.*s - 2*r.*sx - 2*sx.*u ...
                   - 2*s.*ux - 2*sx.*t - 2*s.*tx    ;
    p_x(:,4)  = tx - 2*rx.*t - 2*r.*tx - 2*sx.*t ...
                   - 2*s.*tx - 2*tx.*u - 2*t.*ux    ;
    p_x(:,5)  = 4*rx.*u + 4*r.*ux;
    p_x(:,6)  = 4*rx.*s + 4*r.*sx;
    p_x(:,7)  = 4*sx.*u + 4*s.*ux;
    p_x(:,8)  = 4*rx.*t + 4*r.*tx;
    p_x(:,9)  = 4*sx.*t + 4*s.*tx;
    p_x(:,10) = 4*tx.*u + 4*t.*ux;

    p_y = zeros(rule,n);
    p_y(:,1)  = uy - 2*ry.*u - 2*r.*uy - 2*sy.*u ...
                   - 2*s.*uy - 2*ty.*u - 2*t.*uy    ;
    p_y(:,2)  = ry - 2*ry.*u - 2*r.*uy - 2*ry.*s ...
                   - 2*r.*sy - 2*ry.*t - 2*r.*ty    ;
    p_y(:,3)  = sy - 2*ry.*s - 2*r.*sy - 2*sy.*u ...
                   - 2*s.*uy - 2*sy.*t - 2*s.*ty    ;
    p_y(:,4)  = ty - 2*ry.*t - 2*r.*ty - 2*sy.*t ...
                   - 2*s.*ty - 2*ty.*u - 2*t.*uy    ;
    p_y(:,5)  = 4*ry.*u + 4*r.*uy;
    p_y(:,6)  = 4*ry.*s + 4*r.*sy;
    p_y(:,7)  = 4*sy.*u + 4*s.*uy;
    p_y(:,8)  = 4*ry.*t + 4*r.*ty;
    p_y(:,9)  = 4*sy.*t + 4*s.*ty;
    p_y(:,10) = 4*ty.*u + 4*t.*uy;

    p_z = zeros(rule,n);
    p_z(:,1)  = uz - 2*rz.*u - 2*r.*uz - 2*sz.*u ...
                   - 2*s.*uz - 2*tz.*u - 2*t.*uz    ;
    p_z(:,2)  = rz - 2*rz.*u - 2*r.*uz - 2*rz.*s ...
                   - 2*r.*sz - 2*rz.*t - 2*r.*tz    ;
    p_z(:,3)  = sz - 2*rz.*s - 2*r.*sz - 2*sz.*u ...
                   - 2*s.*uz - 2*sz.*t - 2*s.*tz    ;
    p_z(:,4)  = tz - 2*rz.*t - 2*r.*tz - 2*sz.*t ...
                   - 2*s.*tz - 2*tz.*u - 2*t.*uz    ;
    p_z(:,5)  = 4*rz.*u + 4*r.*uz;
    p_z(:,6)  = 4*rz.*s + 4*r.*sz;
    p_z(:,7)  = 4*sz.*u + 4*s.*uz;
    p_z(:,8)  = 4*rz.*t + 4*r.*tz;
    p_z(:,9)  = 4*sz.*t + 4*s.*tz;
    p_z(:,10) = 4*tz.*u + 4*t.*uz;

  else
    error('Only quadratic isoparametric elements are supported\n')
  end

author: johannesgerer

repo: jburkardt-m


text:
function test_suite = test_get_absolute_path
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
            test_functions=localfunctions();
        catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_get_absolute_path_basics()
    aeq=@(a,b)assertEqual(mocov_get_absolute_path(a),b);

    aeq('/','/');
    aeq('/foo/../','/');
    aeq('/foo/..//','/');
    aeq('/foo/..','/');
    aeq('/foo/../.','/');
    aeq('/foo/.././','/');

    orig_pwd=pwd();
    cleaner=onCleanup(@()cd(orig_pwd));
    p=fileparts(mfilename('fullpath'));
    cd(p);
    aeq('',p);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [F,vars] = sort_internal(t,X,data);%ii,D,V);

% Hack to figure out all the sorted variables, not just this index.
% Sort is implemented in a slighlt different way (general feature in
% future versions) that allows one element in an operator to modell all
% elements. Reduces the number of calls to the operator code.

ii = data.i;
D = data.D;
V = data.V;

var_start = getvariables(t)-ii+1;
n = length(X);
var_end   = getvariables(t)-ii+1 + n-1;
vars = var_start:var_end;

% Is this a location variable instead of the actuial sort variable. If so,
% shift everything back to get the indicies of the sorted variables.
if data.isthisloc == 1
    vars = vars - n;
end

t   = recover(vars);
loc = recover(vars+n);

[M,m] = derivebounds(X);
X = reshape(X,1,n);

% Standard model
F = (sum(D,1) == 1) + (sum(D,2) == 1);
F = F + (t == sum(V,2));
F = F + (diff(t) >= 0);
for i = 1:n
   di = D(i,:);
   vi = V(i,:);
   F = F + (-(-m)'.*(1-di) <= X-vi <= (M)'.*(1-di));
   F = F + (m'.*di <=  vi <= M'.*di);
end

% Cuts
F = F + (X == sum(V,1));
F = F + (sum(t) == sum(X));

% Definition of location
F = F + (loc == D*[(1:n)']) + (1 <= loc <= n);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%GRIDMERGE Merge two grids.
%
%   [ GRID ] = GRIDMERGE( GRID1, IND1, GRID2, IND2, I_DELETE, I_CHECK ) Merges
%   the matching boundaries edges and faces specified by boundary numbers IND1
%   in GRID1 and IND2 in GRID2. The algorithm tries to match up the boundary
%   with the fewest boundary edges/faces and nodes to the closest corresponding
%   ones on GRID2. The optional argument I_DELETE selects which grid to delete
%   duplicate nodes from (I_DELETE is 2 per default). I_CHECK toggles error
%   checking on/off.
%
%   Can optionally be called as GRIDMERGE( GRID1, GRID2 ) where IND1 and IND2
%   are determined by any shared boundary numbers (also in this case boundaries
%   and subdomains are not renumbered).
%
%   Examples:
%
%      1) Create a 2D rivet by joining a rectangle and half circle.
%
%      grid1 = rectgrid( 4, 10, [0.3 0.7;0 1] );
%      grid2 = circgrid( 4, 3, 0.5, [0.5;1] );
%      grid2 = delcells( grid2, 'y<=1' );
%      grid  = gridmerge( grid1, 3, grid2, 4 );
%
%      2) Create a flow over cylinder benchmark grid by merging three grids.
%
%      grid1 = ringgrid( [0.05 0.06 0.08 0.11 0.15], 32, [], [], [0.2;0.2] );
%      grid2 = holegrid( 8, 1, [0 0.41;0 0.41], 0.15, [0.2;0.2] );
%      grid2 = gridmerge( grid1, 5:8, grid2, 1:4 );
%      grid3 = rectgrid( [0.41 0.5 0.7 1 1.4 1.8 2.2], 8, [0.41 2.2;0 0.41] );
%      grid  = gridmerge( grid3, 4, grid2, 6 );
%
%      3) Create a 3D grid with two brackets attached to an I-beam section.
%
%      grid01 = ringgrid( 1, 20, 0.03, 0.06, [0;0] );
%      indc01 = selcells( grid01, 'y<=sqrt(eps)' );
%      grid01 = delcells( grid01, indc01 );
%
%      grid02 = holegrid( 5, 1, .06*[-1 1;-1 1], .03, [0;0] );
%      indc02 = selcells( grid02, 'y>=-sqrt(eps)' );
%      grid02 = delcells( grid02, indc02 );
%      grid2d = gridmerge( grid01, [5 6], grid02, [7 8] );
%
%      grid1 = gridextrude( grid2d, 1, 0.02 );
%      grid1 = gridrotate( grid1, pi/2, 1 );
%      grid2 = grid1;
%      grid1.p(2,:) = grid1.p(2,:) + 0.03;
%      grid2.p(2,:) = grid2.p(2,:) - 0.01;
%
%      x_coord = [ -0.08 linspace(-0.06,0.06,6) 0.08 ];
%      y_coord = [ -0.2 -0.15 -0.1 -0.05 -0.03 -0.01 ...
%                   0.01  0.03  0.05  0.1  0.15  0.2 ];
%      grid3 = blockgrid( x_coord, y_coord, 1, ...
%                         [-0.08 0.08;-0.2 0.2;-0.08 -0.06] );
%      grid4 = blockgrid( 1, y_coord, 5, ...
%                         [-0.01 0.01;-0.2 0.2;-0.18 -0.08] );
%      grid5 = grid3;
%      grid5.p(3,:) = grid5.p(3,:) - 0.12;
%
%      grid = gridmerge( grid1, 8, grid3, 6 );
%      grid = gridmerge( grid2, 8, grid, 19 );
%      grid = gridmerge( grid4, 6, grid, 24, 1 );
%      grid = gridmerge( grid5, 6, grid, 33, 2 );
%
%   See also GRIDEXTRUDE, GRIDREVOLVE, GRIDROTATE, GRIDSCALE

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function bbox = mask2box( mask )
    tmp = regionprops(double(mask),'BoundingBox'); % Double to force a single bbox
    bbox(1) = tmp.BoundingBox(2)+0.5;
    bbox(2) = tmp.BoundingBox(1)+0.5;
    bbox(3) = tmp.BoundingBox(2)+tmp.BoundingBox(4)-0.5;
    bbox(4) = tmp.BoundingBox(1)+tmp.BoundingBox(3)-0.5;
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
clc
clear
close all

% This example shows how to interpolate the fMRI timecourse if the stimulus
% duration is not equal to the duration of the TR. Imagine that there are
% 100 TRs and each lasts for 1s, while your stimulus duration is only 0.5s
% The stimulus starts with the TR onset. To correctly code the stimulus
% onset in your design matrix you need 200 rows in your design matrix
% (but you only have 100 fMRI volumes). Interpolation of the timecourse is
% necessary to match the size of rows in your design matrix.
%%

% Simulated data properties.
TRs = 100;
TR = 1;
stimdur = 0.5;
TRs_after_resampling = TR/stimdur*TRs;

% Use an example hrf and create an example fMRI time series.
cond1 = zeros(TRs,1);
cond1(1:20:end) = 1;

hrf = getcanonicalhrf(0.5,1);
tcs = conv(cond1,hrf);
tcs = tcs(1:TRs);
figure(1);clf
plot(0:TR:TRs-TR,tcs,'-','LineWidth',2); hold on
ylabel('%BOLD')
xlabel('TRs')

% The following line resmaples the timecourse so that each timepoint will
% correspond to 0.5 s instead of 1 s. The tseriesinterp is not a function
% available in GLMsingle but you can download it form github
% https://github.com/cvnlab/knkutils/blob/master/timeseries/tseriesinterp.m

tcs_interp = tseriesinterp(tcs,1,0.5);

% plot the interpolated timeseries
plot(0:stimdur:TRs-stimdur,tcs_interp,'o','MarkerSize',3,'LineWidth',2)
legend box off
set(gca,'FontSize',15)

%%
whos tcs tcs_interp

% Notice that the lenght of the tcs_interp is double the length of tcs.

% With an interpolated timecourse now you can code your design matrix
% correctly. The design matrix is going to consist of 200 columns were 1 will
% specify the stimulus onset. Remember to specify the stimduration for
% GLMsingle as 0.5 s instead of 1 s.

%% show the dm as stem plot with correct length

dm = repelem(cond1,2);
find_rep = diff(dm);
dm(find_rep==-1) = 0;
dm(dm==0) = NaN;
stem(0:stimdur:TRs-stimdur,dm,'filled','LineWidth',2)
legend({'Original tcs';'Interpolated tcs';'Stimulus onset'},'Location','EastOutside')
set(gcf,'Position',[1000        1090        1234         247])

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function result = quickif(condition,ontrue,onfalse)
% Select one of two possible values dependent on a condition
% Result = quickif(Condition,OnTrue,OnFalse)
%
% In:
%   Condition : conditional expression
%
%   OnTrue : value to return if condition evaluates to true
%
%   OnFalse : value to return if condition evaluates to false
%
% Out:
%   Result : either OnTrue or OnFalse
%
% Notes:
%   Note that both branches will always be evaluated, regardless of the value of Condition, since
%   quickif() is a regular function call. This is purely a convenience function that allows to use
%   conditional expressions inside other expressions.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-10-16

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

if condition
    result = ontrue;
else
    result = onfalse;
end

repo: REST


text: %PLOT Mapping overload (generates error)

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [au] = m2au(m)
% Convert length from meters to astronomical units.
% Chad A. Greene 2012
au = m*6.684587122671e-12;

author: Sable

repo: mcbench-benchmarks


text:
%STEREORECTIFYUNCALIBRATED  Computes a rectification transform for an uncalibrated stereo camera
%
%     [H1,H2,success] = cv.stereoRectifyUncalibrated(points1, points2, F, imageSize)
%     [...] = cv.stereoRectifyUncalibrated(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __points1__ Array of feature points in the first image as a cell array of
%   2-element vectors: `{[x1, y1], [x2, y2], ...}` or an Nx2/Nx1x2/1xNx2
%   numeric array. The same formats as in cv.findFundamentalMat are supported.
% * __points2__ The corresponding points in the second image, same size and
%   type as `points1`.
% * __F__ Input 3x3 fundamental matrix. It can be computed from the same set
%   of point pairs using cv.findFundamentalMat.
% * __imageSize__ Size of the image `[w,h]`.
%
% ## Output
% * __H1__ 3x3 rectification homography matrix for the first image.
% * __H2__ 3x3 rectification homography matrix for the second image.
% * __success__ success flag. Returns true if successfull, false otherwise.
%
% ## Options
% * __Threshold__ Optional threshold used to filter out the outliers. If the
%   parameter is greater than zero, all the point pairs that do not comply
%   with the epipolar geometry (that is, the points for which
%   `|points2{i}' * F * points1{i}| > Threshold`) are rejected prior to
%   computing the homographies. Otherwise, all the points are considered
%   inliers. default 5
%
% The function computes the rectification transformations without knowing
% intrinsic parameters of the cameras and their relative position in the
% space, which explains the suffix "uncalibrated". Another related difference
% from cv.stereoRectify is that the function outputs not the rectification
% transformations in the object (3D) space, but the planar perspective
% transformations encoded by the homography matrices `H1` and `H2`. The
% function implements the algorithm [Hartley99].
%
% ### Note
% While the algorithm does not need to know the intrinsic parameters of the
% cameras, it heavily depends on the epipolar geometry. Therefore, if the
% camera lenses have a significant distortion, it would be better to correct
% it before computing the fundamental matrix and calling this function. For
% example, distortion coefficients can be estimated for each head of stereo
% camera separately by using cv.calibrateCamera. Then, the images can be
% corrected using cv.undistort, or just the point coordinates can be corrected
% with cv.undistortPoints.
%
% ## References
% [Hartley99]:
% > Richard I Hartley. "Theory and practice of projective rectification".
% > International Journal of Computer Vision, 35(2):115-127, 1999.
%
% See also: cv.stereoCalibrate, cv.stereoRectify, cv.calibrateCamera,
%  cv.undistort, cv.undistortPoints, estimateUncalibratedRectification,
%  rectifyStereoImages
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [LB,UB,PLB,PUB,fixidx] = boundscheck(x0,LB,UB,PLB,PUB)
%BOUNDSCHECK Initial check of bounds and see if there are fixed variables.

nvars = numel(x0);

% Expand scalar inputs to full vectors
if isscalar(LB); LB = LB*ones(1,nvars); end
if isscalar(UB); UB = UB*ones(1,nvars); end
if isscalar(PLB); PLB = PLB*ones(1,nvars); end
if isscalar(PUB); PUB = PUB*ones(1,nvars); end

if isempty(PLB) || isempty(PUB)
    warning('bads:pbUnspecified', 'Plausible lower/upper bounds PLB and/or PUB not specified. Using hard upper/lower bounds instead.');
    if isempty(PLB); PLB = LB; end
    if isempty(PUB); PUB = UB; end
end

% Test that all vectors have the same length
if any([numel(LB),numel(UB),numel(PLB),numel(PUB)] ~= nvars)
    error('All input vectors (X0, LB, UB, PLB, PUB), if specified, need to have the same size.');
end

% Test that all vectors are row vectors
if ~isvector(x0) || ~isvector(LB) || ~isvector(UB) || ~isvector(PLB) || ~isvector(PUB) ...
        || size(x0,1) ~= 1 || size(LB,1) ~= 1 || size(UB,1) ~= 1 || size(PLB,1) ~= 1 || size(PUB,1) ~= 1
    error('All input vectors should be row vectors.');
end

% Test that plausible bounds are finite
if ~all(isfinite(([PLB, PUB])))
    error('Plausible interval bounds PLB and PUB need to be finite.');
end

% Test that all vectors are real-valued
if ~isreal([x0; LB; UB; PLB; PUB])
    error('All input vectors should be real-valued.');
end

% Fixed variables (all bounds equal)
fixidx = (x0 == LB) & (LB == UB) & (UB == PLB) & (PLB == PUB);

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


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

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

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

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

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

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

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

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

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ states ] = BuildStateList
%BuildStateList builds a state list from a state matrix

% state discretization for the mountain car problem
xdiv  = (0.55-(-1.5))   / 50.0;
xpdiv = (0.07-(-0.07)) / 50.0;

x = -1.5:xdiv:0.5;
xp= -0.07:xpdiv:0.07;

N=size(x,2);
M=size(xp,2);

states=[];
index=1;
for i=1:N
    for j=1:M
        states(index,1)=x(i);
        states(index,2)=xp(j);
        index=index+1;
    end
end

author: epfl-lasa


function [batches, masked_batches, batch_padding] = rcnn_extract_regions(im, sp, reg2sp, boxes, rcnn_model)

% convert image to BGR and single
im = single(im(:,:,[3 2 1]));
num_boxes = size(boxes, 1);
batch_size = rcnn_model.cnn.batch_size;
num_batches = ceil(num_boxes / batch_size);
batch_padding = batch_size - mod(num_boxes, batch_size);
if(mod(num_boxes, batch_size)==0) batch_padding=0; end
crop_mode = rcnn_model.detectors.crop_mode;
image_mean = rcnn_model.cnn.image_mean;
crop_size = size(image_mean,1);
crop_padding = rcnn_model.detectors.crop_padding;

batches = cell(num_batches, 1);
masked_batches = cell(num_batches, 1);
for batch = 1:num_batches
%  disp(batch);
%parfor batch = 1:num_batches
  batch_start = (batch-1)*batch_size+1;
  batch_end = min(num_boxes, batch_start+batch_size-1);

  ims = zeros(crop_size, crop_size, 3, batch_size, 'single');
  masked_ims=zeros(crop_size, crop_size, 3, batch_size, 'single');
  for j = batch_start:batch_end
    bbox = boxes(j,:);
    m1=reg2sp(:,j);
    mask=double(m1(sp));
    [crop, mask_crop] = rcnn_im_crop_mask(im, mask, bbox, crop_mode, crop_size, ...
        crop_padding, image_mean);
    % swap dims 1 and 2 to make width the fastest dimension (for caffe)
    ims(:,:,:,j-batch_start+1) = permute(crop, [2 1 3]);
    masked_ims(:,:,:,j-batch_start+1) = permute(mask_crop, [2 1 3]);
  end

  batches{batch} = ims;
  masked_batches{batch} = masked_ims;
end


text:
% Small technical example.
% Demonstrates an expedient for showing an image.
%
disp('Small technical example.');
disp('Demonstrates an expedient for showing an image.');
disp('FOR MORE INFORMATION:  help printshop');
disp(' ');
disp('See also the report http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04178D.pdf');
disp('Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>');
disp(' (C) 1998-2006 Stichting CWI, Amsterdam, The Netherlands');
disp(' ');
%------------------------------------------------------------------------------
load zenithgray; Orig = zenithgray; clear zenithgray;
printshop('original','Nice watch', Orig, 3, 0, []);


text:
function [mm] = ly2mm(ly)
% Convert length from light years to millimeters.
% Chad A. Greene 2012
mm = ly*9460730472581000000;


text:
% the following matlab functions are used to calculate the Rand,
% adjusted Rand, Wallace and other partition comparison coefficients
%
%    Copyright (C) 2009  UMMI@IMM
%
%    This file is part of Comparing Partitions website <http://www.comparingpartitions.info/>.
%    Comparing Partitions website 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/>.
%
%[a,b,c,d,bc,dn,confmat,res]=PartAgreeCoef(c1,c2)
%Outputs:
%ri=rand Index
%AR = adjusted rand
%jac=jaccard
%w1- Wallace c1->c2
%w2 - Wallace c2->c1
%lar1 - Larsen c1->c2
%lar2 - Larsen c2->c1
%MH - Melia & Heckerman
%VI - variation of information
%nVI - normalized variation of information
% Francisco Pinto fpinto@fm.ul.pt

function res=PartAgreeCoef_ARonly(c1,c2)

c1=c1-min(c1)+1;
c2=c2-min(c2)+1;
n=length(c1); % number cass
ng1=max(c1);
ng2=max(c2);
%dn=n*(n-1)/2; %number of possible pairwise comparisons of cases(a+b+c+d)

%y1=pdist(c1(:),'ham');
%y2=pdist(c2(:),'ham');
%size(y1)
%size(y2)
%ad=sum((y1'==y2')); % number of pairwise concordances (matches (a) and mismatches(d))(a+d)
%bc=sum((y1'~=y2')); % number of pairwise discordances(b+c)
%Rand Index
%res.ri=ad/dn;

%a=sum((y1'==0).*(y2'==0));
%b=sum(y1'==0)-a;
%c=sum(y2'==0)-a;

%res.w1=a/sum(y1'==0);% check is =dn
%res.w1a=a/(a+b);
%res.w2a=a/(a+c);
%res.w2=a/sum(y2'==0);
%d=ad-a;

%Jaccard Index
%res.jac=a/(dn-d);

%confmat=crosstab(c1,c2);
confmat=full(sparse(c1,c2,1,ng1,ng2));

coltot=sum(confmat);
rowtot=sum(confmat')';
%summat=repmat(coltot,ng1,1)+repmat(rowtot,1,ng2);
%larsenmat=2*confmat./summat;
%res.lar1=mean(max((larsenmat'))');
%res.lar2=mean(max(larsenmat)');

%todelmat=larsenmat;
%cumval=0;
%for i=1:min([ng1;ng2])
%    [val]=max(max(todelmat)');
%    [rr,cc]=find(todelmat==val);
%        todelmat(rr(1),:)=0;
%        todelmat(:,cc(1))=0;
%        cumval=cumval+confmat(rr(1),cc(1));
%end
%res.MH=cumval/n;

%H1=-sum((rowtot/n).*log2((rowtot/n)));
%H2=-sum((coltot/n)'.*log2((coltot/n)'));
%indmat=(rowtot/n)*(coltot/n);
%nozeromat=(confmat/n)+(confmat==0);
%H12=-sum(sum((confmat/n).*log2(nozeromat)));
%MI=H1+H2-H12;
%res.VI=H1+H2-2*MI;
%res.nVI=res.VI/log2(n);

nis=sum(rowtot.^2);		%sum of squares of sums of rows
njs=sum(coltot.^2);		%sum of squares of sums of columns

t1=nchoosek(n,2);		%total number of pairs of entities
t2=sum(sum(confmat.^2));	%sum over rows & columnns of nij^2
t3=.5*(nis+njs);

%Expected index (for adjustment)
nc=(n*(n^2+1)-(n+1)*nis-(n+1)*njs+2*(nis*njs)/n)/(2*(n-1));

A=t1+t2-t3;		%no. agreements
%D=  -t2+t3;		%no. disagreements

if t1==nc
   res=0;			%avoid division by zero; if k=1, define Rand = 0
else
   res=(A-nc)/(t1-nc);		%adjusted Rand - Hubert & Arabie 1985
   %res.AR2=(ad-nc)/(dn-nc); % (a+d-nc)/(a+b+c+d-nc)
end

repo: gspbox


text:
function [newbnd] = mesh2edge(bnd)

% MESH2EDGE finds the edge lines from a triangulated mesh or the edge
% surfaces from a tetrahedral or hexahedral mesh. An edge is defined as an
% element that does not border any other element. This also implies that a
% closed triangulated surface has no edges.
%
% Use as
%   [edge] = mesh2edge(mesh)
%
% See also POLY2TRI

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

if isfield(bnd, 'tri')
  % make a list of all edges
  edge1 = bnd.tri(:, [1 2]);
  edge2 = bnd.tri(:, [2 3]);
  edge3 = bnd.tri(:, [3 1]);
  edge = cat(1, edge1, edge2, edge3);

elseif isfield(bnd, 'tet')
  % make a list of all triangles that form the tetraheder
  tri1 = bnd.tet(:, [1 2 3]);
  tri2 = bnd.tet(:, [2 3 4]);
  tri3 = bnd.tet(:, [3 4 1]);
  tri4 = bnd.tet(:, [4 1 2]);
  edge = cat(1, tri1, tri2, tri3, tri4);

elseif isfield(bnd, 'hex')
  % make a list of all "squares" that form the cube/hexaheder
  % FIXME should be checked, this is impossible without a drawing
  square1 = bnd.hex(:, [1 2 3 4]);
  square2 = bnd.hex(:, [5 6 7 8]);
  square3 = bnd.hex(:, [1 2 6 5]);
  square4 = bnd.hex(:, [2 3 7 6]);
  square5 = bnd.hex(:, [3 4 8 7]);
  square6 = bnd.hex(:, [4 1 5 8]);
  edge = cat(1, square1, square2, square3, square4, square5, square6);

end % isfield(bnd)

% soort all polygons in the same direction
% keep the original as "edge" and the sorted one as "sedge"
sedge = sort(edge, 2);

% % find the edges that are not shared -> count the number of occurences
% n = size(sedge,1);
% occurences = ones(n,1);
% for i=1:n
%   for j=(i+1):n
%     if all(sedge(i,:)==sedge(j,:))
%       occurences(i) = occurences(i)+1;
%       occurences(j) = occurences(j)+1;
%     end
%   end
% end
%
% % make the selection in the original, not the sorted version of the edges
% % otherwise the orientation of the edges might get flipped
% edge = edge(occurences==1,:);

% find the edges that are not shared
indx = findsingleoccurringrows(sedge);
edge = edge(indx, :);

% replace pnt by pos
bnd = fixpos(bnd);

% the naming of the output edges depends on what they represent
newbnd.pos  = bnd.pos;
if isfield(bnd, 'tri')
  % these have two vertices in each edge element
  newbnd.line = edge;
elseif isfield(bnd, 'tet')
  % these have three vertices in each edge element
  newbnd.tri = edge;
elseif isfield(bnd, 'hex')
  % these have four vertices in each edge element
  newbnd.poly = edge;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1833#c12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function indx = findsingleoccurringrows(X)
[X, indx] = sortrows(X);
sel  = any(diff([X(1,:)-1; X],1),2) & any(diff([X; X(end,:)+1],1),2);
indx = indx(sel);

function indx = finduniquerows(X)
[X, indx] = sortrows(X);
sel  = any(diff([X(1,:)-1; X],1),2);
indx = indx(sel);

author: spm


text:
function test10
%TEST10 test cs_qr
%
% Example:
%   test10
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

rand ('state', 0) ;

% f = 185 ;
% f = 449 ;
clf

for trials = 1:100

    m = fix (100 * rand (1)) ;
    n = fix (100 * rand (1)) ;
    d = 0.1 * rand (1) ;
    A = sprandn (m, n, d) ;
    [m n] = size (A) ;
    if (m < n)
        A = A' ;
    end
    [m n] = size (A) ;
    sp = sprank (A) ;
    % if (sp < n)
    %   continue ;
    % end

    Aorig = A ;

    % A = A (:, colamd (A)) ;

    tic ;
    R = qr (A) ;
    t1 = toc ;

    % tic ;
    % [Q,R] = qr (A) ;
    % t1 = toc ;

    [c,h,parent] = symbfact (A, 'col') ;                                    %#ok
    rnz = sum (c) ;                                                         %#ok
    tic ;
    [V2,Beta2,p,R2] = cs_qr (sparse(A)) ;
    t2 = toc ;

    C = A ;
    m2 = size (V2,1) ;
    if (m2 > m)
        C = [A ; sparse(m2-m, n)] ;
    end
    C = C (p,:) ;

    [H1,R1] = myqr (C) ;
    err1 = norm (R1-R2,1) / norm (R1) ;
    disp ('err1 = ') ;
    disp (err1) ;
    % [svd(A) svd(R1) svd(full(R2))]
    s1 = svd (full (A)) ;
    s2 = svd (full (R2)) ;
    if (n > 0)
        err2 = norm (s1 - s2) / s1 (1)  ;
        disp ('err2 = ') ;
        disp (err2) ;
    else
        err2 = 0 ;
    end
    fprintf ('%10.6f %10.6f  cs speedup %8.3f sprank %d vs %d\n', t1, t2, t1/t2, sp, n) ;

    % H2 = full (H2)
    % R2 = full (R2)

    subplot (2,4,1) ; spy (A) ;         title ('A colamd') ;
    subplot (2,4,4) ; spy (Aorig) ;     title ('Aorig') ;
    subplot (2,4,2) ; spy (C) ;         title ('A rperm') ;
    subplot (2,4,5) ; spy (abs(R2)>0) ; title ('spqr R, no zeros') ;
    subplot (2,4,6) ; spy (R) ;         title ('matlab R') ;
    subplot (2,4,7) ; spy (R2) ;        title ('spqr R') ;
    subplot (2,4,8) ; spy (V2) ;        title ('spqr H') ;
    drawnow

    if (err2 > 1e-9)
        error ('!') ;
    end

    if (m2 > m)
        fprintf ('added %d rows, sprank %d n %d\n', m2-m, sp, n) ;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
%  Compute the factorization of a multivariate polynomial
%  by Ruppert matrix, projection, Newton polygon, generalized eigenvalue
%  and numerical GCD
%
%  Syntax:   >> fac = NIF(F,tol)
%
%    Input:    F -- (matrix) multivariate polynomial in coeff. matrices
%                            !!!must be squarefree!!!
%            tol -- (numeric) coefficient error tolerence
%
%    Output: fac -- (cell) irreducible factors of G in coeff. matrices
%
%  Example:
%
%        >> F
%
%        F =
%
%             0     1     1     2
%             0     1     0     1
%             0     0     2     2
%            15    10     3     2
%
%        >> fac = NIF(F,1e-10);
%        >> fac{:}
%
%        ans =
%
%                0             1.0000
%                0             1.0000
%                0                  0
%           0.8320 - 0.0035i   0.5547 - 0.0024i
%
%
%        ans =
%
%                0             1.0000
%                0                  0
%                0             2.0000
%           0.9729 + 0.1223i   0.1946 + 0.0245i
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%krf_logpolar 'log-polar transform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros rf_logpolar.pane file
%
% Parameters:
% InputFile: i 'Input (rect)', required: 'input image, rectangular coordinates'
% OutputFile: o 'Output (log-polar)', required: 'output image, log-polar coordinates'
% Double: rmin 'min', default: 1: 'minimum radius'
% Double: rmax 'max', default: 64: 'maximum radius'
% Double: amin 'min', default: -90: 'minimum angle'
% Double: amax 'max', default: 90: 'maximum angle'
%
% Example: o = krf_logpolar(i, {'i','';'o','';'rmin',1;'rmax',64;'amin',-90;'amax',90})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% rf_logpolar - log-polar transform
%
%  DESCRIPTION
% This is a spacial transform (warping) of an image.
% Specifically, we map the specified section of an
% annulus in the input image to a rectangle in the
% output image.  This is done by mapping annulus pixels
% (m,n) to (rho, theta) pixels, where rho is the log of
% the radius and theta is the angle from m-axis.
%
% The center of the annulus is also specified as the
% upper left corner of the image (0,0), or the center
% (N/2, N/2).  In both cases, angles are measured from
% the m (horizontal) axis, and positive in the
% clockwise direction (to remain consistant with the
% left-hand coordinate system typically used in image
% processing).
%
% The dimensions of the output image are set to either
% the minimum size that will retain the resolution of
% the input image, or that size rounded up to the
% nearest power of 2 depending on the -logradsz and -anglesz flags.
%
%
%
%  EXAMPLES
% See the html tutorial in $REGISTER/examples/html/README_FIRST.html
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
% Cideciyan, Artur V.  "Registration of Ocular Fundus
% Images" IEEE Engineering in Medicine and Biology.
% Jan/Feb 1995, pp. 52-58.
%
%  COPYRIGHT
% Copyright (C) 1996 - 1997, University of New Mexico.  All rights reserved.
%

function varargout = krf_logpolar(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = krf_logpolar(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'rmin', 1;'rmax', 64;'amin', -90;'amax', 90};
maxval={0,0,2,2,360,360};
minval={0,0,2,2,-360,-360};
istoggle=[0,0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Double','Double','Double','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'rf_logpolar"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function lh = lighthandles(axishandle)
%lh = lighthandles(axishandle)
% get all light objects associated with the current axis
% see lightRestoreSingle

lh = findobj(axishandle, 'Type', 'Light');

hans = get(axishandle,'Children');
lh2 = find(strcmp(get(hans,'Type'),'light'));
lh = [lh hans(lh2)];

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [struct_irf_record, D_record, gamma_record]=irfsignrespanel(beta_gibbs,sigma_gibbs,It,Bu,IRFperiods,n,p,m,k,signrestable,signresperiods)

% function [struct_irf_record D_record gamma_record]=irfsignrespanel(sigma_gibbs,irf_record,It,Bu,IRFperiods,n,signrestable,signresperiods,checkalgo,checkiter)
% runs the gibbs sampler to obtain draws from the posterior distribution of
% IRFs, orthogonalised with sign restrictions
% inputs:  - matrix'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)
%          - cell 'irf_record': record of the gibbs sampler draws for the IRFs
%          - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
% outputs: - cell 'struct_irf_record': record of the gibbs sampler draws for the orthogonalised IRFs
%          - matrix 'D_record': record of the gibbs sampler draws for the structural matrix D
%          - matrix 'gamma_record': record of the gibbs sampler draws for the structural disturbances variance-covariance matrix gamma

% this function implements the sign restrictions approach for the panel

% preliminary tasks
% create first the cell that will store the results from the simulations
struct_irf_record=cell(n,n);
% storage cell
storage1=cell(It-Bu,1);
storage2=cell(It-Bu,1);

% now identify all the periods concerned with restrictions
% first expand the non-empty entries in signresperiods since they are only expressed in intervals: transform into list
% for instance, translate [1 4] into [1 2 3 4]; I don't think this can done without a loop
temp=cell2mat(signresperiods(~cellfun(@isempty,signresperiods)));
periods=[];
for ii=1:size(temp,1)
    periods=[periods temp(ii,1):temp(ii,2)];
end
% suppress duplicates and sort
periods=sort(unique(periods))';
% count the total number of restriction periods (required for IRF matrix)
nperiods=size(periods,1);

% Identify the restriction matrices
% create five cells, corresponding to the three possible restrictions:
% one cell for sign restrictions, three cells for magnitude restrictions, one cell for zero restrictions
Scell=cell(1,n);
Mcell=cell(1,n);
Mlcell=cell(1,n);
Mucell=cell(1,n);
Zcell=cell(1,n);

% loop over rows and columns of the period matrix
for ii=1:n
    for jj=1:n
        % if entry (ii,jj) of the period matrix is not empty...
        if ~isempty(signresperiods{ii,jj})
            % ... then there is a restriction over one (or several) periods
            % loop overt those periods
            for kk=signresperiods{ii,jj}(1,1):signresperiods{ii,jj}(1,2)
                % identify the position of the considered period within the list of all periods (required to build the matrix)
                position=find(periods==kk);
                % now create the restriction matrix: this will depend on the type of restriction
                % if it is a positive sign restriction...
                if strcmp(signrestable{ii,jj},'+')
                    % ... then input a 1 entry in the corresponding S matrix
                    Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
                    Scell{1,jj}(end,(position-1)*n+ii)=1;
                    % if it is a negative sign restriction...
                elseif strcmp(signrestable{ii,jj},'-')
                    % ... then input a -1 entry in the corresponding S matrix
                    Scell{1,jj}=[Scell{1,jj};zeros(1,n*nperiods)];
                    Scell{1,jj}(end,(position-1)*n+ii)=-1;
                    % if it is a zero restriction...
                elseif strcmp(signrestable{ii,jj},'0')
                    % ... then input a 1 entry in the corresponding Z matrix
                    Zcell{1,jj}=[Zcell{1,jj};zeros(1,n*nperiods)];
                    Zcell{1,jj}(end,(position-1)*n+ii)=1;
                    % else, a non-empty entry being neither a sign nor a zero restriction has to be a magnitude restriction
                else
                    % fill the corresponding M matrices:
                    % input a 1 in M
                    Mcell{1,jj}=[Mcell{1,jj};zeros(1,n*nperiods)];
                    Mcell{1,jj}(end,(position-1)*n+ii)=1;
                    % input the lower value of the interval in Ml
                    temp=str2num(signrestable{ii,jj});
                    Mlcell{1,jj}=[Mlcell{1,jj};temp(1,1)];
                    % input the upper value of the interval in Mu
                    Mucell{1,jj}=[Mucell{1,jj};temp(1,2)];
                end
            end
        end
    end
end

% now check what kind of restrictions apply among sign, zero and magnitude restrictions
% check for sign restrictions: if there are any, at least one entry in the cell Scell is non-empty
if sum(~cellfun(@isempty,Scell))~=0
    signres=1;
else
    signres=0;
end
% similarly check for zero restrictions
if sum(~cellfun(@isempty,Zcell))~=0
    zerores=1;
else
    zerores=0;
end
% and finally, check for magnitude restrictions
if sum(~cellfun(@isempty,Mcell))~=0
    magnres=1;
else
    magnres=0;
end

hbar = bear.parfor_progressbar(It-Bu,'Progress');  %create the progress bar

% step 1: repeat simulations a number of times equal to the number of simulations retained from Gibbs sampling
parfor ii=1:It-Bu
    % initiate the variable 'success'; this variable will be used to check whether the restrictions are satisfied
    % if there are only zero restrictions, they will be satisfied by construction, and 'success' will simply be ignored
    success=0;
    % how the algorithm will be conducted will depend on the types of restrictions implemented

    % if there are only zero restrictions, the algorithm is simple as no checking is required: the conditions are satisfied by construction
    if zerores==1 && signres==0 && magnres==0
        % draw beta and sigma
        beta=beta_gibbs(:,ii);
        sigma=reshape(sigma_gibbs(:,ii),n,n);
        hsigma=chol(bear.nspd(sigma),'lower');
        % obtain orthogonalised IRFs
        [irfmatrix, ortirfmatrix]=bear.irfsim(beta,hsigma,n,m,p,k,max(IRFperiods,max(periods)));
        % generate the stacked IRF matrix
        stackedirfmat=[];
        for jj=1:numel(periods)
            stackedirfmat=[stackedirfmat;ortirfmatrix(:,:,periods(jj,1)+1)];
        end
        % draw an entire random matrix Q satisfying the zero restrictions
        [Q]=bear.qzerores(n,Zcell,stackedirfmat);
        % there is no need to verify the restrictions: there are satisfied by construction

        % if there are sign/magnitude restrictions, possibly associated with zero restrictions
    else
        % the algorithm becomes a bit more complicated as conditions now need to be checked
        % to maintain efficiency, the algorithm proceeds recursively shock by shock, and stops as soon as a condition on the considered shock fails
        % repeat algorithm for the iteration as long as not all conditions are satisfied
        while success==0
            % switch 'success' to 1; it will be turned back to zero if at any time Q is detected as a candidate not satisfying the restrictions
            success=1;
            % draw randomly the vector of VAR coefficients: draw a random index
            index=floor(rand*(It-Bu))+1;
            % then draw a random set of beta and sigma corresponding to this index (this is done to make it possible to draw, if required, an infinite number of values from the gibbs sampler record, with equal probability on each value)
            beta=beta_gibbs(:,index);
            sigma=reshape(sigma_gibbs(:,index),n,n);
            hsigma=chol(bear.nspd(sigma),'lower');
            % obtain orthogonalised IRFs
            [irfmatrix, ortirfmatrix]=bear.irfsim(beta,hsigma,n,m,p,k,max(IRFperiods,max(periods)));
            % generate the stacked IRF matrix
            stackedirfmat=[];
            for jj=1:numel(periods)
                stackedirfmat=[stackedirfmat;ortirfmatrix(:,:,periods(jj,1)+1)];
            end
            % initiate Qj
            Qj=[];
            % now start looping over the shocks and checking sequentially whether conditions on these shocks hold
            % stop as soon as one restriction fails
            kk=1;
            while success==1 && kk<=n
                % build column j of the random matrix Q
                [qj]=bear.qrandj(n,Zcell{1,kk},stackedirfmat,Qj);
                % obtain the candidate column fj
                fj=stackedirfmat*qj;
                % check restrictions: first sign restrictions
                [success qj]=bear.checksignres(Scell{1,kk},qj,fj);
                % if 'success' is still equal to 1, also check for magnitude restrictions
                if success==1
                    [success]=bear.checkmagres(Mcell{1,kk},Mlcell{1,kk},Mucell{1,kk},fj);
                end
                % also, if 'success' is still equal to 1, update Qj by concatenating qj
                if success==1
                    Qj=[Qj qj];
                end
                kk=kk+1;
            end
            % repeat this loop until a succesful draw is obtained
        end
        % with succesful Qj at hand, eventually set Q as Qj
        Q=Qj;
    end

    % store
    for jj=1:IRFperiods
        storage1{ii,1}(:,:,jj)=ortirfmatrix(:,:,jj)*Q;
    end
    storage2{ii,1}=hsigma*Q;

    hbar.iterate(1);   % update progress by one iteration

end

close(hbar);   %close progress bar

% reorganise storage
% loop over iterations
for ii=1:It-Bu
    % loop over IRF periods
    for jj=1:IRFperiods
        % loop over variables
        for kk=1:n
            % loop over shocks
            for ll=1:n
                struct_irf_record{kk,ll}(ii,jj)=storage1{ii,1}(kk,ll,jj);
            end
        end
    end
    D_record(:,ii)=storage2{ii,1}(:);
    gamma_record(:,ii)=bear.vec(eye(n));
end

author: european-central-bank

repo: BEAR-toolbox


text:
function legendre_associated_normalized_values_test ( )

%*****************************************************************************80
%
%% LEGENDRE_ASSOCIATED_NORMALIZED_VALUES_TEST tests LEGENDRE_ASSOCIATED_NORMALIZED_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_ASSOCIATED_NORMALIZED_VALUES_TEST:\n' );
  fprintf ( 1, '  LEGENDRE_ASSOCIATED_NORMALIZED_VALUES stores values of\n' );
  fprintf ( 1, '  the normalized associated Legendre polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N     M    X             P(N,M)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, m, x, fx ] = legendre_associated_normalized_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %4d  %12f  %24.16f\n', n, m, x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [xc,Vi] = canonical_autocorrelation(TR,n,varargin)
% [xc,Vi] = canonical_autocorrelation(TR,n,[a])
% Tor Wager, 2 / 23 / 04
%
% exponential autocorrelation function of the form xc = e ^ (-ax)
% Optional input "a" is an exponent; higher = more white noise, low =
% colored
%
% TR = sampling rate (time per data point)
% n = number of samples in timeseries
%
% xc = vector containing autocorrelation coefficients
% Vi = autocorrelation matrix

f = inline('A * exp(-a * x)','A','a','x');  % exponential function
xa = 0:TR:60;

if length(varargin) > 0
    xc = f(1,varargin{1},xa);
else
    xc = f(.9475,.5323,xa);    % params at 3T from vnl experiment, n = 10
    xc = xc ./ max(xc);
end

% autocorrelation
Vi = getV('make',xc,n);

return


text:
function f = toFunctionOut(disc, coeffs, cutoff)
%TOFUNCTIONOUT   Convert discrete values of an ULTRAS discretization to a
%   CHEBFUN.
%
%   F = TOFUNCTIONOUT(DISC, COEFFS, CUTOFF) converts the coeffs returned by
%   ULTRAS to a CHEBFUN. The input may be piecewise smooth, as indicated by the
%   dimension property of the discretization.

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

% Check for cutoff
if ( nargin == 3 )
    m = cutoff;
else
    m = inf;
end

dom = disc.domain;         % Domain we're working on, including breakpoints
c = mat2cell(full(coeffs), disc.dimension); % Break into smooth pieces
funs = cell(numel(c), 1);  % One FUN per piece

% adjust size of cutoff if necessary
if ( length(m) ~= numel(c) )
    m = max(m)*ones(numel(c), 1);
end

for k = 1:numel(c)
    % Construct CHEBTECH2 objects from each piece:
    coeffs = c{k};
    c{k} = coeffs(1:min(m(k), size(coeffs, 1)), :);
    ct = chebtech2({[], c{k}});
    % Assign each piece to a subinterval with a BNDFUN:
    funs{k} = bndfun(ct, struct('domain', dom(k: k + 1)));
end
% Conver the FUNS cell-array to a CHEBFUN.
f = chebfun(funs);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function M = sympositivedefinitefactory(n)
% Manifold of n-by-n symmetric positive definite matrices with
% the bi-invariant geometry.
%
% function M = sympositivedefinitefactory(n)
%
% A point X on the manifold is represented as a symmetric positive definite
% matrix X (nxn). Tangent vectors are symmetric matrices of the same size
% (but not necessarily definite).
%
% The Riemannian metric is the bi-invariant metric, described notably in
% Chapter 6 of the 2007 book "Positive definite matrices"
% by Rajendra Bhatia, Princeton University Press.
%
%
% The retraction / exponential map involves expm (the matrix exponential).
% If too large a vector is retracted / exponentiated (e.g., a solver tries
% to make too big a step), this may result in NaN's in the returned point,
% which most likely would lead to NaN's in the cost / gradient / ... and
% will result in failure of the optimization. For trustregions, this can be
% controlled by setting options.Delta0 and options.Delta_bar, to prevent
% too large steps.
%
%
% Note also that many of the functions involve solving linear systems in X
% (a point on the manifold), taking matrix exponential and logarithms, etc.
% It could therefore be beneficial to do some precomputation on X (an
% eigenvalue decomposition for example) and store both X and the
% preprocessing in a structure. This would require modifying the present
% factory to work with such structures to represent both points and tangent
% vectors. We omit this in favor of simplicity, but it may be good to keep
% this in mind if efficiency becomes an issue in your application.

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, August 29, 2013.
% Contributors: Nicolas Boumal
% Change log:
%
%   March 5, 2014 (NB)
%       There were a number of mistakes in the code owing to the tacit
%       assumption that if X and eta are symmetric, then X\eta is
%       symmetric too, which is not the case. See discussion on the Manopt
%       forum started on Jan. 19, 2014. Functions norm, dist, exp and log
%       were modified accordingly. Furthermore, they only require matrix
%       inversion (as well as matrix log or matrix exp), not matrix square
%       roots or their inverse.
%
%   July 28, 2014 (NB)
%       The dim() function returned n*(n-1)/2 instead of n*(n+1)/2.
%       Implemented proper parallel transport from Sra and Hosseini (not
%       used by default).
%       Also added symmetrization in exp and log (to be sure).
%
%   April 3, 2015 (NB):
%       Replaced trace(A*B) by a faster equivalent that does not compute
%       the whole product A*B, for inner product, norm and distance.
%
%   May 23, 2017 (NB):
%       As seen in a talk of Wen Huang at the SIAM Optimization Conference
%       today, replaced the retraction of this factory (which was simply
%       equal to the exponential map) with a simpler, second-order
%       retraction. That this retraction is second order can be verified
%       numerically with checkretraction(sympositivedefinitefactory(5));
%       Notice that, for this retraction, it would be cheap to evaluate for
%       many values of t, that is, it is cheap to retract many points along
%       the same tangent direction. This could in principle be exploited to
%       speed up line-searches.
%
%   Jan. 12, 2022 (NB):
%       Simplified code for ehess2rhess by commenting out the reasoning and
%       computing only the end result.

    symm = @(X) .5*(X+X');

    M.name = @() sprintf('Symmetric positive definite geometry of %dx%d matrices', n, n);

    M.dim = @() n*(n+1)/2;

    % Helpers to avoid computing full matrices simply to extract their trace
    vec  = @(A) A(:);
    trAB = @(A, B) vec(A')'*vec(B);  % = trace(A*B)
    trAA = @(A) sqrt(trAB(A, A));    % = sqrt(trace(A^2))

    % Choice of the metric on the orthonormal space is motivated by the
    % symmetry present in the space. The metric on the positive definite
    % cone is its natural bi-invariant metric.
    % The result is equal to: trace( (X\eta) * (X\zeta) )
    M.inner = @(X, eta, zeta) trAB(X\eta, X\zeta);

    % Notice that X\eta is *not* symmetric in general.
    % The result is equal to: sqrt(trace((X\eta)^2)).
    % Thus, we compute the sum of the squared eigenvalues of X\eta, which
    % in general is different from summing the squared singular values:
    % that is why it is not equivalent to compute the Frobenius norm here.
    % There should be no need to take the real part, but rounding errors
    % may cause a small imaginary part to appear, so we discard it.
    M.norm = @(X, eta) real(trAA(X\eta));

    % Same here: X\Y is not symmetric in general.
    % Same remark about taking the real part.
    M.dist = @(X, Y) real(trAA(real(logm(X\Y))));

    M.typicaldist = @() sqrt(n*(n+1)/2);

    M.egrad2rgrad = @egrad2rgrad;
    function eta = egrad2rgrad(X, eta)
        eta = X*symm(eta)*X;
    end

    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)
        % Directional derivatives of the Riemannian gradient
        %  Hess = X*symm(ehess)*X + 2*symm(eta*symm(egrad)*X);
        % Correction factor for the non-constant metric
        %  Hess = Hess - symm(eta*symm(egrad)*X);
        % Combined:
        Hess = X*symm(ehess)*X + symm(eta*symm(egrad)*X);
    end

    M.proj = @(X, eta) symm(eta);

    M.tangent = M.proj;
    M.tangent2ambient = @(X, eta) eta;

    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            teta = eta;
        else
            teta = t*eta;
        end
        % The symm() call is mathematically unnecessary but numerically
        % necessary.
        Y = symm(X + teta + .5*teta*(X\teta));
    end

    M.exp = @exponential;
    function Y = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        % The symm() and real() calls are mathematically not necessary but
        % are numerically necessary.
        Y = symm(X*real(expm(X\(t*eta))));
    end

    M.log = @logarithm;
    function H = logarithm(X, Y)
        % Same remark regarding the calls to symm() and real().
        H = symm(X*real(logm(X\Y)));
    end

    M.hash = @(X) ['z' hashmd5(X(:))];

    % Generate a random symmetric positive definite matrix following a
    % certain distribution. The particular choice of a distribution is of
    % course arbitrary, and specific applications might require different
    % ones.
    M.rand = @random;
    function X = random()
        D = diag(1+rand(n, 1));
        [Q, R] = qr(randn(n)); %#ok
        X = Q*D*Q';
    end

    % Generate a uniformly random unit-norm tangent vector at X.
    M.randvec = @randomvec;
    function eta = randomvec(X)
        eta = symm(randn(n));
        nrm = M.norm(X, eta);
        eta = eta / nrm;
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(X) zeros(n);

    % Poor man's vector transport: exploit the fact that all tangent spaces
    % are the set of symmetric matrices, so that the identity is a sort of
    % vector transport. It may perform poorly if the origin and target (X1
    % and X2) are far apart though. This should not be the case for typical
    % optimization algorithms, which perform small steps.
    M.transp = @(X1, X2, eta) eta;

    % For reference, a proper vector transport is given here, following
    % work by Sra and Hosseini: "Conic geometric optimisation on the
    % manifold of positive definite matrices", in SIAM J. Optim.
    % in 2015; also available here: http://arxiv.org/abs/1312.1039
    % This will not be used by default. To force the use of this transport,
    % execute "M.transp = M.paralleltransp;" on your M returned by the
    % present factory.
    M.paralleltransp = @parallel_transport;
    function zeta = parallel_transport(X, Y, eta)
        E = sqrtm(Y/X);
        zeta = E*eta*E';
    end

    % vec and mat are not isometries, because of the unusual inner metric.
    M.vec = @(X, U) U(:);
    M.mat = @(X, u) reshape(u, n, n);
    M.vecmatareisometries = @() false;

end

author: NicolasBoumal

repo: manopt


text:
classdef myPoint
    %MYPOINT Summary of this class goes here
    %   Detailed explanation goes here

    properties (GetAccess='public', SetAccess='public')
    x;
    y;
    end

    methods
        function obj = myPoint(x,y)
           obj.x = x;
           obj.y = y;
        end
    end

end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function test_int_2d_test02 ( )

%*****************************************************************************80
%
%% TEST_INT_2D_TEST02 applies a product of composite midpoint rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_INT_2D_TEST02\n' );
  fprintf ( 1, '  Apply a product of composite midpoint rules.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeatedly multiply the number of points by 16.\n' );

  problem_num = p00_problem_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1,'   Problem      Points      Approx          Error\n' );

  for problem = 1 : problem_num

    fprintf ( 1, '\n' );

    nx = 1;
    ny = 1;

    for i = 1 : 12

      n = nx * ny;

      x = zeros ( 2, n );

      [ a, b ] = p00_lim ( problem );

      k = 0;

      for ix = 1 : nx

        xval = ( ( 2 * nx - 2 * ix + 1 ) * a(1)   ...
               + (          2 * ix - 1 ) * b(1) ) ...
               / ( 2 * nx              );

        for iy = 1 : ny

          yval = ( ( 2 * ny - 2 * iy + 1 ) * a(2)   ...
                 + (          2 * iy - 1 ) * b(2) ) ...
                 / ( 2 * ny              );

          k = k + 1;
          x(1,k) = xval;
          x(2,k) = yval;

        end

      end

      volume = prod ( b(1:2) - a(1:2) );

      fx = p00_fun ( problem, n, x );

      quad = volume * sum ( fx(1:n) ) / n;

      exact = p00_exact ( problem );

      error = abs ( quad - exact );

      fprintf ( 1, '  %8d  %10d  %14e  %14e\n', ...
        problem, n, quad, error );

      nx = nx * 2;
      ny = ny * 2;

    end

    fprintf ( 1, '  %8d       Exact  %14e\n', ...
      problem, exact );

  end

  return
end


text:
function [atlasA,atlasE,atlasMask]=atlasAjustHoles(atlasA,atlasE,atlasMask)
% [atlasA,atlasE,atlasMask]=atlasAjustHoles(atlasA,atlasE,atlasMask)
% Author: Schira
% Purpose:
%    Clears little holes in an complete "Atlas"
%    basicly this is just a cosmetic operation, assuming any -1 point
%    completly surrounded by "Atlas" points is a "hole"
%    The tool fills in the mean for eccentricity and angle and the median
%    for the mask.

map=zeros(size(atlasE));
ind=find(atlasE~=-1);
map(ind)=1;
map=imdilate(map,[1 1])-map;
ind=find(map==1);
for ii=1:length(ind)
    i=ind(ii);
    b=floor(i/size(atlasE,1));
    a=i-size(atlasE,1)*b;
    b=b+1;
    vekt=cat(2,atlasA(a-1,b),atlasA(a+1,b),atlasA(a,b-1),atlasA(a,b-1));
    if sum(vekt==-1)==0
        atlasA(i)=mean([atlasA(a-1,b),atlasA(a+1,b),atlasA(a,b-1),atlasA(a,b-1)]);
        atlasE(i)=mean([atlasE(a-1,b),atlasE(a+1,b),atlasE(a,b-1),atlasE(a,b-1)]);
        atlasMask(i)=median([atlasMask(a-1,b),atlasMask(a+1,b),atlasMask(a,b-1),atlasMask(a,b-1)]);
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ x, y, z, w ] = ld1730 ( )

%*****************************************************************************80
%
%% LD1730 computes the 1730 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(1730,1);
  y = zeros(1730,1);
  z = zeros(1730,1);
  w = zeros(1730,1);
  a = 0.0;
  b = 0.0;
  v = 0.6309049437420976E-04;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.6398287705571748E-03;
  [ n, x, y, z, w ] = gen_oh ( 2, n, a, b, v, x, y, z, w );
  v = 0.6357185073530720E-03;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.2860923126194662E-01;
  v = 0.2221207162188168E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.7142556767711522E-01;
  v = 0.3475784022286848E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1209199540995559;
  v = 0.4350742443589804E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1738673106594379;
  v = 0.4978569136522127E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2284645438467734;
  v = 0.5435036221998053E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2834807671701512;
  v = 0.5765913388219542E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3379680145467339;
  v = 0.6001200359226003E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3911355454819537;
  v = 0.6162178172717512E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4422860353001403;
  v = 0.6265218152438485E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4907781568726057;
  v = 0.6323987160974212E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5360006153211468;
  v = 0.6350767851540569E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6142105973596603;
  v = 0.6354362775297107E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6459300387977504;
  v = 0.6352302462706235E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6718056125089225;
  v = 0.6358117881417972E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6910888533186254;
  v = 0.6373101590310117E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.7030467416823252;
  v = 0.6390428961368665E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.8354951166354646E-01;
  v = 0.3186913449946576E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.2050143009099486;
  v = 0.4678028558591711E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.3370208290706637;
  v = 0.5538829697598626E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.4689051484233963;
  v = 0.6044475907190476E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.5939400424557334;
  v = 0.6313575103509012E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.1394983311832261;
  b = 0.4097581162050343E-01;
  v = 0.4078626431855630E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1967999180485014;
  b = 0.8851987391293348E-01;
  v = 0.4759933057812725E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2546183732548967;
  b = 0.1397680182969819;
  v = 0.5268151186413440E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3121281074713875;
  b = 0.1929452542226526;
  v = 0.5643048560507316E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3685981078502492;
  b = 0.2467898337061562;
  v = 0.5914501076613073E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4233760321547856;
  b = 0.3003104124785409;
  v = 0.6104561257874195E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4758671236059246;
  b = 0.3526684328175033;
  v = 0.6230252860707806E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5255178579796463;
  b = 0.4031134861145713;
  v = 0.6305618761760796E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5718025633734589;
  b = 0.4509426448342351;
  v = 0.6343092767597889E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2686927772723415;
  b = 0.4711322502423248E-01;
  v = 0.5176268945737826E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3306006819904809;
  b = 0.9784487303942695E-01;
  v = 0.5564840313313692E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3904906850594983;
  b = 0.1505395810025273;
  v = 0.5856426671038980E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4479957951904390;
  b = 0.2039728156296050;
  v = 0.6066386925777091E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5027076848919780;
  b = 0.2571529941121107;
  v = 0.6208824962234458E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5542087392260217;
  b = 0.3092191375815670;
  v = 0.6296314297822907E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6020850887375187;
  b = 0.3593807506130276;
  v = 0.6340423756791859E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4019851409179594;
  b = 0.5063389934378671E-01;
  v = 0.5829627677107342E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4635614567449800;
  b = 0.1032422269160612;
  v = 0.6048693376081110E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5215860931591575;
  b = 0.1566322094006254;
  v = 0.6202362317732461E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5758202499099271;
  b = 0.2098082827491099;
  v = 0.6299005328403779E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6259893683876795;
  b = 0.2618824114553391;
  v = 0.6347722390609353E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5313795124811891;
  b = 0.5263245019338556E-01;
  v = 0.6203778981238834E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5893317955931995;
  b = 0.1061059730982005;
  v = 0.6308414671239979E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6426246321215801;
  b = 0.1594171564034221;
  v = 0.6362706466959498E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6511904367376113;
  b = 0.5354789536565540E-01;
  v = 0.6375414170333233E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end


text:
classdef Stress3DVoigtTensor < SecondOrder3DVoigtTensor ...
                               & StressDescriptor ...

    methods (Access = public)

        function obj = Stress3DVoigtTensor()
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function metaboliteMassBalancePlot(model, metAbbr, solution, N)
% Plots the top N reactions producing and consuming a metabolite in a flux solution
%
% USAGE:
%
%    metaboliteMassBalancePlot(model, metAbbr, solution, N)
%
% INPUTS:
%    model:       COBRA model structure
%    metAbbr:     metabolite abbreviation
%    solution:    `solveCobraLP` output of a solution to FBA problem
%    N:           Number of reactions to include for production/consumption
%
% .. Author: - Ronan M.T. Fleming

[nMet,nRxn]=size(model.S);

metInd=find(strcmp(metAbbr,model.mets));
if isempty(metInd)
    error('metabolite abbreviation not found')
end

tol = 1e-6;

negS=(model.S(metInd,:)<0)';
posS=(model.S(metInd,:)>0)';

v=solution.full;

consumeBool=( negS & v>tol ) | ( posS & v<-tol );
produceBool=( posS & v>tol ) | ( negS & v<-tol );

vConsume=zeros(nRxn,1);
vProduce=zeros(nRxn,1);

vConsume(consumeBool)=v(consumeBool);
vProduce(produceBool)=v(produceBool);

[sortedConsume,XIMax]=sort(abs(vConsume),'descend');
NConsume=nnz(sortedConsume~=0);

[sortedProduce,XIMin]=sort(abs(vProduce),'descend');
NProduce=nnz(sortedProduce~=0);

figure
h1=subplot(2,1,1);
NConsume=min([NConsume,N]);
YI=(1:NConsume)';
barh(YI,v(XIMax(1:NConsume)));
% title(['Consumers (top) and Producers of ' metAbbr],'FontSize',16)

if NConsume==0
    set(h1,'YTickLabel','NONE','FontSize',8);
else
    for n=1:NConsume
        if strcmp(model.rxn(XIMax(n)).directionality,'forward')
            YTickLabelMax{n}=['*' model.rxns{XIMax(n)}];
        else
            YTickLabelMax{n}=model.rxns{XIMax(n)};
        end
    end
    set(h1,'Ytick', 1:NConsume)
    set(h1,'YTickLabel',YTickLabelMax,'FontSize',8);
    ylim([0.5 NConsume+0.5])
end
xlabel(['Consumer of ' metAbbr],'FontSize',12)

h2=subplot(2,1,2);
NProduce=min([NProduce,N]);
YI=(1:NProduce)';
barh(YI,v(XIMin(1:NProduce)));
if NProduce==0
    set(h1,'YTickLabel','NONE','FontSize',8);
else
    for n=1:NProduce
        if strcmp(model.rxn(XIMin(n)).directionality,'forward')
            YTickLabelMin{n}=['*' model.rxns{XIMin(n)}];
        else
            YTickLabelMin{n}=model.rxns{XIMin(n)};
        end
    end
    set(h2,'Ytick', 1:NProduce)
    set(h2,'YTickLabel',YTickLabelMin,'FontSize',8);
    ylim([0.5 NProduce+0.5])
end
xlabel(['Producer of ' metAbbr],'FontSize',12)

author: opencobra

repo: cobratoolbox


text:
%% housekeeping
clear all
close all
clc
%% RISE the model
m=rise('ireland2004');

%% load the data
mydat=load('ych.dat');
data=ts('1948Q1',log(mydat),{'LY','LC','LH'});
data=pages2struct(data);
data.TREND=ts(data.LY.start,(1:data.LY.NumberOfObservations).');
%%
ms=estimate(m,'data',data,'kf_tol',0);
%% Redo filtering?

mfilt=filter(ms,'data',data,'kf_tol',0);

author: jmaih


text:
function odf = mix2
% mix2 sample ODF

CS = crystalSymmetry('cubic');
SS = specimenSymmetry('222');

psi = SO3vonMisesFisherKernel('halfwidth',17*degree);
ori1 = orientation.byEuler(54.736*degree, 45.0*degree, 0.0*degree,CS,SS);
ori2 = orientation.byEuler(62.968*degree, 57.689*degree, 71.565*degree,CS,SS);
ori3 = orientation.byEuler(50.768*degree, 65.905*degree, 63.435*degree,CS,SS);

odf = 0.3095*uniformODF(CS,SS) + 0.315*unimodalODF(ori1,psi) + ...
  0.315*unimodalODF(ori2,psi) + 0.0605*unimodalODF(ori3,psi);

% 3,7
% MIX2
% 3,3,0
% 0.3095
% 54.736,45.0,0.0
% 1,17.0,0.315,2
% 62.968,57.689,71.565
% 1,17.0,0.315,2
% 50.768,65.905,63.435
% 1,17.0,0.0605,1

author: mtex-toolbox

repo: mtex


text:
function month_cal_republican ( y, m )

%*****************************************************************************80
%
%% MONTH_CAL_REPUBLICAN prints a Republican month calendar.
%
%  Format:
%
%    REPUBLICAN CALENDAR
%    Brumaire 3 ER
%
%     1  2  3  4  5  6  7  8  9 10
%    11 12 13 14 15 16 17 18 19 20
%    21 22 23 24 25 26 27 28 29 30
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, the YM date.
%

%
%  Make local copies of the input.
%
  m2 = m;
  y2 = y;
%
%  Check the month and year.
%
  [ y2, m2, ierror ] = ym_check_republican ( y2, m2 );

  if ( ierror ~= 0 )
    return
  end
%
%  Print out a heading.
%
  s1 = month_to_month_name_republican ( m2 );
  s2 = y_to_s_republican ( y2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Republican Calendar\n' );
  fprintf ( 1, '%s %s\n', s1, s2 );
  fprintf ( 1, '\n' );
%
%  Get the days of the week.
%
  for w = 1 : 10
    fprintf ( 1, '%3d', w );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
%
%  Print out a line of day numbers.
%
  iday = 1;

  while ( iday <= month_length_republican ( y2, m2 ) )

    for w = 1 : 10

      if ( month_length_republican ( y2, m2 ) < iday )
        fprintf ( 1, '   ' );
      else
        fprintf ( 1, '%3d', iday );
      end

      iday = iday + 1;

    end

    fprintf ( 1, '\n' );

  end

  return
end