text:
%kPrintMapVal 'PrintMapVal Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros PrintMapVal.pane file
%
% Parameters:
% Integer: printmapvalWidth 'Width of Display (number of cells across)', default: 7: 'width of printmapval display'
% Integer: printmapvalHeight 'Height of Display (number of cells down) ', default: 9: 'height of printmapval display'
% InputFile: printmapvalFilename 'Input File', required: 'input file'
%
% Example: kPrintMapVal(printmapvalFilename, {'printmapvalWidth',7;'printmapvalHeight',9;'printmapvalFilename',''})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The PrintMapVal Visual Object"
% A printmapval object provides a mechanism with wich to roam around in an image
% and print the values of the map columns indexed by the pixel under the pointer.
% It is only useful with data objects having a map segment.  It is made up of one
% or more (width x height) sets of label objects;  typically, there will be three
% such sets, one each for the Red, Green, and Blue map columns of the colormap.
% The background of each label is the color defined by the values of the map
% columns at that pixel, and the value displayed is the value of the corresponding
% map column.
%
% .item "Width of Display"
% Specifies the number of map column values that should be
% displayed in the horizontal direction on the printmapval display.
%
% .item "Height of Display"
% Specifies the number of map column values that should be
% displayed in the vertical direction on the printmapval display.
%
% .item "Show Color"
% Indicates whether or not the color of the pixel under the
% pointer should be displayed as the background color for the label object in
% which the map values are displayed.  When set to FALSE, the background will be
% white.
%
% .item "Update Mode"
% Indicates whether the printmapval display is to be updated
% continuously as the pointer moves across the image, or not until the button
% is clicked at a particular location in the image.
%
% .item "Display Policy"
% Indicates whether the printmapval object is to display
% the values from the colormap that is being used to display the data object
% (ie, the RGB values that may have been converted and/or normalized), or
% the values from the actual map segment of the data object being visualized
% (ie, the map values before they are converted and/or normalized for use
% in defining the colors that appear on the screen).
%
% .item "Input File"
% The file containing the data object for which the map values are to be
% displayed.
% .end tagged

