text:
%maincheck.m
%check whether an earthquake is followed by an equal or bigger sized
%eq in a certain distance
%Last change    11/95 Alexander Allmann
%

%TODO Delete, it was error-full -CGR
global newcat err derr;
newcat=a;
n=0;               %counter
n1=0;
mcut1=5;           %default for magnitudes of interest
mcut2=5;           %default for magnitudes of eqs following first event
dist=30;           %distance in which programs looks for following events
tcut1=30;          %time intervall in which program searches for other events
err=2;derr=2;
newcat=newcat.Magnitude>mcut1;
backnewcat=newcat;
eqtime=clustime();      %onset time of earthquakes

for i=1:newcat.Count

    tmp=find((eqtime-eqtime(i))>0 & (eqtime - eqtime(i))<tcut1);
    if ~isempty(tmp)
        tmp2=newcat.subset(tmp);
        tmp3=find(tmp2(:,6)>=newcat(i,6));
        if ~isempty(tmp3)
            newcat=[newcat.subset(i);tmp2(tmp3,:)];
                [dist1, dist2] = distance(1,1,2:newcat.Count);
                tmp4=dist1<dist;
                if ~isempty(tmp4)
                n=n+1;
                if length(tmp4)>1
                    n1=n1+1;                %more than one bigger event follows in sequence
                end
                end
                newcat=backnewcat;
        end
    end
end
n
n1

author: CelsoReyes


text:
% batch run full-clustering on all fish, sweep param with 2-fold CV

data_masterdir = GetCurrentDataDir();

% range_fish = [5,6,7];
% M_ClusGroup = [2,2,2,2];
% M_Cluster = [1,1,1,1];
range_fish = 8:9;
% M_ClusGroup = 2;
% M_Cluster = 3;
M_stim = 1;
% M_fish_set = [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2];

%%
M_param = 5:5:100;%0.3:0.1:0.8;

Param_nClus = zeros(length(M_param),length(range_fish));
Param_CVscore = zeros(length(M_param),length(range_fish));
Param_CVscore_raw = cell(length(M_param),length(range_fish));

%     thres_split = M_param(k_param);
%     setappdata(hfig,'thres_split',thres_split);

for k_fish = 1:length(range_fish),
    i_fish = range_fish(k_fish);
    disp(i_fish);
    LoadFullFish(hfig,i_fish,0);
    absIX = getappdata(hfig,'absIX');

    %% partitions for CV
    timelists = getappdata(hfig,'timelists');
    timelists_names = getappdata(hfig,'timelists_names');
    periods = getappdata(hfig,'periods');
    if length(periods)>1,
        timelistsCV = cell(length(M_stim),2);

        k_stim = 1;
%         for k_stim = 1:length(M_stim),
        i_stim = M_stim(k_stim);
        TL = timelists{i_stim};
        period = periods(i_stim);
        nrep = size(TL,2)/periods(i_stim); % integer
        n = floor(nrep/2);
        timelistsCV{k_stim,1} = TL(1):TL(n*period);
        timelistsCV{k_stim,2} = TL(1+n*period):TL(2*n*period);
%           end
    end

    for k_param = 1:length(M_param),
        numK2 = M_param(k_param);

        Score = zeros(1,2);%(length(M_stim),2);
        %%
        NumClus = zeros(1,2);
        CIX = cell(1,2);
        GIX = cell(1,2);
        for k = 1:2, % CV halves
            i_ClusGroup = 2;
            i_Cluster = k+4;
            [cIX,gIX] = LoadCluster_Direct(i_fish,i_ClusGroup,i_Cluster,absIX);

            tIX = timelistsCV{k_stim,k};
            M_0 = GetTimeIndexedData_Default_Direct(hfig,[],tIX,'isAllCells');

            isWkmeans = 0;
            [cIX,gIX] = AutoClustering(cIX,gIX,absIX,i_fish,M_0,isWkmeans,numK2);

            NumClus(k) = length(unique(gIX));
            CIX{k} = cIX;
            GIX{k} = gIX;
        end
        % plot cell-matching figure
        Score(1) = HungarianCV(NumClus(1),NumClus(2),CIX{1},CIX{2},GIX{1},GIX{2});% true,timelists_names{i_stim});
        Score(2) = HungarianCV(NumClus(2),NumClus(1),CIX{2},CIX{1},GIX{2},GIX{1});% true,timelists_names{i_stim});
        Param_CVscore(k_param,k_fish) = mean(Score);
        Param_CVscore_raw{k_param,k_fish} = Score;

        nClus1 = length(unique(GIX{1}));
        nClus2 = length(unique(GIX{2}));
        Param_nClus(k_param,k_fish) = mean([nClus1,nClus2]);

    end
end

%%
figure;
subplot(2,1,1)
plot(M_param*20,Param_nClus)
% legend('Fish8','Fish9');
xlabel('total k for kmeans')
ylabel('# of auto-clusters')
subplot(2,1,2)
plot(M_param*20,Param_CVscore)
ylim([0,1])
legend('Fish8','Fish9');
xlabel('total k for kmeans')
ylabel('CV (overlapping cell %)')

author: xiuyechen

repo: FishExplorer


text:
function relm_PaintCumPlot(rRelmTest, sXLabel, hAxes)
    % Creates the result plots of one of the RELM tests
    %
    % relm_PaintCumPlot(rRelmTest, sXLabel, hAxes)
    %
    %
    % Input parameters:
    %   rRelmTest		Record of results from one of the RELM tests
    %   sXLabel		Label of X-axis (optional)
    %   hAxes		Handle of existing axes. If not specified, a figure is created
    %
    % Copyright (C) 2002-2006 by Danijel Schorlemmer
    %
    % This program is free software; you can redistribute it and/or modify
    % it under the terms of the GNU General Public License as published by
    % the Free Software Foundation; either version 2 of the License, or
    % (at your option) any later version.
    %
    % This program is distributed in the hope that it will be useful,
    % but WITHOUT ANY WARRANTY; without even the implied warranty of
    % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    % GNU General Public License for more details.
    %
    % You should have received a copy of the GNU General Public License
    % along with this program; if not, write to the
    % Free Software Foundation, Inc.,
    % 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

    if ~exist('sXLabel')
        sXLabel = '';
    end
    if ~exist('hAxes')
        figure_w_normalized_uicontrolunits('Name', 'Result plot', 'NumberTitle', 'off');
        hAxes = newplot;
    end

    % Activate the given or newly created axes
    axes(hAxes);

    % Plot the CDF and the line of observed data
    vIndex = [1:rRelmTest.nNumberSimulation]/rRelmTest.nNumberSimulation;
    vObservedX = [rRelmTest.fObservedData, rRelmTest.fObservedData];
    vObservedY = [0,1];
    plot(rRelmTest.vSimValues_H, vIndex, 'g', rRelmTest.vSimValues_N, vIndex, 'r', vObservedX, vObservedY, 'k', 'LineWidth', 1);

    % Add the patches
    set(hAxes, 'NextPlot', 'add');
    vPatch1Y = [1 1 0.975 0.975];
    vPatch2Y = [0.025 0.025 0 0];
    vXLim = xlim;
    vPatch1X = [vXLim(1) rRelmTest.fObservedData rRelmTest.fObservedData vXLim(1)];
    vPatch2X = [rRelmTest.fObservedData vXLim(2) vXLim(2) rRelmTest.fObservedData];
    patch(vPatch1X, vPatch1Y, [0.8 0.8 0.8], 'facealpha', 0.7);
    patch(vPatch2X, vPatch2Y, [0.8 0.8 0.8], 'facealpha', 0.7);

    % Add information
    title(['\alpha = ' num2str(rRelmTest.fAlpha) ', \beta = ' num2str(rRelmTest.fBeta)]);
    ylabel('Fraction of cases');
    xlabel(sXLabel);
end


text:
classdef testTopOptDensityPrinting < testTopOptPrinting

    properties (Access = protected)
        fileOutputName = 'testTopOptDensityPrinting';
        printMode = 'DesignVariable';
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Test the annotation procedure, the annotation data is backed up so
% feel free to experiment.
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

clear all;

addpath('..');

training_files = dir('cootes/*.bmp.mat');

for i=1:numel(training_files)
	load(sprintf('cootes/%s', training_files(i).name));

	app = imread(sprintf('cootes/%s', training_files(i).name(1:end-4)));

	% Map RGB colors to [0,1]
	appearances(:,:,:,i) = double(app) ./ 255;

	shapes(:,:,i) = xy2ij(annotations, size(app,1));
end

load('cootes/triangulation.mat');

test = 23;
template = 22;

AAM = build_model_2d(shapes, appearances, 'triangulation', triangulation);
backup_file = sprintf('cootes/%s.bak', training_files(test).name);
copyfile(sprintf('cootes/%s', training_files(test).name), backup_file);

annotate(sprintf('cootes/%s', training_files(test).name(1:end-4)), sprintf('cootes/%s', training_files(template).name), 'AAM', AAM);

copyfile(backup_file, sprintf('cootes/%s', training_files(test).name));
delete(backup_file);

author: Sable

repo: mcbench-benchmarks


text:
function [mssim, ssim_map, mcs, cs_map] = ssim_index_new(img1, img2, K, win, max_val)

if (nargin < 2 || nargin > 5)
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

if (size(img1) ~= size(img2))
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

[M N] = size(img1);

if (nargin == 2)
   if ((M < 11) || (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   win = fspecial('gaussian', 11, 1.5);	%
   K(1) = 0.01;										% default settings
   K(2) = 0.03;										%
end

if (nargin == 3)
   if ((M < 11) || (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   win = fspecial('gaussian', 11, 1.5);
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 4 || nargin == 5)
   [H W] = size(win);
   if ((H*W) < 4 || (H > M) || (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

C1 = (K(1)*max_val)^2;
C2 = (K(2)*max_val)^2;
win = win/sum(sum(win));

mu1   = filter2(win, img1, 'valid');
mu2   = filter2(win, img2, 'valid');
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(win, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(win, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(win, img1.*img2, 'valid') - mu1_mu2;

if (C1 > 0 & C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
   cs_map = (2*sigma12 + C2)./(sigma1_sq + sigma2_sq + C2);
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
	denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;

   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);

   cs_map = ones(size(mu1));
   index = denominator2 > 0;
   cs_map(index) = numerator2(index)./denominator2(index);
end

mssim = mean2(ssim_map);
mcs = mean2(cs_map);

return

author: sooyekim

repo: Deep-SR-ITM


text:
function [out] = RFT_main(X,options,verbose)
% applies RFT to an input 1D-field X sampled on a regular lattice
% function [out] = RFT_main(X,options,verbose)
% RFT (Random Field Theory) allows one to derive the FWE-corrected p-values
% on some topological features of the field (e.g., local peaks and
% upcrossing thresholds), under the assumption of non-independance
% between neighbouring positions on the lattice.
% IN:
%   - X: Lx1 vector of sampled RF
%   - options: structure containing the set-inducing thresholds, i.e.:
%       .u: height threshold (cluster-level inference, default=2.33)
%       .k: extent threshold (set-level inference, default=6)
%       .R: sampled RF of residuals (for smoothness estimation). If empty,
%       then the smoothness is estimated directly on X.
%       .type: type of RF. Can be set to 'norm' (normal, default), 't'
%       (Student) or 'F' (Fisher).
%       .dof: degrees of freedom (only relevant for 't' or 'F' fields).
%   - verbose: a verbose flag
% OUT:
%   - out: a structure containing the following fields:
%       .peaks: summary stats of peak-level inference, i.e.:
%           .ind: npx1 vector of local maxima indices
%           .val: npx1 vector of local maxima values
%           .prft: npx1 vector of local maxima corrected p-values
%           .punc: npx1 vector of local maxima uncorrected p-values
%       where np is the number of local maxima on the field.
%       .clusters: summary stats of cluster-level inference, i.e.:
%           .ind: ncx1 cell-array of sets of indices
%           .k: ncx1 vector of clusters' spatial extent
%           .prft: ncx1 vector of clusters' corrected p-values
%       where nc is the number of upcrossing clusters on the field.
%       .set: summary stats of set-level inference, i.e.:
%           .c: number of upcrossing threshold with a size bigger than k
%           .prft: set-level p-value
%       .xc: the height threshold that corresponds to the nominal FWER
%       .Em: expected number of clusters (under H0)
%       .En: expected number of voxels per cluster (under H0)
%       .fwhm: the estimated smoothness of the field

try,options;catch,options=[];end
try,verbose;catch,verbose=0;end

OUTSTR = cell(0);
OUTSTR{1} = ['Date: ',datestr(clock)];
if verbose
    disp(' ')
    disp('-- 1D-RFT analysis --')
    disp(OUTSTR{1})
end

if isempty(options)
    options.FWER = 0.05;
    options.u = VBA_spm_invNcdf(0.99,0,1);
    options.k = 6;
    options.R = [];
    options.type = 'norm';
    options.dof = NaN; % only relevant for 't' or 'F' fields
    options.mask = [];
    if verbose
        disp(['Using default set-inducing thresholds (X>',num2str(options.u,'%3.2f'),', k>',num2str(options.k),').'])
    end
else
    if ~isfield(options,'FWER')
        options.FWER = 0.05;
    end
    if ~isfield(options,'k')
        options.k = 6;
        if verbose
            disp(['Using default set-inducing threshold (k>',num2str(options.k),').'])
        end
    end
    if ~isfield(options,'R')
        options.R = [];
        if verbose
            disp(['Estimating smoothness on statistical map.'])
        end
    end
    if ~isfield(options,'type')
        options.type = 'norm';
    end
    if ~isfield(options,'dof')
        if ~isequal(options.type,'norm')
            disp('RFT: error: must provide valid degrees of freedom!')
            out = [];
        else
            options.dof = NaN;
        end
    end
    if ~isfield(options,'u')
        switch options.type
            case 'norm'
                options.u = VBA_spm_invNcdf(0.99,0,1);
            case 't'
                options.u = VBA_spm_invTcdf(0.99,options.dof);
            case 'F'
                options.u = VBA_spm_invFcdf(0.99,options.dof(1),options.dof(2));
        end
        if verbose
            disp(['Using default cluster-inducing threshold (X>',num2str(options.u,'%3.2f'),').'])
        end
    end
    if ~isfield(options,'mask')
        options.mask = [];
    end
end

if ~isempty(options.mask)
    if length(options.mask)~=length(X) || ~isequal(VBA_vec(unique(options.mask)),[0;1])
        disp(['RFT-1D: error: invalid mask provided!'])
        out = [];
        return
    end
end

out.options = options;
out.verbose = verbose;

X = VBA_vec(X);
if isempty(options.mask)
    L = length(X);
else
    L = length(find(options.mask==1));
end

% estimate smoothness
if ~isempty(options.R)
    out.fwhm = RFT_smoothness(options.R);
else
    out.fwhm = RFT_smoothness(X);
end
OUTSTR{2} = ['Search volume: L=',num2str(L)];
OUTSTR{3} = ['Estimated smoothness: FWHM=',num2str(out.fwhm,'%3.1f')];
OUTSTR{4} = ['Number of resels: R=',num2str(L./out.fwhm,'%3.1f')];
if verbose
    disp(OUTSTR{2})
    disp(OUTSTR{3})
    disp(OUTSTR{4})
end

% estimate critical height threshold (peak-level)
gridu = 0:1e-3:10;
pgrid = RFT_Pval(gridu',0,1,out.fwhm,L,options.type,options.dof);
d = abs(options.FWER-pgrid);
out.xc = gridu(find(d==min(d)));
OUTSTR{5} = ['Critical height threshold [peak-level]: X>',num2str(out.xc,'%3.2f'),' (test size: FWER=',num2str(round(options.FWER*100)),'%)'];
if verbose
    disp(OUTSTR{5})
end

% peak- level inference
peaks.ind = RFT_localmax(X);
peaks.val = X(peaks.ind);
peaks.prft = RFT_Pval(peaks.val,0,1,out.fwhm,L,options.type,options.dof);
switch options.type
    case 'norm'
        peaks.punc = 1-VBA_spm_Ncdf(peaks.val,0,1);
    case 't'
        peaks.punc = 1-VBA_spm_Tcdf(peaks.val,options.dof);
    case 'F'
        peaks.punc = 1-VBA_spm_Fcdf(peaks.val,options.dof(1),options.dof(2));
end

% cluster-level inference
[clusters.ind,clusters.imax] = RFT_clusters(X,options.u,0);
nc = length(clusters.ind);
clusters.k = [];
clusters.prft = [];
for i=1:nc
    clusters.k(i) = length(clusters.ind{i});
    clusters.prft(i) = RFT_Pval(options.u,clusters.k(i),1,out.fwhm,L,options.type,options.dof);
end

% set-level inference
if nc>0
    bigC = find(clusters.k>=options.k);
else
    bigC = [];
end
set.c = length(bigC);
set.prft = RFT_Pval(options.u,options.k,set.c,out.fwhm,L,options.type,options.dof);

% E[number of voxel per cluster] and E[number of clusters]
out.Em = RFT_expectedTopo(options.u,L,out.fwhm,1,options.type,options.dof);
switch options.type
    case 'norm'
        P0 = 1-VBA_spm_Ncdf(options.u,0,1);
    case 't'
        P0 = 1-VBA_spm_Tcdf(options.u,options.dof);
    case 'F'
        P0 = 1-VBA_spm_Fcdf(options.u,options.dof(1),options.dof(2));
end
out.En = L.*P0./out.Em;

OUTSTR{6} = ['Expected voxels per cluster [cluster-level]: E[k|H0]=',num2str(out.En,'%3.1f')];
OUTSTR{7} = ['Expected number of clusters [set-level]: E[c|H0]=',num2str(out.Em,'%3.1f')];
if verbose
    disp(OUTSTR{6})
    disp(OUTSTR{7})
end

OUTSTR{8} = ['Number of local peaks =',num2str(length(peaks.prft)),];
OUTSTR{9} = ['Number of upcrossing clusters =',num2str(nc),' (X>',num2str(options.u,'%3.2f'),')'];
OUTSTR{10} = ['RFT [set-level]: p=',num2str(set.prft,'%3.3f'),' (c=',num2str(set.c),')'];
if verbose
    disp(OUTSTR{8})
    disp(OUTSTR{9})
    disp(OUTSTR{10})
end

% wrap-up results
out.peaks = peaks;
out.clusters = clusters;
out.set = set;
out.OUTSTR = OUTSTR;

% order peak-pval of clusters' maxima
inc = ismember(peaks.ind,clusters.imax);
[ps,is] = sort(peaks.prft(inc==1),'ascend');
STR.loc = cell(0);
STR.unc = cell(0);
STR.peak = cell(0);
STR.cluster = cell(0);
for i=1:nc % for all clusters
    % find local maxima that belong to each cluster
    ic = is(i);
    ip = find(ismember(peaks.ind,clusters.ind{ic})==1);
    % sort peak-pvalues in each cluster
    [pval,iop] = sort(peaks.prft(ip),'ascend');
    % aggregate info Re: local maxima per (ordered) cluster
    nt = length(STR.loc);
    for j=1:length(ip)
        STR.loc{nt+j} = num2str(peaks.ind(ip(iop(j))));
        STR.unc{nt+j} = num2str(peaks.punc(ip(iop(j))),'%3.3f');
        STR.peak{nt+j} = [num2str(peaks.prft(ip(iop(j))),'%3.3f'),'   (',num2str(peaks.val(ip(iop(j))),'%3.2f'),')'];
        if j==1
            STR.cluster{nt+j} = [num2str(clusters.prft(ic),'%3.3f'),'   (',num2str(clusters.k(ic)),')'];
        end
    end
end
STR.set = [num2str(set.prft,'%3.3f'),'   (',num2str(set.c),')'];
out.STR = STR;

% display results
if verbose
    [out] = RFT_ReDisplay(X,out);
end

author: MBB-team

repo: VBA-toolbox


text:
function [res,fval,it] = muller (f,Z0,itmax,ztol,ftol,option)
% MULLER find a zero of a real or complex function Y = F(Z)
%
% Syntax:
%
% RES = MULLER (F,Z0) find the zero of a complex or real function
% 'F' (either an anonymous function or .m function) using three initial
% guesses contained in the vector Z0. Muller takes the function F and
% evaluetes it at each initial point using feval. F doesn't need to be
% vectorized.
% The initial guesses can be real or complex numbers close to the zero,
% bracketing the zero is not necessary. Parameters ITMAX, ZTOL and
% FTOL are set by default to 1000, 1e-5 and 1e-5, respectively.
%
% RES = MULLER (F,Z0,ITMAX) the maximum number of iterations is set
% equal to ITMAX. ZTOL and FTOL are set by default with the values mentioned
% above.
%
% RES = MULLER (F,Z0,ITMAX,ZTOL) ZTOL is used as a stopping
% criterion. If the absolute difference between the values of Z found in
% the two latest iterations is less than ZTOL, the program is stopped. FTOL
% is set by default with the value mentioned above.
%
% RES = MULLER (F,Z0,ITMAX,ZTOL,FTOL) FTOL is used as a stopping
% criterion. If the value of the function F at the Z found in the last
% iteration is less than FTOL, the program is stopped.
%
% RES = MULLER (F,Z0,ITMAX,ZTOL,FTOL,'both') indicate that both
% criteria ZTOL and FTOL must be satisfied simultaneously. By default,
% MULLER stops if one of the two criteria is fulfilled.
%
% [RES,FVAL] = MULLER (F,Z0,...) return the value of the function
% F at the Z found in the last iteration.
%
% [RES,FVAL,IT] = MULLER (F,Z0,...) return the number of iterations
% used to find the zero.
%
% Example 1:
% myf = @(x) (x-1)^3;
%
% muller(myf,[0 0.1 0.2],[],[],[],'both')
% ans =
%    1.0000 + 0.0000i
%
% Example 2:
%
% [res,fval,it] = muller2('cosh',[0 0.1 0.2],[],[],[],'both')
%
% res =
%   0.0000 + 1.5708i
%
% fval =
%   5.5845e-012 + 3.0132e-012i
%
% it =
%   5
%
% Method taken from:
% Numerical Recipes: The art of scientific computing
% W.H. Press; B.P. Flannery; S.A. Teukolsky; W.T. Vetterling
% 1986
%
% Thanks to John D'Errico for his helpfull review.
%
% Written by Daniel H. Cortes
% MAE Department, West Virginia University
% March, 2008.
%

%=================================================
% Checking proper values of the input parameters
%=================================================

if nargin > 6
    error ('Too many arguments.')
elseif nargin < 2
    error('Too few arguments.')
end

if nargin < 6
    opt = 1;
elseif ischar(option) == 1
    if size(option,2) == 4
        if sum(option == 'both') == 4
            opt = 2;
        else
            error ('Option parameter must be *both*.')
        end
    else
        error ('Option parameter must be *both*.')
    end
else
    error ('Option parameter must be a character array (string).')
end

if nargin < 5
    ftol = 1e-5;
elseif isnumeric(ftol) ~= 1
    error ('FTOL must be a numeric argument.')
elseif isempty(ftol) == 1
    ftol = 1e-5;
elseif  size(ftol,1) ~= 1 || size(ftol,2) ~= 1
    error ('FTOL cannot be an array')
end

if nargin < 4
    ztol = 1e-5;
elseif isnumeric(ztol) ~= 1
    error ('ZTOL must be a numeric argument.')
elseif isempty(ztol) == 1
    ztol = 1e-5;
elseif  size(ztol,1) ~= 1 || size(ztol,2) ~= 1
    error ('ZTOL cannot be an array.')
end

if nargin < 3
    itmax = 1000;
elseif isnumeric(itmax) ~= 1
    error ('ITMAX must be a numeric argument.')
elseif isempty(itmax) == 1
    itmax = 1000;
elseif  size(itmax,1) ~= 1 || size(itmax,2) ~= 1
    error ('ITMAX cannot be an array.')
end

if isnumeric(Z0) ~= 1
    error ('Z0 must be a vector of three numeric arguments.')
elseif isempty(Z0) == 1 || length(Z0) ~= 3 || min(size(Z0)) ~= 1
    error ('Z0 must be a vector of length 3 of either complex or real arguments.')
end

if Z0(1)==Z0(2) || Z0(1)==Z0(3) || Z0(2)==Z0(3)
    error('The initial guesses must be different')
end

%=============================
% Begining of Muller's method
%=============================

z0 = Z0(1);
z1 = Z0(2);
z2 = Z0(3);

y0 = feval ( f, z0);
y1 = feval ( f, z1);
y2 = feval ( f, z2);

for it = 1:itmax
    q = (z2 - z1)/(z1 - z0);
    A = q*y2 - q*(1+q)*y1 + q^2*y0;
    B = (2*q + 1)*y2 - (1 + q)^2*y1 + q^2*y0;
    C = (1 + q)*y2;

    if ( A ~= 0 )

        disc = B^2 - 4*A*C;

        den1 = ( B + sqrt ( disc ) );
        den2 = ( B - sqrt ( disc ) );

        if ( abs ( den1 ) < abs ( den2 ) )
            z3 = z2 - (z2 - z1)*(2*C/den2);
        else
            z3 = z2 - (z2 - z1)*(2*C/den1);
        end

    elseif ( B ~= 0 )
        z3 = z2 - (z2 - z1)*(2*C/B);
    else
        warning('Muller Method failed to find a root. Last iteration result used as an output. Result may not be accurate')
        res = z2;
        fval = y2;
        return
    end

    y3 = feval ( f, z3);

    if opt == 1
        if ( abs (z3 - z2) < ztol || abs ( y3 ) < ftol )
            res = z3;
            fval = y3;
            return
        end
    else
        if ( abs (z3 - z2) < ztol && abs ( y3 ) < ftol )
            res = z3;
            fval = y3;
            return
        end
    end

    z0 = z1;
    z1 = z2;
    z2 = z3;

    y0 = y1;
    y1 = y2;
    y2 = y3;

end

res = z2;
fval = y2;
warning('Maximum number of iterations reached. Result may not be accurate')


text:
function [xd,xn,option] = denoise(x,h,type,option)
%    [xd,xn,option] = denoise(x,h,type,option);
%
%    DENOISE is a generic program for wavelet based denoising.
%    The program will denoise the signal x using the 2-band wavelet
%    system described by the filter h using either the traditional
%    discrete wavelet transform (DWT) or the linear shift invariant
%    discrete wavelet transform (also known as the undecimated DWT
%    (UDWT)).
%
%    Input:
%       x         : 1D or 2D signal to be denoised
%       h         : Scaling filter to be applied
%       type      : Type of transform (Default: type = 0)
%                   0 --> Discrete wavelet transform (DWT)
%                   1 --> Undecimated DWT (UDWT)
%       option    : Default settings is marked with '*':
%                   *type = 0 --> option = [0 3.0 0 0 0 0]
%                   type = 1 --> option = [0 3.6 0 1 0 0]
%       option(1) : Whether to threshold low-pass part
%                   0 --> Don't threshold low pass component
%                   1 --> Threshold low pass component
%       option(2) : Threshold multiplier, c. The threshold is
%                   computed as:
%                     thld = c*MAD(noise_estimate)).
%                   The default values are:
%                     c = 3.0 for the DWT based denoising
%                     c = 3.6 for the UDWT based denoising
%       option(3) : Type of variance estimator
%                   0 --> MAD (mean absolute deviation)
%                   1 --> STD (classical numerical std estimate)
%       option(4) : Type of thresholding
%                   0 --> Soft thresholding
%                   1 --> Hard thresholding
%       option(5) : Number of levels, L, in wavelet decomposition. By
%                   setting this to the default value '0' a maximal
%                   decomposition is used.
%       option(6) : Actual threshold to use (setting this to
%                   anything but 0 will mean that option(3)
%                   is ignored)
%
%    Output:
%       xd     : Estimate of noise free signal
%       xn     : The estimated noise signal (x-xd)
%       option : A vector of actual parameters used by the
%                program. The vector is configured the same way as
%                the input option vector with one added element
%                option(7) = type.
%
%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
%  Cut-and-paste the example you want to run to a new file
%  called ex.m, for example. Delete out the % at the beginning
%  of each line in ex.m (Can use search-and-replace in your editor
%  to replace it with a space). Type 'ex' in matlab and hit return.
%
%    Example 1:
%       h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N);
%       x = s + n/10;     % (approximately 10dB SNR)
%       figure;plot(x);hold on;plot(s,'r');
%
%       %Denoise x with the default method based on the DWT
%       [xd,xn,opt1] = denoise(x,h);
%       figure;plot(xd);hold on;plot(s,'r');
%
%       %Denoise x using the undecimated (LSI) wavelet transform
%       [yd,yn,opt2] = denoise(x,h,1);
%       figure;plot(yd);hold on;plot(s,'r');
%
% Example 2: (on an image)
%      h = daubcqf(6);  load lena;
%      noisyLena = lena + 25 * randn(size(lena));
%      figure; colormap(gray); imagesc(lena); title('Original Image');
%       figure; colormap(gray); imagesc(noisyLena); title('Noisy Image');
%       Denoise lena with the default method based on the DWT
%      [denoisedLena,xn,opt1] = denoise(noisyLena,h);
%      figure; colormap(gray); imagesc(denoisedLena); title('denoised Image');
%
%
%    See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt
%

%File Name: denoise.m
%Last Modification Date: 04/15/97	10:44:28
%Current Version: denoise.m	2.4
%File Creation Date: Mon Feb 20 08:33:15 1995
%Author: Jan Erik Odegard  <odegard@ece.rice.edu>
%
%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
%Created by Jan Erik Odegard, Department of ECE, Rice University.
%
%This software is distributed and licensed to you on a non-exclusive
%basis, free-of-charge. Redistribution and use in source and binary forms,
%with or without modification, are permitted provided that the following
%conditions are met:
%
%1. Redistribution of source code must retain the above copyright notice,
%   this list of conditions and the following disclaimer.
%2. Redistribution in binary form must reproduce the above copyright notice,
%   this list of conditions and the following disclaimer in the
%   documentation and/or other materials provided with the distribution.
%3. All advertising materials mentioning features or use of this software
%   must display the following acknowledgment: This product includes
%   software developed by Rice University, Houston, Texas and its contributors.
%4. Neither the name of the University nor the names of its contributors
%   may be used to endorse or promote products derived from this software
%   without specific prior written permission.
%
%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%For information on commercial licenses, contact Rice University's Office of
%Technology Transfer at techtran@rice.edu or (713) 348-6173
%
%Change History: Fixed output of function and an error in the computation
%                of the threshold for redundant denoising.
%                <Jan Erik Odegard> <Mon Jul 31, 1995>
%
%                This code is composed of several of our old codes for
%                wavelet based denoising. In an effort to make the mess
%                more manageable we decided to create on code that would
%                handle all the various wavelet based denoising methods.
%                However, only time will show (as we discover new and
%                improved forms of denoising) if we can succeed in our goals.
%                <Jan Erik Odegard> <Thu May 11, 1995>
%

if(nargin < 2)
  error('You need to provide at least 2 inputs: x and h');
end;
if(nargin < 3),
  type = 0;
  option = [];
elseif(nargin < 4)
  option = [];
end;
if(isempty(type)),
  type = 0;
end;
if(type == 0),
  default_opt = [0 3.0 0 0 0 0];
elseif(type == 1),
  default_opt = [0 3.6 0 1 0 0];
else,
  error(['Unknown denoising method',10,...
	  'If it is any good we need to have a serious talk :-)']);
end;
option = setopt(option,default_opt);
[mx,nx] = size(x);
dim = min(mx,nx);
if(dim == 1),
  n = max(mx,nx);
else,
  n = dim;
end;
if(option(5) == 0),
  L = floor(log2(n));
else
  L = option(5);
end;
if(type == 0), 			% Denoising by DWT
  xd = mdwt(x,h,L);
  if (option(6) == 0),
    tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx);
    if(option(3) == 0),
      thld = option(2)*median(abs(tmp(:)))/.67;
    elseif(option(3) == 1),
      thld = option(2)*std(tmp(:));
    else
      error('Unknown threshold estimator, Use either MAD or STD');
    end;
  else,
    thld = option(6);
  end;
  if(dim == 1)
    ix = 1:n/(2^L);
    ykeep = xd(ix);
  else
    ix = 1:mx/(2^L);
    jx = 1:nx/(2^L);
    ykeep = xd(ix,jx);
  end;
  if(option(4) == 0),
    xd = SoftTh(xd,thld);
  elseif(option(4) == 1),
    xd = HardTh(xd,thld);
  else,
    error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
  end;
  if (option(1) == 0),
    if(dim == 1),
      xd(ix) = ykeep;
    else,
      xd(ix,jx) = ykeep;
    end;
  end;
  xd = midwt(xd,h,L);
elseif(type == 1), 			% Denoising by UDWT
  [xl,xh] = mrdwt(x,h,L);
  if(dim == 1),
    c_offset = 1;
  else,
    c_offset = 2*nx + 1;
  end;
  if (option(6) == 0),
    tmp = xh(:,c_offset:c_offset+nx-1);
    if(option(3) == 0),
      thld = option(2)*median(abs(tmp(:)))/.67;
    elseif(option(3) == 1),
      thld = option(2)*std(tmp(:));
    else
      error('Unknown threshold estimator, Use either MAD or STD');
    end;
  else,
    thld = option(6);
  end;
  if(option(4) == 0),
    xh = SoftTh(xh,thld);
    if(option(1) == 1),
      xl = SoftTh(xl,thld);
    end;
  elseif(option(4) == 1),
    xh = HardTh(xh,thld);
    if(option(1) == 1),
      xl = HardTh(xl,thld);
    end;
  else,
    error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
  end;
  xd = mirdwt(xl,xh,h,L);
else, 					% Denoising by unknown method
  error(['Unknown denoising method',10,...
         'If it is any good we need to have a serious talk :-)']);
end;
option(6) = thld;
option(7) = type;
xn = x - xd;

author: thomaskuestner


text:
function [out] = read_tmsi_poly5(filename, varargin)

% READ_TMSI_POLY5
%
% Use as
%   hdr = read_tmri_poly5(filename)
%   dat = read_tmsi_poly5(filename, hdr, begblock, endblock)
%
% This implementation is as closely as possible based on the original "tms_read",
% which contains the comments
%
%   Changed on 08-10-2014 by TL, TMSi
%   - Now supports loading a file from different directory than the script file
%   - Feedback on the validity of arguments and whether a file could be found or not.
%   - Dialogue is opened when no argument was given.
%
% 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.

needhdr = (nargin==1);
needdat = (nargin>1);

if needdat
  signal = varargin{1};
  begblock = varargin{2};
  endblock = varargin{3};
end

fid = fopen_or_error(filename, 'rb');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the header
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if needhdr
  pos = 31;
  fseek(fid,pos,'bof');
  version = fread(fid,1,'int16');
  if version == 203
    frewind(fid);
    signal.header.FID                   = fread(fid,31,'uchar');
    signal.header.VersionNumber         = fread(fid, 1,'int16');
  else % version 204
    frewind(fid);
    signal.header.FID                   = fread(fid,32,'uchar');
    signal.header.VersionNumber         = fread(fid, 1,'int16');
  end
  signal.header.MeasurementName       = fread(fid,81,'uchar');
  signal.header.FS                    = fread(fid, 1,'int16');
  signal.header.StorageRate           = fread(fid, 1,'int16');
  signal.header.StorageType           = fread(fid, 1,'uchar');
  signal.header.NumberOfSignals       = fread(fid, 1,'int16');
  signal.header.NumberOfSamplePeriods = fread(fid, 1,'int32');
  signal.header.EMPTYBYTES            = fread(fid, 4,'uchar');
  signal.header.StartMeasurement      = fread(fid,14,'uchar');
  signal.header.NumberSampleBlocks    = fread(fid, 1,'int32');
  signal.header.SamplePeriodsPerBlock = fread(fid, 1,'uint16');
  signal.header.SizeSignalDataBlock   = fread(fid, 1,'uint16');
  signal.header.DeltaCompressionFlag  = fread(fid, 1,'int16');
  signal.header.TrailingZeros         = fread(fid,64,'uchar');

  %conversion to char of text values
  signal.header.FID               = char(signal.header.FID);
  signal.header.MeasurementName   = char(signal.header.MeasurementName(2:signal.header.MeasurementName(1)+1));
  %    signal.header.StartMeasurement  = typecast(int8(signal.header.StartMeasurement),'int16'); % commented by RS 1feb11

  signal.fs = signal.header.FS;

  %check for right fileversion
  % if signal.header.VersionNumber ~= 203
  %    disp('Wrong file version! Imput file must be a Poly5/TMS version 2.03 file!');
  %    return;
  % end;

  % Signal description
  for g=1:signal.header.NumberOfSignals,
    signal.description(g).SignalName        = fread(fid,41,'uchar');
    signal.description(g).Reserved          = fread(fid, 4,'uchar');
    signal.description(g).UnitName          = fread(fid,11,'uchar');
    signal.description(g).UnitLow           = fread(fid, 1,'float32');
    signal.description(g).UnitHigh          = fread(fid, 1,'float32');
    signal.description(g).ADCLow            = fread(fid, 1,'float32');
    signal.description(g).ADCHigh           = fread(fid, 1,'float32');
    signal.description(g).IndexSignalList   = fread(fid, 1,'int16');
    signal.description(g).CacheOffset       = fread(fid, 1,'int16');
    signal.description(g).Reserved2         = fread(fid,60,'uchar');

    % conversion of char values (to right format)
    signal.description(g).SignalName = char(signal.description(g).SignalName(2:signal.description(g).SignalName(1)+1));
    signal.description(g).UnitName   = char(signal.description(g).UnitName(2:signal.description(g).UnitName(1)+1));
  end

  % return the header details
  out = signal;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%read data blocks
NB = signal.header.NumberSampleBlocks;
SD = signal.header.SizeSignalDataBlock;
NS = signal.header.NumberOfSignals;
NS_32bit = NS/2;

if needdat

  for g=begblock:endblock

    % jump to right position in file
    if signal.header.VersionNumber == 203
      pos = 217 + NS*136 + (g-1) *(86+SD);
    else
      pos = 218 + NS*136 + (g-1) *(86+SD);
    end
    fseek(fid,pos,'bof');

    signal.block(g).PI = fread(fid,1,'int32'); %period index
    fread(fid,4,'uchar'); %reserved for extension of previous field to 8 bytes
    signal.block(g).BT = fread(fid,14/2,'int16'); %dostime
    fread(fid,64,'uchar'); %reserved
    data = single(fread(fid,SD/4,'float32'));

    % Convert data to 32bit values.
    % In case also 16bit values have to be measured, these values are typecasted below:
    % data = fread(fid,SD/2,'int16'); %read data
    % data = int16(data);
    % data = typecast(data,'int32');
    % signal.block(g).DATA = data;
    signal.data{g} = reshape(data,NS_32bit,SD/(NS_32bit*4));
  end % for

  % Converting data to a usable format
  dat = double(cell2mat(signal.data(begblock:endblock)));

  for g = 1:NS_32bit  % represent data in [uV]
    dat(g,:) = (dat(g,:) - signal.description(g*2).ADCLow)./(signal.description(g*2).ADCHigh - signal.description(g*2).ADCLow)  .* (signal.description(g*2).UnitHigh - signal.description(g*2).UnitLow) + signal.description(g*2).UnitLow ; %correction for uV
  end

  % return the data
  out = dat;
end

fclose(fid);

author: spm


text:
function plotXY(x1, x2, y1, y2,h_newfig, fontsize)

% fontsize= 20;
if h_newfig
    figure;
    set(gcf,'position',[500,100,1000,650])
    set(gcf,'color','w');
end

mArrow2(x1,0,x2,0,{'color','k'});
mArrow2(0,y1,0,y2,{'color','k'});
xlim([x1,x2])
ylim([y1,y2])
set(gca,'visible','off');

% (x1+1):10:(x2-1)
for i = unique(round([linspace(x1*0.9, x2*0.9,7),0]))
    line([i i],[-(y2-y1)*0.005 (y2-y1)*0.005],'color','k')
    line([-(x2-x1)*0.005 (x2-x1)*0.005],[i i],'color','k')

    % xticks
    if i<=0
        t= text(i-(y2-y1)*0.03, -(x2-x1)*0.025,num2str(round(i)));
    else
        t= text(i-(y2-y1)*0.02, -(x2-x1)*0.025,num2str(round(i)));
    end
    t.FontSize=fontsize;

    % yticks
    if i<0
        t = text((y2-y1)*0.01, i, num2str(round(i)));
    elseif i>0
        t = text((y2-y1)*0.01, i, num2str(round(i)));
    end
    t.FontSize = fontsize;
end

t = text(x2*0.97,(y2-y1)*0.03,'$$x$$','Interpreter','latex');
t.FontSize = fontsize;
t = text(-(x2-x1)*0.03,y2*1,'$$y$$','Interpreter','latex');
t.FontSize = fontsize;
axis square

end

function [ h ] = mArrow2(x1,y1,x2,y2,props)

h = annotation('arrow');
set(h,'parent', gca, ...
    'position', [x1,y1,x2-x1,y2-y1], ...
    'HeadLength', 10, 'HeadWidth', 10, 'HeadStyle', 'cback1', ...
    props{:} );
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
%APPLYCHANNELGAINS  Implements an efficient fixed-point approximation for applying channel gains, which is the last step of multiple white balance algorithms
%
%     dst = cv.applyChannelGains(src, gains)
%
% ## Input
% * __src__ Input 3-channel image in the BGR color space (either `uint8` or
%   `uint16`).
% * __gains__ gains for the [B,G,R] channels.
%
% ## Output
% * __dst__ Output image of the same size and type as `src`.
%
% See also: cv.GrayworldWB, cv.LearningBasedWB, imapplymatrix
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%% FUNCTION example_LSSMTC.m
%   Example of multi-task clustering.
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Quanquan Gu, Jiayu Zhou, and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 17, 2012.
%
%% RELATED PAPERS
%   [1]Quanquan Gu, Jie Zhou. Learning the shared subspace
%      for multi-task clustering and transductive transfer
%      classification. In Proc. of the International Conference
%      on Data Mining (ICDM), Miami, Florida, USA, 2009.

clear
close all
clc

addpath('../MALSAR/functions/mutli-task clustering/');
addpath('../MALSAR/utils/');

%===========================================
%This is a demo on comp v.s. sci data set

path = '../data/Newsgroup/comp.vs.sci/';

m = 2;
cellX = cell(m,1);
cellgnd = cell(m,1);
for i=1:m
    task_data = load([path 'Task' num2str(i)]);
    cellX{i} = task_data.fea';
    % if the evaluation metrices are not needed
    % then cellgnd can be set to [].
    cellgnd{i} = task_data.gnd;
end

c = length(unique(cellgnd{1})); % cluster number

opts.tFlag = 2; % termination: run maximum iteration.
opts.maxIter = 20; % maximum iteration.

l = 2; % subspace dimension
lambda = 0.75; % multi-task regulariation paramter
[W cellM cellP residue cellAcc cellNMI] = LSSMTC(cellX,cellgnd,c,l,lambda,opts);

figure
for i = 1:m
    subplot(1,m,i)
    plot(cellAcc{i})
    ylabel('accuracy')
    xlabel('iteration')
    title(sprintf('Task %u', i))
end
print('-dpdf', '-r600', 'LSSMTC_acc');

figure
for i = 1:m
    subplot(1,m,i)
    plot(cellNMI{i})
    ylabel('normalized mutual information')
    xlabel('iteration')
    title(sprintf('Task %u', i))
end
print('-dpdf', '-r600', 'LSSMTC_nmi');

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [meanDoseV, stdDoseV, DVHm, volsV, shiftXv, shiftYv, shiftZv] = simulate_organ_motion(planC,structNumV,doseNum, igrtFlag)
%function [DVH_out, DVHm, volsV] = simulate_organ_motion(planC,structNumV,doseNum)
%
%APA, 09/21/2012
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

%Number of Fractions and Trials
numFractions = 1; % 48
numTrials    = 1; % 200
numStructs   = length(structNumV);

shiftXv = zeros(1,numTrials);
shiftYv = shiftXv;
shiftZv = shiftXv;

structurenamesC = {planC{indexS.structures}.structureName};
rectumIndex = getMatchingIndex('rectum_o',lower(structurenamesC),'exact');
rectalCS = calculate_structure_cross_sectional_area(rectumIndex,planC);
bladderVol = getMatchingIndex('bladder_o',lower(structurenamesC),'exact');

% % Systematic Shifts in cm (Inter-treatment)
% XdispSys_Std = 0; % L-R
% YdispSys_Std = 0; % A-P
% ZdispSys_Std = 0; % S-I

% % Random Shifts in cm (inter-fraction)
% XdispRnd_Std = 1; % L-P
% YdispRnd_Std = 1; % A-P
% ZdispRnd_Std = 1; % S-I

%Distances are in cm
XpdfMean = 0; % L-P
YpdfMean = 0; % A-P
ZpdfMean = 0; % S-I

%Get scan associated with doseNum
scanNum = getAssociatedScan(planC{indexS.dose}(doseNum).assocScanUID, planC);

if isempty(scanNum) %Assume dose is associated with this scan
    scanNum = getStructureAssociatedScan(structNumV(1),planC);
end

%Get reference transformation matrix for doseNum
if ~isempty(scanNum) && isempty(planC{indexS.dose}(doseNum).transM)
    referenceTransM = planC{indexS.scan}(scanNum).transM;
else
    referenceTransM = planC{indexS.dose}(doseNum).transM;
end

%Store the reference doseArray
referenceDoseArray = planC{indexS.dose}(doseNum).doseArray;

%Try and get a binWidth from stateS.  If it doesnt exist, get it from
%the CERROptions file (allows this function to be called outside CERR)
global stateS;
if ~isempty(stateS) && isfield(stateS, 'optS') && isfield(stateS.optS, 'DVHBinWidth') && ~isempty(stateS.optS.DVHBinWidth)
    binWidth = stateS.optS.DVHBinWidth;
else
    optS = CERROptions;
    binWidth = optS.DVHBinWidth;
end

%Compute DVH at given dose
structCount = 0;
for structNum = structNumV
    structCount = structCount + 1;
    [dosesCurrentV{structCount}, volsV{structCount}] = getDVH(structNum, doseNum, planC);
end

%Divide doseArray into fractions
planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray / numFractions;

%Get systematic errors
% deltaX_systematic = randn(1, numTrials) * XdispSys_Std + XpdfMean;
% deltaY_systematic = randn(1, numTrials) * YdispSys_Std + YpdfMean;
% deltaZ_systematic = randn(1, numTrials) * ZdispSys_Std + ZpdfMean;

hWait = waitbar(0,'Computing plan robustness...');

if isempty(referenceTransM)
    referenceTransM_matrix = eye(4);
else
    referenceTransM_matrix = referenceTransM;
end

%Obtain DVH calculation points in blocks
optS    = planC{indexS.CERROptions};

ROIImageSize = [planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension1  planC{indexS.scan}(scanNum).scanInfo(1).sizeOfDimension2];

deltaY = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;

%Get raster segments for structure.
structCount = 0;
for structNum = structNumV

    structCount = structCount + 1;

    [segmentsM, planC, isError] = getRasterSegments(structNum, planC);

    if isempty(segmentsM)
        isError = 1;
    end
    numSegs = size(segmentsM,1);

    %Relative sampling of ROI voxels in this place, compared to CT spacing.
    %Set when rasterSegments are generated (usually on import).
    sampleRate = optS.ROISampleRate;

    %Sample the rows
    indFullV =  1 : numSegs;
    if sampleRate ~= 1
        rV = 1 : length(indFullV);
        rV([rem(rV+sampleRate-1,sampleRate)~=0]) = [];
        indFullV = rV;
    end

    %Block process to avoid swamping on large structures
    if isfield(optS, 'DVHBlockSize') & ~isempty(optS.DVHBlockSize)
        DVHBlockSize = optS.DVHBlockSize;
    else
        DVHBlockSize = 5000;
    end

    blocks{structCount} = ceil(length(indFullV)/DVHBlockSize);

    start = 1;

    for b = 1 : blocks{structCount}

        %Build the interpolation points matrix

        dummy = zeros(1,DVHBlockSize * ROIImageSize(1));
        x1V = dummy;
        y1V = dummy;
        z1V = dummy;
        volsSectionV{structCount} =  dummy;

        if start+DVHBlockSize > length(indFullV)
            stop = length(indFullV);
        else
            stop = start + DVHBlockSize - 1;
        end

        indV = indFullV(start:stop);

        mark = 1;
        for i = indV

            tmpV = segmentsM(i,1:10);
            delta = tmpV(5) * sampleRate;
            xV = tmpV(3): delta : tmpV(4);
            len = length(xV);
            rangeV = ones(1,len);
            yV = tmpV(2) * rangeV;
            zV = tmpV(1) * rangeV;
            sliceThickness = tmpV(10);
            %v = delta^2 * sliceThickness;
            v = delta * (deltaY*sampleRate) * sliceThickness;
            x1V(mark : mark + len - 1) = xV;
            y1V(mark : mark + len - 1) = yV;
            z1V(mark : mark + len - 1) = zV;
            volsSectionV{structCount}(mark : mark + len - 1) = v;
            mark = mark + len;

        end

        %cut unused matrix elements
        x1V = x1V(1:mark-1);
        y1V = y1V(1:mark-1);
        z1V = z1V(1:mark-1);
        volsSectionV{structCount} = volsSectionV{structCount}(1:mark-1);

        %Get transformation matrices for both dose and structure.
        transMDose    = getTransM('dose', doseNum, planC);
        transMStruct  = getTransM('struct', structNum, planC);

        %Forward transform the structure's coordinates.
        if ~isempty(transMStruct)
            [x1V, y1V, z1V] = applyTransM(transMStruct, x1V, y1V, z1V);
        end

        dvhCalcPtsC{structCount}{b} = [x1V', y1V', z1V'];

        %Interpolate.
        %   [dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);

        %   dosesV = [dosesV, dosesSectionV];
        %   volsV  = [volsV, volsSectionV{structCount}];

        start = stop + 1;

    end

end

%Loop over to generate DVH
for iTrial = 1:numTrials

    [deltaX_systematic, deltaY_systematic, deltaZ_systematic, XdispRnd_Std, YdispRnd_Std, ZdispRnd_Std] = getShiftParameters(rectalCS,bladderVol,igrtFlag);
    deltaX_systematic = min(max(deltaX_systematic,-1.0),1.0);
    deltaY_systematic = min(max(deltaY_systematic,-1.5),1.5);
    deltaZ_systematic = min(max(deltaZ_systematic,-1.0),1.0);
    deltaX = deltaX_systematic + randn(1, numFractions) * XdispRnd_Std + XpdfMean;
    deltaY = deltaY_systematic + randn(1, numFractions) * YdispRnd_Std + YpdfMean;
    deltaZ = deltaZ_systematic + randn(1, numFractions) * ZdispRnd_Std + ZpdfMean;

    % To get planned DVH - comment it otherwise
    deltaX = 0;
    deltaY = 0;
    deltaZ = 0;

    shiftXv(1,iTrial) = mean(deltaX);
    shiftYv(1,iTrial) = mean(deltaY);
    shiftZv(1,iTrial) = mean(deltaZ);

    structCount = 0;
    for structNum = structNumV

        structCount = structCount + 1;

        tmpDoseV = zeros(1,length(dosesCurrentV{structCount}));

        for iFraction = 1:numFractions

            waitbar(((iTrial-1)*numStructs + structCount)/(numTrials*numStructs),hWait)

            transM = referenceTransM_matrix;
            transM(1:3,4) = transM(1:3,4) + [deltaX(iFraction); deltaY(iFraction); deltaZ(iFraction)];

            %Apply the new transM to dose
            planC{indexS.dose}(doseNum).transM = transM;

            % Convert Dose to Gy
            if max(planC{indexS.dose}(doseNum).doseArray(:)) > 150
                planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray/100;
            end

            %Get doses and volumes of points in structure.
            %[dosesV, volsV] = getDVH(structNum, doseNum, planC);

            volsV{structCount}  = [];
            dosesV{structCount} = [];

            for b = 1 : blocks{structCount}
                x1V = dvhCalcPtsC{structCount}{b}(:,1);
                y1V = dvhCalcPtsC{structCount}{b}(:,2);
                z1V = dvhCalcPtsC{structCount}{b}(:,3);

                %Back transform the coordinates into the doses' coordinate system.
                [x1V, y1V, z1V] = applyTransM(inv(transM), x1V, y1V, z1V);

                %Interpolate.
                [dosesSectionV] = getDoseAt(doseNum, x1V, y1V, z1V, planC);

                dosesV{structCount} = [dosesV{structCount}, dosesSectionV];
                volsV{structCount}  = [volsV{structCount}, volsSectionV{structCount}];
            end

            tmpDoseV = tmpDoseV + dosesV{structCount};

        end

        DVHm{structCount}(:,iTrial) = single(tmpDoseV(:));

    end

end

close(hWait)

%Reassign reference transM and doseArray to doseNum
planC{indexS.dose}(doseNum).transM = referenceTransM;
planC{indexS.dose}(doseNum).doseArray = referenceDoseArray;

%Compute Mean and Std Dev in blocks
structCount = 0;
for structNum = structNumV

    structCount = structCount + 1;

    numVoxels = length(volsV{structCount});
    DVHBlockSize = 5000;
    blocks = ceil(numVoxels/DVHBlockSize);
    start = 1;

    for b = 1 : blocks

        if start+DVHBlockSize > numVoxels
            stop = numVoxels;
        else
            stop = start + DVHBlockSize - 1;
        end

        DVH_block = DVHm{structCount}(start:stop,:);

        meanDoseV{structCount}(start:stop) = mean(DVH_block,2);
        stdDoseV{structCount}(start:stop)  = std(DVH_block,0,2);

        start = stop + 1;

    end
end

% %Histogram for the expectation of dose over all trials
% [doseBinsV, volsHistV] = doseHist(meanDoseV, volsV, binWidth);
% DVH_out = [doseBinsV(:)'; volsHistV(:)'];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [fMc, fMc95, fMc90] = calc_McBest(mag, fBinning);
    %CALC_MCBEST Calculate best Mc from combining methods

    % First estimation of magnitude of completeness (maximum curvature)
    [vEvents, vMag] = hist(mag, -2:0.1:6);
    nSel = max(find(vEvents == max(vEvents)));
    fMcStart = vMag(nSel);

    % Data container
    mData = [];

    % Magnitude increment
    if ~exist('fBinning')
        fBinning = 0.1;
    end

    for nCnt = (fMcStart - 0.9):fBinning:(fMcStart + 1.5)
        vSel = mag > (nCnt - (fBinning/2));
        %vSel = mag >= nCnt - 0.0499;
        nNumberEvents = length(mag(vSel));
        if nNumberEvents >= 25;
            [fDummy fBValue fDummy fDummy] =  Catalog.bvalue_lib.bmemag(mag(vSel));

            fStartMag = nCnt; % Starting magnitude (hypothetical Mc)

            % log10(N)=A-B*M
            vMag = [fStartMag:fBinning:15]; % Ending magnitude must be sufficiently high
            vNumber = 10.^(log10(nNumberEvents)-fBValue*(vMag - fStartMag));
            vNumber = round(vNumber);

            % Find the last bin with an event
            nLastEventBin = min(find(vNumber == 0)) - 1;
            if isempty(nLastEventBin)
                nLastEventBin = length(vNumber);
            end

            %    ctM=vMag(ct1);

            % Determine set of all magnitude bins with number of events > 0
            ct = round((vMag(nLastEventBin)-fStartMag)*(1/fBinning) + 1);
            %     ct=0;
            %     for I=fStartMag:fBinning:ctM;
            %       ct=ct+1;
            %     end;

            PM=vMag(1:ct);
            vNumber = vNumber(1:ct);
            [bval, vDummy] = hist(mag(vSel),PM);
            b3 = fliplr(cumsum(fliplr(bval)));    % N for M >= (counted backwards)
            res2 = sum(abs(b3 - vNumber))/sum(b3)*100;
            mData = [mData; nCnt res2];
        else
            mData = [mData; nCnt NaN];
        end
    end

    % Evaluation of results

    % Is fMc90 available
    nSel = min(find(mData(:,2) < 10));
    if isempty(nSel)
        fMc90 = NaN
    else
        fMc90 = mData(nSel,1)
    end

    % Is fMc95 available
    nSel = min(find(mData(:,2) < 5));
    if isempty(nSel)
        fMc95 = NaN;
    else
        fMc95 = mData(nSel,1);
    end

    % ?????
    j =  min(find(mData(:,2) < 10 ));
    if isempty(j) == 1; j =  min(find(mData(:,2) < 15 )); end
    if isempty(j) == 1; j =  min(find(mData(:,2) < 20 )); end
    if isempty(j) == 1; j =  min(find(mData(:,2) < 25 )); end
    %j2 =  min(find(dat(:,2) == min(dat(:,2)) ));

    fMc = mData(j,1);
    if isempty(fMc)
        fMc = NaN;
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function [ n_data, n, p ] = pi_values ( n_data )

%*****************************************************************************80
%
%% PI_VALUES returns values of the Pi function.
%
%  Discussion:
%
%    Pi[n] is the number of primes less than or equal to n.
%
%    In Mathematica, the function can be evaluated by:
%
%      PrimePi[n]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the argument.
%
%    Output, integer P, the value of the function.
%
  n_max = 17;

  n_vec = [ ...
            10, ...
            20, ...
            30, ...
            40, ...
            50, ...
            60, ...
            70, ...
            80, ...
            90, ...
           100, ...
          1000, ...
         10000, ...
        100000, ...
       1000000, ...
      10000000, ...
     100000000, ...
    1000000000 ];

  p_vec = [ ...
             4, ...
             8, ...
            10, ...
            12, ...
            15, ...
            17, ...
            19, ...
            22, ...
            24, ...
            25, ...
           168, ...
          1229, ...
          9592, ...
         78498, ...
        664579, ...
       5761455, ...
      50847534 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    p = 0;
  else
    n = n_vec(n_data);
    p = p_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function hermite_exactness ( quad_filename, degree_max, option )

%*****************************************************************************80
%
%% MAIN is the main program for HERMITE_EXACTNESS.
%
%  Discussion:
%
%    This program investigates a standard Gauss-Hermite quadrature rule
%    by using it to integrate monomials over (-oo,+oo), and comparing the
%    approximate result to the known exact value.
%
%    The user specifies:
%    * the "root" name of the R, W and X files that specify the rule;
%    * DEGREE_MAX, the maximum monomial degree to be checked.
%    * the OPTION (unweighted/physicist weight/probabilist weight)
%
%    OPTION indicates the weight function and normalization:
%    0, Integral ( -oo < x < +oo ) x^n exp(-x*x)               dx.
%    1, Integral ( -oo < x < +oo ) x^n exp(-x*x)               dx.
%    2, Integral ( -oo < x < +oo ) x^n exp(-x*x/2)             dx.
%    3, Integral ( -oo < x < +oo ) x^n exp(-x*x)   / sqrt (pi) dx.
%    4, Integral ( -oo < x < +oo ) x^n exp(-x*x/2) / sqrt(2pi) dx.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_EXACTNESS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Investigate the polynomial exactness of a Gauss-Hermite\n' );
  fprintf ( 1, '  quadrature rule by integrating exponentially weighted\n' );
  fprintf ( 1, '  monomials up to a given degree over the (-oo,+oo) interval.\n' );
%
%  Get the quadrature file root name:
%
  if ( 1 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS:\n' );

    quad_filename = input ( '  Enter the "root" name of the quadrature files: ' );

  end
%
%  Create the names of:
%    the quadrature X file;
%    the quadrature W file;
%    the quadrature R file;
%
  quad_x_filename = strcat ( quad_filename, '_x.txt' );
  quad_w_filename = strcat ( quad_filename, '_w.txt' );
  quad_r_filename = strcat ( quad_filename, '_r.txt' );
%
%  The second command line argument is the maximum degree.
%
  if ( 2 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS:\n' );

    degree_max = input ( '  Please enter the maximum degree to check: ' );

  end
%
%  The third command line argument is OPTION.
%
  if ( 3 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS:\n' );
    fprintf ( 1, '  OPTION specifies the weight function w(x):\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  0, rho(x) = 1\n' );
    fprintf ( 1, '  1, rho(x) = exp ( -x*x  )\n' );
    fprintf ( 1, '  2, rho(x) = exp ( -x*x/2)\n' );
    fprintf ( 1, '  3, rho(x) = exp ( -x*x  ) / sqrt( pi)\n' );
    fprintf ( 1, '  4, rho(x) = exp ( -x*x/2) / sqrt(2pi)\n' );
    fprintf ( 1, '\n' );
    option = input ( '  Please enter OPTION: ' );

  end

  if ( option < 0 || 4 < option )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of OPTION.\n' );
    error ( 'HERMITE_EXACTNESS - Fatal error!' );
  end
%
%  Summarize the input.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_EXACTNESS: User input:\n' );
  fprintf ( 1, '  Quadrature rule X file = "%s".\n', quad_x_filename );
  fprintf ( 1, '  Quadrature rule W file = "%s".\n', quad_w_filename );
  fprintf ( 1, '  Quadrature rule R file = "%s".\n', quad_r_filename );
  fprintf ( 1, '  Maximum degree to check = %d\n', degree_max );
%
%  Read the X file.
%
  [ dim_num, order ] = r8mat_header_read ( quad_x_filename );

  if ( dim_num ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS - Fatal error!\n' );
    fprintf ( 1, '  The spatial dimension should be 1.\n');
    fprintf ( 1, '  The spatial dimension in the X file is %d\n', dim_num );
    error ( 'HERMITE_EXACTNESS - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension = %d\n', dim_num );
  fprintf ( 1, '  Number of points  = %d\n', order );

  x = r8mat_data_read ( quad_x_filename, dim_num, order );
%
%  Read the W file.
%
  [ dim_num2, point_num ] = r8mat_header_read ( quad_w_filename );

  if ( dim_num2 ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS - Fatal error!\n' );
    fprintf ( 1, '  The quadrature weight file should have exactly\n');
    fprintf ( 1, '  one value on each line.\n' );
    error ( 'HERMITE_EXACTNESS - Fatal error!' );
  end

  if ( point_num ~= order )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS - Fatal error!\n' );
    fprintf ( 1, '  The quadrature weight file should have exactly\n' );
    fprintf ( 1, '  the same number of lines as the abscissa file.\n' );
    error ( 'HERMITE_EXACTNESS - Fatal error!' );
  end

  w = r8mat_data_read ( quad_w_filename, 1, order );
%
%  Read the R file.
%
  [ dim_num2, point_num ] = r8mat_header_read ( quad_r_filename );

  if ( dim_num2 ~= dim_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'v - Fatal error!\n' );
    fprintf ( 1, '  The quadrature region file should have the same\n' );
    fprintf ( 1, '  number of values on each line as the abscissa file\n' );
    fprintf ( 1, '  does.\n' );
    error ( 'HERMITE_EXACTNESS - Fatal error!' );
  end

  if ( point_num ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HERMITE_EXACTNESS - Fatal error!\n' );
    fprintf ( 1, '  The quadrature region file should have two lines.\n' );
    error ( 'HERMITE_EXACTNESS - Fatal error!' );
  end

  r = r8mat_data_read ( quad_r_filename, dim_num, 2 );
%
%  Print the input quadrature rule.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test a Gauss-Hermite quadrature rule of\n' );
  fprintf ( 1, '  ORDER = %d\n', order );
  fprintf ( 1, '\n' );
  if ( option == 0 )
    fprintf ( 1, '  OPTION = 0, the unweighted rule for:\n' );
    fprintf ( 1, '  Integral ( -oo < x < +oo ) f(x) dx\n' );
  elseif ( option == 1 )
    fprintf ( 1, '  OPTION = 1, the physicist weighted rule for:\n' );
    fprintf ( 1, '  Integral ( -oo < x < +oo ) f(x) * exp(-x*x) dx\n' );
  elseif ( option == 2 )
    fprintf ( 1, '  OPTION = 2, the probabilist weighted rule for:\n' );
    fprintf ( 1, '  Integral ( -oo < x < +oo ) f(x) * exp(-x*x/2) dx\n' );
  elseif ( option == 3 )
    fprintf ( 1, '  OPTION = 3, the physicist normalized weighted rule for:\n' );
    fprintf ( 1, '  Integral ( -oo < x < +oo ) f(x) * exp(-x*x) / sqrt(pi) dx\n' );
  elseif ( option == 4 )
    fprintf ( 1, '  OPTION = 4, the probabilist normalized weighted rule for:\n' );
    fprintf ( 1, '  Integral ( -oo < x < +oo ) f(x) * exp(-x*x/2) / sqrt(2 pi) dx\n' );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weights W:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : order
    fprintf ( 1, '  w(%d) = %24.16f\n', i, w(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Abscissas X:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : order
    fprintf ( 1, '  x(%d) = %24.16f\n', i, x(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Region R:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 2
    fprintf ( 1, '  r(%d) = %e\n', i, r(i) );
  end
%
%  Explore the monomials.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A Gauss-Hermite rule would be able to exactly\n' );
  fprintf ( 1, '  integrate monomials up to and including \n' );
  fprintf ( 1, '  degree = %d\n', 2 * order - 1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Degree      Error\n' );
  fprintf ( 1, '\n' );

  for degree = 0 : degree_max

    quad_error = hermite_monomial_quadrature ( degree, order, option, w, x );

    fprintf ( 1, '  %2d  %24.16f\n', degree, quad_error );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_EXACTNESS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function value = hermite_integral ( n, option )

%*****************************************************************************80
%
%% HERMITE_INTEGRAL evaluates a monomial Hermite integral.
%
%  Discussion:
%
%    H(n,1) = Integral ( -oo < x < +oo ) x^n exp(-x^2) dx
%    H(n,1) is 0 for n odd.
%    H(n,1) = (n-1)!! * sqrt(pi) / 2^(n/2) for n even.
%
%    H(n,2) = Integral ( -oo < x < +oo ) x^n exp(-x^2/2) dx
%    H(n,2) is 0 for n odd.
%    H(n,2) = (n-1)!! * sqrt(2*pi) for n even.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the integral.
%    0 <= N.
%
%    Input, integer OPTION, the integral has the form:
%    0, Integral ( -oo < x < +oo ) x^n exp(-x*x)               dx.
%    1, Integral ( -oo < x < +oo ) x^n exp(-x*x)               dx.
%    2, Integral ( -oo < x < +oo ) x^n exp(-x*x/2)             dx.
%    3, Integral ( -oo < x < +oo ) x^n exp(-x*x)   / sqrt (pi) dx.
%    4, Integral ( -oo < x < +oo ) x^n exp(-x*x/2) / sqrt(2pi) dx.
%
%    Output, real VALUE, the value of the integral.
%
  if ( n < 0 )

    value = - r8_huge ( );

  elseif ( mod ( n, 2 ) == 1 )

    value = 0.0;

  elseif ( option == 0 )

    value = r8_factorial2 ( n - 1 ) * sqrt ( pi ) / 2.0^( n / 2 );

  elseif ( option == 1 )

    value = r8_factorial2 ( n - 1 ) * sqrt ( pi ) / 2.0^( n / 2 );

  elseif ( option == 2 )

    value = r8_factorial2 ( n - 1 ) * sqrt ( 2.0 * pi );

  elseif ( option == 3 )

    value = r8_factorial2 ( n - 1 ) / 2.0^( n / 2 );

  elseif ( option == 4 )

    value = r8_factorial2 ( n - 1 );

  end

  return
end
function quad_error = hermite_monomial_quadrature ( expon, order, option, w, x )

%*****************************************************************************80
%
%% HERMITE_MONOMIAL_QUADRATURE applies a quadrature rule to a monomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer EXPON, the exponent.
%
%    Input, intege ORDER, the number of points in the rule.
%
%    Input, integer OPTION, the integral has the form:
%    0, Integral ( -oo < x < +oo ) x^n exp(-x*x)               dx.
%    1, Integral ( -oo < x < +oo ) x^n exp(-x*x)               dx.
%    2, Integral ( -oo < x < +oo ) x^n exp(-x*x/2)             dx.
%    3, Integral ( -oo < x < +oo ) x^n exp(-x*x)   / sqrt (pi) dx.
%    4, Integral ( -oo < x < +oo ) x^n exp(-x*x/2) / sqrt(2pi) dx.
%
%    Input, real W(ORDER), the quadrature weights.
%
%    Input, real X(ORDER), the quadrature points.
%
%    Output, real QUAD_ERROR, the quadrature error.
%

%
%  Get the exact value of the integral of the unscaled monomial.
%
  exact = hermite_integral ( expon, option );
%
%  Evaluate the unweighted monomial at the quadrature points.
%
  if ( option == 0 )
    value(1:order) = exp ( - x(1:order).^2 ) .* x(1:order).^expon;
  elseif ( option == 1 )
    value(1:order) = x(1:order).^expon;
  elseif ( option == 2 )
    value(1:order) = x(1:order).^expon;
  elseif ( option == 3 )
    value(1:order) = x(1:order).^expon;
  elseif ( option == 4 )
    value(1:order) = x(1:order).^expon;
  end
%
%  Compute the weighted sum.
%
  quad = ( w(1:order) * value(1:order)' );
%
%  Absolute error for cases where exact integral is zero,
%  Relative error otherwise.
%
  if ( exact == 0.0 )
    quad_error = abs ( quad );
  else
    quad_error = abs ( quad - exact ) / abs ( exact );
  end

  return
end
function value = r8_factorial2 ( n )

%*****************************************************************************80
%
%% R8_FACTORIAL2 computes the double factorial function.
%
%  Discussion:
%
%    FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 )  (N even)
%                    = Product ( N * (N-2) * (N-4) * ... * 1 )  (N odd)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the double factorial function.
%    If N is less than 1, VALUE is returned as 1.
%
%    Output, real VALUE, the value of N!!.
%
  value = 1;

  if ( n < 1 )
    return
  end

  while ( 1 < n )
    value = value * n;
    n = n - 2;
  end

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function rbfard2KernDisplay(kern, spacing)

% RBFARD2KERNDISPLAY Display parameters of the RBFARD2 kernel.
% FORMAT
% DESC displays the parameters of the automatic relevance determination radial basis function
% kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : rbfard2KernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Michalis K. Titsias, 2009

% KERN

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('RBF ARD Variance: %2.4f\n', kern.variance)
fprintf(spacing);
for i = 1:kern.inputDimension
  fprintf(spacing);
  fprintf('RBF ARD Input %d scale: %2.4f\n', i, kern.inputScales(i))
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function C = nextperm(N,K)
%NEXTPERM Loop through permutations without replacement.
% NEXTPERM(N,K), when first called, returns a function handle.  This
% function handle, when called, returns the next permutation without
% replacement of K elements taken from the set 1:N. This can be useful
% when the number of such permutations is too large to hold in memory at
% once.  If the number of such permutations is not too large, use
% COMBINATOR(N,K,'p') (on the FEX) instead.
%
% The number of permutations without replacement is:  N!/(N-K)!
%   where N >= 1, N >= K >= 0
%
% Examples:
%
%     % To use each permutation one at a time, put it in a loop.
%     N = 4;  % Length of the set.
%     K = 3;  % Number of samples taken for each sampling.
%     H = nextperm(N,K);
%     for ii = 1:((prod(1:N)/(prod(1:(N-K)))))
%         A = H();
%         % Do stuff with A: use it as an index, etc.
%     end
%
%
%     % To build all of the permutations, do this (See note below):
%     ROWS = (prod(1:N)/(prod(1:(N-K))));
%     C = ones(ROWS,K);
%     for ii = 1:ROWS
%         C(ii,:) = H();
%     end
%     %Note this is a lot slower than using combinator(N,K,'p')
%
% The function handle will cycle through when the final permutation is
% returned.
%
% See also,  nchoosek, perms, combinator, npermutek (both on the FEX)
%
% Author:   Matt Fig
% Contact:  popkenai@yahoo.com
% Date: 6/9/2009
% Reference:  http://mathworld.wolfram.com/BallPicking.html

% Arg checking.
if K>N
    error('K must be less than or equal to N.')
end

if isempty(N) || K == 0
   C = [];
   return
elseif numel(N)~=1 || N<=0 || ~isreal(N) || floor(N) ~= N
    error('N should be one real, positive integer. See help.')
elseif numel(K)~=1 || K<0 || ~isreal(K) || floor(K) ~= K
    error('K should be one real non-negative integer. See help.')
end

CNT = 0;  % Initializations for the nested function.
G = [];
TMP = [];
op = 1; % These are the variables which are passed to subfunc.
blk = 1;  % Keeps track of the permutation blocks in the subfunc.
idx = 1:N/2; % Index vectors for subfunc.
idx2 = idx;
WV = 1:K;  % Initial WV.
lim = 0;  % Sets the limit for working index.
inc = K;  % Controls which element of WV is being worked on.
cnt = 1;  % Keeps track of which block we are in.  See loop below.
BC = prod(1:N)/(prod(1:(N-K))); % Number of permutations.
L = prod(1:K);   % Size of the blocks.
CNT = 0;  % Counter for blocks.
Floop = BC - prod(1:K);  % Length of blocks.
A2 = (N-K+1):N;  % Seed for the final block.
B2 = 1:K;
ii = 1;

if K==1
    C = @nestfunc2;  % Return argument is a func handle.
else
    C = @nestfunc;
end

    function B = nestfunc2
    % A handle to this function is passed back from the main func if K=1.
        if WV > N
            B = 1;
            WV = 2;
            return
        end
       B = WV(1);
       WV = WV + 1;
    end

    function B = nestfunc
    % A handle to this function is passed back from the main func.
            if ii<=Floop
                if CNT == 0
                    for jj = 1:inc
                        WV(K + jj - inc) = lim + jj;
                    end
                    % This is the first combination.  We will permute it
                    % below for the rest of the calls in this block.
                    B = WV;
                    cnt = cnt + L;

                    if lim<(N-inc)
                        inc = 0;  % Reset this guy.
                    end

                    TMP = WV;  % This serves as seed for perm index.
                    inc = inc+1;  % Increment the counter.
                    lim = WV(K+1-inc);  % Limit for working index.
                    CNT = 1;
                    G = 1:K; % Seed for nextp subfunc
                    op = 1;
                    blk = 1;
                    idx = 1:N/2;
                    idx2 = idx;
                    ii = ii + 1; % "Loop" index.
                else
                    % Permute the seed.
                    [G,op,blk,idx,idx2] = nextp(G,op,blk,idx,idx2,K);
                    B = TMP(G);  % Index into current combination.
                    CNT = CNT + 1;

                    if CNT==(L) % Goes back to if for next seed (combin).
                        CNT = 0;
                    end

                    ii = ii + 1;
                end
            else
                if ii == Floop+1  % We are at the last block
                    op = 1;  % Re-initialize for subfunc.
                    blk = 1;
                    idx = 1:N/2;
                    idx2 = idx;
                    B = A2;  % Seed for this last block.
                    cnt = cnt + 1;
                    ii = ii + 1;
                    return
                elseif ii == BC + 1  % Time to start over.
                    WV = 1:K;  % Seed for first block.
                    op = 1;  % Re-initialize for subfunc.
                    blk = 1;
                    idx = 1:N/2;
                    idx2 = idx;
                    inc = 1; % Reset the incrementer.
                    lim = WV(K+1-inc); % And the lim.
                    cnt = L + 1;  % Reset block counter.
                    CNT = 1;
                    ii = 2;  % Reset "Loop" index.
                    B = WV;
                    TMP = WV;
                    B2 = 1:K;
                    G = 1:K;
                    return
                end
                % Permute seed the seed.
                [B2,op,blk,idx,idx2] = nextp(B2,op,blk,idx,idx2,K);
                B = A2(B2);
                cnt = cnt + 1;
                ii = ii + 1;
            end
    end
end

function [x,op,blk,idx,idx2] = nextp(x,op,blk,idx,idx2,K)
% Delivers one permutation at a time.  This is a modification of an
% algorithm attributed to H. F. Trotter.
% x = 1:4;
% op = 1;
% blk = 1;
% idx = 1:length(x)/2;
% idx2 = idx;
% C(1,:) = x;
% for jj = 2:factorial(length(x))
%    [x,op,blk,idx,idx2] = nextp(x,op,blk,idx,idx2,length(x));
%    C(jj,:) = x;
% end

if op<K
    np = op + 1; % Index of where the 1 goes.
    x(op) = x(np);  % Here we are just doing the switcheroo on adjacents.
    x(np) = 1;  % np is the current position of the 1.
    op = np;  % op is the old position of the 1, for next time.
    return
else
    x(K) = x(1);  % Here we are switching the endpoints of the vect.
    x(1) = 1;
    op = 1;  % Reset to the first position.
    blk = blk + 1;  % Keep track of this block. Reset every N*(N-1)th call.

    if blk<K % Not through with this block yet.
        return
    end

    blk = 1;  % If here, we need to mix internal elements of the vect.
    low = 2;  % Start with the second element.
    upp = K - 1;  % And the next to last element.

    while true  % We always break this loop.
        % In here, on every N*(N-1)th call, we are mixing internal elems.
        % The number of times through the while loop is usually very small,
        % even for large N most calls will go through less than 3 times.
        cur = idx(low);  % Holds the current index to switch.

        if cur==upp
            np = low;
            idx2(low) = idx2(low) + 1;
        else
            np = cur + 1;  % For next iter.
        end

        tmp = x(cur);  % Switcheroo.
        x(cur) = x(np);
        x(np) = tmp;
        idx(low) = np;

        if idx2(low)<upp
            break  % Out of while loop, we've mixed them enough.
        end

        idx2(low) = low;
        low = low + 1; % These two march towards each other.
        upp = upp - 1;
    end
end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function envihdrwrite(filename,dims,datatype,resolu,UL,interleave,zone)
%  write envi format file of unit8 and int16
%  for example: enviwrite('testzz',data,'int16',[30,30],jiUL,'bsq',ZC,'this is for test only');
%  i,j => x, y => rows and cols
%  Input  1)  'filename' is the name of the file
%  Input  2)  'dims' is the dimenstion of the image (rows,cols,bands)
%  Input  3)  'dataype' is the datatype to write 'unit8' or 'int16'
%  Input  4)  'resolu' is the resolution of the pixel
%  Input  5)  'UL' is the UpperLeftPoint X Y of the UL pixel (not center of UL pixel)
%  Input  6)  'interleave' is the bsq bip bil format
%  Input  7)  'zone' is the UTM zone

% if strcmp(datatype,'uint8')
% %     envi_data=uint8(envi_data);
%     dt=1;
% elseif strcmp(datatype,'int16')
% %     envi_data=int16(envi_data);
%     dt=2;
% elseif strcmp(datatype,'uint16')
% %     envi_data=uint16(envi_data);
%     dt=12;
% else
%     fprintf('Invalid write data type!\n');
%     return;
% end

if strcmp(datatype,'uint8')
%    envi_data=uint8(envi_data);
    dt=1;
elseif strcmp(datatype,'int16')
%    envi_data=int16(envi_data);
    dt=2;
elseif strcmp(datatype,'int32')
%    envi_data=int32(envi_data);
    dt=3;
elseif strcmp(datatype,'uint16')
%    envi_data=uint16(envi_data);
    dt=12;
elseif strcmp(datatype,'single')
%    envi_data=single(envi_data);
    dt=4;
else
    fprintf('Invalid write data type!\n');
    return;
end

% n_dims=size(envi_data);

nrows=dims(1);
ncols=dims(2);
bands=1;

if length(dims)==3
    bands=dims(3);
end

% multibandwrite(envi_data,filename,interleave);

filename_hdr=[filename,'.hdr'];

fid_out=fopen(filename_hdr,'wt');

fprintf(fid_out,'ENVI\n');
fprintf(fid_out,'description = {Landsat Scientific Data}\n');
fprintf(fid_out,'samples = %d\n',ncols); % samples is for j
fprintf(fid_out,'lines   = %d\n',nrows); % lines is for i
fprintf(fid_out,'bands   = %d\n',bands);
fprintf(fid_out,'header offset = 0\n');
fprintf(fid_out,'file type = ENVI Standard\n');
fprintf(fid_out,'data type = %d\n',dt);
fprintf(fid_out,'interleave = %s\n',interleave);
fprintf(fid_out,'sensor type = Landsat\n');
fprintf(fid_out,'byte order = 0\n');

if isempty(zone)
    %label as Landsat ARD
    fprintf(fid_out, 'map info = {Albers Conical Equal Area, 1, 1, %d, %d, %d, %d, WGS-84}\n',UL(1),UL(2),resolu(1),resolu(2));
    fprintf(fid_out, 'projection info = {9, 6378140, 6356755.288157528, 23, -96, 0, 0, 29.5, 45.5,WGS-84, Albers Conical Equal Area}\n');
    fprintf(fid_out, 'coordinate system string = {PROJCS["WGS_1984_Albers",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378140.0,298.257]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-96.0],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_origin",23.0],UNIT["Meter",1.0]]}');
else
    if (zone > 0)
        fprintf(fid_out, 'map info = {UTM, 1.000, 1.000, %d, %d, %d, %d, %d, North, WGS-84, units=Meters}',UL(1),UL(2),resolu(1),resolu(2),zone);
    else
        fprintf(fid_out, 'map info = {UTM, 1.000, 1.000, %d, %d, %d, %d, %d, South, WGS-84, units=Meters}',UL(1),UL(2),resolu(1),resolu(2),-zone);
    end
end

fclose(fid_out);
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_KUKA_KR6_2(robot, T)
%   Solves the inverse kinematic problem for the KUKA KR6_2 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC__KUKA_KR6_2 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('kuka', 'KR6_2');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%   Author: C.Escoto, E.Leon, V. Martinez & L. Mijares
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_kuka_kr6_2(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=abs(d(6));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));
%q1=atan2(Py, Px);

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist(robot, q(:,i), T, 1,'geometric'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist(robot, q(:,i), T, -1, 'geometric'); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=abs(a(2));
L3=abs(d(4));
A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%See geometry
L2=abs(a(2));
L3=abs(d(4));

A2 = abs(a(3));

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) =  atan(L3/ A2)-eta;
q3(2) =  atan(L3/ A2)+ eta;

repo: ARTE