text:
%%This Matlab script generates Figure 5 in the paper:
%
%Emil Björnson, Özgecan Özdogan, Erik G. Larsson, “Intelligent Reflecting
%Surface vs. Decode-and-Forward: How Large Surfaces Are Needed to Beat
%Relaying?,” IEEE Wireless Communications Letters, To appear
%
%Download article: https://arxiv.org/pdf/1906.03949
%
%This is version 1.0 (Last edited: 2019-10-28)
%
%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
%paper as described above.

close all;
clear;

%% Set simulation parameters

%Carrier frequency (in GHz)
fc = 3;

%Bandwidth
B = 10e6;

%Noise figure (in dB)
noiseFiguredB = 10;

%Compute the noise power in dBm
sigma2dBm = -174 + 10*log10(B) + noiseFiguredB;
sigma2 = db2pow(sigma2dBm);

%Define the channel gain functions based on the 3GPP Urban Micro in
%"Further advancements for E-UTRA physical layer aspects (Release 9)."
%3GPP TS 36.814, Mar. 2010. Note that x is measured in m and that the
%antenna gains are included later in the code
pathloss_3GPP_LOS = @(x) db2pow(-28-20*log10(fc)-22*log10(x));
pathloss_3GPP_NLOS = @(x) db2pow(-22.7-26*log10(fc)-36.7*log10(x));

%Define the antenna gains at the source, relay/IRS, and destination. The
%numbers are in linear scale
antennaGainS = db2pow(5);
antennaGainR = db2pow(5);
antennaGainD = db2pow(0);

%Set the amplitude reflection coefficient
alpha = 1;

%Set the range of rate values
Rbar = [0.01 0.1:0.1:10];

%Set parameters related to circuit power consumption
Ps = 100; %Power dissipation in the transceiver hardware of the source
Pd = 100; %Power dissipation in the transceiver hardware of the destination
Pe = 5;   %Power dissipation per element in the IRS (mW)
Pr = 100; %Power dissipation in the transceiver hardware of the relay
nu = 0.5; %Efficiency of the power amplifier at the source

%Define distances in simulation setup

d_SR = 80; %Distance between the source and IRS/relay
dv = 10; %Minimum distance between destination and the IRS/relay

%Define the range of d1 values in the simulation setup
d1 = 70;

%Copmute distance between the source and destination
d_SD = sqrt(d1^2+dv^2);

%Compute distance between the IRS/relay and destination
d_RD = sqrt((d1-d_SR)^2+dv^2);

%Compute the channel gains using the 3GPP models and antenna gains
betaSR = pathloss_3GPP_LOS(d_SR)*antennaGainS*antennaGainR;
betaRD = pathloss_3GPP_LOS(d_RD)*antennaGainR*antennaGainD;
betaSD = pathloss_3GPP_NLOS(d_SD)*antennaGainS*antennaGainD;

%Prepare to save simulation results
EE_SISO = zeros(length(Rbar),1);
EE_IRS = zeros(length(Rbar),1);
EE_DF = zeros(length(Rbar),1);
Nopt = zeros(length(Rbar),1);

%% Go through all rate values
for ind = 1:length(Rbar)

    %Compute required SINR values
    SINR = 2^(Rbar(ind))-1; %SISO and IRS
    SINR_DF = 2^(2*Rbar(ind))-1; %DF relaying

    %Compute the transmit power in the SISO case, using Eq. (17)
    P_SISO = SINR*sigma2/betaSD;

    %Compute the energy efficiency in the SISO case
    %(the factor 1000 is used to convert mW to W)
    EE_SISO(ind) = 1000*B*Rbar(ind)/(P_SISO/nu + Ps + Pd);

    %Compute the transmit power in the DF relaying case, using Eq. (19)
    P_DF = SINR_DF*sigma2*(betaSR+betaRD-betaSD)/(2*betaRD*betaSR);

    %Compute the energy efficiency in the DF relaying case
    %(the factor 1000 is used to convert mW to W)
    EE_DF(ind) = 1000*B*Rbar(ind)/(P_DF/nu + Ps/2 + Pd + Pr);

    %Compute the power-minimizing number of reflecting elements
    Nopt(ind) = (2*SINR*sigma2/(alpha^2*betaSR*betaRD*Pe))^(1/3) - sqrt(betaSD/(betaSR*betaRD))/alpha;

    if Nopt(ind)<0
        Nopt(ind) = 0;
    end

    %Compute the transmit power in the IRS case, using Eq. (18)
    P_IRS = SINR*sigma2./(sqrt(betaSD) + Nopt(ind)*alpha*sqrt(betaSR*betaRD)).^2;

    %Compute the energy efficiency in the IRS case
    %(the factor 1000 is used to convert mW to W)
    EE_IRS(ind) = 1000*B*Rbar(ind)/(P_IRS/nu + Ps + Pd + Nopt(ind)*Pe);

end

%Plot simulation results
figure;
hold on; box on;
plot(Rbar,EE_DF/1e6,'b-.','LineWidth',2);
plot(Rbar,EE_IRS/1e6,'r-','LineWidth',2);
plot(Rbar,EE_SISO/1e6,'k--','LineWidth',2);
xlabel('Achievable rate [bit/s/Hz]','Interpreter','Latex');
ylabel('Energy efficiency [Mbit/Joule]','Interpreter','Latex');
legend('DF relay','IRS','SISO','Location','NorthWest');
set(gca,'fontsize',18);

author: emilbjornson

repo: IRS-relaying


text:
function [Y, bias, scale] = scaleData(Y, scaleMethod, scaleVal, bias)
% SCALEDATA Scale and center data
% VARGPLVM

if nargin < 4, bias = []; end
if nargin < 3, scaleVal = [];   end
if nargin < 2, scaleMethod = []; end

d = size(Y,2);
if isempty(bias)
    bias = mean(Y);
end

% Remove bias
m = Y;
for i = 1:d
  m(:, i) = m(:, i) - bias(i);
end

scale = ones(1, d);

if ~isempty(scaleMethod) && scaleMethod ~= 0
    if ~isempty(scaleVal)
        warning('Both scale2var1 and scaleVal set for GP');
    end
    if(scaleMethod == 1) % Scale to variance 1
        scale = std(Y);
    elseif scaleMethod == 2 % Scale so that maximum is 1
        scale = max(max(abs(m)));
    else
        error('Unknown scale option')
    end
    scale(find(scale==0)) = 1;
end

if isscalar(scaleVal)
    if scaleVal
        scale = repmat(scaleVal, 1, d);
    end
elseif ~isempty(scaleVal) && ~sum(scaleVal==0)
    scale = scaleVal;
end

% Apply scale
for i = 1:d
  if scale(i)
    m(:, i) = m(:, i)/scale(i);
  end
end

Y = m;

author: SheffieldML

repo: vargplvm


text:
function [y,O,X,S] = filterAdjust(O)
% :Usage:
% ::
%
%     [y,O,X,S] = filterAdjust(OPTIONS)
%
% :O:
%
%   **O.y:**
%        signal
%
%   **O.HP:**
%        high pass freq. cutoff
%
%   **O.TR:**
%        sampling rate in s
%
%   **O.doHP:**
%        [0 or 1] - do HP filter (spm), default is 0
%
%   **O.doLP:**
%        [0, 1, 2] - do LP filter (spm), default is 0
%           2 = Gaussian filter with TR*2 s length
%
%   **O.firstimg:**
%        sets values for first image in each run to mean of the
%        remaining values.  Good for removing first-image artifacts present in
%        some scanner sequences.  Default is 1.
%
%   **O.cyclecorrection:**
%        checks for unimodal (normally distributed) data
%        within each session, because some bimodal data that cycles between two
%        mean scanner values has been observed in some data.  if a high proportion
%        of outliers are found in non-normal data, subtracts mean of higher mode
%        to adjust data. Sorry-not clear.  Check the code. Default is 0
%
%   **O.cyclecorrection2:**
%        removes large transitions from data, as in
%        detransition.m.  Default is 0.  Artifacts may be acquisition or
%        motion-correction/resampling related.
%
%        We do this after filtering, but if there's trouble, we re-do the
%        scanadjust and filtering, because 'cycling' can affect these.
%
%   **O.scanadjust:**
%        [0 or 1] - adjust to scan means, default is 0
%           - If O.X is entered, assumes this is the session mean matrix,
%           instead of recomputing.
%
%   **O.percent:**
%        [0 or 1] - adjust to percent change from baseline, default is 0
%
%   **O.filtertype:**
%        filter style, default is 'none'
%         - 'spm', use spm's filtering
%         - if O.S is entered, uses this instead of recomputing
%         - 'fourier', Doug's fourier filter
%         - 'fouriernotch', Omit frequencies between HP(1) and HP(2)
%         - 'cheby', chebyshev
%         - 'Luis', Luis' custom filter
%         - 'none', no filtering (or leave field out).
%
%   **O.HP:**
%        for SPM, the filter cutoff in s
%        for fourier, the HP value or the [HP LP] values
%        notches out everything slower than HP and faster than LP, in s
%        (1/s = Hz).
%
%   **O.nruns:**
%        number of runs (scanadjust), default is 1
%
%   **O.adjustmatrix:**
%        custom adjustment matrix to regress out (e.g., movement params)
%
%   **O.plot [0 or 1]:**
%        plots intermediate products, default is 0
%
%   **O.verbose [0 or 1]:**
%        verbose output
%
%   **O.trimts [0 or std]:**
%        trim overall timseries values to std, 0 for no trimming
%
%   **O.lindetrend:**
%        specify linear detrending of timeseries.
%        occurs after adjustment and filtering and windsorizing
%            - detrending option -> what to enter in this field:
%            - no detrending  -> empty, missing field, or 0
%            - detrending every n elements -> single number (n)
%            - piecewise linear detrend -> ROW vector of breakpoints
%              (do not specify 1 as the start of the 1st segment.)
%
% ..
%    by Tor Wager, 09/24/01
%    modified 10/15/02 to add trial baseline adjustment and linear detrending.
% ..

if ~isfield(O,'verbose'), O.verbose = 0;,end
if ~isfield(O,'trimts'), O.trimts = 0;,end
if ~isfield(O,'lindetrend'), O.lindetrend = 0;, end
if ~isfield(O,'percent'), O.percent = 0;, end

% out = voistat('adjusty',y,HChoice,TR,nruns, [opt] adjustmatrix);
% disp('*	voistat adjusty')

y = O.y;
TR = O.TR;
npoints = size(y,1);

if size(y,2) > size(y,1), y = y';, end

nruns = 1;
if isfield(O,'scanadjust'), if O.scanadjust, nruns = O.nruns;,end,else, O.scanadjust = 0;,end
if isfield(O,'adjustmatrix'), docustomadjust = 1;,adjustmatrix = O.adjustmatrix;,else, docustomadjust = 0;,adjustmatrix = [];,end