function varargout = kPrintMapVal(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kPrintMapVal(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'printmapvalWidth', 7;'printmapvalHeight', 9;'printmapvalFilename', '__input'};
maxval={2,2,0};
minval={2,2,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'Integer','Integer','InputFile'};
% 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 'PrintMapVal" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%IM_GRAY Fixed mapping converting multi-band images to gray images
%
%  B = IM_GRAY(A,V)
%  B = A*IM_GRAY([],V)
%  B = A*IM_GRAY(V)
%
% INPUT
%   A    Multiband image or dataset with multi-band images as objects
%   V    Weight vector, one weight per band. Default: equal weights.
%        (weights will be normalized to sum to one)
%
% OUTPUT
%   B    Output image or dataset.
%
% DESCRIPTION
% The multi-band components in the image A (3rd dimension) or in the
% objects in the dataset A are weigthed (default: equal weights) and
% averaged.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, DATAFILES

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = im_gray(varargin)

	argin = shiftargin(varargin,'vector');
  argin = setdefaults(argin,[],[]);
  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');
    b = setname(b,'Color-to-gray conversion');
  else
    [a,v] = deal(argin{:});
    if isa(a,'prdataset') % allows datafiles too
      isobjim(a);
      b = filtm(a,mfilename,v);
      imsize = getfeatsize(a);
      b = setfeatsize(b,imsize(1:2));
    else
      a = double(a);
      imsize = size(a);
      if isempty(v)
        if length(imsize) == 3
          b = mean(a,3);
          b = squeeze(b);
        elseif length(imsize) == 2
          b = a;
        else
          error('Illegal image size')
        end
      else
        if length(imsize) == 2
          b = a;
        else
          if length(v) ~= imsize(3)
            error('Number of weights should be equal to number of bands')
          end
          v = v/sum(v);
          b = zeros(imsize(1),imsize(2));
          for j=1:size(a,3)
            b = b + a(:,:,j)*v(j);
          end
        end
      end
    end
  end
return

author: marianux

repo: ecg-kit

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


text:
function viewNeighborsAuto(inputImages, inputSignals, neighborsCell, varargin)
	% View the neighboring cells, their traces and trace correlations.
	% Biafra Ahanonu
	% started 2013.11.01
	% inputs
		% inputImages - [x y nCells] matrices containing each set of filters
		% inputSignals - [nFilters frames] matrices containing each set of filter traces
		% neighborsCell - {nCells 1} cell array of [nNeighborId 1] vectors of neighbor IDs for each cell matching indices in inputImages
	% outputs
		%
	% changelog
		% 2018 - updated interface, speed improvements.
		% 2019.05.28 [13:49:31] - No longer use montage to create montage, use loops plus cat and cell arrays to improve speed, reduce use of montage(), and make outlines crisp (on 2018b they would become expanded with montage).
		% 2019.05.28 [15:15:43] - General improvements to GUI clarity.
		% 2019.06.18 [14:54:34] - Added support for spatial correlations
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	options.plottingOn = 0;
	options.overlapradius = 10;
	% overlap radius in pixels (so convert properly)
	options.overlapDistance = 10;
	% 'centroid', 'imdilate'
	options.neighborMethod = 'centroid';
	% Input pre-computed x,y coordinates for objects in images
	options.xCoords = [];
	options.yCoords = [];
	%
	options.inputImagesThres = [];
	%
	options.startCellNo = 1;
	%
	options.cropSizeLength = 20;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	%     eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	if isempty(neighborsCell)
		copts = options;
		neighborsCell = identifyNeighborsAuto(inputImages, inputSignals,'options',copts);
	end

	% cell vector info
	cellvec(1) = 1;
	cellvec(2) = size(inputImages,3);
	valid = zeros(1,cellvec(2));

	openFigure(21,'full');
		clf

	exitLoop = 0;
	nCells = cellvec(2);
	cellnum = options.startCellNo;
	cropSizeLength = options.cropSizeLength;
	usePadding = 1;
	xCoords = [];
	yCoords = [];

	rowP = 3;
	colP = 3;

	% instructions
	instructionStr =  ['Controls: left/right arrow:forward/back | e:exit | g:(goto cell #)' 10 ' cell #s low-high from top to bottom | top trace = this cell' 10];

	set(21, 'KeyPressFcn', @(source,eventdata) figure(21));
	while exitLoop==0
		directionOfNextChoice = 0;
		thisFilt = squeeze(inputImages(:,:,cellnum));
		if issparse(thisFilt)
			thisFilt = full(thisFilt);
		end
		nList = [cellnum neighborsCell{cellnum,1}(:)'];
		% nList
		% montage of nearby cells
		subplot(rowP,colP,1)
			rgbImage = [];
			nImgs = inputImages(:,:,nList);
			if issparse(nImgs)
				nImgs = full(nImgs);
			end
			nMatchCells = size(nImgs,3);

			objCutMovie = getObjCutMovie(nImgs,thisFilt,'createMontage',0,'extendedCrosshairs',0,'crossHairVal',NaN,'outlines',1,'waitbarOn',0,'cropSize',cropSizeLength,'addPadding',usePadding,'xCoords',xCoords,'yCoords',yCoords,'outlineVal',NaN);
			objCutMovie = vertcat(objCutMovie{:});
			[thresholdedImages boundaryIndices] = thresholdImages(objCutMovie(:,:,1),'binary',1,'getBoundaryIndex',1,'threshold',0.35,'imageFilter','','waitbarOn',0,'normalizationType','zeroToOne','removeUnconnected',1,'boundaryHoles','noholes');

			for zz = 1:size(nImgs,3)
				tmpImg = squeeze(objCutMovie(:,:,zz));
				tmpImg([boundaryIndices{:}]) = NaN;

				% tmpImg = squeeze(nanmean(...
				%     insertText(tmpImg,[0 0],num2str(nList(zz)),...
				%     'BoxColor',[0 0 0],...
				%     'TextColor',[1e10 1e10 1e10],...
				%     'AnchorPoint','LeftTop',...
				%     'FontSize',8,...
				%     'BoxOpacity',0)...
				% ,3));
				% tmpImg(tmpImg>1e7) = NaN;
				% figure;imagesc(tmpImg)
				objCutMovie(:,:,zz) = tmpImg;
			end
			clear objCutMovie2;
			objCutMovie2(:,:,1,:) = objCutMovie;
			% montage(objCutMovie2)

			croppedPeakImagesCell = {};
			[xPlot yPlot] = getSubplotDimensions(size(objCutMovie,3));
			for iii = 1:xPlot*yPlot
				if iii>size(objCutMovie,3)
					croppedPeakImagesCell{iii} = NaN(size(objCutMovie(:,:,1)));
				else
					croppedPeakImagesCell{iii} = objCutMovie(:,:,iii);
				end
			end
			mNum = 1;
			g = cell([xPlot 1]);
			rowNum = 1;

			xyLoc = ones([nMatchCells 2]);
			xyLoc(:,2) = [3];
			xyLoc(:,1) = [2];

			cellNoHere = 1;
			for xNo = 1:xPlot
				for yNo = 1:yPlot
					g{rowNum} = cat(2,g{rowNum},croppedPeakImagesCell{mNum});

					if sum(isnan(croppedPeakImagesCell{mNum}(:)))~=numel(croppedPeakImagesCell{mNum}(:))&&cellNoHere>1&&yNo>1
						xyLoc(cellNoHere,:) = [xyLoc(cellNoHere-1,1)+options.cropSizeLength*2+1 xyLoc(cellNoHere,2)];
					end
					cellNoHere = cellNoHere+1;
					mNum = mNum + 1;
				end

				xyLoc(cellNoHere:end,2) = [xyLoc(cellNoHere:end,2)+options.cropSizeLength*2+1];
				xyLoc(cellNoHere:end,1) = [2];
				rowNum = rowNum + 1;
			end
			croppedPeakImages2 = cat(1,g{:});

			% croppedPeakImages2 = getimage;
			imAlpha = ones(size(croppedPeakImages2));
			imAlpha(isnan(croppedPeakImages2))=0;
			imagesc(croppedPeakImages2,'AlphaData',imAlpha);
			set(gca,'color',[0 0 0]);
			colormap(gca,customColormap([]))
			title(['neighboring cell ID #s (white)' 10 'top-left = selected cell' 10 '#s increase L->R and T->B'])
			% title(['neighboring cells | top-left = selected cell, #s increase L->R and T->B'])
			axis equal tight
			box off

			% xyLoc
			for zz = 1:size(nImgs,3)
				text(xyLoc(zz,1),xyLoc(zz,2),num2str(nList(zz)),'Color',[1 1 1],'FontSize',options.cropSizeLength/2)
			end

		% plot the traces
		subplot(rowP,colP,[2:3 5:6 8:9])
			% nList = [cellnum neighborsCell{cellnum,1}(:)'];
			% plotSignalsGraph(inputSignals, 'plotList', flipdim(nList(:),1),'newAxisColorOrder','','minAdd',0,'maxIncrementPercent',1.1);
			plotSignalsGraph(inputSignals, 'plotList', flipdim(nList(:),1),'newAxisColorOrder','lines','minAdd',0,'maxIncrementPercent',1.1);
			axis tight
			xlabel('Time (frames)')
			ylabel('Cell activity')
			title('Activity traces for selected and neighboring cells (Zoom enabled)')
			legend(strsplit(num2str(nList),' '),'Location','southoutside','Orientation','horizontal')
			% ciapkg.overloaded.suptitle(instructionStr);
			zoom on

		% plot temporal correlations
		subplot(rowP,colP,4)
			z0 = corr(inputSignals(nList(:),:)');
			imagesc(z0)
			colorbar
			colormap(gca,customColormap([]))
			% colorbar; colormap jet;
			% colormap hot

			% imagesc(l);
			xlabel('Cell ID #s');
			ylabel('Cell ID #s');
			xticks(1:length(nList))
			xticklabels(nList)
			yticks(1:length(nList))
			yticklabels(nList)
			title(['cell-cell activity trace correlations' 10 'cell #' num2str(nList(1)) ' = this cell']);
			box off; axis equal tight
			caxis([0 1])
			% datacursormode on

		% plot spatial correlations
		subplot(rowP,colP,7)

			nCellsH = length(nList);
			fracMaxVal = 0.2;
			tmpImgs = inputImages(:,:,nList);
			if issparse(tmpImgs)
				tmpImgs = full(tmpImgs);
			end
			for cInd1 = 1:nCellsH
				tmpImgs1 = tmpImgs(:,:,cInd1);
				tmpImgs1(tmpImgs1<(nanmax(tmpImgs1(:))*fracMaxVal)) = 0;
				tmpImgs(:,:,cInd1) = tmpImgs1;
			end

			cEst = permute(tmpImgs, [3 1 2]);
			cEst = reshape(cEst, [size(cEst,1), numel(inputImages(:,:,1))]);
			% figure;imagesc(cEst);pause
			corrImgMatrix = 1-squareform(pdist(single(cEst>0),'jaccard'));

			spatialCorr = z0*0;
			spatialCorr2 = z0*0;
			for cInd1 = 1:nCellsH
				cImg1 = tmpImgs(:,:,cInd1);
				cImg1_2 = inputImages(:,:,nList(cInd1));
				for cInd2 = 1:nCellsH
					if cInd2==cInd1
						spatialCorr(cInd1,cInd2) = 1;
					end
					if cInd2>cInd1
						continue;
					end
					cImg2 = tmpImgs(:,:,cInd2);
					cImg2_2 = inputImages(:,:,nList(cInd2));

					% cImg1(cImg1<(nanmax(cImg1(:))*fracMaxVal)) = 0;
					% cImg2(cImg2<(nanmax(cImg2(:))*fracMaxVal)) = 0;

					c1Pixels = find(cImg1>0);
					c2Pixels = find(cImg2>0);
					thisOverlap=length(intersect(c1Pixels,c2Pixels));
					thisOverlap1=thisOverlap/length(c1Pixels);
					thisOverlap2=thisOverlap/length(c2Pixels);

					spatialCorr(cInd1,cInd2) = thisOverlap1;
					spatialCorr(cInd2,cInd1) = thisOverlap2;

					c1Pixels = find(cImg1_2>0);
					c2Pixels = find(cImg2_2>0);
					thisOverlap=length(intersect(c1Pixels,c2Pixels));
					thisOverlap1=thisOverlap/length(c1Pixels);
					thisOverlap2=thisOverlap/length(c2Pixels);

					spatialCorr2(cInd1,cInd2) = thisOverlap1;
					spatialCorr2(cInd2,cInd1) = thisOverlap2;
				end
			end
			imagesc(cat(2,spatialCorr,corrImgMatrix,spatialCorr2));
			colorbar
			colormap(gca,customColormap([]))
			% colorbar; colormap jet;
			% colormap hot

			% imagesc(l);
			xlabel('Cell ID #s');
			ylabel('Cell ID #s');
			nListLen = length(nList);
			xticks(1:(nListLen*3))
			xticklabels([nList; nList; nList])
			yticks(1:length(nList))
			yticklabels(nList)
			title(['cell-cell image correlations']);
			box off; axis equal tight
			hold on;
			% plot([nListLen+0.5 nListLen+0.5],get(gca,'YLim'),'k-','LineWidth',1);
			% plot([nListLen*2+0.5 nListLen*2+0.5],get(gca,'YLim'),'k-','LineWidth',1);
			hold off;
			% datacursormode on
			caxis([0 1])

		ciapkg.overloaded.suptitle([sprintf('cell # %d/%d',cellnum,nCells) 10 instructionStr],'plotregion',0.9);
		set(findall(gcf,'-property','FontSize'),'FontSize',13);
		% [x,y,reply]=ginput(1);

		set(gcf,'currentch','3');
		keyIn = get(gcf,'CurrentCharacter');
		figure(21)

		while strcmp(keyIn,'3')
			keyIn = get(gcf,'CurrentCharacter');
			pause(0.05);
		end
		reply = double(keyIn);
		set(gcf,'currentch','3');

		[valid directionOfNextChoice exitLoop cellnum] = respondToUserInput(reply,cellnum,valid,directionOfNextChoice,exitLoop,nCells);

		cellnum=cellnum+directionOfNextChoice;
		% loop if user gets to either end
		if cellnum<=0
			cellnum = nCells;
		elseif i>nCells;
			cellnum = 1;
		end
	end
	close(21)

end

function [valid directionOfNextChoice exitLoop i] = respondToUserInput(reply,i,valid,directionOfNextChoice,exitLoop,nCells)
	% decide what to do based on input (not a switch due to multiple comparisons)
	if isequal(reply, 3)|isequal(reply, 110)|isequal(reply, 31)
		% n key or right click
		directionOfNextChoice=1;
		% display('invalid IC');
		% set(fig1,'Color',[0.8 0 0]);
		valid(i) = 0;
	elseif isequal(reply, 28)
		% go back, left
		directionOfNextChoice=-1;
	elseif isequal(reply, 29)
		% go forward, right
		directionOfNextChoice=1;
	elseif isequal(reply, 101)
		exitLoop=1;
		% user clicked 'f' for finished, exit loop
		% i=nCells+1;
	elseif isequal(reply, 103)
		% if user clicks 'g' for goto, ask for which IC they want to see
		icChange = inputdlg('enter signal #'); icChange = str2num(icChange{1});
		if icChange>nCells|icChange<1
			% do nothing, invalid command
		else
			i = icChange;
			directionOfNextChoice = 0;
		end
	elseif isequal(reply, 115)
		% 's' if user wants to get ride of the rest of the ICs
		display(['classifying the following signals as bad: ' num2str(i) ':' num2str(nCells)])
		valid(i:nCells) = 0;
		exitLoop=1;
	elseif isequal(reply, 121)|isequal(reply, 1)|isequal(reply, 30)
		% y key or left click
		directionOfNextChoice=1;
		% display('valid IC');
		% set(fig1,'Color',[0 0.8 0]);
		valid(i) = 1;
	else
		% forward=1;
		% valid(i) = 1;
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function SE = functionPowerOptimization_prodSINR(signal,interference,Pmax,prelogFactor)
%Compute DL power allocation that solves the max product SINR problem,
%using the algorithm in Theorem 7.2.
%
%This function require additional software packages to be used, which
%need to be downloaded and installed separately. These packages are
%developed independently and are delivered with separate licenses.
%The implementation uses CVX (http://cvxr.com/cvx) and has been tested
%using CVX version 2.1. We recommend the use of the Mosek solver (we
%have tested using version 7.1.0.12).
%
%INPUT:
%signal       = K x L matrix where element (k,j) is a_jk in (7.2)
%interference = K x L x K x L matrix where (l,i,j,k) is b_lijk in (7.3)
%Pmax         = Maximum transmit power per BS
%prelogFactor = Prelog factor
%
%OUTPUT:
%SE = K x L matrix where element (k,j) is the downlink SE of UE k in cell j
%     using the max product power allocation solution
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.1 (Last edited: 2018-06-08)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Extract number of UEs
K = size(signal,1);

%Extract number of cells
L = size(signal,2);

%% Solve geometric program in (7.8) using CVX
cvx_begin gp
cvx_quiet(true); % This suppresses screen output from the solver

variable rho(K,L);
variable c(K,L);

maximize prod(prod(c))

subject to

for j = 1:L

    for k = 1:K

        if signal(k,j)>0
            %SINR constraints of UE k in cell j
            c(k,j)*(sum(sum(rho.*interference(:,:,k,j))) + 1) <= (rho(k,j)*signal(k,j));

            rho(k,j) >= 0;

        else
            %This applies if UE k in cell j is inactive
            c(k,j) == 1;
            rho(k,j) >= 0;

        end

    end

    sum(rho(:,j)) <= Pmax;

end

cvx_end

%% Analyze the CVX output and prepare the output variables
if isempty(strfind(cvx_status,'Solved')) %The problem was not solved by CVX, for some reason, and we then consider equal power allocation
    rhoSolution = (Pmax/K)*ones(K,L);
else %The problem was solved by CVX
    rhoSolution = rho;
end

%Refine the solution obtained from CVX using the Matlab command fmincon.
%This is particularly important in case CVX fails to solve the problem
A = kron(eye(L),ones(1,K));
B = Pmax*ones(L,1);
options = optimoptions('fmincon','Algorithm','interior-point','MaxFunEvals',50000,'MaxIter',5000);
xend = fmincon(@(x) -functionComputesumSE_DL_poweralloc(x,signal,interference,prelogFactor),rhoSolution(:),A,B,[],[],zeros(K*L,1),[],[],options);
rhoBest = reshape(xend,[K L]);

%Compute the SEs using Theorem 4.6
SE = functionComputeSE_DL_poweralloc(rhoBest,signal,interference,prelogFactor);

function sumSE = functionComputesumSE_DL_poweralloc(rho,signal,interference,prelogFactor)

%Extract number of UEs
K = size(signal,1);

%Extract number of cells
L = size(signal,2);

%Reshape power variable since fmincon optimizes vectors
rho = reshape(rho,[K L]);

%Prepare to save results
SE = zeros(K,L);

for j = 1:L

    for k = 1:K

        if signal(k,j) > 0 %Check if the UE k in cell j is active

            %Compute the SE in Theorem 4.6 using the formulation in (7.1)
            SE(k,j) = prelogFactor*log2((rho(k,j)*signal(k,j)) / (sum(sum(rho.*interference(:,:,k,j))) + 1));

        else %If the UE is inactive

            SE(k,j) = 0;

        end

    end

end

%Compute the sum SE of all cells
sumSE = sum(SE(:));

author: emilbjornson

repo: massivemimobook


text:
function [G]=gsp_spectrum_cdf_approx(G,param)
%GSP_SPECTRUM_CDF_APPROX  Compute an approximation of the cumulative density function of the graph Laplacian eigenvalues
%   Usage:  G=gsp_spectrum_cdf_approx(G);
%           G=gsp_spectrum_cdf_approx(G,param);
%
%   Input parameters:
%         G                     : Graph structure.
%   Output parameters:
%         G                     : Graph structure, including the addition of G.spectral_warp_fn
%         param                 : Structure of additional parameter
%   Additional parameters:
%         param.num_pts         : Number of interpolation points
%         param.use_speedup     : Only perform step 1 (ldlsymbol) of the ldl algorithm once and then repeat the numeric step 2
%         param.use_permutation : Use a sparsity-preserving permutation
%         param.use_ldl_package : If the mex files for Timothy Davis' ldlsparse package are present. Otherwise use MATLAB's 'ldl' function
%         param.ldl_thresh      : Threshold parameter for MATLAB 'ldl' function
%
%
%   'gsp_spectrum_cdf_approx(G)' computes an approximation of the
%   cumulative density function of the graph Laplacian eigenvalues using
%   spectrum slicing techniques. The algorithm is as follows:
%
%   Step 1: Take $Q+1$ (param.num_pts) evenly spaced points on the interval
%   $[0,\lambda_{\max}]$.
%
%   Step 2: For every $q \in \{1,2,\ldots,Q\}$, compute a triangular
%   factorization:
%
%   .. \{\cal L}-\frac{q \lambda_{\max}}{Q}I = L_q \Delta_q L_q^*
%
%   .. math:: {\cal L}-\frac{q \lambda_{\max}}{Q}I = L_q \Delta_q L_q^*
%
%   where $\Delta_q$ is a diagonal matrix and $L_q$ is a lower triangular
%   matrix. By a corollary of Sylverster's law of inertia, the number of
%   negative eigenvalues of the diagonal matrix $\Delta_q$ is equal to the
%   number of eigenvalues of ${\cal L}$ less than $q \lambda_{\max} / Q$.
%
%   Step 3: Use monotonic cubic polynomial interpolation to interpolate the
%   points
%
%   .. { ( q\lambda_{max}/Q , mu_q/(N-1) ) }_q=0,1,\ldots,Q
%
%   .. math:: \left\{\left(\frac{q \lambda_{upper}}{Q},\frac{\mu_q}{N-1} \right)\right\}_{q=0,1,\ldots,Q}
%
%   where $\mu_q$ is the number of diagonal elements of $\Delta_q$ less than zero.
%
%   To perform the triangular factorizations, we use the LDL sparse
%   Cholesky package, written by Timothy Davis. The option
%   param.use_permutation enables a sparsity preserving permutation with
%   the function 'symamd'. If the option param.use_speedup is set to 1, the
%   symbolic analysis step of the LDL algorithm is only performed once for
%   all shifts.
%
%   See also:
%
%   Demos:
%
%   Requires: ldlsymbol_extra, ldlnumeric, ldlsparse
%
%   References: parlett1998symmetric davis2005ldl davis2011ldl
%

%TODO this function need to be cleaned

%   AUTHOR : David I Shuman.
%   TESTING:
%   REFERENCE:

if isfield(G, 'spectral_warp_fn')
    warning('Overwriting spectral warping function');
end

% if isfield(G,'e');
%     [x,y] = gsp_point2dcdf(G.e);
%     G.spectrum_cdf_approx = @(s) gsp_mono_cubic_warp_fn(x,y,s);
%     return
% end

if nargin<2
    param = struct;
end

if ~isfield(param, 'num_pts')
    num_pts=25;
else
    num_pts = param.num_pts;
end

if ~isfield(param, 'use_speedup')
    if ( exist('ldlsymbol_extra','file')==3 && exist('ldlnumeric','file')==3 )
        use_speedup=1;
    else
        use_speedup=0;
    end
else
    use_speedup = param.use_speedup;
end

if ~isfield(param, 'use_ldl_package')
    if exist('ldlsparse','file')==3
        use_ldl_package=1;
    else
        use_ldl_package=0;
    end
else
    use_ldl_package = param.use_ldl_package;
end

if ~isfield(param, 'use_permutation')
    use_permutation=1;
else
    use_permutation = param.use_permutation;
end

if ~isfield(G, 'lmax')
    warning(['GSP_SPECTRUM_CDF_APPROX: The variable lmax is not ',...
            'available. The function will compute it for you. ',...
            'However, if you apply many time this function, you ',...
            'should precompute it using the function: ',...
            'gsp_estimate_lmax']);
    G = gsp_estimate_lmax(G);
end

if ~isfield(param, 'ldl_thresh')
    ldl_thresh=.001;
else
    ldl_thresh = param.ldl_thresh;
end

counts=zeros(num_pts,1);
counts(num_pts)=G.N-1;

interp_x=(0:num_pts-1)*G.lmax/(num_pts-1);
interp_x=interp_x';

identity=speye(G.N);

if use_speedup
    if use_permutation
        P=symamd(G.L);
        [Parent, Lp, PO, PIn, flopcount] = ldlsymbol_extra(G.L,P);
        for i=1:num_pts-2
            mat=G.L-interp_x(i+1)*identity;
            [~, HD]=ldlnumeric(mat,Lp,Parent,PO,PIn);
            counts(i+1)=sum(diag(HD)<0);
        end
    else
        [Parent, Lp, flopcount] = ldlsymbol_extra(G.L);
        for i=1:num_pts-2
            mat=G.L-interp_x(i+1)*identity;
            [~, HD]=ldlnumeric(mat,Lp,Parent);
            counts(i+1)=sum(diag(HD)<0);
        end
    end
else
    if use_permutation
        P=symamd(G.L);
        for i=2:num_pts-1
            mat=G.L-interp_x(i)*identity;
            if use_ldl_package
                [~,HD]=ldlsparse(mat,P);
            else
                [~,HD,~]=ldl(mat(P,P),ldl_thresh);
            end
            counts(i)=sum(diag(HD)<0);
        end
    else
        for i=2:num_pts-1
            mat=G.L-interp_x(i)*identity;
            if use_ldl_package
                [~,HD]=ldlsparse(mat);
            else
                [~,HD,~]=ldl(mat,ldl_thresh);
            end
            counts(i)=sum(diag(HD)<0);
        end
    end
end
interp_y=counts/(G.N-1);

G.spectrum_cdf_approx = @(s) gsp_mono_cubic_warp_fn(interp_x,interp_y,s);

end

repo: gspbox


text:
function varargout = pivots(varargin)
%PIVOTS   Pivot values of a CHEBFUN2.
%
%   PIVOTS(F) returns the pivots used during the construction by the GE
%   algorithm.
%
%   PIVOTS(F, 'normalize'), returns the normalised pivot values. These numbers
%   are scaled so that the columns and rows have unit 2-norm.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function quadrule_test ( )

%*****************************************************************************80
%
%% QUADRULE_FAST_TEST tests the QUADRULE_FAST library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_FAST_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the QUADRULE_FAST library.\n' );

  quadrule_fast_test01 ( );
  quadrule_fast_test02 ( );
  quadrule_fast_test03 ( );

  quadrule_fast_test04 ( );
  quadrule_fast_test05 ( );
  quadrule_fast_test06 ( );

  quadrule_fast_test07 ( );
  quadrule_fast_test08 ( );
  quadrule_fast_test09 ( );

  quadrule_fast_test10 ( );
  quadrule_fast_test11 ( );
  quadrule_fast_test12 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_FAST_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [bnet, names] = mk_bat_dbn()
% MK_BAT_DBN Make the BAT DBN
% [bnet, names] = mk_bat_dbn()
% See
% - Forbes, Huang, Kanazawa and Russell, "The BATmobile: Towards a Bayesian Automated Taxi", IJCAI 95
% - Boyen and Koller, "Tractable Inference for Complex Stochastic Processes", UAI98.

names = {'LeftClr', 'RightClr', 'LatAct', 'Xdot', 'InLane', 'FwdAct', ...
      'Ydot', 'Stopped', 'EngStatus', 'FBStatus', ...
      'LeftClrSens', 'RightClrSens', 'TurnSignalSens', 'XdotSens', 'YdotSens', ...
      'FYdotDiffSens', 'FclrSens', 'BXdotSens', 'BclrSens', 'BYdotDiffSens', ...
      'SensorValid', 'FYdotDiff', 'FcloseSlow', 'Fclr', 'BXdot', 'BcloseFast', 'Bclr', 'BYdotDiff'};
ss = length(names);

intrac = {...
      'LeftClr', 'LeftClrSens';
  'RightClr', 'RightClrSens';
  'LatAct', 'TurnSignalSens'; 'LatAct', 'Xdot';
  'Xdot', 'XdotSens';
  'FwdAct', 'Ydot';
  'Ydot', 'YdotSens'; 'Ydot', 'Stopped';
  'EngStatus', 'Ydot'; 'EngStatus', 'FYdotDiff'; 'EngStatus', 'Fclr'; 'EngStatus', 'BXdot';
  'SensorValid', 'XdotSens';   'SensorValid', 'YdotSens';
  'FYdotDiff', 'FYdotDiffSens'; 'FYdotDiff', 'FcloseSlow';
  'FcloseSlow', 'FBStatus';
  'Fclr', 'FclrSens'; 'Fclr', 'FcloseSlow';
  'BXdot', 'BXdotSens';
  'Bclr', 'BclrSens'; 'Bclr', 'BXdot'; 'Bclr', 'BcloseFast';
  'BcloseFast', 'FBStatus';
  'BYdotDiff', 'BYdotDiffSens'; 'BYdotDiff', 'BcloseFast'};
[intra, names] = mk_adj_mat(intrac, names, 1);

interc = {...
      'LeftClr', 'LeftClr'; 'LeftClr', 'LatAct';
  'RightClr', 'RightClr'; 'RightClr', 'LatAct';
  'LatAct', 'LatAct'; 'LatAct', 'FwdAct';
  'Xdot', 'Xdot'; 'Xdot', 'InLane';
  'InLane', 'InLane'; 'InLane', 'LatAct';
  'FwdAct', 'FwdAct';
  'Ydot', 'Ydot';
  'Stopped', 'Stopped';
  'EngStatus', 'EngStatus';
  'FBStatus', 'FwdAct'; 'FBStatus', 'LatAct'};
inter = mk_adj_mat(interc, names, 0);

obs = {'LeftClrSens', 'RightClrSens', 'TurnSignalSens', 'XdotSens', 'YdotSens', 'FYdotDiffSens', ...
      'FclrSens', 'BXdotSens', 'BclrSens', 'BYdotDiffSens'};

for i=1:length(obs)
  onodes(i) = strmatch(obs{i}, names); %stringmatch(obs{i}, names);
end
onodes = sort(onodes);

dnodes = 1:ss;
ns = 2*ones(1,ss); % binary nodes
bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'observed', onodes, 'eclass2', (1:ss)+ss);

% make rnd params
for i=1:2*ss
  bnet.CPD{i} = tabular_CPD(bnet, i);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% std_renamestudyfiles() - rename files for design 1 if necessary. In design
%                          1, for backward compatibility, files could have
%                          legacy names. For consistency these files now
%                          need to be renamed. Note that the STUDY is
%                          automatically resave on disk to avoid any potential
%                          inconsistency.
%
% Usage:
%   >> STUDY = std_renamestudyfiles(STUDY, ALLEEG)
%
% Inputs:
%   STUDY      - EEGLAB STUDY set
%   ALLEEG     - vector of the EEG datasets included in the STUDY structure
%
% Output:
%      STUDY - The input STUDY with new design files.
%
% Author: Arnaud Delorme, Institute for Neural Computation UCSD, 2013-

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

function STUDY = std_renamestudyfiles(STUDY, ALLEEG)

if nargin < 2
    help std_renamestudyfiles;
    return;
end;

STUDY2 = std_makedesign(STUDY, ALLEEG, 1, STUDY.design(1), 'defaultdesign', 'forceoff', 'verbose', 'off');
allCell1 = {  STUDY.design(1).cell.filebase };
allCell2 = { STUDY2.design(1).cell.filebase };

fileExtensions = { 'daterp' 'datspec' 'datersp' 'daterpim' 'dattimef' 'datitc' 'daterpim' ...
    'icaerp' 'icaspec' 'icaersp' 'icaerpim' 'icatimef' 'icaitc' 'icaerpim' };

if ~isequal(allCell1, allCell2)

    thereIsAFileNotDesign = false;
    for index = 1:length(allCell1), if length(allCell1{index}) < 6 || all(allCell1{index}(1:6) == 'design'), thereIsAFileNotDesign = true; end; end;
    for index = 1:length(allCell1), if length(allCell1{index}) < 6 || all(allCell1{index}(1:6) == 'design'), thereIsAFileNotDesign = true; end; end;

    if thereIsAFileNotDesign
        res = questdlg2(['Old STUDY design data files have been detected.' 10 ...
            'EEGLAB wants to rename these files to improve consistency' 10 ...
            'and stability. No dataset will be renamed, only preprocessed' 10 ...
            'STUDY data files.' ], 'Rename STUDY data files', 'Cancel', ...
            'Rename', 'Rename');

        if strcmpi(res, 'rename')
            STUDY = pop_savestudy(STUDY, ALLEEG, 'savemode', 'resave');
            for iCell = 1:length(allCell1)

                % scan file extensions
                for iExt = 1:length(fileExtensions)
                    files = dir( [ allCell1{iCell} '.' fileExtensions{iExt} ]);
                    if ~isempty(files) && ~strcmpi(allCell1{iCell}, allCell2{iCell})
                        movefile( [ allCell1{iCell} '.' fileExtensions{iExt} ], ...
                            [ allCell2{iCell} '.' fileExtensions{iExt} ]);
                        disp([ 'Moving ' [ allCell1{iCell} '.' fileExtensions{iExt} ] ' to ' [ allCell2{iCell} '.' fileExtensions{iExt} ] ]);
                    end;
                end;
            end;
            STUDY = STUDY2;
            STUDY = pop_savestudy(STUDY, ALLEEG, 'savemode', 'resave');
        end;
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [g,info] = gabwin(g,a,M,varargin);
%GABWIN  Compute a Gabor window from text or cell array
%   Usage: [g,info] = gabwin(g,a,M,L);
%
%   `[g,info]=gabwin(g,a,M,L)` computes a window that fits well with the
%   specified number of channels *M*, time shift *a* and transform length
%   *L*. The window itself is specified by a text description or a cell array
%   containing additional parameters.
%
%   The window can be specified directly as a vector of numerical
%   values. In this case, `gabwin` only checks assumptions about transform
%   sizes etc.
%
%   `[g,info]=gabwin(g,a,M)` does the same, but the window must be a FIR
%   window, as the transform length is unspecified.
%
%   `gabwin(g,a,M,L,lt)` or `gabwin(g,a,M,[],lt)` does as above but for a
%   non-separable lattice specified by *lt*. Please see the help of
%   |matrix2latticetype| for a precise description of the parameter *lt*.
%
%   The window can be specified as one of the following text strings:
%
%     'gauss'      Gaussian window fitted to the lattice,
%                  i.e. $tfr=a\cdot M/L$.
%
%     'dualgauss'  Canonical dual of Gaussian window.
%
%     'tight'      Tight window generated from a Gaussian.
%
%   In these cases, a long window is generated with a length of *L*.
%
%   It is also possible to specify one of the window names from |firwin|. In
%   such a case, `gabwin` will generate the specified FIR window with a length
%   of *M*.
%
%   The window can also be specified as cell array. The possibilities are:
%
%     `{'gauss',...}`
%         Additional parameters are passed to |pgauss|. When no additional
%         parameters are passed, the window is generated according to the
%         defaults in |pgauss|.
%
%     `{'dual',...}`
%         Canonical dual window of whatever follows. See the examples below.
%
%     `{'tight',...}`
%         Canonical tight window of whatever follows.
%
%   It is also possible to specify one of the window names from |firwin| as
%   the first field in the cell array. In this case, the remaining
%   entries of the cell array are passed directly to |firwin|.
%
%   Some examples: To compute a Gaussian window of length *L* fitted for a
%   system with time-shift *a* and *M* channels use::
%
%     g=gabwin('gauss',a,M,L);
%
%   To compute Gaussian window with equal time and frequency support
%   irrespective of *a* and *M*::
%
%     g=gabwin({'gauss'},a,M,L);
%
%   To compute the canonical dual of a Gaussian window fitted for a
%   system with time-shift *a* and *M* channels::
%
%     gd=gabwin('gaussdual',a,M,L);
%
%   To compute the canonical tight window of the Gaussian window fitted
%   for the system::
%
%     gd=gabwin({'tight','gauss'},a,M,L);
%
%   To compute the dual of a Hann window of length 20::
%
%     g=gabwin({'dual',{'hann',20}},a,M,L);
%
%   The structure *info* provides some information about the computed
%   window:
%
%     `info.gauss`
%        True if the window is a Gaussian.
%
%     `info.tfr`
%        Time/frequency support ratio of the window. Set whenever it makes sense.
%
%     `info.wasrow`
%        Input was a row window
%
%     `info.isfir`
%        Input is an FIR window
%
%     `info.isdual`
%        Output is the dual window of the auxiliary window.
%
%     `info.istight`
%        Output is known to be a tight window.
%
%     `info.auxinfo`
%        Info about auxiliary window.
%
%     `info.gl`
%        Length of window.
%
%   See also: pgauss, firwin, wilwin

% Assert correct input.
if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.keyvals.L=[];
definput.keyvals.lt=[0 1];
definput.keyvals.callfun='GABWIN';
definput.flags.phase={'freqinv','timeinv'};
[~,kv,L,lt]=ltfatarghelper({'L','lt'},definput,varargin,'gabwin');

[g,info] = comp_window(g,a,M,L,lt,kv.callfun);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [yn,en,S] = LMSadapt_dec(un,S)

% LMSadapt_dec      LMS Algorithm for Decision-Directed Channel Equalization
%
%                   The history of output, squared error and coefficients
%                   of FIR filter are passed to calling function.
%                   Only for decision-directed mode channel equalization
%
% Arguments:
% un                Input signal
% S                 Adptive filter parameters as defined in LMSinit.m
% yn                History of output signal
% en                History of error signal
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

M = length(S.coeffs);               % Length of FIR filter
mu = S.step;                        % Step size mu
leak = S.leakage;                   % Leaky factor
w = S.coeffs;                       % Weight vector of FIR filter
u = zeros(M,1);                     % Input signal vector

ITER = length(un);                  % Length of input sequence
yn = zeros(1,ITER);                 % Initialize output sequence to zero
en = zeros(1,ITER);                 % Initialize error sequence to zero

if isfield(S,'unknownsys')
    b = S.unknownsys;
    norm_b = norm(b);
    eml = zeros(1,ITER);
    ComputeEML = 1;
else
    ComputeEML = 0;
end

for n = 1:ITER
    u = [un(n); u(1:end-1)];         % Input signal vector [u(n),u(n-1),...,u(n-M+1)]'
    yn(n) = w'*u;                    % Output signal
    dn(n) = sign(yn(n));             % Only for decision directed mode channel equalizer
    en(n) = dn(n) - yn(n);           % Estimation error
    w = (1-mu*leak)*w + (mu*en(n))*u;% Tap-weight adaptation
    S.iter = S.iter + 1;
end

S.coeffs = w;                        % Coefficient values at the final iteration
if ComputeEML == 1;
    S.eml = eml;
end

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function varargout = getCobraSolverParams(solverType, paramNames, paramStructure)
% This function gets the specified paramStructure in `paramNames` from
% paramStructure, the global cobra paramters variable or default values set within
% this script.
%
% It will use values with the following priority
% paramStructure > solver type paramStructure > default paramStructure
%
% The specified paramStructure will be delt to the specified output arguements.
% See examples below.
%
% USAGE:
%
%    varargout = getCobraSolverParams(solverType, paramNames, paramStructure)
%
% INPUTS:
%    solverType:    Type of solver used: 'LP', 'MILP', 'QP', 'MIQP'
%    paramNames:    Cell array of strings containing parameter names OR one
%                   parameter name as string
%
% OPTIONAL INPUTS:
%    paramStructure:    Structure with fields pertaining to parameter values that
%                   should be used in place of global or default paramStructure.
%                   paramStructure can be set to 'default' to use the default
%                   values set within this script.
%
% OUTPUTS:
%    varargout:     Variables which each value corresponding to paramNames
%                   is output to.
%                   In addition, if paramStructure is passed as an input, then
%                   the final output argument is paramStructureOut, with fields
%                   removed corresponding to paramNames
%
% EXAMPLE:
%    paramStructure.saveInput = 'LPproblem.mat';
%    paramStructure.printLevel = 1;
%    [printLevel, saveInput] = getCobraSolverParams('LP', {'printLevel', 'saveInput'}, paramStructure);
%
%    %Example using default values
%    [printLevel, saveInput] = getCobraSolverParams('LP', {'printLevel','saveInput'}, 'default');
%
% .. Authors:
%       - Richard Que (12/01/2009)
%       - Ronan (16/07/2013) default MPS paramStructure are no longer global variables

if nargin < 2
    error('getCobraSolverParams: No paramStructure specified')
end

if nargin < 3
    paramStructure = [];
end
if ~isempty(paramStructure)
    paramStructureOut = paramStructure;
end
% Persistence will make names specific to one type of solver.
% Default Values
% For descriptions of the different settings please have a look at
% getCobraSolverParamsOptionsForType
valDef.minNorm = 0;

%These default tolerances are based on the default values for the Gurobi LP solver
%https://www.gurobi.com/documentation/9.0/refman/paramStructure.html
valDef.feasTol = 1e-6; % (primal) feasibility tolerance
valDef.optTol = 1e-6;  % (dual) optimality tolerance

valDef.objTol = 1e-6; % this should be used only when comparing the values of two objective functions

valDef.minNorm = [];

valDef.printLevel = 0;
valDef.verify = 0;

valDef.primalOnly = 0;

valDef.timeLimit = 1e36;
valDef.iterationLimit = 1000;

%valDef.logFile = ['Cobra' solverType 'Solver.log'];
valDef.logFile = []; %log file should be empty to avoid creating it by default
valDef.saveInput = [];
valDef.PbName = [solverType 'problem'];
valDef.debug = 0;
valDef.lifting = 0;

valDef.method = -1;

% CPLEX paramStructure
valDef.DATACHECK = 1;
valDef.DEPIND = 1;
valDef.checkNaN = 0;
valDef.warning = 0;

% tolerances
valDef.intTol = 1e-12;
valDef.relMipGapTol = 1e-12;
valDef.absMipGapTol = 1e-12;
valDef.NUMERICALEMPHASIS = 1;

if ~iscell(paramNames)
    paramNames = {paramNames};
end

switch solverType
    case 'LP'
        global CBT_LP_PARAMS
        parametersGlobal = CBT_LP_PARAMS;
    case 'QP'
        global CBT_QP_PARAMS
        parametersGlobal = CBT_QP_PARAMS;
    case 'MILP'
        global CBT_MILP_PARAMS
        parametersGlobal = CBT_MILP_PARAMS;
    case 'EP'
        global CBT_EP_PARAMS
        parametersGlobal = CBT_EP_PARAMS;
        valDef.feasTol = 1e-6; % (primal) feasibility tolerance
        valDef.optTol = 1e-6;  % (dual) optimality tolerance
    case 'MIQP'
        global CBT_MIQP_PARAMS
        parametersGlobal = CBT_MIQP_PARAMS;
    case 'NLP'
        global CBT_NLP_PARAMS
        parametersGlobal = CBT_NLP_PARAMS;
    otherwise
        display('Unrecognized solver type')
        return;
end

varargout = cell(1, numel(paramNames));
paramNames = columnVector(paramNames);
for i=1:length(paramNames)
    % set values to default
    if isfield(valDef,paramNames{i})
        varargout{i} = valDef.(paramNames{i});
    end
    if ~strcmp(paramStructure,'default') % skip of using default values
        % set values to global values
        if isfield(parametersGlobal,paramNames{i})
           varargout{i} = parametersGlobal.(paramNames{i});
        end
        % set values to specified values
        if isfield(paramStructure,paramNames{i})
           varargout{i} = paramStructure.(paramNames{i});
           %remove this field so paramStructureOut can be passed directly to the solver, without extra fields
           paramStructureOut = rmfield(paramStructureOut,paramNames{i});
        end
    end
end
if ~isempty(paramStructure)
    varargout{length(paramNames)+1}=paramStructureOut;
end

author: opencobra

repo: cobratoolbox


text:
function [ x, det ] = r8mat_solve_3d ( a, b, det, x )

%*****************************************************************************80
%
%% R8MAT_SOLVE_3D solves a 3 by 3 linear system using Cramer's rule.
%
%  Discussion:
%
%    If the determinant DET is returned as zero, then the matrix A is
%    singular, and does not have an inverse.  In that case, X is
%    returned as the [] vector.
%
%    If DET is nonzero, then its value is roughly an estimate
%    of how nonsingular the matrix A is.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(3,3), the matrix.
%
%    Input, real B(3), the right hand side.
%
%    Output, real DET, the determinant of the matrix A.
%
%    Output, real X(3), the solution of the system,
%    if DET is nonzero.
%

%
%  Compute the determinant.
%
  det =  a(1,1) * ( a(2,2) * a(3,3) - a(2,3) * a(3,2) ) ...
       + a(1,2) * ( a(2,3) * a(3,1) - a(2,1) * a(3,3) ) ...
       + a(1,3) * ( a(2,1) * a(3,2) - a(2,2) * a(3,1) );
%
%  If the determinant is zero, bail out.
%
  if ( det == 0.0 )
    x = [];
    return
  end
%
%  Compute the solution.
%
  x(1) = (   ( a(2,2) * a(3,3) - a(2,3) * a(3,2) ) * b(1) ...
           - ( a(1,2) * a(3,3) - a(1,3) * a(3,2) ) * b(2) ...
           + ( a(1,2) * a(2,3) - a(1,3) * a(2,2) ) * b(3) ) / det;

  x(2) = ( - ( a(2,1) * a(3,3) - a(2,3) * a(3,1) ) * b(1) ...
           + ( a(1,1) * a(3,3) - a(1,3) * a(3,1) ) * b(2) ...
           - ( a(1,1) * a(2,3) - a(1,3) * a(2,1) ) * b(3) ) / det;

  x(3) = (   ( a(2,1) * a(3,2) - a(2,2) * a(3,1) ) * b(1) ...
           - ( a(1,1) * a(3,2) - a(1,2) * a(3,1) ) * b(2) ...
           + ( a(1,1) * a(2,2) - a(1,2) * a(2,1) ) * b(3) ) / det;

  return
end


text:
 function [alpha, beta, ok] = nufft_best_alpha(J, L, K_N)
%function [alpha, beta, ok] = nufft_best_alpha(J, L, K_N)
%
%	return previously numerically optimized alpha and beta
%	use L=0 to return best available choice of any L
%
%	Copyright 2001-12-17	Jeff Fessler	The University of Michigan

if nargin < 1, help(mfilename), error(mfilename), end
if nargin < 2, L=2; end
if nargin < 3, K_N=2; end

Jlist1 = 6;		% list of which J's
alpha1 = [...		% last colum is best beta
	1 -0.46		0.19;	% J=6
];

Jlist2 = 2:10;		% list of which J's
alpha2 = [...		% last colum is best beta
	1 -0.200 -0.04	0.34;	% J=2
	1 -0.485 0.090	0.48;
	1 -0.470 0.085	0.56;	% J=4
	1 -0.4825 0.12	0.495;
	1 -0.57 0.14	0.43;	% J=6
	1 -0.465 0.07	0.65;
	1 -0.540 0.16	0.47;	% J=8
	1 -0.625 0.14	0.325;	% J=9
	1 -0.57 0.185	0.43;	% J=10 5.9707e-07
];

Jlist3 = [4 6];		% list of which J's
alpha3 = [...		% last colum is best beta
	1 -0.5319 0.1522 -0.0199	0.6339;	% J=4 2.5953e-04
	1 -0.6903 0.2138 -0.0191	0.2254;	% J=6 1.0097e-04
];

if K_N == 2

	if L==0
		if any(J == Jlist3)
			L = 3;
		else
			L = 2;
		end	% current best
	end

	if L==1
		alpha = alpha1;
		Jlist = Jlist1;

	elseif L==2
		alpha = alpha2;
		Jlist = Jlist2;

	elseif L==3
		alpha = alpha3;
		Jlist = Jlist3;

	else
		warning 'L not done'
		alpha = nan;
		beta = nan;
		ok = 0;
		return
	end

else
	warning 'K_N not done'
	alpha = 1;
	beta = 0.5;
	ok = logical(1);
	return
end

if any(J == Jlist)
	j = find(J == Jlist);
	beta = alpha(j,end);
	alpha = alpha(j,1:(end-1));
	ok = logical(1);
else
	ok = logical(0);
	alpha = nan;
	beta = nan;
end

author: thomaskuestner


text:
function [flag] = VBA_issymmetric (X)

% // VBA toolbox //////////////////////////////////////////////////////////
%
% [flag] = VBA_issymmetric (X)
% check if X is a symmetric matrix
%
% IN:
%   - X: matrix
% OUT:
%   - flag: - true if X is symmetric
%
% /////////////////////////////////////////////////////////////////////////

try
    flag = issymmetric (X);
catch
    flag = isequal (X, X');
end

author: MBB-team

repo: VBA-toolbox


text:
%DEMO_MEMORYSAVE  Demonstration of memory save option in GPstuff
%
%  Description
%    GP = GP_SET(..., 'savememory', 'on');
%
%    This demo consists of various combinations of covariance
%    functions, likelihoods and full/sparse approximations. This demo
%    is intended for showing how to use memory save option in
%    GPstuff. Unfortunately MATLAB doesn't have memory usage
%    monitoring, so the purpose of the demo is only to test that
%    results are the same whether you use memory save or not and to
%    show that the running times with memory saving are little bit
%    longer (more overhead). Memory saving can be enabled with the
%    following command
%
%    Usually GPstuff computes covariance matrix derivatives with
%    respect to all covariance function parameters at once which
%    avoids repeating computations common to each derivatibe, but
%    requires storage of Q NxN matrices, where Q is the number of
%    covariance function parameters. Q can be large, for example, if
%    there are many covariates each having own lengthsacle
%    parameter. With memory save option on, each matrix derivative is
%    computed separately, increasing the computation time little, but
%    reducing the memory usage. In an application, which motivated the
%    implementation of the option, the memory requirement of GPstuff
%    was reduced from over 60GB to less than 6GB.
%
%  See also
%    GP_SET
%
% Copyright (c) 2012 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% Classification data with full GP's

S = which('demo_classific');
L = strrep(S,'demo_classific.m','demodata/synth.tr');
x=load(L);
y=x(:,end);
y = 2.*y-1;
x(:,end)=[];

lik = lik_logit();

pl = prior_t();
pm = prior_sqrtunif();

gpcf = gpcf_neuralnetwork('weightSigma2', [0.9 0.9], 'biasSigma2', 10);
gpcf = gpcf_neuralnetwork(gpcf, 'weightSigma2_prior', pl,'biasSigma2_prior', pm); %

gpcf2 = gpcf_sexp('lengthScale', [0.9 0.9], 'magnSigma2', 10);
gpcf2 = gpcf_sexp(gpcf2, 'lengthScale_prior', pl,'magnSigma2_prior', pm); %

% GP without memory saving
gp = gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-4, 'latent_method', 'EP');
% GP with memory saving option enabled
gp2 = gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-4, 'latent_method', 'EP', 'savememory', 'on');

fprintf('Classification (EP), Neuralnetwork covariance function without and with memory saving (optimization and prediction)\n');
opt=optimset('TolX',1e-3,'TolFun',1e-3, 'Display', 'off');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'opt',opt);
[Eft,Varft, lpyt]=gp_pred(gp,x,y,x, 'yt', y);toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'opt',opt);
[Eft2,Varft2, lpyt2]=gp_pred(gp2,x,y,x, 'yt', y);toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft==Eft2)); assert(all(Varft==Varft2)); assert(all(lpyt==lpyt2));

gp = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-6, 'latent_method', 'Laplace');
gp2 = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-6, 'latent_method', 'Laplace', 'savememory', 'on');

fprintf('Classification (Laplace), Squared-Exponential covariance function without and with memory saving (optimization and prediction)\n');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'opt',opt);
[Eft,Varft, lpyt]=gp_pred(gp,x,y,x, 'yt', y);toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'opt',opt);
[Eft2,Varft2, lpyt2]=gp_pred(gp2,x,y,x, 'yt', y);toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft==Eft2)); assert(all(Varft==Varft2)); assert(all(lpyt==lpyt2));

% Regression data with sparse approximations

prevstream=setrandstream();
xx=linspace(1,10,901);
x1=logspace(0,1,100);
x1=round(x1*100)-99;
x=xx(x1)';
y=2*sin(4*x)+0.2*randn(size(x));
xt=[1:0.01:14]';
[n,nin] = size(x);

pl = prior_t('s2', 1);
pm = prior_logunif();
pn = prior_logunif();
gpcfse = gpcf_matern32('lengthScale',0.5,'magnSigma2',1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
lik = lik_gaussian('sigma2', 0.1, 'sigma2_prior', pn);
gp = gp_set('lik', lik, 'cf', gpcfse, 'jitterSigma2', 1e-6);
opt=optimset('TolFun',1e-4,'TolX',1e-4);

fprintf('Regression, FIC GP, Matern-3/2 covariance function, without and with memory saving (optimization and prediction)\n')
Xu=round(10+90*rand(18,1))/10; % Random placement

gp_fic = gp_set(gp, 'type','FIC','X_u',Xu,'infer_params','covariance+likelihood+inducing');
gp_fic2 = gp_set(gp, 'type','FIC','X_u',Xu,'infer_params','covariance+likelihood+inducing', 'savememory', 'on');

opt=optimset('TolFun',1e-4,'TolX',1e-4, 'Display', 'off');

% Optimization and prediction without memory saving
tic,gp_fic=gp_optim(gp_fic,x,y,'opt',opt);
[Eft_fic, Varft_fic] = gp_pred(gp_fic, x, y, xt);
Varft_fic = Varft_fic + gp_fic.lik.sigma2;toc

% Optimization and prediction with memory saving
tic,gp_fic2=gp_optim(gp_fic2,x,y,'opt',opt);
[Eft_fic2, Varft_fic2] = gp_pred(gp_fic2, x, y, xt);
Varft_fic2 = Varft_fic2 + gp_fic.lik.sigma2;toc
setrandstream([],prevstream);

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft_fic==Eft_fic2)); assert(all(Varft_fic==Varft_fic2));

gpcfse = gpcf_sexp('lengthScale',0.5,'magnSigma2',1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gp = gp_set('lik', lik, 'cf', gpcfse, 'jitterSigma2', 1e-6);

gp_fic = gp_set(gp, 'type','VAR','X_u',Xu,'infer_params','covariance+likelihood+inducing');
gp_fic2 = gp_set(gp, 'type','VAR','X_u',Xu,'infer_params','covariance+likelihood+inducing', 'savememory', 'on');

fprintf('Regression, VAR GP, Squared-Exponential covariance function, without and with memory saving (optimization and prediction)\n')

% Optimization and prediction without memory saving
tic,gp_fic=gp_optim(gp_fic,x,y,'opt',opt);
[Eft_fic, Varft_fic] = gp_pred(gp_fic, x, y, xt);
Varft_fic = Varft_fic + gp_fic.lik.sigma2;toc

% Optimization and prediction with memory saving
tic,gp_fic2=gp_optim(gp_fic2,x,y,'opt',opt);
[Eft_fic2, Varft_fic2] = gp_pred(gp_fic2, x, y, xt);
Varft_fic2 = Varft_fic2 + gp_fic.lik.sigma2;toc
setrandstream([],prevstream);

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft_fic==Eft_fic2)); assert(all(Varft_fic==Varft_fic2));

% Spatial data with sparse approximations

S = which('demo_spatial1');
data = load(strrep(S,'demo_spatial1.m','demodata/spatial1.txt'));
x = data(1:200,1:2);
ye = data(1:200,3);
y = data(1:200,4);
dims = [1    60     1    35];
[trindex, Xu] = set_PIC(x, dims, 5, 'corners', 0);
[n,nin] = size(x);

pl = prior_t('s2',10);
pm = prior_sqrtunif();

gpcf1 = gpcf_ppcs3('nin',nin,'lengthScale', 1, 'magnSigma2', 0.1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gpcf2 = gpcf_sexp('lengthScale', 5, 'magnSigma2', 0.05, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gpcf3 = gpcf_neuralnetwork('weightSigma2', [1 1], 'biasSigma2', 0.05, 'weightSigma2_prior', pl, 'biasSigma2_prior', pm);

lik = lik_negbin();
% GP without memory save
gp = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu, ...
            'jitterSigma2', 1e-4, 'infer_params', 'covariance+inducing');
% GP with memory saving option enabled
gp2 = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu, ...
            'jitterSigma2', 1e-4, 'infer_params', 'covariance+inducing', 'savememory', 'on');
opt=optimset('TolFun',1e-2,'TolX',1e-2, 'Display', 'off');

fprintf('Spatial process (Laplace), FIC GP, PPCS3 covariance function, without and with memory saving (optimization and prediction)\n');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'z',ye,'opt',opt);
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'z',ye,'opt',opt);
[Ef2, Varf2] = gp_pred(gp2, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Ef==Ef2)); assert(all(Varf==Varf2));

gp = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-4, 'latent_method', 'EP');
gp2 = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-4, 'latent_method', 'EP', 'savememory', 'on');

fprintf('Spatial process (EP), FIC GP, Squared-Exponential covariance function, without and with memory saving (optimization and prediction)\n');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'z',ye,'opt',opt);
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'z',ye,'opt',opt);
[Ef2, Varf2] = gp_pred(gp2, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Ef==Ef2)); assert(all(Varf==Varf2));

gp = gp_set('lik', lik, 'cf', gpcf3, 'jitterSigma2', 1e-4);
gp2 = gp_set('lik', lik, 'cf', gpcf3, 'jitterSigma2', 1e-4, 'savememory', 'on');

fprintf('Spatial process (Laplace), FIC GP, Neuralnetwork covariance function, without and with memory saving (optimization and prediction)\n');
tic,gp=gp_optim(gp,x,y,'z',ye,'opt',opt);
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye, 'tstind', [1:n]); toc

tic,gp2=gp_optim(gp2,x,y,'z',ye,'opt',opt);
[Ef2, Varf2] = gp_pred(gp2, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Ef==Ef2)); assert(all(Varf==Varf2));

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function ind_cell=run_movie(Y, A, C, Cn, min_max, Coor, ctr,  Ncell, tskip, save_avi, avi_name, S)
%% play movies together with selected neurons' calcium traces. It can be used for checking results
%% Input:
%
%% output

%% Author: Pengcheng Zhou, Carnegie Mellon University, 2016
%% parameters
[d1,d2] = size(Cn);
if ismatrix(Y); Y=reshape(Y, d1, d2, []); T=size(Y, 3); else T=size(Y, ndims(Y)); end
if ~exist('save_avi', 'var');    save_avi=false;
elseif ~exist('avi_name', 'var')
    k=0; avi_name = sprintf('example_movie_%d', k);
    while exist(avi_name, file);  k=k+1; avi_name = sprintf('example_movie_%d', k); end
end
if ~exist('tskip', 'var')||isempty(tskip);    tskip=1; end
if ~exist('Ncell', 'var')||isempty(Ncell);    Ncell=5; end
if ~exist('Coor', 'var')||isempty(Coor); figure; Coor=plot_contours(A, Cn, 0.9, 0); close; end
if ~exist('ctr', 'var')|| isempty(ctr); ctr=com(A, d1, d2); end
if ~exist('min_max', 'var') || (isempty(min_max));
    temp = Y(:, :, randi(T, min(100, T), 1));
    min_max = quantile(temp(:), [0.2, 0.9999]);
    min_max(1) = max(min_max(1), 0);
end
ncol = 3;
%% select neurons by clicking the image.
figure('position', [100, 500, 1800, 420]);

% display correlation image for selecting neurons
subplot(Ncell, ncol, 1:ncol:(ncol*Ncell));
imagesc(Cn, [0, 1]);hold on;
axis equal; axis off;
ind_cell = zeros(Ncell,1);
m = 1;
while true
    subplot(Ncell, ncol, 1:ncol:(Ncell*ncol));
    drawnow;
    [x, y] = ginput(1);
    h_dot = plot(x, y, '*m');
    dx = x-ctr(:, 2);
    dy = y-ctr(:, 1);
    dist_xy = sqrt(dx.^2+dy.^2);
    [~, ind_cell(m)] = min(dist_xy);
    temp = Coor{ind_cell(m)};
    h_con = plot(temp(1, 4:end), temp(2, 4:end), 'r');
    h_txt = text(ctr(ind_cell(m), 2), ctr(ind_cell(m), 1), num2str(m), 'fontsize', 10, 'fontweight', 'bold');
    disp(ctr(ind_cell(m), :));
    subplot(Ncell, ncol, ncol*(m-1)+[2, ncol]);cla;
    ylabel(num2str(m));
    plot(C(ind_cell(m), :)); hold on;
    ylabel(sprintf('ID %d', ind_cell(m)));
    if (exist('S', 'var'))&&  (~isempty(S))
        tmp = find(S(ind_cell(m), :)>0);
        plot(tmp, S(ind_cell(m), tmp), '*g');
    end
    set(gca, 'xticklabel', []); axis tight;
    set(gca, 'yticklabel', []);
    xlim([1, T]);

    temp = input('keep the neuron? (y/n)   ', 's');
    if strcmp(temp, 'n')
        delete(h_dot); delete(h_con); delete(h_txt);
        continue;
    else
        ylabel(num2str(m));
        m =  m+1;
    end

    if m>Ncell; break; end
end
subplot(Ncell, ncol, ncol*Ncell+[-1,0]);
set(gca, 'xticklabel', get(gca, 'xtick'));
xlabel('Frame');

%% play movies
% create avi file if sve_avi==true
if save_avi
    avi_file = VideoWriter(avi_name);
    open(avi_file);
end
% draw selected neurons' contours and latex them
subplot(Ncell,ncol,1:ncol:(ncol*Ncell)); cla;  hold on; axis equal; axis off;
for m=1:Ncell
    temp = Coor{ind_cell(m)};
    plot(temp(1, 4:end), temp(2, 4:end), 'r');
    text(ctr(ind_cell(m), 2)+3*randn(1), ctr(ind_cell(m), 1)+3, num2str(m), 'fontsize', ...
        15, 'fontweight', 'bold', 'color', 'g');
end
% play movie
for t=1:tskip:T
    subplot(Ncell,ncol,1:ncol:(ncol*Ncell));
    h_img = imagesc(Y(:, :, t), min_max);
%     colormap gray;
    hold on; axis equal; axis off;
    set(gca, 'children', flipud(get(gca, 'children')));

    % shift the red line indicating the current frame
    for m=1:Ncell
        subplot(Ncell, ncol, ncol*(m-1)+[2, ncol]);
        eval(sprintf('aa%d=plot([t,t], get(gca, ''ylim''), ''r''); ',m));
    end

    % save frame
    if save_avi
        frame = getframe(gcf);
        writeVideo(avi_file, frame);
    end
    drawnow();

    % remove red lines
    for m=1:Ncell
        eval(sprintf('delete(aa%d); ', m));
    end
    delete(h_img);
end
subplot(Ncell,ncol,1:ncol:(ncol*Ncell)); hold on; axis equal; axis off;
imagesc(Cn, [0, 1]);
set(gca, 'children', flipud(get(gca, 'children')));

if save_avi;  close(avi_file);end

author: zhoupc


text:
function cleanVolume()
%
% function cleanVolume()
%
% Deletes:
%   Volume/coords.mat
%   Volume/dataType/corAnal.mat (for each dataType)
%   Volume/dataType/*.mat (all other parameter maps)
%
% If you change this function make parallel changes in:
%     cleanFlat, cleanDataType
%
% djh, 2/2001

global HOMEDIR

curVolumeDir = fullfile(HOMEDIR,'Volume');
backupVolumeDir = fullfile(HOMEDIR,sprintf('deletedVolume_%s', ...
    datestr(now, 'yyyy-dd-mm-hh_MM-ss')));

if exist(curVolumeDir, 'dir')
    movefile(curVolumeDir, backupVolumeDir);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [CVA] = spm_cva_prob (X1,X2,m)
% Probabilistic Canonical Variates Analysis
% FORMAT [CVA] = spm_cva_prob (X1,X2,m)
%
% X1           [d1 x N] matrix of dependent variables
% X2           [d2 x N] matrix of independent variables
% m            dimension of latent variable (min([d1,d2]) by default)
%
% Returns fields:
%
% .U1,.U2      Canonical vectors
% .W1,.W2      Factor matrices
% .L           Log-Likelihood
% .bic         Bayesian Information Criterion
% .aic         Akaike's Information Criterion
%
% Fits probabilistic model
%
% x1 = W1 z + e1
% x2 = W2 z + e2
%
% This algorithm is described in:
%
% F. Bach and M. Jordan (2005) A probabilistic interpretation of canonical
% correlation analysis. Dept. Stats, Univ California, Berkeley CA.
% Tech Rep 688.
%
%___________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_cva_prob.m 4687 2012-03-14 18:15:49Z will $

[d1,N1]=size(X1);
[d2,N2]=size(X2);
d=d1+d2;

if ~N1==N2
    disp('Error in spm_cva_prob: unequal number of samples');
    return
else
    N=N1;
end

if nargin < 3 | isempty(m)
    m=min([d1,d2]);
else
    if m > min([d1,d2]);
        disp('m too large');
        return
    end
end

X=[X1;X2];
Sigma=cov(X')+(10^-8*eye(d1+d2));
Sigma11=Sigma(1:d1,1:d1);

if m==0
    S=diag(diag(Sigma));
    iS=inv(S);
    CVA.L=-0.5*N*(d1+d2)*log(2*pi)-0.5*N*spm_logdet(S)-0.5*N*trace(iS*Sigma);
    CVA.bic=CVA.L;
    CVA.aic=CVA.L;
    CVA.W1=[];CVA.W2=[];CVA.U1=[];CVA.U2=[];
    return
end

Sigma12=Sigma(1:d1,d1+1:d1+d2);
Sigma22=Sigma(d1+1:d,d1+1:d);
Sigma21=Sigma12';

R1=inv(sqrtm(Sigma11));
R2=inv(sqrtm(Sigma22));
A=R1*Sigma12*R2;

[V1,P,V2]=svd(A,0);
sP=diag(P);
rP=sP/max(sP);
U1=R1*V1(:,1:m);
U2=R2*max(sP)*V2(:,1:m);

M1=diag(sqrt(rP(1:m)));
M2=M1;

W1=Sigma11*U1*M1;
W2=Sigma22*U2*M2;
Psi1=diag(diag(Sigma11-W1*W1'));
Psi2=diag(diag(Sigma22-W2*W2'));

W=[W1;W2];
S=W*W'+blkdiag(Psi1,Psi2);

CVA.W1=W1;
CVA.W2=W2;
CVA.U1=U1;
CVA.U2=U2;

iS=inv(S);
CVA.L=-0.5*N*(d1+d2)*log(2*pi)-0.5*N*spm_logdet(S)-0.5*N*trace(iS*Sigma);
k=2*m*(d1+d2); % W1, W2 plus diag Psi's

CVA.bic=CVA.L-0.5*k*log(N);
CVA.aic=CVA.L-k;

author: spm


text:
% Copyright (C) 2014-2017 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/>.

addpath(pwd)
base = 'tests_matlab';
files = dir(base);
cd(base);

% usually I'd just use cputime(), but some of our IPC mechanisms are
% light on CPU and heavy on IO.
totaltime = clock();
totalcputime = cputime();
num_tests = 0;
% do tests in random order:
%rng('shuffle')
%for i=randperm(length(files))
for i=1:length(files)
  mfile = files(i).name;
  % detect tests b/c directory contains other stuff (e.g., surdirs and
  % helper files)
  if ( (~files(i).isdir) && strncmp(mfile, 'test', 4) && mfile(end) ~= '~')
    testtime = clock();
    str = mfile(1:end-2);
    num_tests = num_tests + 1;
    fprintf(['>>> Running test(s) in: ' mfile '  ']);  % no newline
    eval(str)
    testtime = etime(clock(), testtime);
  end
end

totaltime = etime(clock(),totaltime);
totalcputime = cputime() - totalcputime;
fprintf('\n***** Ran tests from %d files, %g seconds (%gs CPU) *****\n', ...
        num_tests, totaltime, totalcputime);
cd('..')

repo: octsympy


text:
function [OV,OF] = combine(IV,IF,EV,EF)
  % COMBINE Two meshes where the second mesh was constructed as the surface of
  % the union of the volume of the first mesh and a third mesh. The order of
  % vertices is such that the first chunk of vertices in the second mesh is a
  % subsequence of vertices in the first mesh (the remaining chunks are a
  % subsequence of vertices from the third mesh, then new vertices). The goal
  % is to create a new mesh where the vertices are: all the vertices from the
  % first mesh then other vertices and faces from the second mesh. This means
  % in the resulting mesh the unreferenced vertices correspond to vertices that
  % were in the first mesh but not in the second mesh.
  %
  % Inputs:
  %   IV  vertex list of the first mesh
  %   IF  face list of the first mesh
  %   EV  vertex list of the second mesh
  %   EF  face list of the second mesh
  % Outputs:
  %   OV  vertex list of the output mesh
  %   OF  face list of the output mesh
  %

  % create new list of vertices which is IV followed by vertices in EV not in
  % IV
  OV = [IV; setdiff(EV,IV,'rows')];
  % find location in new list for each vertex in EV
  [TF,LOC] = ismember(EV,OV,'rows');
  assert(min(LOC) > 0);
  OF = LOC(EF);
end

author: alecjacobson

repo: gptoolbox


text:
function tests = test_imEquivalentEllipsoid
% Test suite for the file imEquivalentEllipsoid.
%
%   Test suite for the file imEquivalentEllipsoid
%
%   Example
%   test_imEquivalentEllipsoid
%
%   See also
%     imEquivalentEllipsoid

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2021-11-03,    using Matlab 9.10.0.1684407 (R2021a) Update 3
% Copyright 2021 INRAE - BIA-BIBS.

tests = functiontests(localfunctions);

function test_Simple(testCase) %#ok<*DEFNU>
% Test call of function without argument.

%
img = readstack(fullfile('files', 'ellipsoid_Center30x27x25_Size20x12x8_Orient40x30x20.tif'));

elli = imEquivalentEllipsoid(img > 0);

assertEqual(testCase, size(elli), [1 9]);
assertEqual(testCase, [30 27 25], elli(1:3), 'AbsTol', 0.5);
assertEqual(testCase, [20 12  8], elli(4:6), 'AbsTol', 0.5);
assertEqual(testCase, [40 30 20], elli(7:9), 'AbsTol', 0.5);

function test_Calibrated(testCase) %#ok<*DEFNU>
% Test call of function without argument.

img = readstack(fullfile('files', 'ellipsoid_Center30x27x25_Size20x12x8_Orient40x30x20.tif'));

elli = imEquivalentEllipsoid(img > 0, 'Spacing', [0.5 0.5 0.5], 'Origin', [0.5 0.5 0.5]);

assertEqual(testCase, size(elli), [1 9]);
assertEqual(testCase, [15 13.5 12.5], elli(1:3), 'AbsTol', 0.5);
assertEqual(testCase, [10 6 4], elli(4:6), 'AbsTol', 0.5);
assertEqual(testCase, [40 30 20], elli(7:9), 'AbsTol', 0.5);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function y = maxstar(x, w)
% maxstar   Log of a sum of exponentials.
%   For vectors, maxstar(x) is equivalent to log(sum(exp(x))).
%   For matrices, maxstar(x) is a row vector and maxstar operates on
%   each column of x. For N-D arrays, maxstar(x) operates along the
%   first non-singleton dimension.
%
%   maxstar(x,w) is the log of a weighted sum of exponentials,
%   equivalent to log(sum(w.*exp(x))). Vectors w and x must be
%   the same length. For matrix x, the weights w can be input as
%   a matrix the same size as x, or as a vector of the same length
%   as columns of x. Weights may be zero or negative, but the result
%   sum(w.*exp(x)) must be greater than zero.
%
%   Acts on first dimensions

if nargin<2
    w = [];
else
    if ~isvector(w)
        error('maxstar: w must be a vector')
    end
    if length(w) ~= size(x,1)
        error('maxstar: weight does not match x')
    end
end
%%
w = w(:);
szx = size(x);
if isempty(w)
    % no weight
    m = max(x);
    y = m + log(sum(exp(bsxfun(@minux,x,m))));
else
    % Move the weight into the exponent xw and find
    % m = max(xw) over terms with positive weights
    wpos = w>0;
    xw = bsxfun(@plus, x(wpos,:), log(w(wpos)));
    m = max(xw);
    exwp = exp( bsxfun(@minus, xw, m) );
    wneg = w<0;
    exwn = exp( x(wneg,:) + bsxfun(@minus,log(-w(wneg)), m) );
    y = m + log(sum(exwp,1) - sum(exwn,1));
end
y = reshape(y,[szx(2:end) 1]);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%BASE64ENCODE encodes a vector of uint8s as a base64 string
% BASE64ENCODE(BYTES)
%    Encodes every tree BYTES in four printable ASCII characters.
%    If LENGTH(BYTES) is not a multiple of three, it is padded with
%    zeros and unused characters are replaced with '='.

% (c) 2014 Bastian Bechtold
% This code is licensed under the BSD 3-clause license

function base64 = base64encode(bytes)
    % pad the base64 string to a multiple of 3
    if mod(length(bytes), 3) ~= 0
        padding = 3-mod(length(bytes), 3);
        bytes = [bytes; zeros(3-mod(length(bytes), 3), 1, 'uint8')];
    else
        padding = 0;
    end

    % convert every three uint8 bytes into four base64 bytes
    base64 = zeros(length(bytes)/3*4, 1, 'uint8');
    base64(1:4:end) = bitshift(bytes(1:3:end), -2);
    base64(2:4:end) = bitor(bitshift(bitand(bytes(1:3:end),  3),  4), ... % two LSB
                      bitshift(bytes(2:3:end), -4));
    base64(3:4:end) = bitor(bitshift(bitand(bytes(2:3:end),  15),  2), ... % four LSB
                            bitshift(bytes(3:3:end), -6));
    base64(4:4:end) = bitand(bytes(3:3:end), 63); % six LSB

    % convert from base64 bytes to string representation
    table = ['A':'Z' 'a':'z' '0':'9' '+' '/'];
    base64 = table(base64+1);
    base64(end-padding+1:end) = '=';
end

author: bastibe

repo: transplant

save_path: github-repos/MATLAB/bastibe-transplant


text:
function [Nv, VX, K, EToV] = MeshGenDistMesh1D()

% function [VX, K, EToV] = MeshGenDistMesh1D()
% Purpose  : Generate 1D mesh using DistMesh;

% distance function for a circle about xc=0
fd=inline('abs(p)-1','p');

% distribution weight function
fh=inline('abs(p)*0.075+0.0125','p');

% generate non-uniform example mesh in 1D using DistMesh
h0 = 0.025;            % chosen element spacing
[p,t]=distmeshnd(fd,fh,h0,[-1;1],[]);
K = size(t,1); Nv = K+1;

% Sort elements in acending order
[x,i]=sort(p(t));
t=t(i,:);
Met=sort(t','ascend')';
EToV = t; VX = p(:,1);
return

repo: Aero-matlab


text:
clear;
addpath(genpath('../'));

output_dir = './CECLM_res_menpo/';
if(~exist(output_dir, 'dir'))
    mkdir(output_dir)
end
%% select database and load bb initializations
db_root = 'E:\datasets\300VW\300VW_Dataset_2015_12_14/';
bb_root = './300VW_dets_mtcnn/';
extra_dir = 'E:\datasets\300VW\300VW_Dataset_2015_12_14\extra';
[ vid_locs, bboxes, gts_all, invalid_frames ] = CollectTestData(db_root, bb_root, extra_dir);

%% loading the patch experts and the PDM
[patches, pdm, clmParams, early_term_params] = Load_CECLM_menpo();
multi_view = true;

%% Select video
for i=1:numel(vid_locs)

    vid = VideoReader(vid_locs{i});

    bounding_boxes = bboxes{i};

    n_frames = size(bounding_boxes,1);
    preds = zeros(68,2,n_frames);
    for f=1:n_frames
        input_image = readFrame(vid);

        %% Initialize from detected bounding boxes every 30 frames
        if (mod(f-1, 30) == 0)
            ind = min(f, size(bounding_boxes,1));
            bb = bounding_boxes(ind, :);
            % If no face detected use the closest detected BB
            if(bb(3) == 0)
               ind_next = ind + find(bounding_boxes(ind+1:end,3), 1);
               if(isempty(bounding_boxes(ind_next)) || bounding_boxes(ind_next,3)==0)
                   ind_next = find(bounding_boxes(1:ind,3), 1, 'last');
               end
               bb = bounding_boxes(ind_next, :);
            end
            bb(3) = bb(1) + bb(3);
            bb(4) = bb(2) + bb(4);
            reset = true;
        else
            reset = false;
        end

        % have a multi-view version for initialization, otherwise use
        % previous shape
        if(reset && multi_view)
            clmParams.window_size = [25,25; 23,23; 21,21; 21,21];
            clmParams.numPatchIters = 4;
            clmParams.startScale = 1;

            views = [0,0,0; 0,-30,0; 0,30,0; 0,-55,0; 0,55,0; 0,0,30; 0,0,-30; 0,-90,0; 0,90,0; 0,-70,40; 0,70,-40];
            views = views * pi/180;

            [shape,g_param,l_param,lhood,lmark_lhood,view_used] =...
                Fitting_from_bb_multi_hyp(input_image, [], bb, pdm, patches, clmParams, views, early_term_params);

        else
            clmParams.window_size = [23,23; 21,21; 19,19; 17,17];
            clmParams.numPatchIters = 3;
            clmParams.startScale = 2;
            [shape,g_param,l_param,lhood,lmark_lhood,view_used] = Fitting_from_bb(input_image, [], bb, pdm, patches, clmParams, 'gparam', g_param, 'lparam', l_param);
        end

        preds(:,:,f) = shape;

        %% plot the result
        imshow(input_image);
        hold on;
        plot(shape(:,1), shape(:,2), '.r');
%         rectangle('Position', [bb(2), bb(1), bb(4), bb(3)]);
        hold off;
        drawnow expose

    end

    %% Grab the ground truth
    gt_landmarks = gts_all{i};

    % Remove unreliable frames
    if(~isempty(invalid_frames{i}))
        preds(:,:,int32(invalid_frames{i}))=[];
    end

    if(size(gt_landmarks,3) ~= size(preds,3))
        fprintf('something went wrong with vid %d\n', i);
    end

    [vid_name,~,~] = fileparts(vid_locs{i});
    [~,vid_name,~] = fileparts(vid_name);
    save([output_dir, '/', vid_name], 'preds', 'gt_landmarks');
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function pv=computePositionVariance(view)
% computePostionVariance - compute position variance across
% cortical surface
%
% pv=computePositionVariance(view)
%
% 2008/05 SD & KA: wrote it

if ~exist('view','var') || isempty(view), error('Need view struct'); end

model = viewGet(view,'rmCurModel');

x = rmGet(model,'x0');
y = rmGet(model,'y0');
w = rmGet(model,'varexp');

edges        = double(view.edges);
numNeighbors = double(view.nodes(4,:));
edgeOffsets  = double(view.nodes(5,:));

pv = varOfNeighbors(x,y,w,edges,edgeOffsets,numNeighbors);

return


text:
function xyz = sphere_cubed_ijk_to_xyz ( n, i, j, k )

%*****************************************************************************80
%
%% SPHERE_CUBED_IJK_TO_XYZ: cubed sphere IJK to XYZ coordinates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sections into which each face of
%    the cube is to be divided.
%
%    Input, integer I, J, K, indices between 0 and N.  Normally,
%    at least one of the indices should have the value 0 or N.
%
%    Output, real XYZ(3,1), coordinates of the point.
%
  if ( i == 0 )
    xc = -1.0;
  elseif ( i == n )
    xc = +1.0;
  else
    xc = tan ( ( 2 * i - n ) * 0.25 * pi / n );
  end

  if ( j == 0 )
    yc = -1.0;
  elseif ( j == n )
    yc = +1.0;
  else
    yc = tan ( ( 2 * j - n ) * 0.25 * pi / n );
  end

  if ( k == 0 )
    zc = -1.0;
  elseif ( k == n )
    zc = +1.0;
  else
    zc = tan ( ( 2 * k - n ) * 0.25 * pi / n );
  end

  xyzn = sqrt ( xc^2 + yc^2 + zc^2 );

  xyz = zeros ( 3, 1 );

  xyz(1) = xc / xyzn;
  xyz(2) = yc / xyzn;
  xyz(3) = zc / xyzn;

  return
end

author: johannesgerer


text:
function geometry_test20701 ( )

%*****************************************************************************80
%
%% TEST20701 tests TRIANGLE_INRADIUS_2D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;

  t = [ ...
    0.0, 1.0; ...
    0.0, 0.0; ...
    1.0, 0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST20701\n' );
  fprintf ( 1, '  For a triangle in 2D,\n' );
  fprintf ( 1, '  TRIANGLE_INRADIUS_2D computes the inradius.\n' );

  r8mat_transpose_print ( dim_num, 3, t, '  Triangle vertices:' );

  r = triangle_inradius_2d ( t );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Incircle radius is %f\n', r );

  return
end


text:
function plot2DPlane(x1,x2,y1,y2,h_newfig)

if h_newfig
    figure;
    set(gcf,'position',[500,100,1000,650])
    set(gcf,'color','w');
end
set(gca,'visible','off')
mArrow2(x1,0,x2,0,{'color','k'});
mArrow2(0,y1,0,y2,{'color','k'});
xlim([x1 x2])
ylim([y1 y2])

for i = (x1+1):(x2-1)
    line([i i],[-0.05 0.05],'color','k')
    if i<=0
        t= text(i-0.2, -0.3,num2str(i));
    else
        t= text(i-0.05, -0.3,num2str(i));
    end
    t.FontSize=15;

end

for i = (y1+1):(y2-1)
    line([-0.05 0.05],[i i],'color','k')

    if i<0
        t = text(-0.4, i,num2str(i));
    elseif i>0
        t = text(-0.3, i,num2str(i));
    end
    t.FontSize = 15;
end
t = text(x2*0.9,0.3,'x');
t.FontSize = 15;
t = text(0.2,y2*0.9,'y');
t.FontSize = 15;
axis square

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols