text:
function z = sppurge(z, options)
% SPPURGE   Purge sparse grid data
%    Z = SPPURGE(Z)  Marks indices that have corresp.
%    hierarchical surplus values larger than the default
%    drop tolerance [0, 100*eps]. The SPPURGE function
%    returns the same sparse grid interpolant data z, but
%    enhanced by a field purgeData that is used by
%    SPINTERP to only consider the marked indices in the
%    interpolation process, thus saving computing time.
%
%    Z = SPPURGE(Z, OPTIONS)  OPTIONS must be an options
%    structure generated with SPSET. Only the value of
%    the DropTol property is used, which enables the user to
%    set any absolute and relative drop tolerance to be used
%    by the purging algorithm.
%
%    See also: SPSET

% Author : Andreas Klimke
% Version: 1.1
% Date   : December 4, 2007

% Change log:
% V1.0   : February 1, 2006
%          Initial version
% V1.1   : December 4, 2007
%          Added Gauss-Patterson grid.

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

if nargin > 1
  % New user-demanded drop tolerance
  dropTol = spget(options, 'DropTol', 'auto');
  % Check if purge data is already present
  if isfield(z, 'dropTol')
    if dropTol == z.dropTol
      % do nothing
    else
      % Remove previous data; it is now invalid.
      if isfield(z,'purgeData')
        z = rmfield(z,'purgeData');
      end
    end
  end
elseif isfield(z, 'dropTol')
  % Update of previously computed sppurge data
  dropTol = z.dropTol;
else
  % Use default drop tolerance.
  dropTol = 'auto';
end
if isa(dropTol, 'char')
	if strcmpi(dropTol, 'off')
		dropTol = [0, 0];
	elseif strcmpi(dropTol, 'auto')
		dropTol = [0, 100*eps];
	else
		error('MATLAB:spinterp:badopt', ...
					['Unknown string value to DropTol property. Did you ' ...
           'accidentally quote the 1x2 drop tolerance vector?']);
	end
elseif isa(dropTol, 'double')
	if length(dropTol) < 2
		error('MATLAB:spinterp:badopt', ...
					['DropTol must be either a string or a 1x2 double ' ...
					 'vector. See ''help spset'' for details.']);
	end
else
	error('MATLAB:spinterp:badopt', ...
				['DropTol must be either a string or a 1x2 double ' ...
				 'vector. See ''help spset'' for details.']);
end

absdroptol = dropTol(1);
reldroptol = dropTol(2);

% Do nothing if drop tolerance is less than or equal to zero
if absdroptol <= 0 && reldroptol <= 0
 % Remove purge data if present
 if isfield(z,'purgeData')
   z = rmfield(z,'purgeData');
   z = rmfield(z,'dropTol');
 end
 return;
end

gridtype = z.gridType;
d = z.d;

if isfield(z, 'indices');
	sparseIndices = 'on';
	indices = z.indices;
else
	sparseIndices = 'off';
end

% Set the currently used drop tolerance
z.dropTol = dropTol;

nout = size(z.vals,1);

if strcmpi(sparseIndices, 'off')
	switch lower(gridtype)
	 case 'clenshaw-curtis'
		getpointsmethod = @spgetnpointscc;
	 case 'maximum'
		getpointsmethod = @spgetnpointsm;
	 case 'noboundary'
		getpointsmethod = @spgetnpointsnb;
	 case 'chebyshev'
		getpointsmethod = @spgetnpointscc;
	 case 'gauss-patterson'
		getpointsmethod = @spgetnpointsnb;
	 otherwise
		error('MATLAB:spinterp:badopt', ...
					['Unknown grid type ''' gridtype '''.']);
	end
	if isfield(z, 'purgeData')
		nfrom = size(z.purgeData,2);
	else
		nfrom = 1;
		z.purgeData{1} = ones(1,nout,'uint8');
		% Omit level 0, since it is clear that there is some data
		% to be considered unless the objective function is f(x) = 0
	end
	nto = size(z.vals,2) - 1;
	% Do this with respect to all outupts
	for k = nfrom:nto
		seq = spgetseq(k,d);
		seqlength = size(seq,1);
		[totalpoints, npoints] = feval(getpointsmethod,seq);
		purgedata = zeros(seqlength,nout,'uint8');
		for l = 1:nout
			droptol = max((z.fevalRange(l,2)- ...
				z.fevalRange(l,1))*reldroptol, absdroptol);
			vals = z.vals{l,k+1};
			zid = uint32(1);
			for m = 1:seqlength
				endzid = zid + npoints(m);
				while zid < endzid
					if abs(vals(zid)) >= droptol
						purgedata(m,l) = 1;
						zid = endzid;
						break;
					else
						zid = zid + 1;
					end
				end
			end
		end
		z.purgeData{k+1} = purgedata;
	end
else
	npoints = z.indices.subGridPoints;
	addr    = z.indices.subGridAddr;
	if isfield(z, 'purgeData')
		idfrom = uint32(size(z.purgeData,1))+1;
	else
		idfrom = uint32(2);
		z.purgeData = ones(1,nout,'uint8');
		% Omit first index set, since it is clear that there is some data
		% to be considered unless the objective function is f(x) = 0
	end
	idto = uint32(size(npoints,1));
	if idto >= idfrom
		purgedata = [z.purgeData; zeros(idto-idfrom+1,nout,'uint8')];
		for l = 1:nout
			droptol = max((z.fevalRange(l,2)- ...
			z.fevalRange(l,1))*reldroptol, absdroptol);
			vals = z.vals{l};
			for m = idfrom:idto
				zid = addr(m);
				endzid = addr(m) + npoints(m);
				while zid < endzid
					if abs(vals(zid)) >= droptol
						purgedata(m,l) = 1;
						break;
					else
						zid = zid + 1;
					end
				end
			end
		end
		z.purgeData = purgedata;
	end
end

author: johannesgerer


text:
function [cameraCBparameters] = RecalculateCBcalibrationParameters(cameraCBparameters,optStruct)
%% function for re-calculating the distortion parameters in STEP0.
% This function is called only in case the user selected a repeated
% analysis.
%
% INPUTS:
% * cameraCBparameters: structure previously created in STEP0
% * optStruct: a structure containing the new parameters of the distortion
%   model selected in the repeated analysis.
%
% OUTPUTS:
% * a figure of all the checkerboard images

%%

%  extract the points already detected in previous analysis
imagePoints=cameraCBparameters.imagePoints;
worldPoints=cameraCBparameters.cameraParameters.WorldPoints;

% Re-calculate camera parameters with the new distortion model
[params,~,estimationErrors] = estimateCameraParameters(imagePoints,worldPoints,...
    'NumRadialDistortionCoefficients',optStruct.NumRadialDistortionCoefficients,'EstimateTangentialDistortion',optStruct.EstimateTangentialDistortion,'EstimateSkew',optStruct.EstimateSkew);

% feed new results into cameraCBparameters
cameraCBparameters.cameraParameters=params;
cameraCBparameters.estimationErrors=estimationErrors;

%parameters after undistortion
imagePointsUndistorted=zeros(size(cameraCBparameters.imagePoints));
for ii=1:size(cameraCBparameters.imagePoints,3)
    imagePointNow=cameraCBparameters.imagePoints(:,:,ii);
    [imagePointsUndistorted(:,:,ii)] = undistortPoints(imagePointNow,cameraCBparameters.cameraParameters);
end
[paramsJ,~,estimationErrorsJ] = estimateCameraParameters(imagePointsUndistorted,worldPoints,...
    'NumRadialDistortionCoefficients',optStruct.NumRadialDistortionCoefficients,'EstimateTangentialDistortion',optStruct.EstimateTangentialDistortion,'EstimateSkew',optStruct.EstimateSkew);

% feed results after undistortion into cameraCBparameters and return
cameraCBparameters.cameraParametersAUD=paramsJ;
cameraCBparameters.estimationErrorsAUD=estimationErrorsJ;
cameraCBparameters.imagePointsAUD=imagePointsUndistorted;

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function cooccurM = calcCooccur(quantizedM, offsetsM, nL, cooccurType)
% function cooccurM = calcCooccur(quantizedM, offsetsM, nL, cooccurType)
%
% This function calculates the cooccurrence matrix for the passed quantized
% image.
%
% INPUTS:
%       quantizedM: quantized 3d matrix obtained, for example, by
%       imquantize_cerr.m
%       offsetsM: Offsets for directionality/neighbors, obtained by
%       getOffsets.m
%       nL: Number of gray levels.
%       cooccurType: flag, 1 or 2.
%                   1: returns a single cooccurrence matrix by combining
%                   contributions from all offsets into one cooccurrence
%                   matrix.
%                   2: returns cooccurM with each column containing
%                   cooccurrence matrix for the row of offsetsM.
% OUTPUT:
%       cooccurM: cooccurrence matrix of size (nL*nL) x 1 for cooccurType =
%       1, or (nL*nL) x size(offsetsM,1) for cooccurType = 2.
%       cooccurM can be passed to cooccurToScalarFeatures.m to get texture
%       features.
%
%
% APA, 05/23/2016

% Default to building cooccurrence by combining all offsets
if ~exist('cooccurType','var')
    cooccurType = 1;
end

% Apply pading of 1 row/col/slc. This assumes offsets are 1. Need to
% parameterize this in case of offsets other than 2. Rarely used for
% medical images.
numColsPad = 1;
numRowsPad = 1;
numSlcsPad = 1;

% Get number of voxels per slice
[numRows, numCols, numSlices] = size(quantizedM);

% Pad quantizedM
if exist('padarray.m','file')
    q = padarray(quantizedM,[numRowsPad numColsPad numSlcsPad],NaN,'both');
else
    q = padarray_oct(quantizedM,[numRowsPad numColsPad numSlcsPad],NaN,'both');
end

% Add level for NaN
lq = nL + 1;
q(isnan(q)) = lq;

q = uint16(q); % q is the quantized image

% Number of offsets
numOffsets = size(offsetsM,1);

% Indices of last level to filter out
nanIndV = false([lq*lq,1]);
nanIndV([lq:lq:lq*lq-lq, lq*lq-lq:lq*lq]) = true;

% Build linear indices column/row-wise for Symmetry
indRowV = zeros(1,lq*lq);
for i=1:lq
    indRowV((i-1)*lq+1:(i-1)*lq+lq) = i:lq:lq*lq;
end

tic
% Initialize cooccurrence matrix (vectorized for speed)
if cooccurType == 1
    %cooccurM = zeros(lq*lq,1,'single');
    cooccurM = sparse(lq*lq,1);
else
    %cooccurM = zeros(lq*lq,numOffsets,'single');
    cooccurM = sparse(lq*lq,numOffsets);
end
for off = 1:numOffsets

    offset = offsetsM(off,:);
    slc1M = q(numRowsPad+(1:numRows),numColsPad+(1:numCols),...
        numSlcsPad+(1:numSlices));
    slc2M = circshift(q,offset);
    slc2M = slc2M(numRowsPad+(1:numRows),numColsPad+(1:numCols),numSlcsPad+(1:numSlices))...
        + (slc1M-1)*lq;
    if cooccurType == 1
        cooccurM = cooccurM + accumarray(slc2M(:),1, [lq*lq,1]);
    else
        cooccurM(:,off) = accumarray(slc2M(:),1, [lq*lq,1]);
    end

end

cooccurM = cooccurM + cooccurM(indRowV,:); % for symmetry
cooccurM(nanIndV,:) = [];
cooccurM = bsxfun(@rdivide,cooccurM,sum(cooccurM,1)+eps);

return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%  Figure 10.71      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig10_71.m is a script to generate Fig. 10.71 the linear RTP response
% LQG method with internal model
% RTP chamber demo Example
clf;
a=[-0.068208813989939 0.014929776357245 0.000000065782442;...
   0.045809672136430 -0.118134773528570 0.021802006306129;...
   0.000000433637498 0.046839194867347 -0.100884399149872];
b3=[0.37873508304055 0.110575403544586 0.022912893467038;...
   0.000000002974222 0.449046982554739 0.073572900808161;...
   0.000000027324116 0.000702342292121 0.417797228783230];
c3=eye(3,3);
d=0*eye(3,3);
sysp=ss(a,b3,c3,d);
[eol]=eig(a);
[zol]=tzero(sysp);
[zol22]=tzero(a,b3(:,2),c3(2,:),d(2,2));
% Combine 3 lamps into a single actuator
b=b3(:,1)+b3(:,2)+b3(:,3);
% Select center temperature
c=c3(2,:);
aa=[zeros(1,1) c;zeros(3,1) a];
bb=[zeros(1,1);b];
%
% weight temperature differences
q1hat=[eye(1,1) zeros(1,3);0 2 -1 -1; 0 -1 2 -1; 0 -1 -1 2]+1e-6*eye(4,4);
q1hat=diag([1 10 10 10])*q1hat;
q2hat=eye(1,1);
[kk,sric,ee]=lqr(aa,bb,q1hat,q2hat);
k1=kk(1,1);
ko=kk(:,2:4);
acl=[a-b*ko b;-k1*c zeros(1,1)];
bcl=[zeros(3,1);k1];
ccl=[c zeros(1,1)];
dcl=zeros(1,1);
[ecl]=eig(acl);
[zcl]=tzero(acl,bcl,ccl,dcl);
%CL DC gain
[cldcgain]=dcl-ccl*inv(acl)*bcl;
%CL Step Response
%[y,t]=step(acl,bcl,ccl,dcl);
%s1y=plot(y(:,:,1));
%grid;
%pause;
%Control effort
%cclu=[-ko eye(1,1)];
%dclu=[0];
%[uu,t]=step(acl,bcl,cclu,dclu);
%su=plot(uu(:,:,1));
%grid;
%pause
%Step in all channels
%t=0:.1:100;
%u=[25*ones(1,251) 25*ones(1,500)0*ones(1,250)];
%u10=[u'];
%sysc1=ss(acl,bcl,ccl,dcl);
%[yy1,t]=lsim(sysc1,u10,t);
%plot(t,yy1)
%grid
%pause;
%cclu=[-ko eye(1,1)];
%dclu=[0];
%syscl2=ss(acl,bcl,cclu,dclu);
%[uu1,t]=lsim(syscl2,u10,t);
%su=plot(uu1(:,:,1));
%grid;
%pause
% Estimator design
qe=eye(1,1);
re=0.001*eye(1,1);
[ll,pp,el]=lqe(a,b,c,qe,re);
ac=zeros(1,1);
bc=-k1;
cc=eye(1,1);
dc=-ko;
acle=[a, b*cc, -b*ko;
   bc*c, ac, zeros(1,3);
   ll*c, b*cc, a-ll*c-b*ko];
bcle=[zeros(3,1);-bc;zeros(3,1)];
ccle=[c, zeros(1,4)];
dcle=zeros(1,1);
[ecle]=eig(acle);
[zcle]=tzero(acle,bcle,ccle,dcle);
dcgain=dcle-ccle*inv(acle)*bcle;
t=0:.1:100;
R=[0:.1:25, 25*ones(1,500), 0*ones(1,250)];
R11=[R'];
R11=[R11 R11 R11 R11 R11 R11];
syscl=ss(acle,bcle,ccle,dcle);
[yy,t]=lsim(syscl,R,t);
figure(1);
plot(t,yy,'--','LineWidth',2);
grid on;
hold on;
plot(t,R11,'-');
xlabel('Time (sec)');
ylabel('Temperature (K)');
title('Fig. 10.71 (a) Internal model controller: temperature tracking response');
%grid
h_axes = findobj(get(gcf,'Children'),'Type','axes');
grey = [0.7,0.7,0.7];
set(h_axes,'xcolor',grey,'ycolor',grey, ...
    'GridLineStyle','-','MinorGridLineStyle','-', ...
    'Units','pixels');
grid on
c11=copyobj(h_axes,gcf);
set(c11,'color','none','xcolor','k', ...
    'xgrid','off','ycolor','k', ...
    'ygrid','off');
%legend('y')
%pause;
hold off;
cclu=[zeros(1,3), eye(1,1), -ko];
dclu=zeros(1,1);
syscu=ss(acle,bcle,cclu,dclu);
[uuu,t]=lsim(syscu,R,t);
figure(2);
plot(t(1:753,:),uuu(1:753,:),'LineWidth',2);
hold on;
plot(t(752:818,:),0*ones(67,3),'LineWidth',2);
hold on;
plot(t(753:818,:),uuu(753:818,:),'--','LineWidth',2);
hold on;
plot(t(819:1001,:),uuu(819:1001,:),'LineWidth',2);
xlabel('Time (sec)');
ylabel('Lamp voltage');
title('Fig. 10.71 (b) Internal model controller: control effort');
legend('u');
nicegrid;

author: Sable

repo: mcbench-benchmarks


text:
function [ n_data_new, n, c ] = sigma_values ( n_data )

%*****************************************************************************80
%
%% SIGMA_VALUES returns some values of the Sigma function.
%
%  Definition:
%
%    SIGMA(N) is the sum of the distinct divisors of N, including 1 and N.
%
%  First values:
%
%     N  SIGMA(N)
%
%     1    1
%     2    3
%     3    4
%     4    7
%     5    6
%     6   12
%     7    8
%     8   15
%     9   13
%    10   18
%    11   12
%    12   28
%    13   14
%    14   24
%    15   24
%    16   31
%    17   18
%    18   39
%    19   20
%    20   42
%
%  Formula:
%
%    SIGMA(U*V) = SIGMA(U) * SIGMA(V) if U and V are relatively prime.
%
%    SIGMA(P**K) = ( P**(K+1) - 1 ) / ( P - 1 ) if P is prime.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA_NEW
%    from the previous call.
%
%    Output, integer N_DATA_NEW, the index of the test data.
%
%    Output, integer N, the argument of the Sigma function.
%
%    Output, integer C, the value of the Sigma function.
%
  n_max = 20;
  c_vec = [ ...
     1,    3,    4,    7,    6,   12,    8,   15,   13,   18, ...
    72,  128,  255,  176,  576, 1170,  618,  984, 2232, 2340 ];
  n_vec = [ ...
      1,   2,   3,   4,   5,   6,   7,   8,   9,  10, ...
     30, 127, 128, 129, 210, 360, 617, 815, 816,1000 ];

  n_data_new = n_data;

  if ( n_data_new < 0 )
    n_data_new = 0;
  end

  n_data_new = n_data_new + 1;

  if ( n_max < n_data_new )
    n_data_new = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data_new);
    c = c_vec(n_data_new);
  end

  return
end

repo: jburkardt-m


text:
function plotsid(t,iflaws,k);
%PLOTSID Schematic interference diagram of FM signals.
%	PLOTSID(T,IFLAWS,K) plots the schematic interference diagram of
%	(analytic) FM signals.
%
%	T : time instants,
%	IFLAWS : matrix of instantaneous frequencies,
%	         with as may columns as signal components.
%	K : distribution		(default : 2):
%	  K = 2     : Wigner-Ville
%	  K = 1/2   : D-Flandrin
%	  K = 0     : Bertrand (unitary)
%	  K = -1    : Unterberger (active)
%	  K = inf   : Margenhau-Hill-Rihaczek
%
%	Example :
%	 Nt=90; [y,iflaw]=fmlin(Nt,0.05,0.25);
%	 [y2,iflaw2]=fmconst(50,0.4);
%	 iflaw(:,2)=[NaN*ones(10,1);iflaw2;NaN*ones(Nt-60,1)];
%	 plotsid(1:Nt,iflaw,0);
%
%	See also PLOTIFL, MIDPOINT.

%	P. Flandrin, September 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin==1),
 error('at least 2 parameters : t and iflaws');
elseif (nargin==2),
 k=2; % Wigner-Ville
end;

indices=find(1-isnan(iflaws));
if (min(iflaws(indices))<0)|(max(iflaws(indices))>0.5),
 error ('each element of IFLAWS must be between 0 and 0.5');
end;

[iflawrow,iflawcol]=size(iflaws);
tcol=length(t);
clf; figure(gcf); plotifl(t,iflaws);
hold on;
col=['b','m','c','r'];

% auto-terms
for j=1:iflawcol,
 indices=find(1-isnan(iflaws(:,j)));
 Nbpoints=length(indices);
 for i=1:Nbpoints-1,
  ta=       t(indices(i))  *ones(1,Nbpoints-i);
  fa=  iflaws(indices(i),j)*ones(1,Nbpoints-i);
  tb= t(indices(i+1:Nbpoints));
  fb=iflaws(indices(i+1:Nbpoints),j)';
  [ti,fi]=midscomp(ta,fa,tb,fb,k);
  plot(ti,fi,['.',num2str(col(rem(j-1,4)+1))]);
 end;
end;

% cross-terms
for j1=1:iflawcol,
 indices1=find(1-isnan(iflaws(:,j1)));
 Nbpoints1=length(indices1);
 for j2=j1+1:iflawcol,
  indices2=find(1-isnan(iflaws(:,j2)));
  Nbpoints2=length(indices2);
  for i=1:Nbpoints1,
   ta=       t(indices1(i))   *ones(1,Nbpoints2);
   fa=  iflaws(indices1(i),j1)*ones(1,Nbpoints2);
   tb= t(indices2);
   fb=iflaws(indices2,j2)';
   [ti,fi]=midscomp(ta,fa,tb,fb,k);
   plot(ti,fi,'.g')
  end;
 end;
end;

hold off
axis([t(1) t(tcol) 0 0.5]);
grid
if k==2,
 dist=' of the Wigner-Ville distribution';
elseif k==1/2,
 dist=' of the D-Flandrin distribution';
elseif k==0,
 dist=' of the (unitary) Bertrand distribution';
elseif k==-1,
 dist=' of the (active) Unterberger distribution';
elseif k>1/sqrt(eps),
 dist=' of the Margenhau-Hill-Rihaczek distribution';
else
 dist='';
end

title(['Interference diagram',dist,' (k = ',num2str(k),')']);

author: HeLiangHIT


text:
function [uHessMat,uTriMat,Q,Z]=HessenbergTriangRed(A,B)
%%HESSENBERGTRIANGRED Perform a Hessenberg triangular reduction of the nXn
%               matrices A and B. This finds Q and Z such that Q'*A*Z= an
%               upper Hessenberg matrix and Q'*B*Z is an upper triangular
%               matrix.
%
%INPUTS: A, B Two real nXn matrices.
%
%OUTPUTS: uHessMat The upper Hessenberg matrix obtained from Q'*A*Z.
%          uTriMat The upper triangular matrix obtained with Q'*B*Z.
%
%This function implements algorithm 7.7.1 in [1].
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%September 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(A,1);

if(nargout>2)
    Z=eye(n,n);
    %B=Q*R. Replace B with Q'*B and A with Q'*A
    [Q,R]=qr(B);
    B=R;%Q'*B=R.
    A=Q'*A;

    for j=1:(n-2)
        for i=n:-1:(j+2)
            %Zeros A(i,j)
            [c,s]=GivensCS(A(i-1,j),A(i,j));
            rotMat=[c, s;
                   -s, c];
            A((i-1):i,j:n)=rotMat'*A((i-1):i,j:n);
            B((i-1):i,(i-1):n)=rotMat'*B((i-1):i,(i-1):n);

            Q(:,(i-1):i) = Q(:,(i-1):i)*rotMat;

            %Zeros B(i,i-1).
            [c,s]=GivensCS(-B(i,i),B(i,i-1));
            rotMat=[c, s;
                   -s, c];
            B(1:i,(i-1):i) = B(1:i,(i-1):i)*rotMat;
            A(1:n,(i-1):i) = A(1:n,(i-1):i)*rotMat;

            Z(:,(i-1):i) = Z(:,(i-1):i)*rotMat;
        end
    end
    uHessMat=A;
    uTriMat=B;
else
    %If Q and Z are not desired, don't compute them.

    %B=Q*R. Replace B with Q'*B and A with Q'*A
    [Q,R]=qr(B);
    B=R;%Q'*B=R.
    A=Q'*A;
    for j=1:(n-2)
        for i=n:-1:(j+2)
            %Zeros A(i,j)
            [c,s]=GivensCS(A(i-1,j),A(i,j));
            rotMat=[c, s;
                   -s, c];
            A((i-1):i,j:n)=rotMat'*A((i-1):i,j:n);
            B((i-1):i,(i-1):n)=rotMat'*B((i-1):i,(i-1):n);

            %Zeros B(i,i-1).
            [c,s]=GivensCS(-B(i,i),B(i,i-1));
            rotMat=[c, s;
                   -s, c];
            B(1:i,(i-1):i) = B(1:i,(i-1):i)*rotMat;
            A(1:n,(i-1):i) = A(1:n,(i-1):i)*rotMat;
        end
    end
    uHessMat=A;
    uTriMat=B;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function RegMeanSquareVersor(handles, metric)
%
% 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/>.

global planC stateS;
indexS = planC{end};

    [originF, spacingF, centerF] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationBaseDataset));

    [originM, spacingM, centerM] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationMovDataset));

    Opt_S_X = str2double(get(handles.versor_s_x, 'string'));
    Opt_S_Y = str2double(get(handles.versor_s_y, 'string'));
    Opt_S_Z = str2double(get(handles.versor_s_z, 'string'));
    tscale  = str2double(get(handles.versor_s_t, 'string'));
    Min_Step_Len = str2double(get(handles.versor_i_minstep, 'string'));
    Max_Step_Len = str2double(get(handles.versor_i_maxstep, 'string'));
    Iteration_Num = str2double(get(handles.versor_i_num, 'string'));

    FImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
    MImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray;

    dimF = size(FImg);
    dimM = size(MImg);
    %generate cliped base dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseTrans'));
    if ~isempty(clipBox)
        FImg = FImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originF(1) = originF(1) + spacingF(1)*double(clipBox(1)-1);
        originF(2) = originF(2) + spacingF(2)*double(uint16(dimF(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseSag'));
    if ~isempty(clipBox)
        FImg = FImg(:, :, clipBox(2):clipBox(4));
        originF(3) = originF(3) + spacingF(3)*double(clipBox(2)-1);
    end

    %generate cliped move dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movTrans'));
    if ~isempty(clipBox)
        MImg = MImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originM(1) = originM(1) + spacingM(1)*double(clipBox(1)-1);
        originM(2) = originM(2) + spacingM(2)*double(uint16(dimM(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movSag'));
    if ~isempty(clipBox)
        MImg = MImg(:, :, clipBox(2):clipBox(4));
        originM(3) = originM(3) + spacingM(3)*double(clipBox(2)-1);
    end

%downsample the input datasets
    dsampled = 0;
    if (get(handles.dsampleCheck,'Value') == get(handles.dsampleCheck,'Max'))

        tic;
        xyScale =  2; zScale = 2;
        spacingF = [spacingF(1)*xyScale spacingF(2)*xyScale spacingF(3)*zScale];
        spacingM = [spacingM(1)*xyScale spacingM(2)*xyScale spacingM(3)*zScale];
        disp('downSampling ...');
%         FImg = imdesample3d(FImg,xyScale,zScale);
%         MImg = imdesample3d(MImg,xyScale,zScale);
%         FImg=GPReduce2D(FImg,1);
%         MImg=GPReduce2(MImg,1);
        FImg = FImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);
        MImg = MImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);

        toc;
        dsampled = 1;
    end

%     call registration method
    switch metric
        case 'mean squares'
            metricIndex = 0;
        case 'normalized correlation'
            metricIndex = 1;
    end

    tic;
tic;

%flip the source datasets on X for itk coordinate system
    fdim = 1;
    FImg = flipdim(FImg, fdim);
    MImg = flipdim(MImg, fdim);

% transform initializing modes 0:MomentON 1:GeometryOn 2:initail transform On
    initMode = get(handles.InitTrans, 'value');

%prepare the initial transform matrix
    transMB = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).transM;
    transMM = planC{indexS.scan}(stateS.imageRegistrationMovDataset).transM;
    if isempty(transMM), transMM = eye(4); end;
    if isempty(transMB), transMB = eye(4); end;
    transM = inv(transMB) * transMM;

%flip the moving dataset to match the initial transM
    Tf = eye(4);
    flipM = 0;
    if (get(handles.flipX, 'value'))
        MImg = flipdim(MImg, 2);
        Tf = [-1 0 0 2*centerM(1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
        flipM = 2;
    end
    if (get(handles.flipY, 'value'))
        MImg = flipdim(MImg, 1);
        Tf = [1 0 0 0; 0 -1 0 2*centerM(2); 0 0 1 0; 0 0 0 1];
        flipM = 1;
    end
    if (get(handles.flipZ, 'value'))
        MImg = flipdim(MImg, 3);
        Tf = [1 0 0 0; 0 1 0 0; 0 0 -1 2*centerM(3); 0 0 0 1];
        flipM = 3;
    end
    transM = transM * inv(Tf);
    rotM = transM(1:3, 1:3); transV = transM(1:3, 4);
    rotM = rotM';
    transV = -transV;

    try
        [im, Rotation, Offset] = VersorRigid3D(int16(FImg), originF, spacingF, ...
                                           int16(MImg), originM, spacingM, ...
                                           Opt_S_X, Opt_S_Y, Opt_S_Z, tscale, Min_Step_Len, Max_Step_Len, ...
                                           Iteration_Num, metricIndex, rotM, transV, initMode);
    catch
        [im, Rotation, Offset] = VersorRigid3D_64(int16(FImg), originF, spacingF, ...
                                           int16(MImg), originM, spacingM, ...
                                           Opt_S_X, Opt_S_Y, Opt_S_Z, tscale, Min_Step_Len, Max_Step_Len, ...
                                           Iteration_Num, metricIndex, rotM, transV, initMode);
    end

    toc;

    im = flipdim(im, fdim);

    output = cell(1, 16);

    output{1} = ['Translation X  = ' num2str(Offset(1))];
    output{2} = ['Translation Y  = ' num2str(Offset(2))];
    output{3} = ['Translation Z  = ' num2str(Offset(3))];
    output{4} = ['Rotation Angle = ' num2str(asin(Rotation(2))*180/pi)];
    output{5} = ['Iterations     = ' num2str(Offset(4))];
    output{6} = ['Metric value   = ' num2str(Offset(5))];
    output{7} = ['RCenter X   = ' num2str(Offset(7))];
    output{8} = ['RCenter Y   = ' num2str(Offset(8))];
    output{9} = ['RCenter Z   = ' num2str(Offset(9))];
    output{10} = ['Offset X   = ' num2str(Offset(10))];
    output{11} = ['Offset Y   = ' num2str(Offset(11))];
    output{12} = ['Offset Z   = ' num2str(Offset(12))];

    set(handles.OutputList, 'string', output);

%update the transM;

    Tv = [Offset(1) Offset(2) Offset(3)];
    Cv = [Offset(7) Offset(8) Offset(9)];
    Cv(3) = -Cv(3); %CERR z axis is opposite to the ITK.

    rot = reshape(Rotation, 3,3);
    offset = [Offset(10) Offset(11) Offset(12)]';

%     rot(2,1) = -rot(2,1); rot(1,2) = -rot(1,2);
%     rot(1,3) = -rot(1,3); rot(3,1) = -rot(3,1);
%     rot(2,3) = -rot(2,3); rot(3,2) = -rot(3,2);
%
%     offset(1) = -offset(1);
%     offset(2) = -offset(2);
%     offset(3) = -offset(3);

    TM = eye(4);
    TM(:,4) = [offset; 1];

    RM = eye(4);
    RM(1:3, 1:3) = rot;

    newTransform = inv(TM*RM);
    newTransform =  transMB * newTransform * Tf;

    scanSetM = stateS.imageRegistrationMovDataset;
    scanSetF = stateS.imageRegistrationBaseDataset;
    planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = newTransform;

% save the resampled dataset
    if (get(handles.saveCheck, 'value'))
        if (~dsampled)
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = im;
        else
            %imReg = resample(scanSetF, scanSetM); % need to up-sample im ????
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = imReg;
        end
        controlFrame('refresh');
    end

    sliceCallBack('refresh');

end


text:
function a = double(s)
%DOUBLE Converts a sparse tensor to a dense multidimensional array.
%
%  See also SPTENSOR, SPTENSOR/FULL.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

a = zeros([size(s) 1 1]);
if nnz(s) > 0
    a(tt_sub2ind(size(s),s.subs)) = s.vals;
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% The COBRAToolbox: testCheckObjective.m
%
% Purpose:
%     - Tests the changeObjective function
%

global CBTDIR

currentDir = pwd; % save the current path

% initialize the test
fileDir = fileparts(which('testCheckObjective'));
cd(fileDir);

% load the ecoli core model
model = getDistributedModel('ecoli_core_model.mat');

abbr = checkObjective(model);

assert(strcmp(abbr{1}, 'Biomass_Ecoli_core_w_GAM'));

% test if there are more than 1 objective
model.c(1) = 1;
abbr = checkObjective(model);
assert(strcmp(abbr{1}, 'ACALD'));
assert(strcmp(abbr{2}, 'Biomass_Ecoli_core_w_GAM'));

% test if there is no objective - will throw a warning message
w = warning ('off','all');
model.c(1) = 0;
model.c(13) = 0;
try
    abbr = checkObjective(model);
catch ME
    assert(length(ME.message) > 0)
end
w = warning ('on','all');

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function tc = tc_visualizeCorAnal(tc, params, parent);
%
% tc = tc_visualizeCorAnal(tc, params, <parent=tc.ui.plot panel>);
%
% Apply a corAnal to an ROI time course, and visualize the
% results in the TCUI.
%
% params: struct w/ params for the corAnal:
%   nCycles: # of cycles per scan for the fitted sinusoid.
%   frames: frames to use from the wholeTc. Default: use all frames.
% Pops up a dialog if omitted.
%
% ras, 09/2005.
if notDefined('tc'), tc = get(gcf,'UserData'); end

if notDefined('params')
    if isfield(tc, 'corAnal')
        % we can use the existing corAnal
    else
        tc = tc_applyCorAnal(tc);
    end
else
    tc = tc_applyCorAnal(tc, params);
end

if notDefined('parent'),  parent = tc.ui.plot;                end
if parent==gcf | parent==get(gcf, 'CurrentAxes')
    % make a uipanel to fit on the target
    parent = uipanel('Parent', parent, ...
        'Units', 'normalized', ...
        'BackgroundColor', get(gcf, 'Color'), ...
        'Position', [0 0 1 1])
end

nCycles = tc.corAnal.nCycles;
frames = tc.corAnal.frames;

%%%%%get the cor anal data
C = tc.corAnal;

%%%%%visualize the results
% delete existing objects in display
delete( findobj('Parent', tc.ui.plot) );

% show FFT, highlighting corAnal peak
tc = tc_plotFFT(tc,C.nCycles,[.12 .58 .3 .3]);
set(gca,'Box','off')

%%%%%plot amp, ph, and co values
% (we also check if there are polar angle retinotopy params defined,
% and if so, use it to match the angle to the estimated true polar angle)
subplot(222);
x = C.co .* cos(C.ph);
y = C.co .* sin(C.ph);

% check if retinotopy params have been set for this scan,
% and if so, and it's a polar angle scan, make it plot
% the expected polar angle:
try
    V = getCurView;
    p = retinoGetParams(V);
    if ~isempty(p) & isequal(p.type, 'polar_angle')
        disp('Setting Phase to match expected polar angle...')
        theta = polarAngle(C.ph, p)
        theta = -deg2rad( theta - 90 );
        x = C.co .* cos(theta);
        y = C.co .* sin(theta);
    end

catch
    % don't worry
end

params.grid = 'on';
params.line = 'off';
params.gridColor = [.4 .4 .4];
params.fontSize = 12;
params.symbol = 'o';
params.size = 5;
params.color = 'w';
params.fillColor = 'w';
params.maxAmp = 1;
params.ringTicks = [0:0.2:1];
polarPlot(0, params);
h = plot(x, y, 'ro', 'MarkerSize', params.size);
set(h,'MarkerFaceColor','r')
title('Coherence Vs. Phase Polar Plot', 'FontSize', 14);
text(1, .5, sprintf('Co: %2.2f',C.co), 'FontWeight', 'bold');
text(1, .25, sprintf('Amp: %2.2f',C.amp), 'FontWeight', 'bold');
text(1, 0, sprintf('Ph: %2.2f rad (%2.2f deg)',C.ph, rad2deg(C.ph)), ...
    'FontWeight', 'bold');
if exist('p', 'var') & ~isempty(p) & isequal(p.type, 'polar_angle')
    text(1, -.25, '(Expected Polar Angle)', 'FontWeight', 'bold');
end

%%%%%plot mean 2 cycles (for clarity) and fitted predictor
subplot(223);
framesPerCycle = length(frames)/nCycles;
t = [1:framesPerCycle] .* tc.TR;
cycles = reshape(tc.wholeTc(frames),[framesPerCycle nCycles]);
meanCycle = mean(cycles,2);
semCycle = std(cycles,1,2) / sqrt(nCycles);
predCycle = mean(reshape(C.predictor,[framesPerCycle nCycles]),2);
hold on
errorbar(t,meanCycle,semCycle,'k-','LineWidth',2.5);
plot(t, predCycle, 'r', 'LineWidth', 1.5);
xlabel('Cycle time, sec','FontSize',14);
ylabel('% Signal','FontSize',14);
set(gca,'Box','off')
title('Mean Cycle + Predictor', 'FontSize', 14);

%%%%%%show time course + selected predictors
subplot(224);
t = [1:length(frames)] .* tc.TR;
hold on
plot(t, tc.wholeTc(frames), 'k-.', 'LineWidth', 2.5);
plot(t, C.predictor, 'r', 'LineWidth', 1.5);
xlabel('Time, sec','FontSize',14)
ylabel('% Signal','FontSize',14)
title('Time course + Scaled Predictor', 'FontSize', 14)
set(gca,'Box','off')

% for longer time courses, may want a UI control
scale = 300;  % max seconds to nicely plot TC
scrollbar(gca, scale);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Copyright (C) 2014-2016, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym acosd (@var{x})
%% Symbolic inverse cos function with output in degrees.
%%
%% Example:
%% @example
%% @group
%% acosd (sqrt (sym (2))/2)
%%   @result{} (sym) 45
%%
%% syms x
%% y = acosd (x)
%%   @result{} y = (sym)
%%       180⋅acos(x)
%%       ───────────
%%            π
%% @end group
%% @end example
%%
%% @seealso{@@sym/cosd, @@sym/acos}
%% @end defmethod

function y = acosd(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('lambda a: deg(acos(a))', x);
end

%!error acosd (sym(1), 2)
%!assert (isequaln (acosd (sym(nan)), sym(nan)))

%!test
%! f1 = acosd (sym(1)/2);
%! f2 = acosd (1/2);
%! assert (double (f1), f2, -eps)

%!test
%! D = [1 2; 3 4]/4;
%! A = sym([1 2; 3 4])/4;
%! f1 = acosd (A);
%! f2 = acosd (D);
%! assert (double (f1), f2, -eps)

repo: octsympy


text:
function B = nDshift(B, sig, varargin)

idims = 2;
if ~isempty(varargin)
    idims = varargin{1};
end
if numel(idims)>1 && numel(sig)>1
    sigall = sig;
else
    sigall = repmat(sig, numel(idims), 1);
end

for i = 1:length(idims)
    sig = sigall(i);

    idim = idims(i);
    Nd = ndims(B);

    B = permute(B, [idim 1:idim-1 idim+1:Nd]);

    dsnew = size(B);

    B = reshape(B, dsnew(1), []);
    dsnew2 = size(B);

    if sig>0
        irange = sig+1:dsnew2(1);
        B(irange, :) = B(irange-sig, :);
        B(1:sig, :) = 0;
    else
        sig = -sig;
        irange = 1:dsnew2(1)-sig;
        B(irange, :) = B(irange+sig, :);
        B([1:sig] + dsnew2(1)-sig, :) = 0;
    end

     B = reshape(B, dsnew);

    B = permute(B, [2:idim 1 idim+1:Nd]);
end

author: cortex-lab


text:
function [ft2] = mm22ft2(mm2)
% Convert area from square millimeters to square feet.
% Chad A. Greene 2012
ft2 = mm2*0.00001076391041671;


text:
%kkman 'Displays and prints man pages for software objects'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kman.pane file
%
% Parameters:
% String: tb 'Toolbox', default: ' ': 'toolbox in which to search for the manpage'
% String: k 'Keyword', default: ' ': 'keyword to search the whatis file for'
% String: f 'Roff Formatter', default: ' ': 'roff formater command for printing man pages'
% Toggle: r 'Reformat', default: 0: 'reformats the man page'
% Toggle: man 'man', default: 0: 'format using manpage macros'
% Toggle: hlp 'help', default: 0: 'format using helppage macros'
% Toggle: manual 'manual', default: 0: 'format using manual macros'
%
% Example: kkman( {'tb',' ';'k',' ';'f',' ';'r',0;'man',0;'hlp',0;'manual',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kman - Displays and prints man pages for software objects
%
%  DESCRIPTION
% .I kman
% is used to display and print out man pages for Khoros.
%
% .I kman
% will search all toolboxes specified in the user's KHOROS_TOOLBOX file
% for the man page specified.  If a section is specified, kman will
% search all toolboxes in the section's directory for the specified man page.
%
% For a printed hardcopy of the man page, you must specify both the -f and
% -P options. It is recommended that you use groff if available for the formatter.
%
%
%
%  EXAMPLES
% % kman -f groff -Plp kman
%
% Prints the man page for kman to the printer 'lp' using the roff
% formatter 'groff'.
% % kman guise
%
% Displays the man page for guise to the screen.
%
%  "SEE ALSO"
% kgenwhatis(1)
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkman(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kkman(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'tb', ' ';'k', ' ';'f', ' ';'r', 0;'man', 0;'hlp', 0;'manual', 0};
maxval={0,0,0,0,0,0,0};
minval={0,0,0,0,0,0,0};
istoggle=[1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'String','String','String','Toggle','Toggle','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; 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(0);
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
callKhoros([w 'kman"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function nurbs = nrbtform(nurbs,tmat)
%
% Function Name:
%
%   nrbtform - Apply transformation matrix to the NURBS.
%
% Calling Sequence:
%
%   tnurbs = nrbtform(nurbs,tmatrix);
%
% Parameters:
%
%   nurbs	: NURBS data structure (see nrbmak for details).
%
%   tmatrix     : Transformation matrix, a matrix of size (4,4) defining
%                 a single or multiple transformations.
%
%   tnurbs	: The return transformed NURBS data structure.
%
% Description:
%
%   The NURBS is transform as defined a transformation matrix of size (4,4),
%   such as a rotation, translation or change in scale. The transformation
%   matrix can define a single transformation or multiple series of
%   transformations. The matrix can be simple constructed by the functions
%   vecscale, vectrans, vecrotx, vecroty, and vecrotz.
%
% Examples:
%
%   Rotate a square by 45 degrees about the z axis.
%
%   rsqr = nrbtform(nrbrect(), vecrotz(deg2rad(45)));
%   nrbplot(rsqr, 10);
%
% See:
%
%   vecscale, vectrans, vecrotx, vecroty, vecrotz

%  D.M. Spink
%  Copyright (c) 2000

if nargin < 2
  error('Not enough input arguments!');
end;

if iscell(nurbs.knots)
  % NURBS is a surface
  [dim,nu,nv] = size(nurbs.coefs);
  nurbs.coefs = reshape(tmat*reshape(nurbs.coefs,dim,nu*nv),[dim nu nv]);
else
  % NURBS is a curve
  nurbs.coefs = tmat*nurbs.coefs;
end


text:
function mean = dirichlet_mix_mean ( comp_num, elem_num, a, comp_weight )

%*****************************************************************************80
%
%% DIRICHLET_MIX_MEAN returns the means of a Dirichlet mixture PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer COMP_NUM, the number of components in the Dirichlet
%    mixture density, that is, the number of distinct Dirichlet PDF's
%    that are mixed together.
%
%    Input, integer ELEM_NUM, the number of elements of an observation.
%
%    Input, real A(ELEM_NUM,COMP_NUM), the probabilities for
%    element ELEM_NUM in component COMP_NUM.
%    Each A(I,J) should be positive.
%
%    Input, real COMP_WEIGHT(COMP_NUM), the mixture weights of the densities.
%    These do not need to be normalized.  The weight of a given component is
%    the relative probability that that component will be used to generate
%    the sample.
%
%    Output, real MEAN(ELEM_NUM), the means for each element.
%
  comp_weight_sum = sum ( comp_weight );

  mean(1:elem_num) = 0.0;

  for comp_i = 1 : comp_num
    comp_mean(1:elem_num) = dirichlet_mean ( elem_num, a(1:elem_num,comp_i) );
    mean(1:elem_num) = mean(1:elem_num) ...
      + comp_weight(comp_i) * comp_mean(1:elem_num);
  end

  mean(1:elem_num) = mean(1:elem_num) / comp_weight_sum;

  return
end


text:
function b = power(a,n)
% function B=power(A,N)
%
% DESCRIPTION
%   Element-by-element powers
%
% INPUTS
%   A: polynomial
%   N: matrix of natural numbers
%
% OUTPUTS
%   B: polynomial, result of element-by-element power A.^N.
%
% SYNTAX
%   B = A.^N
%     Element-by-element powers. A and B must have the same dimensions
%     unless one is a scalar. A scalar can operate into anything.
%   B=power(A,N)
%     Function-call form of power.

% 6/7/2002 PJS  Initial Coding
% 12/12/2010 PJS Call matrix.^matrix code for scalar expansion cases

% Check number of inputs
error(nargchk(2,2,nargin));
sza=size(a);
szn=size(n);

if isempty(a) || isempty(n)
    if isempty(a) && all(szn==[1 1])
        % empty.^scalar = empty(sza)
        b=polynomial(zeros(sza));
        return;
    elseif isempty(n) && all(sza==[1 1])
        % scalar.^empty = empty(szn)
        b=polynomial(zeros(szn));
        return;
    elseif all(sza==szn)
        b=polynomial(zeros(sza));
        return;
    else
        error('Matrix dimensions must agree.');
    end

elseif sza==szn
    % Matrix .^ Matrix
    b = zeros(sza);

    % exponent = 0
    idx = find( n==0 );
    b(idx) = 1;
    b = polynomial(b);

    % exponent = 1
    idx = find( n==1 );
    L.type = '()';
    L.subs = {idx};
    b = subsasgn(b,L,subsref(a,L));

    % exponent = 2;
    idx = find( n>1 );
    idx = idx(:)';
    for i1=idx
        L.subs = {i1};
        ai = subsref(a,L);
        Nt = size(ai.coefficient,1);
        if Nt==1
            bi = ai;
            bi.coefficient = bi.coefficient^n(i1);
            bi.degmat = bi.degmat*n(i1);
        else
            bi = mpower(ai,n(i1));
        end
        b = subsasgn(b,L,bi);
    end

    %         % XXX This is pretty slow.
    %         b = a;
    %         L.type = '()';
    %         for i1 = 1:sza(1);
    %             for i2 = 1:sza(2);
    %                 L.subs = {i1,i2};
    %                 aij = subsref(a,L);
    %                 bij = mpower(aij,n(i1,i2));
    %                 b = subsasgn(b,L,bij);
    %             end
    %         end
    %     end
elseif all(szn==[1 1])
    % Matrix.^Scalar:
    n = repmat(n,sza);
    b = power(a,n);

    %     b = polynomial(1);
    %     for i1 = 1:n
    %         b = b.*a;
    %     end
elseif all(sza==[1 1])
    % Scalar .^ Matrix:
    a = repmat(a,szn);
    b = power(a,n);

    %     % XXX This is pretty slow.
    %     b = polynomial(zeros(szn));
    %     L.type = '()';
    %     for i1 = 1:szn(1);
    %         for i2 = 1:szn(2);
    %             L.subs = {i1,i2};
    %             bij = mpower(a,n(i1,i2));
    %             b = subsasgn(b,L,bij);
    %         end
    %     end
else
    error('Matrix dimensions must agree');
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [M_V] = vox2ras_ksolve(M_R, Vc_Ps, varargin)
%%
%% NAME
%%
%%     vox2ras_ksolve.m (vox2ras_k{-space column}solve)
%%
%% AUTHOR
%%
%%	Rudolph Pienaar
%%
%% SYNOPSIS
%%
%%     [M_V] = vox2ras_ksolve(M_R, Vc_Ps, Vr_logicalSpace)
%%
%% ARGUMENTS
%%
%%	M_R		in	4x4 vox2ras candidate matrix with correct
%%					direction cosines, but invalid
%%					k-space column
%%      Vc_Ps		in      column vector defining the slice position
%%					(typically read from meas.asc)
%%	Vr_logicalSpace	in/opt	row vector that defines the logical dimension
%%					size. If omitted, this is assumed
%%					to be [256 256 128]
%%	M_V		out	complete 4x4 vox2ras matrix
%%
%% DESCRIPTION
%%
%%	"vox2ras_ksolve" uses the given rotation components of M_R and the
%%	starting point, Vc_Ps, to determine the center of k-space and thus
%%	create a fully-qualified vox2ras matrix.
%%
%% PRECONDITIONS
%%
%%	o M_R is 4x4 where the first 3x3 submatrix contains the direction
%%		cosines. The 4th column is ignored, as is the 4th row.
%%	o The vector Vc_Ps is typically read from a Siemens meas.asc file such
%%	  that
%%		Vc_Ps(1)	= sSliceArray.asSlice[0].sPosition.dSag
%%		Vc_Ps(2)	= sSliceArray.asSlice[0].sPosition.dCor
%%		Vc_ps(3)	= sSliceArray.asSlice[0].sPosition.dTra
%%
%% POSTCONDITIONS
%%
%%	o M_V(1:3, 1:3) = M_R(1:3, 1:3)
%%	o M_V(:,4)	= k-space center
%%
%% SEE ALSO
%%
%%	vox2ras_rsolveAA- determine the rotational component of a vox2ras matrix
%%				using Siemens reference orientations directly
%%	vox2ras_rsolve	- determine the rotational component of a vox2ras matrix
%%				using Siemens reference orientations indirectly
%%	vox2ras_dfmeas	- main function: determines the vox2ras matrix from a
%%			  Siemens meas.asc file.
%%
%% HISTORY
%%
%% 26 May 2004
%% o Initial design and coding.
%%

%
% vox2ras_ksolve.m
%
% Original Author: Rudolph Pienaar
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

M_V		= M_R;
M_V(4,:)	= M_R(4,:)*0;
M_V(:,4)	= M_R(:,4)*0;
M_V(4,4)	= 1;

Vr_logicalSpace	= [ 256 256 128 ];
if length(varargin)
	Vr_logicalSpace = varargin{1};
end

%% First read the direction cosines
Vc_x 		= M_R(1:3,1);
Vc_y 		= M_R(1:3,2);
Vc_z 		= M_R(1:3,3);

xoff		= Vr_logicalSpace(1)/2;
yoff		= Vr_logicalSpace(2)/2;
zoff		= Vr_logicalSpace(3)/2;

%% Solve for the k-space center:
%% The first two components are found by moving along the direction cosines
Vc_Pe1		=  - (Vc_Ps + (xoff*Vc_x + yoff*Vc_y + zoff*Vc_z));
%% and the last component is found by moving against the direction cosines
Vc_Pe2		= (Vc_Ps - (xoff*Vc_x + yoff*Vc_y + zoff*Vc_z));
Vc_Pe(1)	= Vc_Pe1(1) -0.67;		%% a strange correction?
Vc_Pe(2)	= Vc_Pe1(2);
Vc_Pe(3)	= Vc_Pe2(3);
M_V(1:3, 4)	= Vc_Pe';

%% All done!

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef GeneralizedHoughGuil < handle
    %GENERALIZEDHOUGHGUIL  Generalized Hough transform
    %
    % Finds arbitrary template in the grayscale image using
    % Generalized Hough Transform.
    %
    % Detects position, translation and rotation.
    %
    % ## References
    % > Guil, N., Gonzalez-Linares, J.M. and Zapata, E.L. (1999).
    % > Bidimensional shape detection using an invariant approach.
    % > Pattern Recognition 32 (6): 1025-1038.
    %
    % See also: cv.GeneralizedHoughBallard
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % Canny low threshold. default 50
        CannyLowThresh
        % Canny high threshold. default 100
        CannyHighThresh
        % Minimum distance between the centers of the detected objects.
        % default 1.0
        MinDist
        % Inverse ratio of the accumulator resolution to the image resolution.
        % default 1.0
        Dp
        % Maximal size of inner buffers. default 1000
        MaxBufferSize

        % Angle difference in degrees between two points in feature.
        % default 90.0
        Xi
        % Maximal difference between angles that treated as equal. default 1.0
        AngleEpsilon
        % Feature table levels. default 360
        Levels
        % Minimal rotation angle to detect in degrees. default 0.0
        MinAngle
        % Maximal rotation angle to detect in degrees. default 360
        MaxAngle
        % Angle step in degrees. default 1.0
        AngleStep
        % Angle votes threshold. default 15000
        AngleThresh
        % Minimal scale to detect. default 0.5
        MinScale
        % Maximal scale to detect. default 2.0
        MaxScale
        % Scale step. default 0.05
        ScaleStep
        % Scale votes threshold. default 1000
        ScaleThresh
        % Position votes threshold. default 100
        PosThresh
    end

    %% GeneralizedHoughGuil
    methods
        function this = GeneralizedHoughGuil()
            %GENERALIZEDHOUGHGUIL  Constructor
            %
            %     obj = cv.GeneralizedHoughGuil()
            %
            % See also: cv.GeneralizedHoughGuil.detect
            %
            this.id = GeneralizedHoughGuil_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.GeneralizedHoughGuil
            %
            if isempty(this.id), return; end
            GeneralizedHoughGuil_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.GeneralizedHoughGuil.empty,
            %  cv.GeneralizedHoughGuil.load
            %
            GeneralizedHoughGuil_(this.id, 'clear');
        end

        function load(this, filename, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.GeneralizedHoughGuil.save
            %
            GeneralizedHoughGuil_(this.id, 'load', filename, varargin{:});
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.GeneralizedHoughGuil.load
            %
            GeneralizedHoughGuil_(this.id, 'save', filename);
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.GeneralizedHoughGuil.clear,
            %  cv.GeneralizedHoughGuil.load
            %
            b = GeneralizedHoughGuil_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.GeneralizedHoughGuil.save,
            %  cv.GeneralizedHoughGuil.load
            %
            name = GeneralizedHoughGuil_(this.id, 'getDefaultName');
        end
    end

    %% GeneralizedHough
    methods
        function [positions,votes] = detect(this, varargin)
            %DETECT  Find template on image
            %
            %     [positions,votes] = hough.detect(image)
            %     [positions,votes] = hough.detect(edges, dx, dy)
            %
            % ## Input
            % * __image__ input image, 8-bit 1-channel image
            % * __edges__ image edges
            % * __dx__ image x-derivative of the same size as `edges` and
            %   single-precision floating type.
            % * __dy__ image y-derivate of the same size and type as `dx`.
            %
            % ## Output
            % * __positions__ Cell array of 4-element vectors, each of the
            %   form: `[posx, posy, scale, angle]`
            % * __votes__ Cell array of 3-element vectors, of the same length
            %   as `positions`.
            %
            [positions,votes] = GeneralizedHoughGuil_(this.id, 'detect', varargin{:});
        end

        function setTemplate(this, varargin)
            %SETTEMPLATE  Set template to search
            %
            %     hough.setTemplate(templ)
            %     hough.setTemplate(edges, dx, dy)
            %     hough.setTemplate(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __templ__ template, 8-bit 1-channel image
            % * __edges__ template edges
            % * __dx__ template x-derivative of the same size as `edges` and
            %   single-precision floating type.
            % * __dy__ template y-derivate of the same size and type as `dx`.
            %
            % ## Options
            % * __Center__ Template center `[x,y]`. The default `[-1,-1]`
            %   will use `[size(templ,2) size(templ,1)]./2` as center.
            %
            % In the first variant of the function (with the `templ` input),
            % the `edges` are internally calculated using the cv.Canny filter,
            % and the derivatives `dx` and `dy` using cv.Sobel operator.
            %
            GeneralizedHoughGuil_(this.id, 'setTemplate', varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.CannyHighThresh(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'CannyHighThresh');
        end
        function set.CannyHighThresh(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'CannyHighThresh', value);
        end

        function value = get.CannyLowThresh(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'CannyLowThresh');
        end
        function set.CannyLowThresh(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'CannyLowThresh', value);
        end

        function value = get.Dp(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'Dp');
        end
        function set.Dp(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'Dp', value);
        end

        function value = get.MaxBufferSize(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'MaxBufferSize');
        end
        function set.MaxBufferSize(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'MaxBufferSize', value);
        end

        function value = get.MinDist(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'MinDist');
        end
        function set.MinDist(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'MinDist', value);
        end

        function value = get.AngleEpsilon(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'AngleEpsilon');
        end
        function set.AngleEpsilon(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'AngleEpsilon', value);
        end

        function value = get.AngleStep(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'AngleStep');
        end
        function set.AngleStep(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'AngleStep', value);
        end

        function value = get.AngleThresh(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'AngleThresh');
        end
        function set.AngleThresh(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'AngleThresh', value);
        end

        function value = get.Levels(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'Levels');
        end
        function set.Levels(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'Levels', value);
        end

        function value = get.MaxAngle(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'MaxAngle');
        end
        function set.MaxAngle(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'MaxAngle', value);
        end

        function value = get.MaxScale(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'MaxScale');
        end
        function set.MaxScale(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'MaxScale', value);
        end

        function value = get.MinAngle(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'MinAngle');
        end
        function set.MinAngle(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'MinAngle', value);
        end

        function value = get.MinScale(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'MinScale');
        end
        function set.MinScale(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'MinScale', value);
        end

        function value = get.PosThresh(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'PosThresh');
        end
        function set.PosThresh(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'PosThresh', value);
        end

        function value = get.ScaleStep(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'ScaleStep');
        end
        function set.ScaleStep(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'ScaleStep', value);
        end

        function value = get.ScaleThresh(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'ScaleThresh');
        end
        function set.ScaleThresh(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'ScaleThresh', value);
        end

        function value = get.Xi(this)
            value = GeneralizedHoughGuil_(this.id, 'get', 'Xi');
        end
        function set.Xi(this, value)
            GeneralizedHoughGuil_(this.id, 'set', 'Xi', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [d,a] =  training(a,d)

%if alpha empty in child | alpha full in rss model, train child (first train/retrain).
%if alpha full in child, empty in rss, only train rss.

if isempty(a.child.alpha) | ~isempty(a.alpha) % train algorithm first
   [r a.child]=train(a.child,d);
end

alpha=a.child.alpha;
svs=find(sum(abs(alpha)',1)>1e-5);

if(isfield(struct(a.child),'X'))
   a.Xsv=get(a.child.X,svs);
else
    a.Xsv=get(a.child.Xsv,svs);
end

K=calc(a.child.child,a.Xsv);
newalpha=alpha(svs,:)*0; alpha=alpha(svs,:); origalpha=alpha;
w2=[]; for i=1:size(alpha,2) w2(i)=alpha(:,i)'*K*alpha(:,i); end;
worig=w2; origsvs=svs;
loops=1;
disp('compressing..');
svs=[];

while max(w2)>a.tolerance & loops<a.max_loops
   %w is determined by alpha ONLY

   wx=alpha(:,1)*0;
   for i=1:size(alpha,2)
    f=find(abs(alpha(:,i))>0);
    tmp= w2(i) - (((K(:,f)*alpha(f,i)).^2) ./ diag(K)) ;  % calc all (w.x_i)^2/||x_i||^2 for each w
    if a.bal_w
     wx=wx+ (tmp/w2(i));  % was just tmp, attempt to normalize problems
    else
     wx=wx+tmp;
    end
   end
   if a.dont_revisit
     wx(svs)=Inf;    % ensure do not revisit existing basis function
   end
   [val,I]=min(wx); % find arg max w.x_i
   if ~isempty(intersect(I,svs)) & a.dont_revisit
        break; % seen everything
    end;
   oldsvs=svs;
   if length(oldsvs)<length(union(svs,I))
    svs=[oldsvs;I];
   end

   if loops>0 & (mod(loops,a.backfit)==0  | (loops<a.backfit_at_start) | ...
		 (sum(loops==a.test_on)==1 & a.backfit>0) )
       % find optimal alphas by backfitting

       %replace this by finvupdate
      if ~strcmp(a.optimizer,'iterative')
       %inv=pinv(K(svs,svs)); % only need to compute inverse once
       minv=inv(K(svs,svs)+eye(length(svs))*1e-6);
                             % only need to compute inverse once
       else
                newsv=I;
                 if length(oldsvs)==0
                   R=1/sqrt(K(newsv,newsv)); minv=1/(K(newsv,newsv));
               else
                  newsv=I;
                  [R,minv]=finvupdate(a,R,K(oldsvs,newsv),K(newsv,newsv));
         end
       end

     alpha=origalpha;
     for i=1:size(alpha,2)
      beta=minv*K(svs,:)*origalpha(:,i);
      newalpha(svs,i)=beta;  alpha(svs,i)=alpha(svs,i)-beta;
     end

	if ~isempty(a.dtst) & sum(loops==a.test_on)==1
	                 % run on separate test / validation set
                         %  classification only right now

	  a.alpha=newalpha;a.Xsv=get(a.child.Xsv,origsvs);
	  fin=find(sum(abs(a.alpha)',1)>a.alpha_cutoff);
	  a.alpha=a.alpha(fin,:); a.Xsv=get(a.Xsv,fin);
          a.b0=a.child.b0;

         if a.reoptimize_b % reoptimize b // only for pat. rec. right now
         a.b0=0;
	   r=test(a,d); a.b0=[];
	  for i=1:size(newalpha,2)
	     rr=r.X; rr=rr(:,i); [x s2]=sort(rr); y=r.Y(s2,i);
            xs=cumsum(y(end:-1:1)==1);
            [m1 m2]=max(cumsum(y==-1)+xs(end:-1:1));
            a.b0=[a.b0 -x(m2)];
            %if m2+1<=size(x,1) a.b0=[a.b0 -(x(m2)+x(m2+1))/2]; end;
         end
         end

        [r2]=test(a,a.dtst);
        [m1 c]=max(a.dtst.Y');
        [m1 c1]=max(r2.X');
        err=sum(c1~=c)/length(c);
        a.res=[a.res; loops err sum(sum(abs(newalpha)'>0)>0) ];
	 a.res(:,1)'
	 a.res(:,2)'
        disp(sprintf('svs: %d    err: %f ',sum(sum(abs(newalpha)'>0)>0),err));
     end

   else % only update one alpha
    for i=1:size(alpha,2)
     ai= (K(I,:)*alpha(:,i)) / K(I,I);   %   a_i = (x_i,w) / (x_i,x_i)  , optimal alpha
     alpha(I,i)=alpha(I,i)-ai;          %   w <- w - a_i x_i, project out learnt part of w
     newalpha(I,i)=newalpha(I,i)+ai;    %   wnew <- wnew + a_i x_i, add to new w
    end
   end

    w2=[]; for i=1:size(alpha,2) f=find(abs(alpha(:,i))>0); w2(i)=alpha(f,i)'*K(f,f)*alpha(f,i); end;
    a.w2=w2; % store w2
    %% this line could be sped up

   if mod(loops,25)==0
   txt='iteration %d : max||w_orig-w_new||^2=%1.3f svs=%d';
   disp(sprintf(txt,[loops max(w2) sum(sum(abs(newalpha)',1)>0)]));
   end;
   loops=loops+1;
end

a.alpha=newalpha;

if(isfield(struct(a.child),'X'))
    a.Xsv=get(a.child.X,origsvs);
else
    a.Xsv=get(a.child.Xsv,origsvs);
end

fin=find(sum(abs(a.alpha)',1)>a.alpha_cutoff);
a.alpha=a.alpha(fin,:);
a.Xsv=get(a.Xsv,fin);
a.b0=a.child.b0;

if a.reoptimize_b % reoptimize b // only for pat. rec. right now
   a.b0=0; r=test(a,d); a.b0=[];
   for i=1:size(newalpha,2)
    rr=r.X; rr=rr(:,i); [x s2]=sort(rr); y=r.Y(s2,i);
    xs=cumsum(y(end:-1:1)==1);
    [m1 m2]=max(cumsum(y==-1)+xs(end:-1:1));
    a.b0=[a.b0 -x(m2)];
    %if m2+1<=size(x,1) a.b0=[a.b0 -(x(m2)+x(m2+1))/2]; end;
   end
end

if a.algorithm.do_not_evaluate_training_error==1
  d=set_x(d,get_y(d));
else
  d=test(a,d);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [sortedBodyNames, sortedBodyInfo] = ma_getSortedBodyNames(celBodyData, varargin)
%ma_getSortedBodyNames Summary of this function goes here
%   Detailed explanation goes here

    if(isempty(varargin))
        curLevel = 1;
        maxLevel = Inf;
    elseif(length(varargin)==1)
        curLevel = 1;
        maxLevel = varargin{1};
    else
        curLevel = varargin{2};
        maxLevel = varargin{1};
    end

    if(curLevel == 1)
        topLevelBodyInfo = getTopLevelCentralBody(celBodyData);
        curParent = topLevelBodyInfo;
        sortedBodyNames = {curParent.name};
    else
        curParent = varargin{3};
        sortedBodyNames = {};
    end

    if(curLevel < maxLevel)
        numIndents = curLevel;
        childBodies = getChildrenOfParentInfo(celBodyData, curParent.name);
        childBodies = sortChildrenBySma(childBodies);
        for(i=1:length(childBodies)) %#ok<*NO4LP>
            childBody = childBodies{i};
            sortedBodyNames{end+1} = paddStrLeft(childBody.name, length(childBody.name)+numIndents); %#ok<AGROW>

            sortedChildBodyNames = ma_getSortedBodyNames(celBodyData, maxLevel, curLevel+1, childBody);
            for(j=1:length(sortedChildBodyNames))
                sortedBodyNames{end+1} = sortedChildBodyNames{j}; %#ok<AGROW>
            end
        end
    end

    sortedBodyNames = sortedBodyNames';
    sortedBodyInfo = cell(size(sortedBodyNames));
    for(i=1:length(sortedBodyNames))
        sortedBodyInfo{i} = celBodyData.(strtrim(lower(sortedBodyNames{i})));
    end
end

function smaSortedChildren = sortChildrenBySma(childBodies)
    [~,I] = sort(cellfun(@sortSmaCellFun, childBodies));
    smaSortedChildren = childBodies(I);
end

function sma = sortSmaCellFun(bodyInfo)
    sma = bodyInfo.sma;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ c, discrepancy ] = partition_brute ( n, w )

%*****************************************************************************80
%
%% PARTITION_BRUTE approaches the partition problem using brute force.
%
%  Discussion:
%
%    We are given a set of N integers W.
%
%    We seek to partition W into subsets W0 and W1, such that the subsets
%    have equal sums.
%
%    The "discrepancy" is the absolute value of the difference between the
%    two sums, and will be zero if we have solved the problem.
%
%    For a given set of integers, there may be zero, one, or many solutions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the set.
%
%    Input, integer W(N), the integers.
%
%    Output, integer C(N), indicates the proposed solution.
%    C(I) is 0 for items in set W0 and 1 for items in set W1.
%
%    Output, integer DISCREPANCY, the discrepancy.
%
  w = w(:);
  w_sum = sum ( w(1:n) );
  discrepancy = w_sum;

  d = [];
  rank = -1;

  while ( 1 )

    [ d, rank ] = subset_next ( n, d, rank );

    if ( rank == -1 )
      break
    end

    d_discrepancy = abs ( w_sum - 2 * d' * w );

    if ( d_discrepancy < discrepancy )
      discrepancy = d_discrepancy;
      c(1:n) = d(1:n);
    end

    if ( discrepancy == 0 )
      break
    end

  end

  return
end


text:
function result = t_most_3b_3_1_2(quiet)
%T_MOST_3B_1_1_2  Tests for MOST.

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

if nargin < 1
    quiet = 0;
end

n_tests = 51;

t_begin(n_tests, quiet);

casename = 't_case3_most';
fudging = struct( ...       %% paramters for fudging reserve contract for sopf2
    'fudge',    0.05, ...   %% initial value (MW)
    'step',     0.01, ...   %% if necessary, increase by this amount and retry (MW)
    'lim',      0.1);       %% upper limit (MW), give up if no convergence
                            %% with fudge equal to this limit

%% options
algs.dc     = {'MIPS'}; %% opf.dc.solver sequence to try for c3sopf (DC run)
algs.ac     = {'MIPS'}; %% opf.ac.solver sequence to try for c3sopf (AC run)
mpopt = mpoption('verbose', 0, 'out.all', 0);
mpopt = mpoption(mpopt, 'opf.violation', 5e-7, 'mips.comptol', 5e-8);
mpopt = mpoption(mpopt, 'sopf.force_Pc_eq_P0', 0);  %% don't constrain contracted == base case dispatch
mpoptac = mpoption(mpopt, 'model', 'AC');
mpoptdc = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'most.solver', algs.dc{1});
if have_feature('octave') && have_feature('octave', 'vnum') > 7
    mpopt = mpoption(mpopt, 'mips.linsolver', 'LU');
end

%% turn off warnings
if have_feature('octave')
    s = warning('query', 'Octave:nearly-singular-matrix');
    warning('off', 'Octave:nearly-singular-matrix');
else
    s7 = warning('query', 'MATLAB:nearlySingularMatrix');
    s6 = warning('query', 'MATLAB:nearlySingularMatrixUMFPACK');
    warning('off', 'MATLAB:nearlySingularMatrix');
    warning('off', 'MATLAB:nearlySingularMatrixUMFPACK');
end

%% define named indices into data matrices
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

%% reserve and delta offers
xgd_table.colnames = {
    'PositiveActiveReservePrice', ...
            'PositiveActiveReserveQuantity', ...
                    'NegativeActiveReservePrice', ...
                            'NegativeActiveReserveQuantity', ...
                                    'PositiveActiveDeltaPrice', ...
                                            'NegativeActiveDeltaPrice', ...
};
xgd_table.data = [
    1       400     2       400     0.01    0.01;
    3       300     4       300     0.01    0.01;
    0.001   450     0.002   450     0       0;
];

%% contingency table
% label probty  type        row column      chgtype newvalue
contab = [
    1   0.01    CT_TBRCH    1   BR_STATUS   CT_REP  0;      %% line 1-2
    2   0.01    CT_TGEN     1   GEN_STATUS  CT_REP  0;      %% gen 1 at bus 1
%   2   0.00049621965   CT_TGEN     1   GEN_STATUS  CT_REP  0;      %% gen 1 at bus 1
];
clist = unique(contab(:, CT_LABEL));
nc = length(clist);

%% load the case
mpc = loadcase(casename);
gbus = mpc.gen(:, GEN_BUS);

%%-----  get c3sopf results  -----
rdc = c3sopf_retry(algs.dc, mpc, xgd_table.data, contab, mpoptdc);
% rac = c3sopf_retry(algs.ac, mpc, xgd_table.data, contab, mpoptac);
% save t_most3_soln rdc rac -v6
% s = load('t_most3_soln');
s.rdc = rdc;
% s.rac = rac;

%%-----  set up data for DC run (most)  -----
ng = size(mpc.gen, 1);      %% number of gens
nt = 3;
xgd = loadxgendata(xgd_table, mpc);
md = loadmd(mpc, nt, xgd, [], contab);

%%-----  do DC run (most)  -----
r = most(md, mpopt);

%%-----  test the results  -----
t = 'success1';
t_ok(s.rdc.opf_results.success, t);
t = 'success2';
t_is(r.QP.exitflag, 1, 12, t);

t = 'f';
t_is(r.results.f/sum(r.StepProb), s.rdc.opf_results.f, 4, t);

for tt = 1:nt

t = sprintf('(t=%d) Pg : base', tt);
t_is(r.flow(tt,1,1).mpc.gen(:, PG), s.rdc.base.gen(:, PG), 5, t);
t = sprintf('(t=%d) Pg : cont ', tt);
for k = 1:nc
    t_is(r.flow(tt,1,k+1).mpc.gen(:, PG), s.rdc.cont(k).gen(:, PG), 5, sprintf('%s %d', t, k));
end

% t = sprintf('(t=%d) gen : base', tt);
% t_is(r.flow(tt,1,1).mpc.gen(:,1:MU_PMIN), s.rdc.base.gen(:,1:MU_PMIN), 3, t);
% t = sprintf('(t=%d) gen : cont ', tt);
% for k = 1:nc
%     t_is(r.flow(tt,1,k+1).mpc.gen(:,1:MU_PMIN), s.rdc.cont(k).gen(:,1:MU_PMIN), 3, sprintf('%s %d', t, k));
% end

t = sprintf('(t=%d) energy prices', tt);
t_is(r.results.GenPrices(:,tt)/r.StepProb(tt), s.rdc.energy.prc.sum_bus_lam_p(gbus), 6, t);

t = sprintf('(t=%d) Pc', tt);
t_is(r.results.Pc(:,tt), s.rdc.energy.Pc, 4, t);

t = sprintf('(t=%d) Gmin', tt);
t_is(r.results.Pc(:,tt) - r.results.Rpm(:,tt), s.rdc.energy.Gmin, 4, t);

t = sprintf('(t=%d) Gmax', tt);
t_is(r.results.Pc(:,tt) + r.results.Rpp(:,tt), s.rdc.energy.Gmax, 4, t);

t = sprintf('(t=%d) upward contingency reserve quantities', tt);
t_is(r.results.Rpp(:,tt), s.rdc.reserve.qty.Rp_pos, 4, t);

t = sprintf('(t=%d) downward contingency reserve quantities', tt);
t_is(r.results.Rpm(:,tt), s.rdc.reserve.qty.Rp_neg, 4, t);

t = sprintf('(t=%d) upward contingency reserve prices', tt);
t_is(r.results.RppPrices(:,tt)/r.StepProb(tt), s.rdc.reserve.prc.Rp_pos, 6, t);

t = sprintf('(t=%d) downward contingency reserve prices', tt);
t_is(r.results.RpmPrices(:,tt)/r.StepProb(tt), s.rdc.reserve.prc.Rp_neg, 6, t);

t = sprintf('(t=%d) contingency physical ramp price', tt);
[vv, ll] = r.om.get_idx();
Ramp_P_max = zeros(ng, nc);
sum_muPmax = zeros(ng, 1);
sum_muPmin = zeros(ng, 1);
for k = 1:nc+1
    ii = find(r.flow(tt,1,k).mpc.gen(:, GEN_STATUS) > 0);
    if k > 1
        Ramp_P_max(ii,k-1) = (r.QP.lambda.mu_u(ll.i1.rampcont(tt,1,k):ll.iN.rampcont(tt,1,k)) - r.QP.lambda.mu_l(ll.i1.rampcont(tt,1,k):ll.iN.rampcont(tt,1,k))) / mpc.baseMVA;
    end
    sum_muPmax(ii) = sum_muPmax(ii) + r.flow(tt,1,k).mpc.gen(ii, MU_PMAX);
    sum_muPmin(ii) = sum_muPmin(ii) + r.flow(tt,1,k).mpc.gen(ii, MU_PMIN);
end
t_is(Ramp_P_max/r.StepProb(tt), s.rdc.energy.mu.Ramp_P_max, 6, t);

t = sprintf('(t=%d) sum_muPmax', tt);
t_is(sum_muPmax/r.StepProb(tt), s.rdc.energy.sum_muPmax, 2, t);

t = sprintf('(t=%d) sum_muPmin', tt);
t_is(sum_muPmin/r.StepProb(tt), s.rdc.energy.sum_muPmin, 2, t);

t = sprintf('(t=%d) Rpmax_pos', tt);
Rpmax_pos = (r.QP.lambda.upper(vv.i1.Rpp(1):vv.iN.Rpp(1)) - r.QP.lambda.lower(vv.i1.Rpp(1):vv.iN.Rpp(1))) / mpc.baseMVA;
t_is(Rpmax_pos, s.rdc.reserve.mu.Rpmax_pos, 6, t);

t = sprintf('(t=%d) Rpmax_neg', tt);
Rpmax_neg = (r.QP.lambda.upper(vv.i1.Rpm(1):vv.iN.Rpm(1)) - r.QP.lambda.lower(vv.i1.Rpm(1):vv.iN.Rpm(1))) / mpc.baseMVA;
t_is(Rpmax_neg, s.rdc.reserve.mu.Rpmax_neg, 6, t);

end

% g1 = s.rdc.base.gen(:, PG);
% g2 = r.flow(1,1,1).mpc.gen(:, PG);
% for k = 1:nc
%     g1 = [ g1 s.rdc.cont(k).gen(:, PG) ];
%     g2 = [ g2 r.flow(1,1,k+1).mpc.gen(:, PG) ];
% end
% [m,n] = size(g1);
% for j = 1:n
%     fprintf('\n');
%     for i = 1:m
%         fprintf('%9.2f  %9.2f\n', g1(i,j), g2(i,j));
%     end
% end

%%-----  do AC run (most)  -----
%mostac;

%% turn warnings back on
if have_feature('octave')
    warning(s.state, 'Octave:nearly-singular-matrix');
else
    warning(s7.state, 'MATLAB:nearlySingularMatrix');
    warning(s6.state, 'MATLAB:nearlySingularMatrixUMFPACK');
end

t_end;

if nargout
    result = r;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [metabolite_structure] = queryLipidMaps(metabolite_structure,startSearch,endSearch)
%
% the function will search for metabolite names
%https://www.lipidmaps.org/search/quicksearch.php?Name=2-methyl-dodecanedioic+acid
%
%
% INPUT
% metabolite_structure  metabolite structure
% startSearch           specify where the search should start in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
% endSearch             specify where the search should end in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
%
% OUTPUT
% metabolite_structure  updated metabolite structure
%
%
% Ines Thiele, 09/2021

Mets = fieldnames(metabolite_structure);

if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(Mets);
end

annotationSource = 'Lipid Maps website by name search';
annotationType = 'automatic';

for i = startSearch : endSearch
    if isempty(metabolite_structure.(Mets{i}).lipidmaps) || length(find(isnan(metabolite_structure.(Mets{i}).lipidmaps)))>0
        metN = metabolite_structure.(Mets{i}).metNames;
        % make name fit for internet
        metN = regexprep(metN,'"','');
        metInt = regexprep(metN,' ','+') ;

        url=['https://www.lipidmaps.org/search/quicksearch.php?Name=' metInt];
        syst = urlread(url);
        lipidmaps = '';
        if exist('syst','var') && ~isempty(strfind(syst,'1 matches'))
            tok = split(syst,'/data/LMSDRecord.php');
            tok2 = split(tok{2},'>');
            tok3 = split(tok2{2},'<');
            lipidmaps = tok3{1};
            % no real information on there - could be useful for getting mol
            % files
            %url = ['https://www.lipidmaps.org/databases/lmsd/' lipidmaps];
            if ~isempty(lipidmaps)
                metabolite_structure.(Mets{i}).lipidmaps = lipidmaps;
                metabolite_structure.(Mets{i}).lipidmaps_source = [annotationSource,':',annotationType,':',':',datestr(now)];
            end
        end
    end
end


text:
function op = linop_TV3D( sz, variation, action )

%LINOP_TV3D   3D Total-Variation (TV) linear operator.
%    OP = LINOP_TV3D( SZ ) returns a handle to a TFOCS linear operator that
%      implements the total variation linear operator on an M x N x P grid;
%      that is, to be applied to volume stacks of size [M,N,P].
%      By default, it expects to operate on M*N*P x 1 vectors
%      but if SZ = {M,N,P}, then expects to operate on M x N x P matrices
%
%    TV = LINOP_TV3D( X ) returns ||X||_TV  if X is bigger than 2 x 2 x 2
%
%    OP = LINOP_TV3D( SZ, VARIANT )
%       if VARIANT is 'regular' (default),
%               ... TODO
%       if VARIANT is 'circular',
%               ... TODO
%
%    [...] = LINOP_TV3D(SZ, VARIATION, ACTION )
%       if ACTION is 'handle', returns a TFOCS function handle (default)
%       if ACTION is 'cvx', returns a function handle suitable for CVX
%       if ACTION is 'matrix', returns the explicit TV matrix
%           (real part corresponds to horizontal differences,
%            imaginary part correspond to vertical differences)
%       if ACTION is 'norm', returns an estimate of the norm
%
%   Contributed by  Mahdi Hosseini (mahdi.hosseini@mail.utoronto.ca)
%   Has not been extensively tested

error(nargchk(1,3,nargin));
if nargin < 2 || isempty(variation), variation = 'regular'; end
if nargin < 3 || isempty(action), action = 'handle'; end

CALCULATE_TV = false;
if numel(sz) > 6
    CALCULATE_TV = true;
    X   = sz;
    sz  = size(X);
end
nDim = numel(sz);

if iscell(sz)
    n1 = sz{1};
    n2 = sz{2};
    n3 = sz{3};
else
    n1 = sz(1);
    n2 = sz(2);
    n3 = sz(3);
end

% Setup the Total-Variation operators
mat = @(x) reshape(x,n1,n2,n3);

if strcmpi(action,'matrix') || strcmpi(action,'cvx')
    I1 = eye(n1);
    I2 = eye(n2);
    I3 = eye(n3);
    switch lower(variation)
        case 'regular'
            e = ones(max([n1,n2,n3]),1);
            e2 = e;
            e2(n2:end) = 0;
            J = spdiags([-e2,e], 0:1,n2,n2);
            Dh = kron(I3,kron(J,I1));  % horizontal differences, sparse matrix
            % see also blkdiag

            e2 = e;
            e2(n1:end) = 0;
            J = spdiags([-e2,e], 0:1,n1,n1);
            Dv = kron(I3,kron(I2,J));  % vertical differences, sparse matrix

            e2 = e;
            e2(n3:end) = 0;
            J = spdiags([-e2,e], 0:1,n3,n3);
            Dd = kron(J,kron(I2,I1));  % Depth differences, sparse matrix
        case 'circular'
            e = ones(max([n1,n2,n3]),1);
            e2 = e;
            %             e2(n2:end) = 0;
            J = spdiags([-e2,e], 0:1,n2,n2);
            J(end,1) = 1;
            Dh = kron(I3, kron(J,I1));  % horizontal differences, sparse matrix
            % see also blkdiag

            e2 = e;
            %             e2(n1:end) = 0;
            J = spdiags([-e2,e], 0:1,n1,n1);
            J(end,1) = 1;
            Dv = kron(I3, kron(I2,J));  % vertical differences, sparse matrix

            e2 = e;
            %             e2(n1:end) = 0;
            J = spdiags([-e2,e], 0:1,n3,n3);
            J(end,1) = 1;
            Dd = kron(J, kron(I2,I1));  % vertical differences, sparse matrix
    end
    if strcmpi(action,'matrix')
        op = [Dh;Dv;Dd];
    else
        % "norms" is a CVX function
        op = @(X) sum( norms( [Dh*X(:), Dv*X(:), Dd*X(:)]' ) );
    end
    return;
end

switch lower(variation)
    case 'regular'
        Dh     = @(X) vec( [diff(X,1,2), zeros(n1,1,n3)] );
        Dv     = @(X) vec( [diff(X,1,1); zeros(1,n2,n3)] );
        Dd     = @(X) [vec(diff(X,1,3)); vec(zeros(n1,n2,1))];

        diff_h = @(X) [zeros(n1,1,n3),X(:,1:end-1,:)] - ...
            [X(:,1:end-1,:),zeros(n1,1,n3)];
        diff_v = @(X) [zeros(1,n2,n3);X(1:end-1,:,:)] - ...
            [X(1:end-1,:,:);zeros(1,n2,n3)];
        diff_d = @(X) mat([vec(zeros(n1,n2,1));vec(X(:,:,1:end-1))] - ...
            [vec(X(:,:,1:end-1));vec(zeros(n1,n2,1))]);
    case 'circular'
        % For circular version, 2 x 2 case is special.
        %         error('not yet implemented');
        Dh     = @(X) vec( [diff(X,1,2), X(:,1,:) - X(:,end,:)] );
        Dv     = @(X) vec( [diff(X,1,1); X(1,:,:) - X(end,:,:)] );
        Dd     = @(X) [vec(diff(X,1,3)); vec(X(:,:,1) - X(:,:,end))];
        % diff_h needs to be checked
        diff_h = @(X) [X(:,end,:),X(:,1:end-1,:)] - X;
        % diff_v needs to be checked
        diff_v = @(X) [X(end,:,:);X(1:end-1,:,:)] - X;
        % diff_d needs to be checked
        diff_d = @(X) mat([vec(X(:,:,end));vec(X(:,:,1:end-1))]) - X;
    otherwise
        error('Bad variation parameter');
end
if iscell(sz)
    Dh_transpose = @(X)      diff_h(mat(X))  ;
    Dv_transpose = @(X)      diff_v(mat(X))  ;
    Dd_transpose = @(X)      diff_d(mat(X))  ;
else
    Dh_transpose = @(X) vec( diff_h(mat(X)) );
    Dv_transpose = @(X) vec( diff_v(mat(X)) );
    Dd_transpose = @(X) vec( diff_d(mat(X)) );
end

%%  TV & TVt Definitions
TV  = @(x) [Dh(mat(x)); Dv(mat(x)); Dd(mat(x))];

firstThird = @(x) x(1: n1*n2*n3);
secondThird= @(x) x(n1*n2*n3+1: 2*n1*n2*n3);
thirdThird= @(x) x(2*n1*n2*n3+1: 3*n1*n2*n3);
TVt = @(z) ( Dh_transpose(firstThird(z)) +...
    Dv_transpose(secondThird(z)) +...
    Dd_transpose(thirdThird(z)));

%%
if CALCULATE_TV
    op = norm( TV(X), 1 );
    return;
end

if strcmpi(action,'norm')
    % to compute max eigenvalue, I use a vector
    % that is very likely to be the max eigenvector:
    %  matrix with every entry alternating -1 and 1
    even = @(n) ~( n - 2*round(n/2) );  % returns 1 if even, 0 if odd
    Y = zeros( n1 + even(n1), n2 + even(n2), n3 + even(n3) );
    nn = numel(Y);
    Y(:) = (-1).^(1:nn);
    Y = Y(1:n1,1:n2,1:n3);
    op = norm( TV(Y) )/norm(Y(:));

    % Nearly equivalent to:
    % norm(full( [real(tv); imag(tv)] ) )
    % where tv is the matrix form
else
    if iscell(sz)
        szW = { [n1,n2,n3], [n1*n2*n3,1] };
    else
        szW = prod(sz); % n1 * n2
        szW = [nDim*szW, nDim*szW];
    end
    op = @(x,mode)linop_tv_r2c(szW,TV,TVt,x,mode);
end

function y = linop_tv_r2c( sz, TV, TVt, x, mode )
switch mode,
    case 0, y = sz;
    case 1, y = TV( realcheck( x ) );
    case 2, y = realcheck( TVt( x ) );
end

function y = realcheck( y )
if ~isreal( y ),
    error( 'Unexpected complex value in linear operation.' );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
%--------------------------------------------------------------------------
%   p3(signal,mark)
%--------------------------------------------------------------------------
%   功能：
%   复信号绘制画图工具，用处较多
%   常见于查看时域复信号，傅里叶变换后信号的幅度和相位信息
%--------------------------------------------------------------------------
%   输入：
%           signal              按照列排列复信号
%           mark                绘制标记，例如'r--','k-.'
%--------------------------------------------------------------------------
%   例子：
%   暂无
%--------------------------------------------------------------------------
function p3(signal,mark)
if nargin==2
    plot3(1:size(signal,1),real(signal),imag(signal),mark);grid on
    xlabel('点数');ylabel('实部');zlabel('虚部')
elseif nargin ==1
    plot3(1:size(signal,1),real(signal),imag(signal));grid on
    xlabel('点数');ylabel('实部');zlabel('虚部')
end

repo: SP