scanlen = npoints./nruns;
if scanlen ~= round(scanlen),warning(['Scan length is not an integer! scanlen = ' num2str(scanlen) ', rounded = ' num2str(round(scanlen))]),end
scanlen = round(scanlen);

wholeyavg = mean(y(~isnan(y)));
X = [];

if ~isfield(O,'doHP'), O.doHP = 0;,HChoice = [];,else, HChoice = O.HP;end
if O.doHP,if isempty(HChoice),error('filterAdjust: specify HP filter length in input options!'),end,end
if ~isfield(O,'doLP'), O.doLP = 0;,end
if ~isfield(O,'filtertype'), O.filtertype = 'none';,end
if ~isfield(O,'plot'), O.plot = 0;,end
if ~isfield(O,'cyclecorrection'), O.cyclecorrection = 0;,end
if ~isfield(O,'cyclecorrection2'), O.cyclecorrection2 = 0;,end
if ~isfield(O,'firstimg'), O.firstimg = 1;,end

% plot raw
if O.plot, f1 = figure('Color','w'); hold on; set(gca,'FontSize',16); plot(scale(y),'k');,title('Raw timeseries (standardized for display)'); pause(1);,end

if ((O.doHP | O.doLP) & strcmp(O.filtertype,'spm'))
	% MAKE FILTER - use_spm_filter.m
	% --------------------------------------------------------
	% Ran out of memory doing this with whole ts, so trying to make more efficient...do it scan by scan.
   if O.doHP | O.doLP
       if isfield(O,'S')
           if O.verbose, disp(	'		...voistat adjusty: Using input S filter: filter shape plot is final S matrix only'), end
           S = O.S; KL = O.S; KH = O.S;
       else

           if O.verbose, disp(	'		...voistat adjusty: making S filter'), end

           if O.doHP & O.doLP

               if O.doLP == 2

                   [S,KL,KH] = use_spm_filter(TR,scanlen,'Gaussian','specify',HChoice,TR*2);

               else

                   [S,KL,KH] = use_spm_filter(TR,scanlen,'hrf','specify',HChoice);

               end

           elseif O.doHP

               [S,KL,KH] = use_spm_filter(TR,scanlen,'none','specify',HChoice);

           elseif O.doLP == 2

               [S,KL,KH] = use_spm_filter(TR,scanlen,'Gaussian','none',[],TR*2);

           elseif O.doLP

               [S,KL,KH] = use_spm_filter(TR,scanlen,'hrf','none',[]);

           else error('Script bug.  Check filterAdjust.m')

           end

       end  % if S

   else warning('No filtering specified, though spm is filtertype.'),
      S = [];
   end

	if O.plot, figure;
   		subplot(1,3,1);imagesc(S); title('Smoothing matrix (spm) *used')
   		subplot(1,3,2);imagesc(KL); title('Low-pass filter(spm)')
      subplot(1,3,3);imagesc(KH); title('High-pass filter(spm)')
      pause(1)
      close
    end

end

if (O.doHP | O.doLP) & strcmp(O.filtertype,'none')
    error('For filtering to be used, enter ''spm'' ''cheby'' ''doug'', etc. in O.filtertype.')
end

if O.cyclecorrection2
    if O.verbose, disp('		...removing large transitions from data (presumably artifacts).'),end

    for startimg = 1:scanlen:npoints
        tmp = y(startimg:startimg+scanlen-1);
        tmp = detransition(tmp,O.plot);
        y(startimg:startimg+scanlen-1) = tmp;
    end

    if O.plot, figure(f1); hold off; plot(scale(y),'k'); hold on; title('Raw (detransitioned, std for display)');,pause(1);end
end

if O.scanadjust
	% adjust scan means
	% --------------------------------------------------------
	if O.verbose, disp('		...adjusting scan means to session mean.'),end

    if isfield(O,'X')
        X = O.X;
    else
	    index = 1;
	    for startimg = 1:scanlen:npoints
		    X(startimg:startimg+scanlen-1,index) = 1;
		    index = index + 1;
	    end
    end

	if size(X,1) > size(y,1)
		disp(['	* * * voistat.m adjusty: WARNING: X is longer than y.  truncating.'])
		X = X(1:length(y),:);
	end

   %y = y - (X * (X \ y));	% y-xB to adjust to mean.
   y = y - X * pinv(X) * y;

   if O.plot, figure;
      imagesc(X); title('Model matrix for effects to regress out.')
      pause(1),close

      figure(f1); hold off; plot(scale(y),'k'); hold on; title('Session means removed, std for display)');,pause(1);
   end

end

if O.firstimg
    if O.verbose, disp('		...Setting values of first 2 images in each run to timeseries mean value.'),end
    tmp = 1:scanlen:length(y);
    tmp = [tmp 2:scanlen:length(y)];

    % first adjust to session means - rough estimate
    tmpn = 1:length(y); tmpn(tmp) = [];
    y(tmp) = mean(y(tmpn));

    % then adjust to expected
    tmpx = (1:length(y))';
    [P] = polyfit(tmpx,y,3);
    yy = polyval(P,tmpx);
    y(tmp) = yy(tmp);

    if O.plot, figure(f1); plot(scale(y),'m'); hold on;
        plot(scale(yy));
        title('Raw (1-2 removed, std for display)');,
    end
end

switch O.filtertype

case 'spm'
% filter each scan
% --------------------------------------------------------
if ~isempty(S)
   filty = [];
	if O.verbose, disp('		...applying to each session: high-pass filtering.'),end

	% whos y

	for startimg = 1:scanlen:npoints
		try
			sessy = y(startimg:startimg+scanlen-1,:);
		catch
			startimg
			whos y
			error(['	voistat adjusty: ERROR.  y is too short, npoints is wrong, or numscans is wrong? Check GUI values.'])
		end

		try
			sessy = S * sessy;
		catch
			try
				sessy = (S * sessy')';
			catch
				disp('		can''t multiply S and sessy: matrix dimensions not correct. Try transposing y?')
				whos S
				whos sessy
				disp(['Scanlength = ' num2str(scanlen)])
				disp(['total time pts = ' num2str(npoints)])
				disp(['Num Sessions = ' num2str(nruns)])
				error('Exiting now...')
			end
		end

		filty = [filty; sessy];
   end

% whos filty

   if O.plot
      figure(f1);hold off; plot(scale(y),'b','LineWidth',1);, hold on; plot(scale(filty),'r','LineWidth',1);
      legend({'Before filtering' 'After filtering'})

      %pause(1),close
   end

   y = filty;

end

   % Old custom filter stuff.
       %Wn = .35;
       %[B,A] = cheby2(3,40,Wn);
       %y = filter(B,A,y);

% Doug's custom fourier LP filter.
%len = length(y) /2;
%ind = [1:len] ./ len;
%ff = 1./(1+exp((ind-.30)./0.05));
%ff2 = [ff ff(len:-1:1)];
%y = real(ifft(fft(y).*ff2'));

case 'cheby'
if ~isempty(HChoice)
	nyquist = TR/2;	% in seconds, not frequencies TR/nyquist = .5
	Wn = [TR/HChoice .35];
	if O.verbose, disp([		'voistat adjusty: Chebyshev filter freqency is ' num2str(Wn)]),end
	for startimg = 1:scanlen:npoints
		sessy = y(startimg:startimg+scanlen-1);
 		[B,A] = cheby2(1,20,Wn);
		sessy = filter(B,A,y);
	end
end

case 'fourier'
if ~isempty(HChoice)
	filty = []; nyquist = (1/TR)/2; lowerlim = 1/HChoice(1);
	hzfreqs = (1:scanlen) ./ (scanlen * TR);
	if length(HChoice) == 1
        ftopass = hzfreqs >= lowerlim & hzfreqs <= nyquist;
    else
        ftopass = hzfreqs >= lowerlim & hzfreqs <= 1./(HChoice(2));
    end
	if O.verbose, disp(['		...applying to each session: notch filtering: ' num2str(lowerlim) ' to ' num2str(nyquist) ' Hz.']),end
	for startimg = 1:scanlen:npoints
		myfft = fft(y(startimg:startimg+scanlen-1));
		mypass = zeros(1,scanlen);
		mypass(ftopass) = myfft(ftopass);
		% plot(abs(myfft));hold on;plot(abs(mypass),'rx');
		sessy = real(ifft(mypass))';
		filty = [filty; sessy];
	end
	y = filty;
    %hold on; plot(y,'r');
end

case 'fouriernotch'
if ~isempty(HChoice)
	filty = []; nyquist = (1/TR)/2; lowerlim = 1/HChoice(1);
	hzfreqs = (1:scanlen) ./ (scanlen * TR);
    ftopass = ones(size(hzfreqs));
    ftopass(hzfreqs > lowerlim & hzfreqs < 1./(HChoice(2))) = 0;
    ftopass = find(ftopass);

	if O.verbose, disp(['		...applying to each session: notch filtering: ' num2str(lowerlim) ' to ' num2str(nyquist) ' Hz.']),end
	for startimg = 1:scanlen:npoints
		myfft = fft(y(startimg:startimg+scanlen-1));
		mypass = zeros(1,scanlen);
		mypass(ftopass) = myfft(ftopass);
		%figure; plot(abs(myfft));hold on;plot(abs(mypass),'rx');
		sessy = real(ifft(mypass))';
		filty = [filty; sessy];
	end
	y = filty;
    %hold on; plot(y,'r');
end

case 'bandpass'
if ~isempty(HChoice)
	filty = []; nyquist = (1/TR)/2; lowerlim = 1/HChoice(1);
	hzfreqs = (1:scanlen) ./ (scanlen * TR);
    ftopass = zeros(size(hzfreqs));
    ftopass(hzfreqs > lowerlim & hzfreqs < 1./(HChoice(2))) = 1;
    ftopass = find(ftopass);

	if O.verbose, disp(['		...applying to each session: band pass filtering: ' num2str(lowerlim) ' to ' num2str(1./(HChoice(2))) ' Hz.']),end
	for startimg = 1:scanlen:npoints
		myfft = fft(y(startimg:startimg+scanlen-1));
		mypass = zeros(1,scanlen);
		mypass(ftopass) = myfft(ftopass);
		%figure; plot(abs(myfft));hold on;plot(abs(mypass),'rx');
		sessy = real(ifft(mypass))';
		filty = [filty; sessy];
	end
	y = filty;
    %hold on; plot(y,'r');
end

case 'Luis'
   if O.plot
      y = luisFilter(y,O.TR,.35,'v');pause(1),close,pause(1),close
   else
      y = luisFilter(y,O.TR,.35);
   end

case 'none'

end % end switch

if O.cyclecorrection
    % check for bimodal 'cycling' of mean BOLD signal within sessions, and
    % correct if necessary, before removing scan means, etc.
    % if we find problems, then re-do session-centering
	% --------------------------------------------------------
    if O.verbose, disp('		...checking for normally distributed data (to avoid bimodal cycling observed in some data).'),end
	cycor = []; ind = 1; classes = [];
    for startimg = 1:scanlen:npoints
		   tmp = y(startimg:startimg+scanlen-1);

           [tmp,IDX] = mixture_model(tmp,O.plot);

           if max(IDX) > 1 & ~isempty(S)
               tmp = S * tmp;
               if O.plot, subplot(1,3,1),hold on;plot(tmp,'g'),end
           end

           y(startimg:startimg+scanlen-1) = tmp - mean(tmp);
           classes(ind) = max(IDX);

            ind = ind+1;
    end
    cycor = classes > 1;

        if O.verbose,

            fprintf(1,'		...%3.0f Sessions were not normally distributed.\n',sum(cycor)),
            if any(cycor), fprintf(1,['\t\t\tClasses by session: ' num2str(classes) '\n']),end
        else
            if any(cycor), fprintf(1,['\t\t\tClasses by session: ' num2str(classes) '\n']),end
        end
end

% do overall trimming here
% --------------------------------------------------------
if O.trimts
    [y,ntrimmed] = trimts(y,O.trimts,[],1);
    if O.verbose,fprintf(1,['\t\t...Windsorized ' num2str(ntrimmed) ' points from overall timeseries to ' num2str(O.trimts) ' std. deviations.\n']), end

    if O.plot
       figure(f1); plot(scale(y),'k--');, legend({'Before filtering' 'After filtering' 'After filt+trimming'})
    end
    %tpts = abs(y) > O.trimts .* std(y);
    %y(tpts) = NaN;
end

if O.lindetrend
     % linear detrending at specified knot points

    if O.verbose, fprintf(1,['\t\t...piecewise linear detrending requested with ' num2str(length(O.lindetrend)) ' breakpoints.\n']),end
    % set bp to [] for removal of one linear trend from the whole column
    if O.lindetrend == 1, O.lindetrend = [];, end
    % if a single number, interpret as knotpoints every n elements.

    if length(O.lindetrend) == 1, O.lindetrend = O.lindetrend:O.lindetrend:length(y);, end

    if O.plot,
        figure; plot(y);,
        set(gcf,'Position',[144         542        1306         504],'Color','w')
        disp(['         Breakpoints: ' num2str(O.lindetrend)])
    end
    hold on;
    y = detrend(y,'linear',O.lindetrend);
    y = y + wholeyavg;  % preserve mean of the original timeseries

    if O.plot,
        figure(f1); hold on; plot(y,'m');
        plot([O.lindetrend' O.lindetrend']',[ones(length(O.lindetrend),1) * min(y) ones(length(O.lindetrend),1)*max(y)]','k')
        legend({'Original' 'Detrended'})

        mystd = std(y) * 3;
        plot([0 length(y)],[mean(y)+mystd mean(y)+mystd],'k--')
        plot([0 length(y)],[mean(y)-mystd mean(y)-mystd],'k--')
        text(length(y)-10,mean(y)+mystd,'3 std lines')

        drawnow
        pause(1)
        %close
    end
end

if O.percent

	if O.verbose, disp(['		...converting to % change from overall mean:  original mean is ' num2str(wholeyavg)]),end
	y = (y-mean(y))*100 / wholeyavg;		% percent change from 0.

end

if docustomadjust
		if O.verbose, disp('		...adjusting for user-entered covariates -e.g. movement params'),end
 		if sum(sum(isnan(adjustmatrix))) > 0,
			disp('		...WARNING: NaN values in user covariates!  Setting these to 0.')
			adjustmatrix(isnan(adjustmatrix)) = 0;
		end
		try
			X = [adjustmatrix];
            X(:,end+1) = 1;     % add intercept
		catch
			warning('		voistat ''adjusty'': X and adjustmatrix are different sizes.  no covariates entered.')
			whos X
			whos adjustmatrix
        end
        O.custombeta = pinv(X) * y;
        y = y - X * O.custombeta;
 end

% figure; plot(abs(fft(y)),'ro');set(gca,'YLim',[0 1000])

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/smv_phase_transition_noiseless.mat';
N = 256;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function [userCount, licenseCount] = licenseuse(toolbox,verbose)
%LICENSEUSE Counts licenses currently checked out from a license server.
% [USERCOUNT, LICENSECOUNT] = LICENSEUSE(TOOLBOX,VERBOSE)
%  Returns the number of current users of the feature TOOLBOX for a Matlab
%   installation that uses a shared network license. LICENSEUSE can also
%   display the usernames and computer names of the users who have the
%   licenses checked out. LICENSEUSE uses the Flexlm(R) license manager
%   utility to get information about the network license status.
%
% TOOLBOX is the Flexlm license manager's name for each Matlab toolbox.
%   Leave empty to count all available toolboxes plus Matlab base licenses.
%   The license manager's name for a toolbox can be found by inspecting the
%   output from
%     >> licenseuse([],2);
%   Correct capitalization is required.
% VERBOSE (0|1|2) is the level of messages printed on the command line.
%   Use VERBOSE=2 to see a list of current users for each toolbox, or to
%   see a list of all available toolboxes if TOOLBOX=[].
%
% USERCOUNT     The total number of MATLAB licenses in use.
% LICENSECOUNT  The total number of licenses available to you. If TOOLBOX
%   is not included in your license, LICENSECOUNT=-1.
%  (Note: The number of licenses remaining that can be checked out for
%      TOOLBOX is equal to LICENSECOUNT - USERCOUNT)
%
% Examples:
%  To display a list of the current users of the Instrument Control Toolbox,
%   use:
%
%  >> licenseuse('Instr_Control_Toolbox',2);
%
%  To return a count of the number of current users of MATLAB on your
%   license, use
%
%  >> [usersMatlab, licensesMatlab] = licenseuse('MATLAB');
%
%
% Note that the only way to release a license is to exit Matlab.
%
%   See also LICENSEUSER, LICENSE, VER.
%
%
% Matt Hopcroft
%  mhopeng@gmail.com
%

% MH Nov 2012
% v2.11 find file for Windows or *nix
% v2.1  update for new location of lm utilities (R2010b+)
% MH Apr2012
% v2.0  counts user/licenses and returns results
%       changed "all" behaviour- counts MATLAB + toolboxes
%       added VERBOSE parameter, improved results reporting
%       fixed help file/comments
% MH Apr2011
% v1.1  added error message
% MH Oct2010
% v1.0
%

% Note: See the an updated MATLAB help page for FlexLM utilities at:
%  http://www.mathworks.fr/support/tech-notes/1300/1303.html

%% Welcome
if nargin < 2
    if nargin < 1 || isempty(toolbox) || strcmp(toolbox,'MATLAB')
        verbose=2;
    else
        verbose=1;
    end
end
if verbose >= 1, fprintf(1,'licenseuse v2.11: %s\n',datestr(clock)); end

if nargin < 1 || isempty(toolbox)
    if verbose >= 1, fprintf(1,'licenseuse: Checking for users of all MATLAB products on your license...\n'); end
    toolboxStr = '-a';
else
    if verbose >= 1, fprintf(1,'licenseuse: Checking for current users of MATLAB "%s"...\n',toolbox); end
    toolboxStr = ['-f ' toolbox];
end

cdir=cd;

%% Determine path to lmutil
if isunix % v2.11
    lmExe='lmutil';
else
    lmExe='lmutil.exe';
end

lmPath=[matlabroot filesep 'bin' filesep computer('arch') filesep]; % prior to R2010b
if ~exist([lmPath lmExe],'file')
    lmPath=[matlabroot filesep 'etc' filesep computer('arch') filesep]; % R2010b+
    if ~exist([lmPath lmExe],'file')
        fprintf(1,'ERROR: "%s" not found in bin/ or etc/ locations!\n\n',lmExe);
        return
    end
end
if verbose >= 3, fprintf(1,'License manager utilities found at:\n  %s\n',lmPath); end

cd(lmPath);

%% Execute license manager command
[ltstat, lt]=system(['.' filesep 'lmutil lmstat -c "' matlabroot filesep 'licenses' filesep 'network.lic" ' toolboxStr]);
cd(cdir);

if ltstat ~= 0
    fprintf(1,'ERROR: "%s"\n',strtrim(lt));
    if strfind(lt, 'No such file')
        fprintf(1,'On Linux you may need to upgrade the package "lsb"\n  (see %s)\n',...
            'https://www.mathworks.com/support/solutions/en/data/1-GLXUHV/index.html?product=ML&solution=1-GLXUHV');
    end
    error('OS returned error %d when executing lmutil',ltstat);
else
    if verbose >= 2, disp(lt); end % display the license manager status
end

%% Count number of users
tlines=textscan(lt,'%s','delimiter','\n'); tlines=tlines{1};
tIndC=strfind(tlines,'Total of');
licenseCount=0; % number of licenses available to you ("issued")
userCount=0;    % number of users that have checked out a license ("in use")
toolboxFound=0;
for i = 1:length(tIndC)
    if ~isempty(tIndC{i})
        if strcmp(toolboxStr,'-a') % if -a, count users of matlab and all users
            if ~isempty(strfind(tlines{i},'Users of MATLAB:'))
                mCount=textscan(tlines{i},'Users of MATLAB:  (Total of %d %*s issued; Total of %d %*s in use)',1);
                userCountMatlab=mCount{2};
                licenseCountMatlab=mCount{1};
            else
                uCount=textscan(tlines{i},'Users of %*s (Total of %d %*s issued; Total of %d %*s in use)',1);
                userCount=userCount+uCount{2};
                licenseCount=licenseCount+uCount{1};
            end
        else % else, look for only the requested toolbox
            if ~isempty(strfind(tlines{i},toolbox))
                toolboxFound=1;
                uCount=textscan(tlines{i},['Users of ' toolbox ': (Total of %d %*s issued; Total of %d %*s in use)'],1);
                userCount=uCount{2};
                licenseCount=uCount{1};
            end
        end
    end
end

%% Report results
if verbose >= 1
    if verbose >= 2, fprintf(1,'licenseuse: Results:\n'); end
    if strcmp(toolboxStr,'-a')
        fprintf(1,'  There are %d MATLAB licenses in use out of %d available to you,\n',userCountMatlab,licenseCountMatlab);
        fprintf(1,'   and %d additional toolbox licenses in use out of %d available to you.\n\n',userCount,licenseCount);
        userCount=userCount+userCountMatlab;
        licenseCount=licenseCount+licenseCountMatlab;
    else
        if toolboxFound==0
            fprintf(1,'  There do not appear to be any licenses for "%s" available to you.\n',toolbox);
            fprintf(1,'  Check for typos in the toolbox name, and/or try using the ''ver'' command,\n');
            fprintf(1,'   and/or verify that the license server is up (use VERBOSE=2).\n\n');
            licenseCount=-1;
        else
            fprintf(1,'  There are %d %s licenses in use out of %d available to you.\n\n',userCount,toolbox,licenseCount);
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function ds = setfile(ds, desc_string, varargin)
   %SETFILE associates a datasource object with a file (or directory)
   %   ds = setfile(ds, descriptor_string, [var1, var2, ... varN])
   %   where;
   %   - DS is the datasource object
   %   - DESCRIPTOR_STRING is a format statement, like that from sprintf.
   %   - var1...varN are optional variable names that are used to fill in the
   %     information from the format statement.
   %      valid variable names (of type integer, or format string of '%d')
   %         YEAR, DAY, HOUR, MINUTE, SECOND, JDAY (<-- julian day)
   %
   %      valid variable names (of type string, or format string of '%s')
   %         station, channel, location, network
   %
   %   this provides a powerful tool that allows the program to navigate a
   %   wide variety of files and directories based upon date or SCNL
   %   information.
   %
   % ------------------------------------------------------------------------
   % EXAMPLE #1:
   % Assume data files are stored in the following directories:
   %   /DATA/year2006/04/OKCF.dat
   %   /DATA/year2006/04/OKFG.dat
   %   /DATA/year2006/05/OKCF.dat
   %   ...etc...
   %
   %  To find any given peice of data requires the year, month (zero padded!),
   %  and the station name.
   %
   %  You can specify which file to use by setting the filename as follows:
   %  ds = setfile(ds,'/DATA/year%04d/%02d/%s.dat','year','month','station');
   %
   %  * notice that each percent sign (%) has a corresponding term.  That is:
   %    %04d --> 'year',   %02d --> 'month', and %s -->'station'
   %  * The numbers and symbols between the '%' and the first alphabetic
   %    character provide additional formatting details.
   %    * any number greater than zero means "use up this many spaces"
   %    * a leading zero means "fill up leading space with zeros"
   %
   %    ex. %04d -> write the number with 4 spaces, and pad with zeros.
   %
   % ------------------------------------------------------------------------
   % EXAMPLE #2:
   % The method waveform/which_archive can be supplanted by
   % the following setfile call...
   % ds = setfile(ds, ...
   %            '/iwrun/op/db/archive/archive_%04d/archive_%04d_%02d_%02d',...
   %            'year','year','month','day')
   % - - - -
   % For access to a specific, date independent, SCNL independent database:
   % ds = setfile(ds,'mydatabasename');
   %
   % see also sprintf, datasource/getfilename

   if ~exist('desc_string','var') || isempty(desc_string)
      ds.usefile = false;
      ds.file_string = {};
      ds.file_args = {};
      return
   end

   % the number of arguments should equal the number of %somethings
   % ... except where '%' is preceded by a  '\'
   n = find(desc_string == '%');
   z=0; %currently not accepting % in filename!!!
   %z = sum('\'== desc_string(n(n>2)-1));
   expectedArgumentCount = numel(n) - z;
   if numel(varargin) ~= expectedArgumentCount
      error('Datasource:FileArgumentMismatch', ...
         'the number of file arguments provided (%d) doesn''t match the number of arguments requested by the filestring %d',...
         numel(varargin),expectedArgumentCount);

   end
   varargin = lower(varargin);

   %check to be sure arguments are valid
   validVariables = {'year', 'month', 'day', 'hour', 'minute', 'second',...
      'station', 'channel','location','network', 'jday'};
   for i=1:numel(varargin)
      if ~ismember(varargin{i},validVariables)
         error('Datasource:FileArgumentError',...
            'unrecognized file argument [%s].',varargin{i});
      end
   end

   ds.file_string = desc_string;
   ds.file_args = varargin;
   ds.usefile = true;
end

author: geoscience-community-codes

repo: GISMO


text:
%kvnopen 'Opening by N times (with respect to Minkowski addition) the structuring element'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vnopen.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% Integer: n 'Size of Structure El.', default: 5: 'size of structure element'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvnopen({i, str}, {'i','';'str','';'n',5;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vnopen - Opening by N times (with respect to Minkowski addition) the structuring element
%
%  DESCRIPTION
% .I vnopen
% will do N successive openings by the same structuring element.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
%
%
%  EXAMPLES
% vnopen -i feath.xv -str b.str -n 5 -o outimage.xv
%
% Will do 5 successive openings (vnopen) of the input image "feath.xv" by the structuring element "b.str" and the resulting image will be written to  "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vnopen
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% open a connection to the VREP simulator
%
% Notes::
% - the VREP constructor needs to know where V-REP is installed.
% - This can come from the environment variable VREP
% - Or you can edit this file to something like
%    vrep=VREP('/Applications/V-REP_PRO_EDU_V3_1_3_rev2b_Mac');

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% create a connection to a running instance of V-REP
vrep = VREP();

% let's change our viewpoint, move the camera further away
%
% first we get an object to mirror the default camera
cam = vrep.camera('DefaultCamera');

% then get the camera position
p = cam.getpos

% then set the camera be twice as far away
cam.setpos(2*p)

% add the human figure Bill
bill = vrep.loadmodel('people/Walking Bill');

% get the initial pose of Bill
T = bill.getpose()

% now we can change his position
bill.setpos([0.1, 0.2, 0]);

% make him turn to his right
bill.setorient([0 0 -pi/4])

% make him lean forward a bit
bill.setorient([0 pi/8 0])

% now set him back to his initial position and orientation
bill.setpose(T);

% now we will setup some nested loops to make him shuffle around a square
orient = [0 0 0];
for j=1:4
    step = [0.05 0 0];
    bill.setorient(orient);
    % move forward
    for i=1:30
        % move relative to current pose
        bill.setpos(step, bill);
    end
    % turn at right angles
    orient(3) = orient(3) + pi/2;
end

% remove the model from the scene
bill.remove();

% now we will load a scene into VREP that contains Bill, a tree, a mirror
% and a camera
vrep.loadscene(which('scene1.ttt'), 'local');

% we will get a handle to the camera
camera = vrep.camera('Vision_sensor');

% now we start the simulator, it must be running in order for the
% camera to render

vrep.simstart();

% now we loop as the camera creeps forward (in the z-direction) and we
% display the frames
for i=1:20
    camera.setpos([0 0 0.05], camera);
    camera.grab();
    drawnow
end

% stop the simulation
vrep.simstop();

% clear the scene, remove all objects
vrep.clearscene();

% add a model of an ABB IRB140 robot
model = vrep.loadmodel('robots/non-mobile/ABB IRB 140');

% and create a VREP_arm object
arm = vrep.arm(model);

% set the pose of the robot's base
arm.setpose( transl(0.3, 0.4, 0));

% get the joint angles
q = arm.getq();

% now move the robot's joints
arm.setq(q+0.3);

% clear the scene, remove all objects
vrep.clearscene();

% close the connection to the VREP simulator
clear vrep

author: petercorke

repo: robotics-toolbox-matlab


text:
classdef ComparingSuperEllipseVsRectangle < handle

    properties (Access = private)
       mesh
       dataRes
       m1S
       m2S
       m1R
       m2R
    end

    properties (Access = private)
        path
        fileNameMesh
        fileResName
        pNorm
        microCase
    end

    methods (Access = public)

        function obj = ComparingSuperEllipseVsRectangle()
            obj.init()
            obj.wrapResAndMesh();
            obj.computeDesignParams();
            obj.computeStresses();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.microCase = 'SuperEllipse';'Rectangle';
            meshType  = 'Medium';'Small';
            finalIter = '589';'485';'405';'250';'296';
            fCase = [obj.microCase,'Rotation',meshType];
            obj.path = ['/media/alex/My Passport/LatticeResults/StressNorm',fCase,'/'];
            obj.fileResName  = ['ExperimentingPlot',finalIter];
            obj.fileNameMesh = ['CantileverSquare',meshType];
            obj.pNorm = 2:2:32;
        end

        function wrapResAndMesh(obj)
            s.folderPath = obj.path;
            s.fileName   = obj.fileResName;
            w = WrapperMshResFiles(s);
            w.compute();
            obj.mesh    = w.mesh;
            obj.dataRes = w.dataRes;
        end

        function computeDesignParams(obj)
            switch obj.microCase
                case 'Rectangle'
                   obj.computeDesignParamsFromrectangleExample();
                case 'SuperEllipse'
                   obj.computeDesignParamsFromSuperEllipseExample();
            end
        end

        function computeDesignParamsFromSuperEllipseExample(obj)
            obj.m1S = obj.dataRes.DesignVar1;
            obj.m2S = obj.dataRes.DesignVar2;
            [obj.m1R,obj.m2R] = obj.computeRectangleM1M2FromSuperEllipseM1M2(obj.m1S,obj.m2S);
        end

        function computeDesignParamsFromrectangleExample(obj)
            obj.m1S = obj.dataRes.DesignVar1;
            obj.m2S = obj.dataRes.DesignVar2;
            [obj.m1R,obj.m2R] = obj.computeRectangleM1M2FromSuperEllipseM1M2(obj.m1S,obj.m2S);
        end

        function computeStresses(obj)
           valueR = obj.computeStressNormRectangle();
           valueS = obj.computeStressNormSuperEllipse();
           f = figure();
           hold on
           p{1} = plot(obj.pNorm,valueR);
           p{2} = plot(obj.pNorm,valueS);
           legend({'Rectangle','SuperEllipse'},'Location','Best');
           xlabel('p');
           ylabel('||\sigma||_p');
           pP = plotPrinter(f,p);
           fName = '/home/alex/Dropbox/GregMeeting30Octubre/Comparison';
           pP.print(fName)
        end

        function value = computeStressNormSuperEllipse(obj)
            s.m1            = obj.m1S;
            s.m2            = obj.m2S;
            s.alpha         = obj.dataRes.AlphaGauss';
            s.vademecumName = 'SuperEllipseQOptAnalytic';
            s.mesh          = obj.mesh;
            s.fileName      = obj.fileNameMesh;
            s.pNorm         = obj.pNorm;
            sC = StressNormFromVademecumComputer(s);
            value = sC.compute();
        end

        function value = computeStressNormRectangle(obj)
            s.m1            = obj.m1R;
            s.m2            = obj.m2R;
            s.alpha         = obj.dataRes.AlphaGauss';
            s.vademecumName = 'SuperEllipseQMax';
            s.mesh          = obj.mesh;
            s.fileName      = obj.fileNameMesh;
            s.pNorm         = obj.pNorm;
            sC = StressNormFromVademecumComputer(s);
            value = sC.compute();
        end

    end

    methods (Access = private, Static)

        function [m1R,m2R] = computeRectangleM1M2FromSuperEllipseM1M2(m1S,m2S)
            sE   = SuperEllipseParamsRelator;
            s.type  = 'Optimal';
            s.m1    = m1S;
            s.m2    = m2S;
            sM = SmoothingExponentComputer.create(s);
            qS(:,1) = sM.compute();
            rho  = sE.rho(m1S,m2S,qS);
            xi   = sE.xi(m1S,m2S);
            s.type  = 'Given';
            s.q    = 32*ones(size(m1S));
            sM = SmoothingExponentComputer.create(s);
            qR(:,1) = sM.compute();
            m1R = sE.mx(xi,rho,qR);
            m2R = sE.my(xi,rho,qR);
        end

        function [m1S,m2S] = computeSuperEllipseM1M2FromRectangleM1M2(m1R,m2R)
            sE   = SuperEllipseParamsRelator;
            rho  = sE.rho(m1R,m2R,q);
            xi   = sE.xi(m1R,m2R);
            s.type  = 'Optimal';
            s.m1    = [];
            s.m2    = [];
            sM = SmoothingExponentComputer.create(s);
            qS(:,1) = sM.computeQ(rho,xi);
            m1S = sE.mx(xi,rho,qS);
            m2S = sE.my(xi,rho,qS);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function LoadFullFish(hfig,i_fish,isFullData,hdf5_dir,mat_dir)
M_fish_set = getappdata(hfig,'M_fish_set');
fishset = M_fish_set(i_fish);
setappdata(hfig,'fishset',fishset);

%% load data from file
disp(['loading fish ' num2str(i_fish) '...']);
loadStart = tic;
isLoadTS = 1;
if ~exist('hdf5_dir','var') || ~exist('mat_dir','var'),
    data_masterdir = getappdata(hfig,'data_masterdir');
    data_dir = fullfile(data_masterdir,['subject_' num2str(i_fish)]);
    if exist('isFullData','var'),
        if isFullData % true
             hdf5_dir = fullfile(data_dir,'TimeSeries.h5');
        else % % %isFullData == -1,
            hdf5_dir = fullfile(data_dir,'TimeSeries.h5');
            absIX = h5read(hdf5_dir,'/absIX');
            setappdata(hfig,'absIX',absIX);
            isLoadTS = 0; % don't load Time series
%         else % false
%             hdf5_dir = fullfile(data_dir,'TimeSeries_half.h5');
        end
    else % default without isfulldata input is full-data
        hdf5_dir = fullfile(data_dir,'TimeSeries.h5');
    end
    mat_dir = fullfile(data_dir,'data_full.mat');
end

% load 'data'
load(mat_dir,'data'); % struct with many fields
names = fieldnames(data); % cell of strings
for i = 1:length(names),
    setappdata(hfig,names{i},eval(['data.',names{i}]));
end

if isLoadTS,
    % load time series (hdf5 file)
    CellResp = h5read(hdf5_dir,'/CellResp');
    CellRespZ = h5read(hdf5_dir,'/CellRespZ');
    CellRespAvr = h5read(hdf5_dir,'/CellRespAvr');
    CellRespAvrZ = h5read(hdf5_dir,'/CellRespAvrZ');
    absIX = h5read(hdf5_dir,'/absIX');

    setappdata(hfig,'CellResp',CellResp);
    setappdata(hfig,'CellRespZ',CellRespZ);
    setappdata(hfig,'CellRespAvr',CellRespAvr);
    setappdata(hfig,'CellRespAvrZ',CellRespAvrZ);
    setappdata(hfig,'absIX',absIX);
end

loadTime = toc(loadStart);
disp([num2str(loadTime) ' sec']);

end

author: xiuyechen

repo: FishExplorer


text:
function brewermap_plot()
% Simple plot of all ColorBrewer colorscheme nodes in one figure.
%
% (c) 2014-2022 Stephen Cobeldick
%
%%% Syntax:
% brewermap_plot()
%
% See also BREWERMAP BREWERMAP_VIEW CUBEHELIX MAXDISTCOLOR
% LBMAP PARULA LINES RGBPLOT COLORMAP COLORBAR PLOT PLOT3 AXES SET

[mcs,nmn,pyt] = brewermap('list');
%
persistent cbh axh
%
xmx = max(nmn);
ymx = numel(pyt);
%
if ishghandle(cbh)
	figure(cbh);
	delete(axh);
else
	cbh = figure('HandleVisibility','callback', 'IntegerHandle','off',...
		'NumberTitle','off', 'Name',mfilename,'Color','white',...
		'MenuBar','figure', 'Toolbar','none', 'Tag',mfilename);
	set(cbh,'Units','pixels')
	pos = get(cbh,'Position');
	pos(1:2) = pos(1:2) - 123;
	pos(3:4) = max(pos(3:4),[842,532]);
	set(cbh,'Position',pos)
end
%
axh = axes('Parent',cbh, 'Color','none',...
	'XTick',0:xmx, 'YTick',0.5:ymx, 'YTickLabel',mcs, 'YDir','reverse');
title(axh,'ColorBrewer Color Schemes (brewermap.m)', 'Interpreter','none')
xlabel(axh,'Scheme Nodes')
ylabel(axh,'Scheme Name')
axf = get(axh,'FontName');
%
for y = 1:ymx
	N = nmn(y);
	M = brewermap(N,mcs{y});
	for x = 1:N
		patch([x-1,x-1,x,x],[y-1,y,y,y-1],1, 'FaceColor',M(x,:), 'Parent',axh)
	end
	text(xmx+0.1,y-0.5,pyt{y}, 'Parent',axh, 'FontName',axf)
end
%
drawnow()
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%brewermap_plot

author: DrosteEffect

repo: BrewerMap


text:
function gpts = APPcreateGroundPoints(lmap, gndmap)
% gpts = APPcreateGroundPoints(lmap, gndmap)
% Gets the ground-vertical boundary from the label map
% high y --> low in image, get max y for each x
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Permission granted to non-commercial enterprises for
% modification/redistribution under GNU GPL.
% Current Version: 1.0  09/30/2005

height = size(lmap, 1);
width = size(lmap, 2);

count = 0;
[y, x] = find(lmap(1:end-1, :) & gndmap(2:end, :));
[y2, x2] = find(lmap(end, :));
gpts = ([(y(:)+1) x(:) ; y2(:)+(height-1) x2(:)]);
gpts(:, 1) = (gpts(:, 1)-1)/(height-1);
gpts(:, 2) = (gpts(:, 2)-1)/(width-1);

author: Cloud-CV

repo: object-proposals


text:
function linplus_test345 ( )

%*****************************************************************************80
%
%% TEST345 tests R8GE_FSS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  nb = 3;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST345\n' );
  fprintf ( 1, '  For a matrix in general storage,\n' );
  fprintf ( 1, '  R8GE_FSS factors and solves multiple linear systems.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8ge_random ( n, n, seed );
%
%  Set the desired solutions.
%
  x(1:n,1) = 1.0;
  x(1:n,2) = 1:n;
  x(1:n,3) = mod ( 0 : n - 1, 3 ) + 1;

  b(1:n,1:nb) = a(1:n,1:n) * x(1:n,1:nb);
%
%  Factor and solve the system.
%
  x = r8ge_fss ( n, a, nb, b );

  r8ge_print ( n, nb, x, '  Solution:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [RHS, RHSx, RHSy, RHSz] = convectionTvdRHS(u, phi, FL)
% This function uses the TVD scheme to discretize a
% convection term in the form $\grad (u \phi)$ where u is a face vactor
% It also returns the x, y, x parts of the matrix of coefficient.
%
% SYNOPSIS:
%   [RHS, RHSx, RHSy, RHSz] = convectionTvdRHS(u, phi, FL)
%
% PARAMETERS:
%   u  - velocity vector, FaceVariable
%   phi  - value of phi from the previous time step or iteration, CellVariable
%   FL  - Flux Limiter function
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

d = u.domain.dimension;
switch d
    case 1
        RHS = convectionTvdRHS1D(u, phi, FL);
    case 1.5
        RHS = convectionTvdRHSCylindrical1D(u, phi, FL);
    case 1.8
        RHS = convectionTvdRHSSpherical1D(u, phi, FL);
    case 2
        [RHS, RHSx, RHSy] = convectionTvdRHS2D(u, phi, FL);
    case 2.5
        [RHS, RHSx, RHSy] = convectionTvdRHSCylindrical2D(u, phi, FL);
    case 2.8
        [RHS, RHSx, RHSy] = ...
            convectionTvdRHSRadial2D(u, phi, FL);
    case 3
        [RHS, RHSx, RHSy, RHSz] = convectionTvdRHS3D(u, phi, FL);
    case 3.2
        [RHS, RHSx, RHSy, RHSz] = ...
            convectionTvdRHSCylindrical3D(u, phi, FL);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function Corr = tapas_Cov2Corr(Cov)
% Converts a covariance matrix into a correlation matrix
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Check if Cov is symmetric
if any(any(Cov'~=Cov))
    error('tapas:hgf:Cov2Corr:MatNotSymm', 'Input matrix is not symmetric.');
end

% Check if Cov is positive semi-definite
if any(isinf(Cov(:))) || any(isnan(Cov(:))) || any(eig(Cov)<0)
    error('tapas:hgf:Cov2Corr:MatNotPosDef', 'Input matrix is not positive semi-definite.');
end

sdev = sqrt(diag(Cov));
Norm = sdev * sdev';
Corr = Cov./Norm;

return;

author: translationalneuromodeling

repo: tapas


text:
classdef ODE113Integrator < AbstractFirstOrderIntegrator
    %ODE45Integrator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        options(1,1) BuiltInIntegratorOptions = BuiltInIntegratorOptions();
        integratorEnum = IntegratorEnum.ODE113;
    end

    methods
        function obj = ODE113Integrator(options)
            if(nargin > 0)
                obj.options = options;
            else
                obj.options = BuiltInIntegratorOptions();
            end
        end

        function [t,y,te,ye,ie] = integrate(obj, odefun, tspan, y0, evtsFunc, odeOutputFun)
            odeSetOptions = obj.options.getIntegratorOptions();
            optionsToUse = odeset(odeSetOptions, 'Events',evtsFunc, 'OutputFcn',odeOutputFun);

            [t,y,te,ye,ie] = ode113(odefun, tspan, y0, optionsToUse);
        end

        function options = getOptions(obj)
            options = obj.options;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef Rx < adi.AD9361.Rx
    % adi.FMComms2.Rx Receive data from the FMComms2 evaluation platform
    %   The adi.FMComms2.Rx System object is a signal source that can
    %   receive complex data from the FMComms2.
    %
    %   rx = adi.FMComms2.Rx;
    %   rx = adi.FMComms2.Rx('uri','192.168.2.1');
    %
    %   Product Page
    %
    %   See also adi.AD9361.Rx
    methods
        %% Constructor
        function obj = Rx(varargin)
            % Returns the matlabshared.libiio.base object
            coder.allowpcode('plain');
            obj = obj@adi.AD9361.Rx(varargin{:});
        end
    end

end

author: analogdevicesinc


text:
function createTreeView(handles, FT_Width, FT_Height)
%
% 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/>.

%     Para.Rigid_Single_Modality = 'MS';

    Para.Rigid = 'MI';

    Para.Non_Rigid.OpticalFlow = 'OpticalFlow';
    Para.Non_Rigid.Demons3D = 'Demons3D';
%     Para.Non_Rigid.SymmetricDemons3D = 'SymmetricDemons';
%     Para.Non_Rigid.BSplineMethod = 'BSplineMethod';
%     Para.Non_Rigid.LevelSetMethod = 'LevelSetMethod';

    expFcn = createParaTree(Para, handles, [0, 1-FT_Height, FT_Width, FT_Height]);
    handles = guidata(handles.mainframe);
    set(handles.FTreeView, 'Visible', true);

    guidata(handles.mainframe, handles);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


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

function test_mocov_is_absolute_path_basics()
    if ispc()
        assertTrue(mocov_is_absolute_path('C:\\'));
        assertTrue(mocov_is_absolute_path('C:'));
        assertTrue(mocov_is_absolute_path('C:\\aa\\bb'));
        assertTrue(mocov_is_absolute_path('D:/bb/cc'));
        assertFalse(mocov_is_absolute_path('C\\'));
        assertFalse(mocov_is_absolute_path('/bb/cc'));
        assertFalse(mocov_is_absolute_path('./bb/cc'));
        assertFalse(mocov_is_absolute_path('../bb/cc'));
    else
        assertFalse(mocov_is_absolute_path('C:\\'));
        assertFalse(mocov_is_absolute_path('C:'));
        assertTrue(mocov_is_absolute_path('/bb/cc'));
        assertTrue(mocov_is_absolute_path('/'));
        assertFalse(mocov_is_absolute_path('./bb/cc'));
        assertFalse(mocov_is_absolute_path('../bb/cc'));
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [imageSpm4DArray, selectionArray] = split_into_MrImageSpm4D(this, dimLabelsSpm4D)
% Splits high dimensional MrImage (>4) into array of 4D MrImageSpm4D
% objects to be used for SPM interfacing
%
%   Y = MrImage(); % >4 dimensions
%   imageSpm4DArray = Y.split_into_MrImageSpm4D(inputs)
%
% This is a method of class MrImage.
%
% IN
%   dimLabelsSpm4D      cell(1,4) of dimLabels that shall represent the 1st
%                       four dimensions that will be in each element of the
%                       output 4D SPM-compatible image array
%                       default: {'x','y','z','t'}
%                       NOTE: If less than 4 dimensions are specified, the
%                             remaining up to 4th dimension is included in
%                             the 4D images
%                             e.g.
%                               specified: {'x','y','z','t'}
%                               actual dimLabels:  {'x','z','echo', 'y', 'coil'}
%                               => {'x','y','z','echo'} will be in 4D image
%                       NOTE2: If certain labels do not exist, the
%                              existing dimensions later in the array are
%                              moved up in their dimension
%                               e.g.,
%                               specified: {'x','y','z','t'}
%                               actual dimLabels:  {'x','y','coil', 't', 'echo'}
%                               => {'x','y','t','coil'} will be in 4D image
%
% OUT
%   imageSpm4DArray         cell(nElementsSplitDim1, ..., nElementsSplitDimN) of
%                           MrImageSpm4D, split along splitDims,
%                           e.g. 6-D MrImage of 4D time series with
%                           16 coils and 3 echoes will result in a 16x3
%                           cell of MrImageSpm4D
%   selectionArray
%                           cell(nElementsSplitDim1, ..., nElementsSplitDimN) of
%                           selections, defined as cells containing propertyName/value
%                           pairs over split dimensions, e.g.
%                           {'t', 5, 'dr', 3, 'echo', 4}
%
% EXAMPLE
%   split_into_MrImageSpm4D
%
%   See also MrImage

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2018-05-03
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if nargin < 2
    dimLabelsSpm4D = {'x','y','z','t'};
end

%% permute dimLabelsSpm4D as first 4 dimensions
permutedThis = this.permute(this.dimInfo.get_dim_index(dimLabelsSpm4D));

%% split along all higher dimensions, keep chunks of 1st four together
[imageSpm4DArray, selectionArray] = permutedThis.split('splitDims', 5:this.dimInfo.nDims, ...
    'doRemoveDims', false); % dimInfo for singleton dimension kept to facilitate later recombination; MrImageSpm4D is more of an internal class
nImages = numel(imageSpm4DArray);
for iImage = 1:nImages
    imageSpm4DArray{iImage} = imageSpm4DArray{iImage}.recast_as_MrImageSpm4D;
end


text:
function cvx_optval = huber_pos( x, M, t ) %#ok

%HUBER   Internal cvx version.

%
% Check arguments
%

narginchk(1,3);
if ~cvx_isconvex( x ),
    error( 'Disciplined convex programming error:\n    HUBER_POS is convex and nondecreasing in X, so X must be convex.', 1 ); %#ok
end
if nargin < 2,
    M = 1;
elseif isa( M, 'cvx' ),
    error( 'Second argument must be numeric.' );
elseif ~isreal( M ) || any( M( : ) <= 0 ),
    error( 'Second argument must be real and positive.' );
end
if nargin < 3,
    t = 1;
elseif ~isreal( t ),
    error( 'Third argument must be real.' );
elseif cvx_isconstant( t ) && nnz( cvx_constant( t ) <= 0 ),
    error( 'Third argument must be real and positive.' );
elseif ~cvx_isconcave( t ),
    error( 'Disciplined convex programming error:\n    HUBER_POS is convex and nonincreasing in T, so T must be concave.', 1 ); %#ok
end
sz = cvx_size_check( x, M, t );
if isempty( sz ),
    error( 'Sizes are incompatible.' );
end

%
% Compute result
%

v = []; w = [];
cvx_begin separable
    variables v( sz ) w( sz )
    minimize( quad_over_lin( w, t, 0 ) + 2 .* M .* v )
    x <= w + v; %#ok
    w <= M * t; %#ok
    v >= 0; %#ok
cvx_end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%EXnmr_cgls_mrnsd Example script, 2D NMR relaxometry

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% Clear workspace and command window.
clear, clc

% Choose if you would like to see the results displayed in a single figure
% window ('subplots') or in multiple figure windows ('manyplots').
dispres = 'subplots';
% dispres = 'manyplots';

LW = 2;  % Plot line width.
MS = 10; % Size of markers on plots.

rng(0);  % Make sure this test is repeatable.

% Define the test problem.
n = 64;
NoiseLevel = 0.05;
[A, b, x, ProbInfo] = PRnmr(n);
bn = PRnoise(b, NoiseLevel);

% Compute a CGLS solution, specifying the noise level for the discrepancy
% principle, and use true solution to compute error norms.
eta = 1.02;
options = IRset('x_true', x, 'NoiseLevel', NoiseLevel, 'eta', eta, 'NoStop', 'on');
[x_cgls, IterInfo_cgls] = IRcgls(A, bn, 1:500, options);

% Compute MRNSD reconstruction.
K = [1, 100:100:20000];
[x_mrnsd, IterInfo_mrnsd] = IRmrnsd(A, bn, K, options);

% Display the reconstructions;
% uncomment as appropriate to avoid displaying titles and legends.
if strcmp(dispres, 'subplots')
    figure(1), clf
    subplot(3,3,1), PRshowx(x, ProbInfo), colormap hsv
    title('True solution','interpreter','latex','fontsize',24)
    set(gca,'fontsize',10)
    %
    subplot(3,3,4), PRshowx(IterInfo_cgls.BestReg.X, ProbInfo), colormap hsv
    title('Best CGLS solution','interpreter','latex','fontsize',24)
    set(gca,'fontsize',10)
    %
    subplot(3,3,2), semilogy(IterInfo_mrnsd.Enrm,'linewidth',1.5)
    hold on
    hl = legend('{\tt info.Enrm}');
    % hl = legend('IRmrnsd error');
    set(hl,'interpreter','latex','fontsize',12)
    semilogy(IterInfo_mrnsd.BestReg.It, IterInfo_mrnsd.BestReg.Enrm, 'ro', 'LineWidth', 1.5, 'MarkerSize', 6)
    axis([0 max(K) 0.08 1.2])
    semilogy(IterInfo_mrnsd.StopReg.It, IterInfo_mrnsd.StopReg.Enrm, 'ms', 'LineWidth', 1.5, 'MarkerSize', 6)
    set(gca,'fontsize',12)
    title('Error history','interpreter','latex','fontsize',24)
    set(gca,'fontsize',10)
    %
    subplot(3,3,5), PRshowx(IterInfo_mrnsd.BestReg.X, ProbInfo)
    title(['Best MRNSD sol., $k$ = ',num2str(IterInfo_mrnsd.BestReg.It)],...
    'interpreter','latex','fontsize',24)
    set(gca,'fontsize',10)
    %
    subplot(3,3,3)
    semilogy(K,IterInfo_mrnsd.Rnrm(K),'-',K,eta*NoiseLevel*ones(size(K)),'--','linewidth',1.5)
    hl = legend('{\tt info.Rnrm}','{\tt eta*NoiseLevel}','location','northwest');
    % hl = legend('IRmrnsd residual','{\tt eta*NoiseLevel}','location','northwest');
    set(hl,'interpreter','latex')
    axis([0 max(K) 0.045 0.08])
    title('Residula history','interpreter','latex','fontsize',24)
    set(gca,'fontsize',10)
    %
    subplot(3,3,6), PRshowx(IterInfo_mrnsd.StopReg.X, ProbInfo)
    title(['DP MRNSD sol., $k$ = ',num2str(IterInfo_mrnsd.StopReg.It)],...
    'interpreter','latex','fontsize',24)
    set(gca,'fontsize',10)
elseif strcmp(dispres, 'manyplots')
    figure(1), clf
    PRshowx(x,ProbInfo), colormap hsv
    title('True solution','interpreter','latex','fontsize',24)
    set(gca,'fontsize',24)
    %
    figure(2), clf
    PRshowx(IterInfo_cgls.BestReg.X, ProbInfo), colormap hsv
    title('Best CGLS solution','interpreter','latex','fontsize',24)
    set(gca,'fontsize',24)
    %
    figure(3), clf
    semilogy(IterInfo_mrnsd.Enrm,'linewidth',LW)
    hold on
    semilogy(IterInfo_mrnsd.BestReg.It, IterInfo_mrnsd.BestReg.Enrm, 'ro', 'LineWidth', LW, 'MarkerSize', MS)
    semilogy(IterInfo_mrnsd.StopReg.It, IterInfo_mrnsd.StopReg.Enrm, 'ms', 'LineWidth', LW, 'MarkerSize', MS)
    axis([0 max(K) 0.08 1.2])
    % hl = legend('{\tt info.Enrm}','location','North');
    % set(hl,'interpreter','latex','fontsize',24)
    hl = legend('IRmrnsd error','optimal stopping iteration', ...
      'DP stopping iteration','location','North');
    set(hl,'interpreter','latex','fontsize',16)
    title('Error history','interpreter','latex','fontsize',24)
    set(gca,'fontsize',30)
    %
    figure(4), clf
    PRshowx(IterInfo_mrnsd.BestReg.X, ProbInfo), colormap hsv
    title(['Best MRNSD sol., $k$ = ',num2str(IterInfo_mrnsd.BestReg.It)],...
    'interpreter','latex','fontsize',24)
    set(gca,'fontsize',24)
    %
    figure(5), clf
    semilogy(K,IterInfo_mrnsd.Rnrm(K),'-',K,eta*NoiseLevel*ones(size(K)),'--','linewidth',LW)
    % hl = legend('{\tt info.Rnrm}','{\tt eta*NoiseLevel}','location','north');
    hl = legend('IRmrnsd residual','{\tt eta*NoiseLevel}','location','north');
    set(hl,'interpreter','latex','fontsize',24)
    title('Residual history','interpreter','latex','fontsize',24)
    axis([0 max(K) 0.045 0.08])
    set(gca,'fontsize',30)
    %
    figure(6), clf
    PRshowx(IterInfo_mrnsd.StopReg.X,ProbInfo), colormap hsv
    title(['DP MRNSD sol., $k$ = ',num2str(IterInfo_mrnsd.StopReg.It)],...
    'interpreter','latex','fontsize',24)
    set(gca,'fontsize',24)
end

return

% A number of instructions useful to save the displayed figures follow;
% the defualt is not to execute them. If you wish to save the displayed
% figures in the dedicated 'Results' folder, please comment the above
% return statement
oldcd = cd;
if strcmp(dispres, 'subplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    figure(1), print -dpng -r300 EXnmr
elseif strcmp(dispres, 'manyplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    % save as eps
    figure(1), print -depsc -r300 EXnmr_a.eps
    figure(2), print -depsc -r300 EXnmr_d.eps
    figure(3), print -depsc -r300 EXnmr_b.eps
    figure(4), print -depsc -r300 EXnmr_e.eps
    figure(5), print -depsc -r300 EXnmr_c.eps
    figure(6), print -depsc -r300 EXnmr_f.eps
    % save as png
    figure(1), print -dpng -r300 EXnmr_a
    figure(2), print -dpng -r300 EXnmr_d
    figure(3), print -dpng -r300 EXnmr_b
    figure(4), print -dpng -r300 EXnmr_e
    figure(5), print -dpng -r300 EXnmr_c
    figure(6), print -dpng -r300 EXnmr_f
end
cd(oldcd)

% Uncomment the following return statement if you wish to save the
% displayed figures as MATLAB figures

% return

oldcd = cd;
if strcmp(dispres, 'subplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    figure(1), saveas('EXnmr.fig')
elseif strcmp(dispres, 'manyplots')
    try
        cd('Results')
    catch
        mkdir('Results')
        cd('Results')
    end
    saveas(figure(1), 'EXnmr_a.fig')
    saveas(figure(2), 'EXnmr_b.fig')
    saveas(figure(3), 'EXnmr_c.fig')
    saveas(figure(4), 'EXnmr_d.fig')
    saveas(figure(5), 'EXnmr_e.fig')
    saveas(figure(6), 'EXnmr_f.fig')
end
cd(oldcd)

repo: IRtools


text:
function ncc_set_test ( )

%*****************************************************************************80
%
%% NCC_SET_TEST tests NCC_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NCC_SET_TEST\n' );
  fprintf ( 1, '  NCC_SET sets up a Newton-Cotes Closed quadrature rule;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index       X             W\n' );
  fprintf ( 1, '\n' );

  for n = 1 : 10

    [ x, w ] = ncc_set ( n );

    fprintf ( 1, '\n' );

    for i = 1 : n
      fprintf ( 1, '  %2d  %12g  %12g\n', i, x(i), w(i) );
    end

  end

  return
end


text:
%DEMO_SPATIAL2  Demonstration for a disease mapping problem with Gaussian
%               process prior and negative binomial observation
%               model
%
%  Description
%    The disease mapping problem consist of a data with number of
%    death cases, Y, and background population, N, appointed to
%    co-ordinates, X. The goal is to find a relative risk surface,
%    which describes if the number of death cases in certain areas
%    is lower or higher than expected. The data is simulated.
%
%    The model is constructed as follows:
%
%    The number of death cases Y_i in area i is assumed to satisfy
%
%         Y_i ~ Neg-Bin(Y_i| d, E_i * r_i)
%
%    where E_i is the expected number of deaths (see Vanhatalo and
%    Vehtari (2007), how E_i is evaluated) at area i, r_i is the
%    relative risk and d is the dispersion parameter coverning the
%    variance.
%
%    We place a zero mean Gaussian process prior for log(R), R =
%    [r_1, r_2,...,r_n], which implies that at the observed input
%    locations latent values, f, have prior
%
%         f = log(R) ~ N(0, K),
%
%    where K is the covariance matrix, whose elements are given as
%    K_ij = k(x_i, x_j | th). The function k(x_i, x_j | th) is
%    covariance function and th its parameters. We place a prior
%    for parameters, p(th).
%
%    The inference is conducted first with Laplace approximation
%    and then with EP. We use compactly supported covariance
%    function which leads to sparse covariance matrix.
%
%  See also
%    DEMO_REGRESSION1, DEMO_CLASSIFIC1, DEMO_SPATIAL1
%
% Copyright (c) 2008-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

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

% =====================================
% Laplace approximation
% =====================================

% load the data
S = which('demo_spatial2');
data = load(strrep(S,'demo_spatial2.m','demodata/spatial2.txt'));

x = data(:,1:2);
ye = data(:,3);
y = data(:,4);

% Now we have loaded the following parameters
% x = co-ordinates
% y = number of deaths
% ye = the expexted number of deaths

fprintf(['GP with negative-binomial observation model, Laplace\n' ...
         'integration over the latent values and MAP estimate\n' ...
         'for the parameters\n']);

% Create the covariance functions
pl = prior_t();
pm = prior_sqrtt('s2', 0.3);
if ~exist('ldlchol')
  warning('GPstuff:SuiteSparseMissing',...
    ['SuiteSparse is not properly installed. \n' ...
    'Using gpcf_sexp (non-compact support) instead of gpcf_ppcs2 (compact support)']);
  gpcf1 = gpcf_sexp('lengthScale', 5, 'magnSigma2', 0.05);
  gpcf1 = gpcf_sexp(gpcf1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
else
  gpcf1 = gpcf_ppcs2('nin', 2, 'lengthScale', 5, 'magnSigma2', 0.05);
  gpcf1 = gpcf_ppcs2(gpcf1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
end

% Create the likelihood structure
lik = lik_negbin();

% Create the GP structure
gp = gp_set('lik', lik, 'cf', gpcf1, 'jitterSigma2', 1e-4);

% Set the approximate inference method to Laplace
gp = gp_set(gp, 'latent_method', 'Laplace');

% Set the options for the scaled conjugate optimization
opt=optimset('TolFun',1e-2,'TolX',1e-2,'Display','iter');
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'z',ye,'opt',opt);

% Visualize sparsity pattern
figure
C = gp_trcov(gp,x);
fprintf('Proportion of non-zeros is %.4f\n',nnz(C) / prod(size(C)))
p = amd(C);
spy(C(p,p))

% Make prediction to the data points
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye);

% Define help parameters for plotting
xii=sub2ind([120 70],x(:,2),x(:,1));
[X1,X2]=meshgrid(1:70,1:120);

% Plot the figures
figure
G=repmat(NaN,size(X1));
G(xii)=exp(Ef);
pcolor(X1,X2,G),shading flat
colormap(mapcolor(G)),colorbar
set(gca, 'Clim', [0.2    1.5])
axis equal
axis([0 70 0 120])
title('Posterior median of the relative risk (Laplace)')

figure
G=repmat(NaN,size(X1));
G(xii)=(exp(Varf) - 1).*exp(2*Ef+Varf);
pcolor(X1,X2,G),shading flat
colormap(mapcolor(G)),colorbar
%set(gca, 'Clim', [0.005    0.03])
axis equal
axis([0 70 0 120])
title('Posterior variance of the relative risk (Laplace)')

% =====================================
% EP approximation
% =====================================
fprintf(['GP with negative-binomial observation model, EP\n' ...
         'integration over the latent values and MAP estimate\n' ...
         'for the parameters\n']);

% Set the approximate inference method to EP
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the scaled conjugate optimization
opt=optimset('TolFun',1e-2,'TolX',1e-2,'Display','iter');
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'z',ye,'opt',opt);

% Visualize sparsity pattern
figure
C = gp_trcov(gp,x);
fprintf('Proportion of non-zeros is %.4f\n',nnz(C) / prod(size(C)))
p = amd(C);
spy(C(p,p))

% make prediction to the data points
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye);

% Define help parameters for plotting
xii=sub2ind([120 70],x(:,2),x(:,1));
[X1,X2]=meshgrid(1:70,1:120);

% Plot the figures
figure
G=repmat(NaN,size(X1));
G(xii)=exp(Ef);
pcolor(X1,X2,G),shading flat
colormap(mapcolor(G)),colorbar
set(gca, 'Clim', [0.2    1.5])
axis equal
axis([0 70 0 120])
title('Posterior median of the relative risk (EP)')

figure
G=repmat(NaN,size(X1));
G(xii)=(exp(Varf) - 1).*exp(2*Ef+Varf);
pcolor(X1,X2,G),shading flat
colormap(mapcolor(G)),colorbar
%set(gca, 'Clim', [0.005    0.03])
axis equal
axis([0 70 0 120])
title('Posterior variance of the relative risk (EP)')

author: gpstuff-dev

repo: gpstuff

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


text:
function mi = BF_MutualInformation(v1,v2,r1,r2,numBins)
% BF_MutualInformation  Mutual information between two data vectors using bin counting.
%
% Mutual information computed using a histogram-based, bin-counting method.
%
%---INPUTS:
% v1, the first input vector
% v2, the second input vector
% r1, the bin-partitioning method for the first input vector, v1
% r2, the bin-partitioning method for the second input vector, v2
% numBins, the number of bins to partition each vector into.
%
% NB: r1 and r2 can also be two-component vectors, that specify a custom range
%     for binning
%
%---OUTPUT:
% mi, the mutual information computed between v1 and v2

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check inputs and set defaults:
% ------------------------------------------------------------------------------
% by default, take a range equal to the range of the vectors
if nargin < 3 || isempty(r1)
    r1 = 'range';
end
if nargin < 4 || isempty(r2)
    r2 = 'range';
end
% use 10 bins for the histograms (=100 bins in 2d)
if nargin < 5 || isempty(numBins)
    numBins = 10;
end
% ------------------------------------------------------------------------------

if length(v1) ~= length(v2)
    error('input vectors must be the same length');
end

N = length(v1); % length of vectors

% Make sure both column vectors
if size(v1,2) > size(v1,1), v1 = v1'; end
if size(v2,2) > size(v2,1), v2 = v2'; end

% Create histograms:
% (i) in v1
edgesi = SUB_GiveMeEdges(r1,v1,numBins);
[ni, bini] = histc(v1, edgesi);

% (ii) in v2
edgesj = SUB_GiveMeEdges(r2,v2,numBins);
[nj, binj] = histc(v2, edgesj);

% ------------------------------------------------------------------------------
%% Create a joint histogram:
% ------------------------------------------------------------------------------

% We have the edges in each dimension: edgesi, and edgesj
histxy = zeros(numBins);
for i2 = 1:numBins
    for j2 = 1:numBins
        histxy(i2,j2) = sum(bini==i2 & binj==j2);
    end
end

% Normalize counts to probabilities
p_i = ni(1:numBins)/N;
p_j = nj(1:numBins)/N;
p_ij = histxy/N;
p_ixp_j = p_i*p_j';
summe = (p_ixp_j > 0 & p_ij > 0);

% Do a matrix-sum mutual information calculation:
if any(summe(:) == 1)
    mi = sum(p_ij(summe).*log(p_ij(summe)./p_ixp_j(summe)));
else
    fprintf(1,['The histograms aren''t catching any points?? ' ...
            'Perhaps due to an inappropriate custom range for binning the data...\n']);
    mi = NaN; return
end

% ------------------------------------------------------------------------------
% ------------------------------------------------------------------------------
function edges = SUB_GiveMeEdges(r,v,nbins)
    EE = 1E-6; % this small addition gets lost in the last bin
    if strcmp(r,'range')
        edges = linspace(min(v),max(v)+EE,nbins+1);

    elseif strcmp(r,'quantile') % bin edges based on quantiles
        edges = quantile(v,linspace(0,1,nbins+1));
%             edges(1) = edges(1) - 0.1;
        edges(end) = edges(end) + EE;
%             edges = sort(unique(edges)); % in case you have many repeated values -- will bias MI calculation
    elseif length(r)==2 % a two-component vector
        edges = linspace(r(1),r(2)+EE,nbins+1);
    else
        error('Unknown partitioning method ''%s''',r);
    end
end
% ------------------------------------------------------------------------------
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function y = sum_square_pos( x, dim )

%SUM_SQUARE_POS   Sum of squares of the positive parts.
%   For vectors, SUM_SQUARE(X) is the sum of the squares of the positive
%   parts of X; i.e., SUM( MAX(X,0)^2 ). X must be real.
%
%   For matrices, SUM_SQUARE_POS(X) is a row vector containing the
%   application of SUM_SQUARE_POS to each column. For N-D arrays, the
%   SUM_SQUARE_POS operation is applied to the first non-singleton
%   dimension of X.
%
%   SUM_SQUARE_POS(X,DIM) takes the sum along the dimension DIM of X.
%
%   Disciplined convex programming information:
%       SUM_SQUARE_POS(X,...) is convex and nondecreasing in X. Thus, when
%       used in CVX expressions, X must be convex (or affine). DIM must
%       always be constant.

error( nargchk( 1, 2, nargin ) );
if nargin == 2,
    y = sum( square_pos( x ), dim );
else
    y = sum( square_pos( x ) );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
% Selams, getSPLabels returns labels of Superpixels as well as number of
% total pixels in the image. It takes four parameters. One is image_input
% second one is Superpixel count limit. Third one is type. Use 1 for SLIC
% algorithm and other numbers for SLICO algorithm. BEWARE, if you are using
% type 1, which is SLIC, you have to set fourth parameter, compactness too.
% For more information, check SLIC and SLICO algorithms. If you are using
% type 2 or something else, just put random value to compactness, preferably
% ~ symbol for ignoring.
function [labels,numLabels] = getSPLabels(image_input,SPLimit,type,compactness)
% First, lets see if we have that files
conf;

% If 1 == SLIC, if 2 == SLICO
if type == 1
    [labels, numLabels] = slicmex(image_input,SPLimit,compactness);
else
    [labels, numLabels] = slicomex(image_input,SPLimit);
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm

path: github-repos/MATLAB/zhangqianqianQQ-MachineVisionAlgorithm/MachineVisionAlgorithm-683338f6c3b1aab9fa2b80026915fe936aebf0ee/识别算法/Object Recognition based on super pixel/getSPLabels.m


text:
function varargout = dsxy2figxy(varargin)
% dsxy2figxy -- Transform point or position from axis to figure coords
% Transforms [axx axy] or [xypos] from axes hAx (data) coords into coords
% wrt GCF for placing annotation objects that use figure coords into data
% space. The annotation objects this can be used for are
%    arrow, doublearrow, textarrow
%    ellipses (coordinates must be transformed to [x, y, width, height])
% Note that line, text, and rectangle anno objects already are placed
% on a plot using axes coordinates and must be located within an axes.
% Usage: Compute a position and apply to an annotation, e.g.,
%   [axx axy] = ginput(2);
%   [figx figy] = getaxannopos(gca, axx, axy);
%   har = annotation('textarrow',figx,figy);
%   set(har,'String',['(' num2str(axx(2)) ',' num2str(axy(2)) ')'])

%% Obtain arguments (only limited argument checking is performed).
% Determine if axes handle is specified
if length(varargin{1})== 1 && ishandle(varargin{1}) && ...
  strcmp(get(varargin{1},'type'),'axes')
 hAx = varargin{1};
 varargin = varargin(2:end);
else
 hAx = gca;
end;
% Parse either a position vector or two 2-D point tuples
if length(varargin)==1 % Must be a 4-element POS vector
 pos = varargin{1};
else
 [x,y] = deal(varargin{:});  % Two tuples (start & end points)
end
%% Get limits
axun = get(hAx,'Units');
set(hAx,'Units','normalized');  % Need normaized units to do the xform
axpos = get(hAx,'Position');
axlim = axis(hAx);              % Get the axis limits [xlim ylim (zlim)]
axwidth = diff(axlim(1:2));
axheight = diff(axlim(3:4));
yisreversed = strcmpi(get(hAx,'ydir'),'reverse');
xisreversed = strcmpi(get(hAx,'xdir'),'reverse');

%% Transform data from figure space to data space
if exist('x','var')     % Transform a and return pair of points

  xout = (x-axlim(1))/axwidth;
  yout = (y-axlim(3))/axheight;
  if xisreversed,
    xout = 1 - xout;
  end
  if yisreversed
    yout = 1 - yout;
  end

  varargout{1} = xout*axpos(3) + axpos(1);
  varargout{2} = yout*axpos(4) + axpos(2);
else                    % Transform and return a position rectangle
  xout = (pos(1)-axlim(1))/axwidth;
  yout = (pos(2)-axlim(3))/axheight;
  if xisreversed,
    xout = 1 - xout;
  end
  if yisreversed,
    yout = 1 - yout;
  end
  pos(1) = xout*axpos(3) + axpos(1);
  pos(2) = yout*axpos(4) + axpos(2);
  pos(3) = pos(3)*axpos(3)/axwidth;
  pos(4) = pos(4)*axpos(4)/axheight;
  varargout{1} = pos;
end
%% Restore axes units
set(hAx,'Units',axun)

author: kristinbranson

repo: JAABA


text:
function t = logtform(rmin, rmax, nr, nw)
% LOGTFORM makes a log-polar transform structure for imtransform
%     T = LOGTFORM(RMIN, RMAX, NR, NW) returns the transform structure for
%     a system with minimum ring radius RMIN, maximum ring radius RMAX, NR
%     rings and NR wedges. The empty matrix may be given for any one (but
%     only one) of these in which case the circular-samples condition
%     RMAX=RMIN*exp(2*pi*(NR-1)/NW) will be applied (with an adjustment
%     to RMIN if necessary to make NR and NW integers.
%
% See also LOGSAMPLE, LOGSAMPBACK

% Copyright David Young 2010

[rmin, rmax, nr, nw, k] = complete_args(rmin, rmax, nr, nw);
tdata = struct('rmin', rmin, 'rmax', rmax, 'nr', nr, 'nw', nw, 'k', k);
t = maketform('custom', 2, 2, @contorth, @rthtocon, tdata);
end

function x = contorth(u, t)
% Conventional to log-polar. See maketform.
td = t.tdata;
[th, p] = cart2pol(u(:,1), u(:, 2));
p(~p) = td.rmin/2;            % Omit centre point
x = [td.k * log(p/td.rmin),  td.nw*mod(th/(2*pi), 1)];
end

function u = rthtocon(x, t)
% Log-polar to conventional. See maketform.
td = t.tdata;
p = td.rmin * exp(x(:, 1)/td.k);
th = (2*pi/td.nw) * x(:, 2);
[x, y] = pol2cart(th, p);
u = [x, y];
end

function [rmin, rmax, nr, nw, k] = complete_args(rmin, rmax, nr, nw)
% Circular pixels condition
if isempty(rmin)
    k = nw / (2*pi);
    rmin = rmax * exp((1-nr)/k);
elseif isempty(rmax)
    k = nw / (2*pi);
    rmax = rmin * exp((nr-1)/k);
elseif isempty(nw)
    k = (nr-1) / log(rmax/rmin);
    nw = round(2 * pi * k);
    k = nw / (2*pi);
    rmin = rmax * exp((1-nr)/k);
elseif isempty(nr)
    k = nw / (2*pi);
    nr = round(k * log(rmax/rmin) + 1);
    rmin = rmax * exp((1-nr)/k);
else
    k = (nr-1) / log(rmax/rmin);
end
end

author: Sable

repo: mcbench-benchmarks


text:
%compute dycentralhead_mm

function [data,units]=compute_dycentralhead_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dycentralhead_mm=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dycentralhead_mm{1,i}=(trx(larva).ycentralhead_mm(2:end)-trx(larva).ycentralhead_mm(1:end-1))./trx(larva).dt;
end

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


text:
function view = setCothresh(view, value)
%
% view = setCothresh(view, value)
%
% Sets cothreshSlider value
% ras 01/07 -- update for hidden views.
if ~exist('view', 'var') | isempty(view), view = getCurView; end

if ~checkfields(view, 'ui', 'cothresh') | ~isstruct(view.ui.cothresh)
    view.settings.cothresh = value;

else    % slider structure
    setSlider(view,view.ui.cothresh,value);

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft