text:
function Tc = slcompresstext(T, varargin)
%SLCOMPRESSTEXT Compresses a cell array of text
%
% $ Syntax $
%   - Tc = slcompresstext(T, ...)
%
% $ Arguments $
%   - T:        the source text to be compressed
%   - Tc:       the compressed text
%
% $ Description $
%   - Tc = slcompresstext(T, ...) compresses the text T, represented in
%     cell array of lines. You can specify the properties to control the
%     process of compression.
%     \*
%     \t    Table 1. Properties of Text Compression
%     \h    name        &    description
%          'rmempty'    &  whether to remove empty line (default = true)
%          'proc'       &  the method of processing each line
%                          (default='trim')
%                          - 'off': do not process
%                          - 'trim': trim the leading and trailing spaces
%                          - 'deblank': trim only the trailing spaces
%
% $ History $
%   - Created by Dahua Lin, on Aug 9th, 2006
%

%% parse and verify input arguments

if ~iscell(T)
    error('sltoolbox:invalidarg', ...
        'T should be a cell array of strings');
end

opts.rmempty = true;
opts.proc = 'trim';
opts = slparseprops(opts, varargin{:});

switch opts.proc
    case 'off'
        procfunc = [];
    case 'trim'
        procfunc = 'strtrim';
    case 'deblank'
        procfunc = 'deblank';
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid string processing method %s', opts.proc);
end

%% Process

if ~isempty(procfunc)
    Tc = slfiltertext(T, procfunc);
else
    Tc = T;
end

%% Select

if opts.rmempty
    nlines = length(Tc);
    is_effline = true(nlines, 1);
    has_deleted = false;
    for i = 1 : nlines
        if isempty(Tc{i})
            is_effline(i) = false;
            has_deleted = true;
        end
    end

    if has_deleted
        Tc = Tc(is_effline);
    end
end

author: lmthang


text:
 function [C, rj] = penalty2_design(type, varargin)
%|function [C, rj] = penalty2_design(type, ['leak'|'tight'], wang, ang, mask)
%|
%| Design the penalty matrix "C" (and penalty coefficients "rj")
%| for a quadratic penalty R(x) = 1/2 x' C * C * x with:
%| 1st-order differences and a 2nd-order neighborhood (8 neighbors).
%| For 2D parallel-beam tomography with shift-invariant blur.
%| Design based on Fourier method in fessler:03:aat (IEEE MIC, 2003).
%|
%| in
%|	type			'test' or 'quad,d1,n2' or 'quad,d1,n1'
%|				(n2 for usual 2nd order neighborhood)
%|	wang	[na nx ny]	angular weights for each pixel
%|	ang	[na]		angles
%|	mask	[nx ny]		(logical) reconstruction support
%| out:
%|	C	[4*nx*ny nx*ny] "modified" C containing sqrt{r_j} factors
%|	rj	[nx ny 4]	horiz,vert,diag1,diag2 coefficients
%|
%| Copyright 2003-5-23, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

% run a self-test to compare my analytical solution
% to the numerical solution computed using NNLS.
if nargin == 1 && streq(type, 'test')
	penalty2_design_test
return
end

% analytical design of 1st-order difference, 2nd-order neighborhood
if streq(type, 'quad,d1', 7)
	if length(varargin) < 1 || length(varargin) > 3, ir_usage, end

	if ischar(varargin{1})
		Ctype = varargin{1};
		varargin = varargin{2:end};
	else
		Ctype = 'leak';
	end

	wang = varargin{1};
	[na nx ny] = size(wang);
	np = nx*ny;
	wang = reshape(wang, [na np]);

	if length(varargin) >= 2
		ang = col(varargin{2});
		if na ~= length(ang), error 'angle dim', end
	else
		ang = [0:(na-1)]'/na * pi; % [0,pi)
	end

	if length(varargin) == 3
		mask = varargin{3};
		if any(size(mask) ~= [nx ny]), error 'mask dims', end
	else
		mask = true(nx,ny);
	end

	% dot product of each basis with wj
	% the "mean" takes care of normalization
	basis = [ones(size(ang)) cos(2*ang) sin(2*ang)];
	dot_products = zeros(ncol(basis),np);
	for ib=1:ncol(basis)
		b = basis(:,ib);
		b = repmat(basis(:,ib), 1, np);
		dot_products(ib,:) = mean(b .* wang, 1);
	end, clear ib b

	sptmp = @(tmp) spdiag(tmp, 'nowarn'); % for now

	if streq(type, 'quad,d1,n1') % 1st-order neighborhood
		rj = penalty2_design_d1_n1(dot_products); % [2 np]
		rj = permute(reshape(rj, [2 nx ny]), [2 3 1]); % [nx ny 2]
		C = C2sparse('leak', mask, 4, 0, 0);
		i1 = 1:(nx*ny);
		C = [	sptmp(sqrt(col(rj(:,:,1)))) * C(0*nx*ny+i1,:);
			sptmp(sqrt(col(rj(:,:,2)))) * C(1*nx*ny+i1,:)];

	elseif streq(type, 'quad,d1,n2') % 2nd-order neighborhood
		rj = penalty2_design_d1_n2(dot_products); % [4 np]
		rj = permute(reshape(rj, [4 nx ny]), [2 3 1]); % [nx ny 4]
		% fix: this must depend on flip_y !!
		printf('Warn: penalty2_design.m may fail if flip_y = -1')
		rj = rj(:,:,[1 2 4 3]); % re-order 45,135 (empirical)

		C = C2sparse('leak', mask, 8, 0, 0);
		i1 = 1:(nx*ny);
		% analysis assumes 1/sqrt(2) for diagonal 1st differences!
		C = [	sptmp(sqrt(col(rj(:,:,1))))	* C(0*nx*ny+i1,:);
			sptmp(sqrt(col(rj(:,:,2))))	* C(1*nx*ny+i1,:);
			sptmp(sqrt(col(rj(:,:,3))/2))	* C(2*nx*ny+i1,:);
			sptmp(sqrt(col(rj(:,:,4))/2))	* C(3*nx*ny+i1,:)];
	else
		error 'bad type'
	end

else
	error(['unknown type: ' type])
end

%
% in:	iprod	[3 np]		inner products of angular weights with
%					[1 cos(2a) sin(2a)]
% out:	rj	[2 np]		horiz,vert
%
function rj = penalty2_design_d1_n1(iprod)

d1 = iprod(1,:);
d2 = iprod(2,:);

if any(d1 < 0), error 'bad inner products: d1<0', end
if any(abs(d2) > d1), error 'bad inner products: d2>d1', end

rj = [d1+2*d2; d1-2*d2];

ii = d2 > d1/2;
r1 = 4/3 * (d1 + d2);
rj(1,ii) = r1(ii);
rj(2,ii) = 0;

ii = d2 < -d1/2;
r2 = 4/3 * (d1 - d2);
rj(1,ii) = 0;
rj(2,ii) = r2(ii);

%
% in:
%	iprod	[3 np]		inner products of angular weights with
%					[1 cos(2a) sin(2a)]
%					typically np = # of pixels
% out:
%	rj	[4 np]		horiz,vert,diag1,diag2 coefficients
%
% caution: diagonals coefficients are designed to be used with
% diagonal differences that are normalized by 1/\sqrt{2}.
%
function rj = penalty2_design_d1_n2(iprod)

d1 = iprod(1,:);
dp2 = iprod(2,:);
dp3 = iprod(3,:);
d2 = max(abs(dp2), abs(dp3)); % d2 >= d3
d3 = min(abs(dp2), abs(dp3));
if any(d1 < 0), error 'bad inner products: d1<0', end
if any(d2 > d1), error 'bad inner products: d2>d1', end

rj = zeros(4,length(d1));

ii = (d2 >= d1/2) & (d3 <= 2/3*d2 - d1/3); % case 1
r1 = 4/3 * (d1 + d2);
rj(1,ii) = r1(ii);

ii = (d3 >= 2/3*d2 - d1/3) & (d3 + d2 >= d1/2); % case 2
r1 = 8/5 * (d1/2 + 3/2*d2 - d3);
r3 = 12/5 * (d3 - 2/3*d2 + 1/3*d1);
rj(1,ii) = r1(ii);
rj(3,ii) = r3(ii);

ii = (d3 + d2 <= d1/2) & (d2 >= d1/4); % case 3
rj(1,ii) = 4*d2(ii);
r3 = d1 - 2*d2 + 2*d3;
r4 = d1 - 2*d2 - 2*d3;
rj(3,ii) = r3(ii);
rj(4,ii) = r4(ii);

ii = (d2 <= d1/4); % case 4
r1 = d1/2 + 2*d2;
r2 = d1/2 - 2*d2;
r3 = d1/2 + 2*d3;
r4 = d1/2 - 2*d3;
rj(1,ii) = r1(ii);
rj(2,ii) = r2(ii);
rj(3,ii) = r3(ii);
rj(4,ii) = r4(ii);

% symmetries
ii = abs(dp3) > abs(dp2);
rj(:,ii) = rj([3 4 1 2],ii);

ii = dp3 < 0;
rj([3 4],ii) = rj([4 3],ii);

ii = dp2 < 0;
rj([1 2],ii) = rj([2 1],ii);

%
% brute-force NNLS approach
% dot	[3 np]
% rj	[4 np]
%
function rj = penalty2_design_nnls(Amat, dot)

ww = [dot(1,:); sqrt(2)*dot([2 3],:)];
opt = optimset('lsqnonneg');
opt = optimset(opt, 'tolx', 10*eps);

rj = zeros(4, ncol(ww));
for jj=1:ncol(ww)
	rj(:,jj) = lsqnonneg(Amat, ww(:,jj));
	if ~rem(jj,100), printf('%d of %d', jj, ncol(ww)), end
end

% penalty2_design_test()
function penalty2_design_test
Amat = 0.5 * [1 1 1 1;
		1/sqrt(2) -1/sqrt(2) 0 0 ;
		0 0 1/sqrt(2) -1/sqrt(2)];
n2 = 41;
n3 = 43;
d2 = linspace(-1,1,n2);
d3 = linspace(-1,1,n3);
dot = zeros(3,n2*n3);
[t2, t3] = ndgrid(d2, d3);
dot(1,:) = 1;
dot(2,:) = t2(:)';
dot(3,:) = t3(:)';

rj_anal = penalty2_design_d1_n2(dot);
rj_nnls = penalty2_design_nnls(Amat, dot); % [4 n2*n3]

penalty2_design_test_figure(Amat, dot, d2, d3, rj_anal, rj_nnls)

function penalty2_design_test_figure(Amat, dot, d2, d3, ra, rn)

n2 = length(d2);
n3 = length(d3);

ww = [dot(1,:); sqrt(2)*dot([2 3],:)];
erra = permute(reshape(Amat * ra - ww, [3 n2 n3]), [2 3 1]); % [n2 n3 3]
errn = permute(reshape(Amat * rn - ww, [3 n2 n3]), [2 3 1]);

ra = permute(reshape(ra, [4 n2 n3]), [2 3 1]); % [n2 n3 4]
rn = permute(reshape(rn, [4 n2 n3]), [2 3 1]);

if im
	clf, pl=440;
	for ib=1:4
		im(pl+ib+0, d2, d3, ra(:,:,ib)), axis xy, cbar
		title(sprintf('Anal. ib=%d', ib))
		im(pl+ib+4, d2, d3, rn(:,:,ib)), axis xy, cbar
		title NNLS
	end

	subplot(4,4,9)
	im(d2, d3, sum(ra > 0,3), 'Anal. sum(r > 0)'), axis xy, cbar
	subplot(4,4,10)
	im(d2, d3, sum(rn > 0,3), 'NNLS sum(r > 0)'), axis xy, cbar

	% hold on, plot(d2, 2/3*(d2 - 1/sqrt(2)), 'c-'), hold off
	% hold on, plot(d2, sqrt(2)/3 - 2/3*d2, 'r-'), hold off

	tmp = sqrt(sum(erra.^2, 3));
	subplot(4, 4, 11)
	im(d2, d3, tmp, 'Anal. |err|'), axis xy, cbar
	subplot(4, 4, 15)
	im(d2, d3, tmp < 10*eps, 'Anal. |err|=0'), axis xy, cbar

	tmp = sqrt(sum(errn.^2, 3));
	subplot(4, 4, 12)
	im(d2, d3, tmp, 'NNLS |err|'), axis xy, cbar
	subplot(4, 4, 16)
	im(d2, d3, tmp < 10*eps, 'NNLS |err|=0'), axis xy, cbar

	subplot(4, 4, 13)
	%im(d2, d3, sum(erra.^2, 3) < sum(errn.^2, 3)+10*eps, 'Anal<NNLS')
	im(d2, d3, sum(errn.^2, 3) - sum(erra.^2, 3), 'NNLS-Anal')
	axis xy, cbar
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


function A = PairMetricLearning(metric_learn_alg, idxa, idxb, matches, X, A0, params)
% X ... data in rows (hopefully)
% idxa ... idx in x of image A of a pair
% idxb ... idx in x of image B of a pair
% matches ... wheter the pair matches or not
% der koe
%

if (~exist('params')),
    params = struct();
end
params = SetDefaultParams(params);

if (~exist('A0')),
    A0 = eye(size(X,2));
end

[l, u] = ComputeDistanceExtremes(X, 5, 95, A0);

C = zeros(length(idxa),4);

C(:,1) = idxa;
C(:,2) = idxb;
C( matches,3)  = 1;
C(~matches,3) = -1;

C( matches,4)  = l;
C(~matches,4) =  u;

try
    A = feval(metric_learn_alg, C, X, A0, params);
catch ME
    disp('Unable to learn mahal matrix');
    le = lasterror;
    disp(le.message);
    A = zeros(d,d);
end


text:
%GSP_DEMO_GRAPH_TV Reconstruction of missing sample on a graph using TV
%
%   In this demo, we try to reconstruct missing sample of a piece-wise
%   smooth signal on a graph. To do so, we will minimize the well-known TV
%   norm defined on the graph.
%
%   For this example, you need the unlocbox. You can download it here:
%   http://unlocbox.sourceforge.net/download
%
%   We express the recovery problem as a convex optimization problem of the
%   following form:
%
%   ..   argmin   ||grad(x)||_1   s. t. ||Mx-b||_2 < epsilon
%
%   .. math:: arg \min_x  \|\nabla(x)\|_1 \text{ s. t. } \|Mx-b\|_2 \leq \epsilon
%
%   Where b represents the known measurements, M is an operator
%   representing the mask and $\epsilon$ is the radius of the l2 ball.
%
%   We set
%
%   * $f_1(x)=||\nabla x ||_1$
%     We define the prox of $f_1$ as:
%
%     .. prox_{f1,gamma} (z) = argmin_{x} 1/2 ||x-z||_2^2  +  gamma  ||grad(z)||_1
%
%     .. math:: prox_{f1,\gamma} (z) = arg \min_{x} \frac{1}{2} \|x-z\|_2^2 +  \gamma \| \nabla z \|_1
%
%   * $f_2$ is the indicator function of the set S define by $||Mx-b||_2 < \epsilon$
%     We define the prox of $f_2$ as
%
%     .. prox_{f2,gamma} (z) = argmin_{x} 1/2 ||x-z||_2^2  +  gamma i_S( x ),
%
%     .. math:: prox_{f2,\gamma} (z) = arg \min_{x} \frac{1}{2} \|x-z\|_2^2   + i_S(x) ,
%
%     with $i_S(x)$ is zero if x is in the set S and infinity otherwise.
%     This previous problem has an identical solution as:
%
%     .. argmin_{z} ||x - z||_2^2   s.t.  ||b - M z||_2 < epsilon
%
%     .. math:: arg \min_{z} \|x - z\|_2^2   \hspace{1cm} such \hspace{0.25cm} that \hspace{1cm} \|Mz-b\|_2 \leq \epsilon
%
%     It is simply a projection on the B2-ball.
%
%   Results
%   -------
%
%   .. figure::
%
%      Original signal on graph
%
%      This figure shows the original signal on graph.
%
%   .. figure::
%
%      Depleted signal on graph
%
%      This figure shows the signal on graph after the application of the
%      mask and addition of noise. Half of the vertices are set to 0.
%
%   .. figure::
%
%      Reconstructed signal on graph usign TV
%
%      This figure shows the reconstructed signal thanks to the algorithm.
%
%   Comparison with Tikhonov regularization
%   ---------------------------------------
%
%   We can also use the Tikhonov regularizer that will promote smoothness.
%   In this case, we solve:
%
%   ..   argmin   ||grad(x)||_2^2   s. t. ||Mx-b||_2 < epsilon
%
%   .. math:: arg \min_x \tau \|\nabla(x)\|_2^2 \text{ s. t. } \|Mx-b\|_2 \leq \epsilon
%
%   The result is presented in the following figure:
%
%   .. figure::
%
%      Reconstructed signal on graph using Tikhonov
%
%      This figure shows the reconstructed signal thanks to the algorithm.
%

% Author: Nathanael Perraudin
% Date: 4th March 2014

%% Initialisation
clear;
close all;
% Load toolbox
init_unlocbox();

%% Important parameters
% size of the graph for the demo
N = 250;
% probability of having a label on a vertex.
p = 0.1;
verbose = 1;    % verbosity level
sigma = 0.0;

%% Create a graph
%  paramgraph.distribute = 1;
% G = gsp_sensor(N, paramgraph);
G = gsp_sensor(N);
G = gsp_community(N);

% for a path with irregular weights and distances
% G = gsp_path(N);
% G = gsp_update_coordinates(G, G.coords + [rand(G.N, 1), zeros(G.N, 1)]);
% G = gsp_update_weights(G, G.W .* exp(-gsp_distanz(G.coords').^2));

G = gsp_adj2vec(G);
G = gsp_estimate_lmax(G);
G = gsp_compute_fourier_basis(G);

% x_0 = (1 + sign(G.U(:,4))) / 2;
x_0 = round((1 + sign(G.U(:,2))) * 2+(1 + sign(G.U(:,3))) +(1 + sign(G.U(:,4))) / 2);
% x_0 = round(linspace(0, 5, N))';
x_0 = G.info.node_com;

%% Set the labels
% create the mask
% ind_obs = boolean(full(sparse(1:4:G.N, 1, 1, G.N, 1)));
ind_obs = rand(G.N, 1) < p;
ind_unobs = not(ind_obs);

%applying the Mask to the data
x_obs = ind_obs.*(x_0+sigma*randn(G.N,1));

%% Method 1: Graph TV
% setting different parameter for the simulation
param_solver.verbose = verbose;  % display parameter
param_solver.tol = 1e-5;
param_solver.maxit = 4000;

tic;
sol_tv = gsp_regression_tv(G, ind_obs, x_obs, 0, param_solver);
toc

%% Method 2: Tikhonov
tic;
G2 = G;
G2.L = G.L^3;
%sol_tik = gsp_regression_tik(G, ind_obs, x_obs, 0, param_solver);
sol_tik = gsp_regression_tik(G2, ind_obs, x_obs, 0, param_solver);
toc

%% Method 3: My idea: inpaint using graph TV on one variable that is l-2
%% away from another one that is smooth
%%
% x, y are primal variables: z is the dual
%      min_x,y,z ||grad(x_all)||_1 + a/2||x-y||^2 + b/2 ||grad(y_all)||^2
%
% grad(x_all) = grad([x;x_obs]) = A*x + A_obs * x_obs = A*x - b
%
% only rows of gradient corresponding to unknown values used:
%      min_x,y,z ||Ax - b||_1  + a/2||x-y||^2 + b/2 ||Ay - b||^2
% Get rid of operator A by using dual variable:
%      min_x,y,z ||z - b||_1   + a/2||x-y||^2 + b/2 ||Ay - b||^2
% Think of all primal variables as one vector w = [x;y] containing ONLY
% the elements on the unlabeled nodes

% choose parameters for solving the problem:
alpha = .1;
beta = 1;

% compute the Diff matrix of the graph
G = gsp_adj2vec(G);

% operators based on the Diff
A_obs = G.Diff(:, ind_obs);
A_unobs = G.Diff(:, not(ind_obs));
n_unobs = nnz(not(ind_obs));
zeros_x = zeros(n_unobs, 1);
zeros_y = zeros_x;

% graph TV term is the l-1 on the dual:
b = - A_obs * x_obs(ind_obs);
f1_params.y = b;
f1.eval = @(z) norm(z - b, 1);
f1.prox = @(z, gamma) prox_l1(z, gamma, f1_params);
f1.L = @(w) A_unobs * w(1:n_unobs);
f1.Lt = @(z) [A_unobs' * z; zeros_y];
f1.norm_L = normest(A_unobs);
%f1.norm_L = sqrt(G.lmax);

% use prox for the second function
% w = [x;y], w0 = [x0;y0]
% prox_g(x) = ((1+g)x + gy)/(1+2g)
f2.prox = @(w, gamma) [ ((1+alpha*gamma) * w(1:n_unobs) + ...
    alpha*gamma     * w(n_unobs+1:end))/(1+2*alpha*gamma);...
    (alpha*gamma     * w(1:n_unobs) + ...
    (1+alpha*gamma) * w(n_unobs+1:end))/(1+2*alpha*gamma)];
f2.eval = @(w) alpha/2 * norm(w(1:n_unobs) - w(n_unobs+1:end))^2;

% use gradient for the last function
AtA = A_unobs' * A_unobs;
Atb = A_unobs' * b;
f3.grad = @(w) [zeros_x; beta * (AtA*w(n_unobs+1:end) - Atb)];
f3.eval = @(w) beta/2 * norm(A_unobs * w(n_unobs+1:end) - b)^2;
f3.beta = beta * normest(AtA);

%% solve the problem
param_solver.verbose = 1;
param_solver.algo = 'FBF_PRIMAL_DUAL';
param_solver.normalized_timestep = 0.99;
tic;
[sol_mine, info_sol] = solvep([sol_tv(not(ind_obs));  sol_tik(not(ind_obs))], {f1,f2,f3}, param_solver);
toc

sol_mine_x = x_obs;
sol_mine_x(ind_unobs) = sol_mine(1:n_unobs);
sol_mine_y = x_obs;
sol_mine_y(ind_unobs) = sol_mine(n_unobs+1:end);

%% Method 4: L-1 of Lx ||Lx||_1
tic;
param_solver = rmfield(param_solver, 'algo');
sol_Lx_l1 = gsp_regression_Lx_l1(G, ind_obs, x_obs, 0, param_solver);
toc

%% Compute the errors
ans_tv = round(sol_tv);
ans_tik = round(sol_tik);
ans_mine = round(sol_mine_x);
ans_Lxl1 = round(sol_Lx_l1);

err_tv      = nnz(ans_tv ~= x_0) / G.N
err_tik     = nnz(ans_tik ~= x_0) / G.N
err_mine    = nnz(ans_mine ~= x_0) / G.N
err_Lxl1    = nnz(ans_Lxl1 ~= x_0) / G.N

%% Print the result
paramplot.show_edges = 1;
val_lims = lin_map([-.05, 1.05], [min([x_0; sol_tv; sol_tik]), max([x_0; sol_tv; sol_tik])], [0, 1]);

%% Plot the original graph
if not(strcmp(G.type, 'path'))
    figure(1)
    paramplot.vertex_highlight = ind_obs;
    gsp_plot_signal(G, x_0, paramplot)
    caxis(val_lims);
    title('Original signal: highlighted known values')

    % % Let show depleted graph
    % figure(2)
    % gsp_plot_signal(G,depleted_graph_value,paramplot)
    % caxis([-1 1])
    % title('Measurement')
    % Let show the reconstructed graph
    figure(3)
    paramplot.vertex_highlight = (ans_tv ~= x_0);
    gsp_plot_signal(G, ans_tv, paramplot)
    caxis(val_lims)
    title('TV solution: highlighted mistakes')

    % Let show the reconstructed graph
    figure(4)
    paramplot.vertex_highlight = (ans_tik ~= x_0);
    gsp_plot_signal(G, ans_tik, paramplot)
    caxis(val_lims)
    title('Tikhonov solution: highlighted mistakes')

    % Let show the reconstructed graph
    figure(5)
    paramplot.vertex_highlight = (ans_mine ~= x_0);
    gsp_plot_signal(G, ans_mine, paramplot)
    caxis(val_lims)
    title('My solution: highlighted mistakes')

    % Let show the reconstructed graph
    figure(6)
    paramplot.vertex_highlight = (ans_Lxl1 ~= x_0);
    gsp_plot_signal(G, ans_Lxl1, paramplot)
    caxis(val_lims)
    title('||Lx||_1 solution: highlighted mistakes')

    %%
else strcmp(G.type, 'path')
    figure;
    plot(G.coords(:, 1), sol_tv)
    hold on;
    plot(G.coords(:, 1), sol_tik)
    plot(G.coords(:, 1), x_0);
    plot(G.coords(:, 1), sol_Lx_l1)
    plot(G.coords(ind_obs, 1), x_0(ind_obs), 'o');
    legend('TV', 'Tikhonov', '||Lx||_1', 'true');
end

figure; plot(sol_mine_x)
hold on; plot(sol_mine_y)

repo: gspbox


text:
function fn = tapas_physio_get_filename_from_id_philips(fileID, dirRaw, ext)
%searches for matching filename with fileID in dirRaw
%
%   output = get_filename_from_id(input)
%
% IN
%   fileID  e.g. 8
%   dirRaw  folder with raw files
%   ext     {'.nii'} file extension
% OUT
%
% EXAMPLE
%   get_filename_from_id
%
%   See also

% Author: Lars Kasper
% Created: 2013-11-04
% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich.

if nargin < 3
    ext = '.nii';
end

try
    fn = regexprep(regexprep(ls(fullfile(dirRaw ...
        ,['*_' num2str(fileID) '_1_*' ext])), [dirRaw filesep], ''),'\n','');
catch
    fn = regexprep(regexprep(ls(fullfile(dirRaw ...
        ,['*_' num2str(fileID) '_2_*' ext])), [dirRaw filesep], ''),'\n','');
end

author: translationalneuromodeling

repo: tapas


text:
function w = ncoh_compute_weights ( n )

%*****************************************************************************80
%
%% NCOH_COMPUTE_WEIGHTS computes weights for a Newton-Cotes "open half" quadrature rule.
%
%  Discussion:
%
%    The input value N is used to define N equal subintervals of [-1,+1].
%    The I-th abscissa is the center of the I-th subinterval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real W(N), the weights.
%
  x = ncoh_compute_points ( n );

  w = nc_compute ( n, x_min, x_max, x );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function write_surf(fname, vert, face)

% write_surf - FreeSurfer I/O function to write a surface file
%
% write_surf(fname, vert, face)
%
% writes a surface triangulation into a binary file
% fname - name of file to write
% vert  - Nx3 matrix of vertex coordinates
% face  - Mx3 matrix of face triangulation indices
%
% The face matrix here must be matlab compatible
% (no zero indices).  It is converted to FreeSurfer
% indices that start at zero.
%
% See also freesurfer_read_surf, freesurfer_write_curv, freesurfer_write_wfile

if(nargin ~= 3)
  fprintf('USAGE: freesurfer_write_surf(fname, vert, face)\n');
  return;
end

if size(vert,2) ~= 3,
    error('vert must be Nx3 matrix');
end

if size(face,2) ~= 3,
    error('face must be Mx3 matrix');
end

fprintf('...subtracting 1 from face indices for FreeSurfer compatibility.\n');
face = face - 1;

% open it as a big-endian file
fid = fopen(fname, 'wb', 'b') ;

TRIANGLE_FILE_MAGIC_NUMBER = 16777214 ;
fwrite3(fid, TRIANGLE_FILE_MAGIC_NUMBER);

vnum = size(vert,1) ;  % number of vertices
fnum = size(face,1) ;  % number of faces

% Ouput a couple of text lines with creation date
fprintf(fid,'created by %s on %s\n\n',getenv('USER'),datestr(now)); % creation date

fwrite(fid, vnum,'int32');
fwrite(fid, fnum,'int32');

% reshape vert into column array and write
vert = reshape(vert',size(vert,1)*size(vert,2),1);
fwrite(fid, vert,'float32');

% reshape face into column array and write
face = reshape(face',size(face,1)*size(face,2),1);
fwrite(fid, face,'int32');

fclose(fid) ;

return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function lvmClassVisualise(call)

% LVMCLASSVISUALISE Callback function for visualising data.
% FORMAT
% DESC contains the callback functions for visualizing points from the
% latent space in the higher dimension space.
% ARG call : either 'click', 'move', 'toggleDynamics',
% 'dynamicsSliderChange'
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2009
%
% SEEALSO : lvmResultsDynamic

% MLTOOLS

global visualiseInfo

switch call
 case 'click'
  [x, y]  = localCheckPointPosition(visualiseInfo);
  if ~isempty(x)
    visualiseInfo.latentPos(visualiseInfo.dim1) = x;
    visualiseInfo.latentPos(visualiseInfo.dim2) = y;
  end
  visualiseInfo.clicked = ~visualiseInfo.clicked;
  if isfield(visualiseInfo.model, 'dynamics') & ~isempty(visualiseInfo.model.dynamics)
    if visualiseInfo.runDynamics
      visualiseInfo.dynamicsRunning = 1;
      fhandle = str2func([visualiseInfo.model.type 'DynamicsRun']);
      feval(fhandle);
      visualiseInfo.dynamicsRunning = 0;
    end
  else
    visualiseInfo.dynamicsRunning = 0;
  end
 case 'move'
  if visualiseInfo.clicked & ~visualiseInfo.runDynamics
    [x, y]  = localCheckPointPosition(visualiseInfo);
    if ~isempty(x)
      % This should be changed to a model specific visualisation.
      visualiseInfo.latentPos(visualiseInfo.dim1) = x;
      visualiseInfo.latentPos(visualiseInfo.dim2) = y;
      set(visualiseInfo.latentHandle, 'xdata', x, 'ydata', y);
      fhandle = str2func([visualiseInfo.model.type 'PosteriorMeanVar']);
      [mu, varsigma] = fhandle(visualiseInfo.model, visualiseInfo.latentPos);
      if isfield(visualiseInfo.model, 'noise')
        Y = noiseOut(visualiseInfo.model.noise, mu, varsigma);
      else
        Y = mu;
      end
      visualiseInfo.visualiseModify(visualiseInfo.visHandle, ...
                                    Y, visualiseInfo.varargin{:});
    end
  end
 case 'toggleDynamics'
  visualiseInfo.runDynamics = ~visualiseInfo.runDynamics;
  set(visualiseInfo.dynamicsRadio, 'value', visualiseInfo.runDynamics);

 case 'dynamicsSliderChange'
  X = modelOut(visualiseInfo.model.dynamics, get(visualiseInfo.dynamicsSlider, 'value'));
  x = X(1);
  y = X(2);
  visualiseInfo.latentPos(visualiseInfo.dim1) = x;
  visualiseInfo.latentPos(visualiseInfo.dim2) = y;
  set(visualiseInfo.latentHandle, 'xdata', x, 'ydata', y);
  fhandle = str2func([visualiseInfo.model.type 'PosteriorMeanVar']);
  [mu, varsigma] = fhandle(visualiseInfo.model, visualiseInfo.latentPos);
  if isfield(visualiseInfo.model, 'noise')
    Y = noiseOut(visualiseInfo.model.noise, mu, varsigma);
  else
    Y = mu;
  end
  visualiseInfo.visualiseModify(visualiseInfo.visHandle, ...
                                Y, visualiseInfo.varargin{:});

 case 'latentSliderChange'
  counter = 0;
  for i = size(visualiseInfo.latentPos, 2):-1:1
    % Go through setting latent positions to slider positions.
    if i ~= visualiseInfo.dim1 && i ~= visualiseInfo.dim2
      counter = counter + 1;
      visualiseInfo.latentPos(i) = get(visualiseInfo.latentSlider(counter), 'value');
      set(visualiseInfo.sliderTextVal(counter), 'string', num2str(visualiseInfo.latentPos(i)));
    end
  end
  lvmSetPlot;
  fhandle = str2func([visualiseInfo.model.type 'PosteriorMeanVar']);
  [mu, varsigma] = fhandle(visualiseInfo.model, visualiseInfo.latentPos);
  if isfield(visualiseInfo.model, 'noise')
    Y = noiseOut(visualiseInfo.model.noise, mu, varsigma);
  else
    Y = mu;
  end
  visualiseInfo.visualiseModify(visualiseInfo.visHandle, ...
                                    Y, visualiseInfo.varargin{:});
  if verLessThan('matlab', 'R2014a')
    visualiseInfo.latentHandle = line(visualiseInfo.latentPos(visualiseInfo.dim1), ...
                                      visualiseInfo.latentPos(visualiseInfo.dim2), ...
                                      'markersize', 20, 'color', [0.5 0.5 0.5], ...
                                      'marker', '.', 'visible', 'on', ...
                                      'erasemode', 'xor');
  else
    visualiseInfo.latentHandle = line(visualiseInfo.latentPos(visualiseInfo.dim1), ...
                                      visualiseInfo.latentPos(visualiseInfo.dim2), ...
                                      'markersize', 20, 'color', [0.5 0.5 0.5], ...
                                      'marker', '.', 'visible', ...
                                      'on');
  end

 case 'updateLatentRepresentation'
  visualiseInfo.dim1 = get(visualiseInfo.xDimension, 'value');
  visualiseInfo.dim2 = get(visualiseInfo.yDimension, 'value');
  lvmSetPlot;
  if verLessThan('matlab', 'R2014a')
    visualiseInfo.latentHandle = line(visualiseInfo.latentPos(visualiseInfo.dim1), ...
                                      visualiseInfo.latentPos(visualiseInfo.dim2), ...
                                      'markersize', 20, 'color', [0.5 0.5 0.5], ...
                                      'marker', '.', 'visible', 'on', ...
                                      'erasemode', 'xor');
  else
    visualiseInfo.latentHandle = line(visualiseInfo.latentPos(visualiseInfo.dim1), ...
                                      visualiseInfo.latentPos(visualiseInfo.dim2), ...
                                      'markersize', 20, 'color', [0.5 0.5 0.5], ...
                                      'marker', '.', 'visible', 'on');
  end

end

function point = localGetNormCursorPoint(figHandle)

point = get(figHandle, 'currentPoint');
figPos = get(figHandle, 'Position');
% Normalise the point of the curstor
point(1) = point(1)/figPos(3);
point(2) = point(2)/figPos(4);

function [x, y] = localGetNormAxesPoint(point, axesHandle)

position = get(axesHandle, 'Position');
x = (point(1) - position(1))/position(3);
y = (point(2) - position(2))/position(4);
lim = get(axesHandle, 'XLim');
x = x*(lim(2) - lim(1));
x = x + lim(1);
lim = get(axesHandle, 'YLim');
y = y*(lim(2) - lim(1));
y = y + lim(1);

function [x, y] = localCheckPointPosition(visualiseInfo)

% Get the point of the cursor
point = localGetNormCursorPoint(gcf);

% get the position of the axes
position = get(visualiseInfo.plotAxes, 'Position');

% Check if the pointer is in the axes
if point(1) > position(1) ...
      & point(1) < position(1) + position(3) ...
      & point(2) > position(2) ...
      & point(2) < position(2) + position(4);

  % Rescale the point according to the axes
  [x y] = localGetNormAxesPoint(point, visualiseInfo.plotAxes);

  % Find the nearest point
else
  % Return nothing
  x = [];
  y = [];
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function h = gp_plot(gp, x, y, varargin)
%GP_PLOT  Make a plot with Gaussian process
%
%  Description
%    GP_PLOT(GP, X, Y, XT, OPTIONS)
%    takes a GP structure together with matrix X of training
%    inputs and vector Y of training targets and plots predictions
%    evaluated at test inputs XT.
%
%    [EF, VARF, LPY, EY, VARY] = GP_PLOT(GP, X, Y, OPTIONS)
%    plots predictions evaluated at training inputs X.
%
%    The form of the plot depends on the dimensionality of X and options.
%      - 1D with Gaussian: mean and 5% and 95% quantiles of f and y
%      - 1D with non-Gaussian: median and 5% and 95% quantiles of mu or f
%      - 2D: Conditional predictions, see GP_CPRED + 2D plot of mu or f
%      - ND with N>2: Conditional predictions, see GP_CPRED
%
%    OPTIONS is optional parameter-value pair
%      target - option for choosing what is computed 'mu' (default)
%               or 'f'
%      normdata - a structure with fields xmean, xstd, ymean, and ystd
%               to allow plotting in the original data scale (see
%               functions normdata and denormdata)
%      xlabels - a cell array of covariate label strings
%      predcf - an index vector telling which covariance functions are
%               used for prediction. Default is all (1:gpcfn).
%               See additional information below.
%      tstind - a vector/cell array defining, which rows of X belong
%               to which training block in *IC type sparse models.
%               Default is []. In case of PIC, a cell array
%               containing index vectors specifying the blocking
%               structure for test data. In FIC and CS+FIC a
%               vector of length n that points out the test inputs
%               that are also in the training set (if none, set
%               TSTIND = []).
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, the expected
%               value for the ith case.
%      fcorr  - Method used for latent marginal posterior corrections.
%               Default is 'off'. Possible methods are 'fact' for EP
%               and either 'fact' or 'cm2' for Laplace. If method is
%               'on', 'fact' is used for EP and 'cm2' for Laplace.
%
%    NOTE! In case of FIC and PIC sparse approximation the
%    prediction for only some PREDCF covariance functions is just
%    an approximation since the covariance functions are coupled in
%    the approximation and are not strictly speaking additive
%    anymore.
%
%    For example, if you use covariance such as K = K1 + K2 your
%    predictions Ef1 = gp_plot(GP, X, Y, X, 'predcf', 1) and Ef2 =
%    gp_plot(gp, x, y, x, 'predcf', 2) should sum up to Ef =
%    gp_plot(gp, x, y, x). That is Ef = Ef1 + Ef2. With FULL model
%    this is true but with FIC and PIC this is true only
%    approximately. That is Ef \approx Ef1 + Ef2.
%
%    With CS+FIC the predictions are exact if the PREDCF covariance
%    functions are all in the FIC part or if they are CS
%    covariances.
%
%    NOTE! When making predictions with a subset of covariance
%    functions with FIC approximation the predictive variance can
%    in some cases be ill-behaved i.e. negative or
%    unrealistically small. This may happen because of the
%    approximative nature of the prediction.
%
%  See also
%    GP_SET, GP_OPTIM, DEMO_REGRESSION*
%
% Copyright (c) 2014 Aki Vehtari

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

ip=inputParser;
ip.FunctionName = 'GP_PLOT';
ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('xlabels', [], @(x) isempty(x) || iscell(x));
ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                 isvector(x) && isreal(x) && all(isfinite(x)&x>=0))
ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                 (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
ip.addParamValue('fcorr', 'off', @(x) ismember(x, {'off', 'fact', 'cm2', 'on'}));
ip.addParamValue('tr', 0.25, @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('target', 'mu', @(x) ismember(x,{'f','mu'}))
ip.addParamValue('normdata', struct(), @(x) isempty(x) || isstruct(x))
if numel(varargin)==0 || isnumeric(varargin{1})
  % inputParser should handle this, but it doesn't
  ip.parse(gp, x, y, varargin{:});
else
  ip.parse(gp, x, y, [], varargin{:});
end
xt=ip.Results.xt;
zt=ip.Results.zt;
options=struct();
options.predcf=ip.Results.predcf;
options.tstind=ip.Results.tstind;
z=ip.Results.z;
if ~isempty(z)
  options.z=z;
end
if ~isempty(zt)
  options.zt=zt;
end
if isempty(zt)
  options.zt=z;
end
if isempty(xt)
  if size(x,2)==1
    [xt,xi]=sort(x);
    if ~isempty(z)
      zt=z(xi);
    else
      zt=[];
    end
  else
    xt=x;
    zt=z;
  end
end
target = ip.Results.target;
if iscell(gp)
  liktype=gp{1}.lik.type;
else
  liktype=gp.lik.type;
end
if isequal(liktype, 'Coxph') && isequal(target,'mu')
    target='f';
    warning('GP_CPRED: Target ''mu'' not applicable for a Cox-PH model. Switching to target ''f''')
end
tr = ip.Results.tr;
xlabels=ip.Results.xlabels;
% normdata
nd=ip.Results.normdata;
ipnd=inputParser;
ipnd.FunctionName = 'normdata';
ipnd.addParamValue('xmean',0,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('xstd',1,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('xlog',0,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('ymean',0,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('ystd',1,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.addParamValue('ylog',0,@(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))));
ipnd.parse(nd);
nd=ipnd.Results;

if iscell(gp)
  liktype=gp{1}.lik.type;
else
  liktype=gp.lik.type;
end

[n,m]=size(x);

switch m
  case 1
    if isequal(liktype, 'Gaussian')
      [Ef, Varf,~,Ey,Vary] = gp_pred(gp, x, y, xt, options);
      Stdf=sqrt(Varf);
      Stdy=sqrt(Vary);
      xt=denormdata(xt,nd.xmean,nd.xstd);
      Ef=denormdata(Ef,nd.ymean,nd.ystd);
      Stdf=Stdf*nd.ystd;
      Ey=denormdata(Ey,nd.ymean,nd.ystd);
      Stdy=Stdy*nd.ystd;
      hh=plot(xt, Ef, '-b', xt, Ef-1.64*Stdf, '--b', xt, Ef+1.64*Stdf, '--b',xt, Ey-1.64*Stdy, ':b', xt, Ey+1.64*Stdy, ':b');
    else
      switch target
        case 'f'
          [Ef, Varf] = gp_pred(gp, x, y, xt, options);
          Stdf=sqrt(Varf);
          xt=denormdata(xt,nd.xmean,nd.xstd);
          Ef=denormdata(Ef,nd.ymean,nd.ystd);
          Stdf=Stdf*nd.ystd;
          hh=plot(xt, Ef, '-b', xt, Ef-1.64*Stdf, '--b', xt, Ef+1.64*Stdf, '--b');
        case 'mu'
          prctmu = gp_predprctmu(gp, x, y, xt, options);
          xt=denormdata(xt,nd.xmean,nd.xstd);
          prctmu=denormdata(prctmu,nd.ymean,nd.ystd);
          Ef = prctmu; Varf = [];
          hh=plot(xt, prctmu(:,2), '-b', xt, prctmu(:,1), '--b', xt, prctmu(:,3), '--b');
      end
    end
    if ~isempty(xlabels)
        xlabel(xlabels{1})
    else
        xlabel('x')
    end
  case 2
    subplot(2,2,1);
    gp_cpred(gp,x,y,xt,1,'z',z,'zt',zt,'target',target,'plot','on',varargin{:});
    if ~isempty(xlabels)
        xlabel(xlabels{1})
    else
        xlabel('x1')
    end
    subplot(2,2,2);
    gp_cpred(gp,x,y,xt,2,'z',z,'zt',zt,'target',target,'plot','on',varargin{:});
    if ~isempty(xlabels)
        xlabel(xlabels{2})
    else
        xlabel('x2')
    end
    subplot(2,1,2);
    gp_cpred(gp,x,y,xt,[1 2],'z',z,'zt',zt,'target',target,'plot','on','tr',1e9,varargin{:});
    axis square
    if ~isempty(xlabels)
        xlabel(xlabels{1})
        ylabel(xlabels{2})
    else
        xlabel('x1')
        ylabel('x2')
    end
    view(3)
    shading faceted
    hhh=get(gca,'children');
    hh=hhh(1);
  otherwise
    sn1=ceil(sqrt(m));
    sn2=ceil(m/sn1);
    for xi=1:m
      subplot(sn1,sn2,xi);
      gp_cpred(gp,x,y,xt,xi,'z',z,'zt',zt,'target',target,'plot','on','normdata',nd,options);
      if ~isempty(xlabels)
          xlabel(xlabels{xi})
      else
          xlabel(sprintf('x%d',xi))
      end
      drawnow
      hhh=get(gca,'children');
      hh(xi)=hhh(1);
    end
end
if nargout>0
  h=hh;
end

author: gpstuff-dev

repo: gpstuff

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


text:
function [Decoded,Decoder_Chip,Temp_Decoded,intgrl] =  ...
  CDMA_decode(OutSignal,Chipbit,User_to_Decode, ...
  SamplesPerBit,SamplePerChip,TotalDataBit)

% ........................ CDMA Decoding Starts Here .....................
clipto = TotalDataBit*SamplesPerBit;

TotalChips = length(Chipbit(:,User_to_Decode));
Decoder_Chip1 =[];
for l1 = 1:TotalDataBit
    Decoder_Chip1 = vertcat(Decoder_Chip1,Chipbit(:,User_to_Decode));
end
Decoder_Chip = MakeSampled(Decoder_Chip1,length(OutSignal),SamplePerChip);
Temp_Decoded = OutSignal.*Decoder_Chip;
Temp_Decoded = Temp_Decoded(1:TotalDataBit*SamplesPerBit);

intgrl = zeros(length(Temp_Decoded),1);
Avg_Ingl = 0;
Decoded = zeros(TotalDataBit,1);
    bitno = 0;
    for l1 = 1:TotalChips*SamplePerChip:length(Temp_Decoded)
        bitno = bitno+1;
        for l2 = l1:(l1+TotalChips*SamplePerChip-1)
            if (l2==1),intgrl(1)=Temp_Decoded(1);continue;end
            intgrl(l2) = intgrl(l2-1) + Temp_Decoded(l2);
        end
        if (intgrl(l2)>0),
            Decoded(bitno)=1;
        else
            Decoded(bitno)=-1;
        end
        intgrl(l2) = 0;
    end
intgrl=vertcat(intgrl,zeros(length(OutSignal)-length(intgrl),1));
Temp_Decoded = vertcat(Temp_Decoded,zeros(length(OutSignal)-length(Temp_Decoded),1));
% ........................ CDMA Decoding Ends Here .....................

author: Sable

repo: mcbench-benchmarks


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym hilb (@var{n})
%% Return the symbolic Hilbert matrix.
%%
%% Example:
%% @example
%% @group
%% hilb (sym(2))
%%   @result{} ans = (sym 2×2 matrix)
%%       ⎡ 1   1/2⎤
%%       ⎢        ⎥
%%       ⎣1/2  1/3⎦
%% @end group
%% @end example
%%
%% @seealso{@@sym/invhilb}
%% @end defmethod

function y = hilb(x)
  if (nargin ~= 1)
    print_usage ();
  end

  y = pycall_sympy__ ('return Matrix(_ins[0], _ins[0], lambda i,j: 1 / (i + j + 1)),', x);

end

%!test
%! A = hilb (sym(3));
%! B = [sym(1) sym(1)/2 sym(1)/3; sym(1)/2 sym(1)/3 sym(1)/4; sym(1)/3 sym(1)/4 sym(1)/5];
%! assert (isequal (A, B))

repo: octsympy


text:
function excelrecord5fcn(stringdates2, endo, plotdata, n, pref)
% script excelrecord5
% records the information contained in the worksheet 'forecasts' of the excel spreadsheet 'results.xls'

% create the cell that will be saved on excel
forecastcell={};

% build preliminary elements: space between the tables
vertspace=repmat({''},size(stringdates2,1)+3,1);

% loop over variables (horizontal dimension)
for ii=1:n
    % create cell of forecast record for variable ii
    temp=['forecasts: ' endo{ii,1}];
    forecast_i=[temp {''} {''} {''} {''};{''} {''} {''} {''} {''};{''} {'sample'} {'lw. bound'} {'median'} {'up. bound'};stringdates2 num2cell((plotdata{ii,1})')];

    % switch NaN entries to empty
    forecast_i(~cellfun(@any,forecast_i))={[]};
    forecastcell=[forecastcell forecast_i vertspace];
end

% trim
forecastcell=forecastcell(:,1:end-1);

% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),forecastcell,'forecasts','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [res,I_patch,I]=getCNNFeature(obj,theNet,theConf,IsFlip,meanData)
theNet.layers{end}.class=theNet.layers{end}.class(:,:,:,1);
theNet.layers{end}.class(:)=1;
[I_patch,I]=getI(obj,theConf,theNet.meta.normalization.imageSize(1:2),IsFlip);
im_=I_patch-meanData;
if(isa(theNet.layers{1}.weights{1},'gpuArray'))
    res=our_vl_simplenn(theNet,gpuArray(im_));
    for i=1:length(res)
        res(i).x=gather(res(i).x);
    end
else
    res=our_vl_simplenn(theNet,im_);
end
I_patch=uint8(I_patch);
end

repo: interpretableCNN


text:
% VL_SIMDCTRL  Toggle VLFeat SIMD optimizations
%   VL_SIMDCTRL() returns TRUE if VLFeat SIMD optimizations are enabled
%   and FALSE otherwise.
%
%   VL_SIMDCTRL(TRUE) enables VLFeat SIMD optimizations. VL_SIMDCTRL(FALSE)
%   disables them.
%
%   See also: VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
% Test file for @chebfun/subspace.m.

function pass = test_subspace(pref)

if ( nargin == 0 )
   pref = chebfunpref();
end

% Orthonormal array-valued chebfun:
A = chebfun(@(t) [1/sqrt(2)+0*t cos(t) sin(2*t) sin(3*t)]/sqrt(pi), ...
    [0 1 2*pi], pref);

f = chebfun(@(t) sin(10*t)/sqrt(pi), [0 2 2*pi], pref);
alpha = [1e-10 pi/5 pi/2-1e-10];

for k = 1:length(alpha)
    B = cos(alpha(k))*A(:,k) + sin(alpha(k))*f;
    angle = subspace(A, B);
    pass(k) = abs(angle - alpha(k)) < 1e3*eps;
    angle = subspace(B, A);
    pass(k) = pass(k) && (abs(angle - alpha(k)) < 1e3*eps);
end

% Check subspaces with multiple columns.
pass(4) = abs(subspace(A(:,1:2), A(:,3:4)) - pi/2) < 10*vscale(A)*eps;
pass(5) = abs(subspace(A(:,1:3), A(:,4)) - pi/2) < 10*vscale(A)*eps;

% Check operation for row chebfuns.
At = A.';
pass(6) = abs(subspace(At(1:2,:), At(3:4,:)) - pi/2) < ...
    10*vscale(At)*eps;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
        % Ciclu repetitiv pentru calcularea coeficientilor binomiali
for k=1:10
    binom(k,1:k+1)=mfun('binomial',k,0:k);
end
        % Afisarea rezultatului
disp('Coeficientii binomiali')
disp(binom)


text:
        % Lansarea utilitarului Profiler cu optiunea cea mai completa
profile on -detail operator -history
        % Secventa simpla de program de analizat
x=0:pi/1e6:2*pi;
y=funct1(x);
z=funct2(x);
        % Comandarea generarii unui raport si salvarea lui in fisierul raport
profile report raport
        % Comandarea generarii reprezentarii rezultatelor statistice
profile plot
        % Stergerea tuturor inregistrarilor
profile clear


text:
function mtrFigureDOCCCorrelations(machineDir, saveDir)

subjVec = {'ss040804','mho040625','bg040719','md040714'};
threshVec = [1000];
paramNames = {'kLength','kSmooth','kMidSD'};
midP = [0 18 0.175];
bSaveImages = 1;

for pp = 1:3
    figure;
    for ss = 1:length(subjVec)
        subjDir = [machineDir subjVec{ss}];

        fgDir = [subjDir '/conTrack/resamp_LDOCC'];
        disp(['cd ' fgDir]);
        cd(fgDir);
        % Correlation between CTs
        if(pp==1)
            mtrPlotCorr('cc_-2',threshVec,paramNames,midP,pp,'b'); hold on;
        else
            mtrPlotCorr('cc',threshVec,paramNames,midP,pp,'b'); hold on;
        end
        % Correlation between CTs and STTs
        %mtrPlotCorr('cc_STT',threshVec,paramNames,midP,1,'g'); hold on;
        % Save out the figures to the image dir

        fgDir = [subjDir '/conTrack/resamp_RDOCC'];
        disp(['cd ' fgDir]);
        cd(fgDir);
        % Correlation between CTs
        if(pp==1)
            mtrPlotCorr('cc_-2',threshVec,paramNames,midP,pp,'b'); hold on;
        else
            mtrPlotCorr('cc',threshVec,paramNames,midP,pp,'b'); hold on;
        end
        oldaxis = axis;
        if pp == 1
            axis([-4 1 0 1]);
        else
            axis([oldaxis(1) oldaxis(2) 0 1]);
        end
        % Correlation between CTs and STTs
        %mtrPlotCorr('cc_STT',threshVec,paramNames,midP,1,'g');
    end
    % Save out the figures to the image dir
    set(gcf,'Position',[395   447   289   240]);
    if bSaveImages
        figFilename = fullfile(saveDir,['param_' paramNames{pp} '_corr.png']);
        set(gcf,'PaperPositionMode','auto');
        print('-dpng', figFilename);
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [g1L, g2L, g3L, g4L] = sdlfmvXsdlfmKernGradientICBlock(lfmKern1, ...
    lfmKern2, tInit1, tInit2, kyy, kyv, kvy, kvv, i, j, generalConst, ...
    generalConstGrad, gkyy1, gkyy2, gkyy3, gkyy4, gkyv1, gkyv2, ...
    gkyv3, gkyv4, gkvy1, gkvy2, gkvy3, gkvy4, gkvv1, gkvv2, gkvv3, gkvv4)

% SDLFMVXSDLFMKERNGRADIENTICBLOCK Partial derivatives initial conditions
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

g1L = cell(2);
g2L = cell(2);
g3L = cell(2);
g4L = cell(2);

typeParam = {'sdlfmv', 'sdlfm'};

% Compute for tInit1(1) and tInit2(2)

[g1L{1,2}, g2L{1,2}, g3L{1,2}] = sdlfmvXsdlfmKernGradientBlock(lfmKern1, lfmKern2, tInit1(1), ...
    tInit2(2), kyy, kyv, kvy, kvv, i, j, generalConst, generalConstGrad, 1);

g1Plus = sdlfmaXsdlfmKernComputeBlock(lfmKern1, lfmKern2, tInit1(1), ...
    tInit2(2), kyy, kyv, kvy, kvv, i, j, generalConst);

g2Plus = sdlfmvXsdlfmvKernComputeBlock(lfmKern1, lfmKern2, tInit1(1), ...
    tInit2(2), kyy, kyv, kvy, kvv, i, j, generalConst);

if (i==1 && j==1) || (i==1) && (j~=1)

    if i~=1 || j~=1
        [gkyy1IC, gkyy2IC, gkyy3IC, gkyy4IC] = sdlfmXsdlfmKernGradientIC(lfmKern1, lfmKern2, ...
            tInit1(1), tInit2(2), gkyy1{i,j}, gkyy2{i,j}, gkyy3{i,j}, gkyy4{i,j}, gkyv1{i,j}, gkyv2{i,j}, ...
            gkyv3{i,j}, gkyv4{i,j}, gkvy1{i,j}, gkvy2{i,j}, gkvy3{i,j}, gkvy4{i,j}, gkvv1{i,j}, gkvv2{i,j}, ...
            gkvv3{i,j}, gkvv4{i,j}, 1, typeParam);
        for k=1:3
            g1L{1,2}{k} = g1L{1,2}{k} + gkyy1IC{k};
            g2L{1,2}{k} = g2L{1,2}{k} + gkyy2IC{k};
        end
        cLength = length(gkyy3IC);
        g3L{1,2}(1:cLength) = g3L{1,2}(1:cLength) + gkyy3IC;
        g3L{1,2}(i)   = g3L{1,2}(i) + g1Plus;
        g3L{1,2}(j+1) = g2Plus;
        g4L{1,2} = gkyy4IC;
    else
        g3L{1,2}(i)   = g3L{1,2}(i) + g1Plus;
        g3L{1,2}(j+1) = g2Plus;
        g4L{1,2} = sdlfmKernMeanCovPartial(lfmKern1(1), lfmKern2(1), tInit1(1), ...
            tInit2(2), 1, typeParam);
    end
end

% Compute for tInit1(2) and tInit2(1)

[g1L{2,1}, g2L{2,1}, g3L{2,1}] = sdlfmvXsdlfmKernGradientBlock(lfmKern1, lfmKern2, tInit1(2), ...
    tInit2(1), kyy, kyv, kvy, kvv, i, j, generalConst, generalConstGrad, 1);

g1Plus = sdlfmaXsdlfmKernComputeBlock(lfmKern1, lfmKern2, tInit1(2), ...
    tInit2(1), kyy, kyv, kvy, kvv, i, j, generalConst);

g2Plus = sdlfmvXsdlfmvKernComputeBlock(lfmKern1, lfmKern2, tInit1(2), ...
    tInit2(1), kyy, kyv, kvy, kvv, i, j, generalConst);

if (i==1 && j==1) || (i~=1) && (j==1)
    if i~=1 || j~=1
        [gkyy1IC, gkyy2IC, gkyy3IC, gkyy4IC] = sdlfmXsdlfmKernGradientIC(lfmKern1, lfmKern2, ...
            tInit1(2), tInit2(1), gkyy1{i,j}, gkyy2{i,j}, gkyy3{i,j}, gkyy4{i,j}, gkyv1{i,j}, gkyv2{i,j}, ...
            gkyv3{i,j}, gkyv4{i,j}, gkvy1{i,j}, gkvy2{i,j}, gkvy3{i,j}, gkvy4{i,j}, gkvv1{i,j}, gkvv2{i,j}, ...
            gkvv3{i,j}, gkvv4{i,j}, 1, typeParam);
        for k=1:3
            g1L{2,1}{k} = g1L{2,1}{k} + gkyy1IC{k};
            g2L{2,1}{k} = g2L{2,1}{k} + gkyy2IC{k};
        end
        cLength = length(gkyy3IC);
        g3L{2,1}(1:cLength) = g3L{2,1}(1:cLength) + gkyy3IC;
        g3L{2,1}(j)   =  g3L{2,1}(j) + g2Plus;
        g3L{2,1}(i+1) =  g1Plus;
        g4L{2,1} = gkyy4IC;
    else
        g3L{2,1}(j)   =  g3L{2,1}(j) + g2Plus;
        g3L{2,1}(i+1) =  g1Plus;
        g4L{2,1} = sdlfmKernMeanCovPartial(lfmKern1(1), lfmKern2(1), tInit1(2), ...
            tInit2(1), 1, typeParam);
    end
end

% Compute for tInit1(2) and tInit2(2)

[g1L{2,2}, g2L{2,2}, g3Local] = sdlfmvXsdlfmKernGradientBlock(lfmKern1, lfmKern2, tInit1(2), ...
    tInit2(2), kyy, kyv, kvy, kvv, i, j, generalConst, generalConstGrad, 1);

g1Plus = sdlfmaXsdlfmKernComputeBlock(lfmKern1, lfmKern2, tInit1(2), ...
    tInit2(2), kyy, kyv, kvy, kvv, i, j, generalConst);

g2Plus = sdlfmvXsdlfmvKernComputeBlock(lfmKern1, lfmKern2, tInit1(2), ...
    tInit2(2), kyy, kyv, kvy, kvv, i, j, generalConst);

if i~=1 || j~=1
    [gkyy1IC, gkyy2IC, gkyy3IC, gkyy4IC] = sdlfmXsdlfmKernGradientIC(lfmKern1, lfmKern2, ...
        tInit1(2), tInit2(2), gkyy1{i,j}, gkyy2{i,j}, gkyy3{i,j}, gkyy4{i,j}, gkyv1{i,j}, gkyv2{i,j}, ...
        gkyv3{i,j}, gkyv4{i,j}, gkvy1{i,j}, gkvy2{i,j}, gkvy3{i,j}, gkvy4{i,j}, gkvv1{i,j}, gkvv2{i,j}, ...
        gkvv3{i,j}, gkvv4{i,j}, 1, typeParam);
    if i>=2 && j>=2
        g1L{2,2}{1} = g1L{2,2}{1} + gkyy1IC{1};
        g2L{2,2}{1} = g2L{2,2}{1} + gkyy2IC{1};
        g1L{2,2}{2} = [gkyy1IC{2} g1L{2,2}{2}];
        g2L{2,2}{2} = [gkyy2IC{2} g2L{2,2}{2}];
        g1L{2,2}{3} = [gkyy1IC{3} g1L{2,2}{3}];
        g2L{2,2}{3} = [gkyy2IC{3} g2L{2,2}{3}];
    else
        for k=1:3
            g1L{2,2}{k} = g1L{2,2}{k} + gkyy1IC{k};
            g2L{2,2}{k} = g2L{2,2}{k} + gkyy2IC{k};
        end
    end
    cLength = length(g3Local);
    g3L{2,2} = zeros(1, cLength+1);
    g3L{2,2}(1:cLength) = g3L{2,2}(1:cLength) + g3Local;
    cLength = length(gkyy3IC);
    g3L{2,2}(1:cLength) = g3L{2,2}(1:cLength) + gkyy3IC;
    g3L{2,2}(i+1) = g3L{2,2}(i+1) + g1Plus;
    g3L{2,2}(j+1) = g3L{2,2}(j+1) + g2Plus;
    g4L{2,2} = gkyy4IC;
else
    g3L{2,2} = g3Local;
    g3L{2,2}(i+1) = g1Plus;
    g3L{2,2}(j+1) = g3L{2,2}(j+1) + g2Plus;
    g4L{2,2} = sdlfmKernMeanCovPartial(lfmKern1(1), lfmKern2(1), tInit1(2), ...
        tInit2(2), 1, typeParam);
end


text:
function c = char(u)

%tstoolbox/@unit/char
%   gives the unit's label (e.g. V for Volt) back.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

c = u.label;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [Utest, Vtest] = trainMF(X, MFparam)

% Input:
%          X: n*d, n is the number of images
%          ITQparam:
%                           ITQparam.pcaW---PCA of all the database
%                           ITQparam.nbits---encoding length
% Output:
%             ITQparam:
%                              ITQparam.pcaW---PCA of all the database
%                              ITQparam.nbits---encoding length
%                              ITQparam.r---ITQ rotation projection

gamma =0.0001;

lambda = 0.001;

iter = 1;
lastF = 99999999;
threshold = 0.001;

[n, d] = size(X);

%V = rand(bits, ntrain);

% initialize with a orthogonal random rotation

d = MFparam.nbits;

U = randn(n, d);

%V = randn(d, d);

% ITQ to find optimal rotation
while (true)

    V = (U'*U+gamma*eye(d))\(U'*X);

    U = (X*V')/(V*V'+lambda*eye(d));

%     U = (X*V')/(V*V'+lambda*eye(d)); % this is not good
%     V = (U'*U+gamma*eye(d))\(U'*X);

    % compute objective function
    norm1 = norm(X - U * V, 'fro');
    norm2 = lambda * norm(U, 'fro') + gamma * norm(V, 'fro');
    currentF= norm1 + norm2;

    fprintf('\nobj at iteration %d: %.4f\n reconstruction error for matrix factorization: %.4f,\n regularization term: %.4f\n\n', iter, currentF, norm1 , norm2);
    if (lastF - currentF) < threshold
        fprintf('algorithm converges...\n');
        fprintf('final obj: %.4f\n reconstruction error for matrix factorization: %.4f, \n regularization term: %.4f\n\n', currentF,norm1, norm2);

        Utest = U;
        Vtest = V;
        return
    end
    iter = iter + 1;
    lastF = currentF;
end

% make B binary
%B = UX;
%B(B<0) = 0;

%ITQTparam.r = R;

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function [P]=patchSmooth(F,V,IND_V,cPar)

% function [P]=patchSmooth(F,V,IND_V,cPar)
% ------------------------------------------------------------------------
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 2014/06/02
%------------------------------------------------------------------------

%%

%Get/set method
if isfield(cPar,'Method')
    smoothMethod=cPar.Method;
else
    smoothMethod='LAP'; %DEFAULT
end

if ~isfield(cPar,'n')
    cPar.n=1; %DEFAULT
end

if isempty(IND_V)
    [~,IND_V]=patchIND(F,V,2);
%     if isa(F,'cell')
%         IND_V=[];
%         for q=1:1:numel(F)
%             C=patchConnectivity(F{q},V,'vv');
%             IND_V=[IND_V C.vertex.vertex];
%         end
%     else
%         C=patchConnectivity(F,V,'vv');
%         IND_V=C.vertex.vertex;
%     end

end

%Smooth
switch smoothMethod
    case 'LAP' %Laplacian
        [P]=tesSmooth_LAP(F,V,IND_V,cPar);
    case 'HC' %Humphreys Classes
        [P]=tesSmooth_HC(F,V,IND_V,cPar);
    case 'tLAP' %Tangent Laplacian

        %Invert face orientation if required
        [logicFlip]=isGlobalSurfDirOutward(F,V);
        if ~logicFlip
            F=fliplr(F); %Flip faces
        end

        %Set control parameters for Laplacian smoothening iterations
        cPar_t=cPar;
        cPar_t.n=1;
        cPar_t.Tolerance=[];
        P=V;
        for q=1:1:cPar.n
            [Ps]=tesSmooth_LAP(F,P,IND_V,cPar_t);  %The Laplacian smoothened coordinate set
            D=Ps-P; %smoothening intended displacement vectors
            [Dt]=patchVectorTangent(F,P,D,[]); %Tangential component of displacement
            P=P+Dt; %Displace mesh
        end

    case 'tHC' %Tangent HC

        %Invert face orientation if required
        [logicFlip]=isGlobalSurfDirOutward(F,V);
        if ~logicFlip
            F=fliplr(F); %Flip faces
        end

        %Set control parameters for Laplacian smoothening iterations
        cPar_t=cPar;
        cPar_t.n=1;
        cPar_t.Tolerance=[];
        P=V;
        for q=1:1:cPar.n
            [Ps]=tesSmooth_HC(F,P,IND_V,cPar_t);  %The smoothened coordinate set
            D=Ps-P; %smoothening intended displacement vectors
            [Dt]=patchVectorTangent(F,P,D,[]); %Tangential component of displacement
            P=P+Dt; %Displace mesh
        end

    otherwise
        error('Invalid smooth method specified');
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function kts = kmph2kts(kmph)
%KMPH2KTS Convert speed from kilometers per hour to knots
%
%  kts = KMPH2KTS(kmph) converts speeds from kilometers per hour to knots.
%
%  See also KMPH2FTPS, KMPH2MPH, KMPH2MPS, KTS2KMPH.

% Jonathan Sullivan
% Original: May 2011
% jonathan.sullivan@ll.mit.edu

kts = kmph/1.852;


text:
function count = compute_output(plaza)
count = 0;
[a, b] = size(plaza);
for j = 1:b
count = count + (plaza(a,j) > 0);
end

author: Eurus-Holmes


text:
% pop_firpm() - Filter data using Parks-McClellan FIR filter
%
% Usage:
%   >> [EEG, com, b] = pop_firpm(EEG); % pop-up window mode
%   >> [EEG, com, b] = pop_firpm(EEG, 'key1', value1, 'key2', ...
%                                value2, 'keyn', valuen);
%
% Inputs:
%   EEG       - EEGLAB EEG structure
%   'fcutoff' - vector or scalar of cutoff frequency/ies (~-6 dB; Hz)
%   'ftrans'  - scalar transition band width
%   'ftype'   - char array filter type. 'bandpass', 'highpass',
%               'lowpass', or 'bandstop'
%   'forder'  - scalar filter order. Mandatory even
%
% Optional inputs:
%   'wtpass'  - scalar passband weight
%   'wtstop'  - scalar stopband weight
%
% Outputs:
%   EEG       - filtered EEGLAB EEG structure
%   com       - history string
%   b         - filter coefficients
%
% Note:
%   Requires the signal processing toolbox.
%
% Author: Andreas Widmann, University of Leipzig, 2005
%
% See also:
%   firfilt, pop_firpmord, plotfresp, firpm, firpmord

%123456789012345678901234567890123456789012345678901234567890123456789012

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

function [EEG, com, b] = pop_firpm(EEG, varargin)

    if ~(exist('firpm', 'file') == 2 || exist('firpm', 'file') == 6)
       error('Requires the signal processing toolbox.');
    end

    com = '';
    if nargin < 1
        help pop_firpm;
        return;
    end
    if isempty(EEG.data)
        error('Cannot process empty dataset');
    end

    if nargin < 2
        drawnow;
        ftypes = {'bandpass' 'highpass' 'lowpass' 'bandstop'};
        uigeom = {[1 0.75 0.75] [1 0.75 0.75] [1 0.75 0.75] 1 [1 0.75 0.75] [1 0.75 0.75] [1 0.75 0.75] 1 [1 0.75 0.75]};
        uilist = {{'Style' 'text' 'String' 'Cutoff frequency(ies) [hp lp] (~-6 dB; Hz):'} ...
                  {'Style' 'edit' 'String' '' 'Tag' 'fcutoffedit'} {} ...
                  {'Style' 'text' 'String' 'Transition band width:'} ...
                  {'Style' 'edit' 'String' '' 'Tag' 'ftransedit'} {} ...
                  {'Style' 'text' 'String' 'Filter type:'} ...
                  {'Style' 'popupmenu' 'String' ftypes 'Tag' 'ftypepop'} {} ...
                  {} ...
                  {'Style' 'text' 'String' 'Passband weight:'} ...
                  {'Style' 'edit' 'String' '' 'Tag' 'wtpassedit'} {} ...
                  {'Style' 'text' 'String' 'Stopband weight:'} ...
                  {'Style' 'edit' 'String' '' 'Tag' 'wtstopedit'} {} ...
                  {'Style' 'text' 'String' 'Filter order (mandatory even):'} ...
                  {'Style' 'edit' 'String' '' 'Tag' 'forderedit'} ...
                  {'Style' 'pushbutton' 'String' 'Estimate' 'Tag' 'orderpush' 'Callback' {@comcb, ftypes, EEG.srate}} ...
                  {} ...
                  {} {} {'Style' 'pushbutton' 'String', 'Plot filter responses' 'Tag' 'plotpush' 'Callback' {@comcb, ftypes, EEG.srate}}};
        result = inputgui(uigeom, uilist, 'pophelp(''pop_firpm'')', 'Filter the data -- pop_firpm()');
        if isempty(result), return; end

        args = {};
        if ~isempty(result{1})
            args = [args {'fcutoff'} {str2num(result{1})}];
        end
        if ~isempty(result{2})
            args = [args {'ftrans'} {str2double(result{2})}];
        end
        args = [args {'ftype'} ftypes(result{3})];
        if ~isempty(result{4})
            args = [args {'wtpass'} {str2double(result{4})}];
        end
        if ~isempty(result{5})
            args = [args {'wtstop'} {str2double(result{5})}];
        end
        if ~isempty(result{6})
            args = [args {'forder'} {str2double(result{6})}];
        end
    else
        args = varargin;
    end

    % Convert args to structure
    args = struct(args{:});

    c = parseargs(args, EEG.srate);
    if ~isfield(args, 'forder') || isempty(args.forder)
        error('Not enough input arguments');
    end
    b = firpm(args.forder, c{:});

    % Filter
    disp('pop_firpm() - filtering the data');
    EEG = firfilt(EEG, b);

    % History string
    com = sprintf('%s = pop_firpm(%s', inputname(1), inputname(1));
    for c = fieldnames(args)'
        if ischar(args.(c{:}))
            com = [com sprintf(', ''%s'', ''%s''', c{:}, args.(c{:}))];
        else
            com = [com sprintf(', ''%s'', %s', c{:}, mat2str(args.(c{:})))];
        end
    end
    com = [com ');'];

% Convert structure args to cell array firpm parameters
function c = parseargs(args, srate)

    if ~isfield(args, 'fcutoff') || ~isfield(args, 'ftype') || ~isfield(args, 'ftrans') || isempty(args.fcutoff) || isempty(args.ftype) || isempty(args.ftrans)
        error('Not enough input arguments.');
    end

    % Cutoff frequencies
    args.fcutoff = [args.fcutoff - args.ftrans / 2 args.fcutoff + args.ftrans / 2];
    args.fcutoff = sort(args.fcutoff / (srate / 2)); % Sorting and normalization
    if any(args.fcutoff < 0)
        error('Cutoff frequencies - transition band width / 2 must not be < DC');
    elseif any(args.fcutoff > 1)
        error('Cutoff frequencies + transition band width / 2 must not be > Nyquist');
    end
    c = {[0 args.fcutoff 1]};

    % Filter type
    switch args.ftype
        case 'bandpass'
            c = [c {[0 0 1 1 0 0]}];
        case 'bandstop'
            c = [c {[1 1 0 0 1 1]}];
        case 'highpass'
            c = [c {[0 0 1 1]}];
        case 'lowpass'
            c = [c {[1 1 0 0]}];
    end

    %Filter weights
    if all(isfield(args, {'wtpass', 'wtstop'})) && ~isempty(args.wtpass) && ~isempty(args.wtstop)
        w = [args.wtstop args.wtpass];
        c{3} = w(c{2}(1:2:end) + 1);
    end

% Callback
function comcb(obj, evt, ftypes, srate)

    args.fcutoff = str2num(get(findobj(gcbf, 'Tag', 'fcutoffedit'), 'String'));
    args.ftype = ftypes{get(findobj(gcbf, 'Tag', 'ftypepop'), 'Value')};
    args.ftrans = str2double(get(findobj(gcbf, 'Tag', 'ftransedit'), 'String'));
    args.wtpass = str2double(get(findobj(gcbf, 'Tag', 'wtpassedit'), 'String'));
    args.wtstop = str2double(get(findobj(gcbf, 'Tag', 'wtstopedit'), 'String'));
    c = parseargs(args, srate);

    switch get(obj, 'Tag')
        case 'orderpush'
            [args.forder, args.wtpass, args.wtstop] = pop_firpmord(c{1}(2:end - 1), c{2}(1:2:end));
            if ~isempty(args.forder) || ~isempty(args.wtpass) || ~isempty(args.wtstop)
                set(findobj(gcbf, 'Tag', 'forderedit'), 'String', ceil(args.forder / 2) * 2);
                set(findobj(gcbf, 'Tag', 'wtpassedit'), 'String', args.wtpass);
                set(findobj(gcbf, 'Tag', 'wtstopedit'), 'String', args.wtstop);
            end

        case 'plotpush'
            args.forder = str2double(get(findobj(gcbf, 'Tag', 'forderedit'), 'String'));
            if isempty(args.forder)
                error('Not enough input arguments');
            end
            b = firpm(args.forder, c{:});
            H = findobj('Tag', 'filter responses', 'Type', 'figure');
            if ~isempty(H)
                figure(H);
            else
                H = figure;
                set(H, 'color', [.93 .96 1], 'Tag', 'filter responses');
            end
            plotfresp(b, 1, [], srate);
    end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ x, err, iter, flag ] = sor(A, x, b, w, max_it, tol)
%
% [ x, err, iter, flag ] = sor(A, x, b, w, max_it, tol)
%
% SOR Successive Over-Relaxation Method
%   This function solves linear equation systems such as Ax=b using SOR
%   method (Successive Over-Relaxation).
%   When the relaxation scalar w=1, the method used is Gauss-Seidel.
%
% Input:
%   A - input matrix
%   x - inicial vector
%   b - vector b
%   w - relaxation scalar
%   max_it - maximum number of iterations
%   tol - tolerance
%
% Output:
%   x - solution vector
%   err - norm err estimate
%   iter - nu,ber of performed iterations
%   flag - 0 = a solution was found / 1 = no convergence for max_it
%
% Example:
%   [ x, err, iter, flag ] = sor( [.5 .125; .125 .25], zeros(2,1),
%       ones(2,1), .5, 1e4, 1e-4 )
%
% Author:	Tashi Ravach
% Version:	1.0
% Date:     08/06/2006
%

    if nargin == 3
        w = .5;
        max_it = 1e4;
        tol = 1e-4;
    elseif nargin == 4
        max_it = 1e4;
        tol = 1e-4;
    elseif nargin == 5
        tol = 1e-4;
    elseif nargin ~= 6
        error('sor: invalid input parameters');
    end

    flag = 0;
    iter = 0;

    norma2_b = norm(b);
    if (norma2_b == 0.0)
        norma2_b = 1.0;
    end

    r = b - A * x;
    err = norm(r) / norma2_b;
    if (err < tol)
        return
    end

    % separate A into several matrix for SOR/Gauss-Seidel
    [ M, N, b ] = matsep(A, b, w, 2);

    for iter = 1 : max_it
        x_1 = x;
        x = M \ (N * x + b); % adjust the aproximation
        %err = norm(x - x_1) / norm(x); % compute error
        err = norm(x_1 - x, 1); % compute error
        if (err <= tol) % check for convergence
            break
        end
    end
    b = b / w; % vector b

    if (err > tol) % no convergence
        flag = 1;
    end

end

function [ M, N, b ] = matsep(A, b, w, flag)
%
% [ M, N, b ] = matsep(A, b, w, flag)
%
% MATSEP Matrix Separation
%   Input matrix is splitted into several others in diferent ways depending
%   on the method to be used: Jacobi and SOR (Gauss-Seidel when w = 1)
%
% Input:
%   A - input matrix
%   x - inicial vector
%   b - vector b
%   flag - 1 = Jacobi / 2 = SOR
%
% Output:
%   M - matrix M
%   N - matrix N
%   b - vector b (modified for SOR)
%
% Author:	Tashi Ravach
% Version:	1.0
% Date:     08/06/2006
%

    if nargin ~= 4
        error('matsep: invalid input parameters');
    end

    [ m, n ] = size(A);
    if m ~= n
        error('matsep: input matrix A must have dimension nXn');
    end

    [ l, o ] = size(b);
    if l ~= n && o ~= 1
        error('matsep: input matrix b must have dimension nX1');
    end

    if (flag == 1) % separation for Jacobi
        M = diag(diag(A));
        N = diag(diag(A)) - A;
    elseif (flag == 2) % separation for SOR/Gauss-Seidel
        b = w * b;
        M =  w * tril(A, -1) + diag(diag(A));
        N = -w * triu(A,  1) + (1.0 - w) * diag(diag(A));
    end

end


text:
function [pdf,X1,X2]=akde(X,grid,gam)
%% adaptive kernel density estimation in high dimensions;
%  optimal accuracy/speed tradeoff, controlled via parameter "gam";
% INPUTS:   X  - data as a 'n' by 'd' vector;
%
%         grid - 'm' points of dimension 'd' over which pdf is computed;
%                default provided only for 2-dimensional data;
%                see example on how to construct it in higher dimensions
%
%          gam - cost/accuracy tradeoff parameter, where gam<n;
%                default value is gam=ceil(n^(1/2)); larger values
%                may result in better accuracy, but always reduce speed;
%                to speedup the code, reduce the value of "gam";
%
% OUTPUT: pdf   - the value of the estimated density at 'grid'
%         X1,X2 - grid only for 2 dimensional data
%
%%  EXAMPLE in 2 dimensions:
%   L=chol([1,-0.999;-0.999,1],'lower');L1=chol([1,0.999;0.999,1],'lower');
%   data=[(L1*randn(10^3,2)')';(L*randn(10^3,2)')'*2;rand(10^4,2)*5-2.5];
%   [pdf,X1,X2]=akde(data);pdf=reshape(pdf,size(X1));contour(X1,X2,pdf,20)
%
%%  EXAMPLE in 3 dimensions:
%  data=[randn(10^3,3);randn(10^3,3)/2+2]; % three dimensional data
%  [n,d]=size(data); ng=100; % total grid points = ng^d
%  MAX=max(data,[],1); MIN=min(data,[],1); scaling=MAX-MIN;
%  % create meshgrid in 3-dimensions
%  [X1,X2,X3]=meshgrid(MIN(1):scaling(1)/(ng-1):MAX(1),...
%      MIN(2):scaling(2)/(ng-1):MAX(2),MIN(3):scaling(3)/(ng-1):MAX(3));
%  grid=reshape([X1(:),X2(:),X3(:)],ng^d,d); % create points for plotting
%  pdf=akde(data,grid); % run adaptive kde
%  pdf=reshape(pdf,size(X1)); % reshape pdf for use with meshgrid
%  for iso=[0.005:0.005:0.015] % isosurfaces with pdf = 0.005,0.01,0.015
%      isosurface(X1,X2,X3,pdf,iso),view(3),alpha(.3),box on,hold on
%      colormap cool
%  end
%
%%  Reference:
%  Kernel density estimation via diffusion
%  Z. I. Botev, J. F. Grotowski, and D. P. Kroese (2010)
%  Annals of Statistics, Volume 38, Number 5, pages 2916-2957.
[n,d]=size(X);
% begin scaling preprocessing
MAX=max(X,[],1);MIN=min(X,[],1);scaling=MAX-MIN;
MAX=MAX+scaling/10;MIN=MIN-scaling/10;scaling=MAX-MIN;
X=bsxfun(@minus,X,MIN);X=bsxfun(@rdivide,X,scaling);
if (nargin<2)|isempty(grid) % failing to provide grid
    warning('Assuming data is 2 dimensional. For higher dimensions, provide a grid as in example.')
    % create meshgrid in 2-dimensions
    [X1,X2]=meshgrid(MIN(1):scaling(1)/(2^7-1):MAX(1),...
           MIN(2):scaling(2)/(2^7-1):MAX(2));
    grid=reshape([X1(:),X2(:)],2^14,d); % create grid for plotting
end
mesh=bsxfun(@minus,grid,MIN);mesh=bsxfun(@rdivide,mesh,scaling);
if nargin<3 % failing to provide speed/accuracy tradeoff
    gam=ceil(n^(1/2));
end
% end preprocessing
% algorithm initialization
del=.1/n^(d/(d+4));perm=randperm(n);mu=X(perm(1:gam),:);w=rand(1,gam);
w=w/sum(w);Sig=bsxfun(@times,rand(d,d,gam),eye(d)*del);ent=-Inf;
for iter=1:1500 % begin algorithm
    Eold=ent;
    [w,mu,Sig,del,ent]=regEM(w,mu,Sig,del,X); % update parameters
    err=abs((ent-Eold)/ent); % stopping condition
    fprintf('Iter.    Tol.      Bandwidth \n');
    fprintf('%4i    %8.2e   %8.2e\n',iter,err,del);
    fprintf('----------------------------\n');
    if (err<10^-4)|(iter>200), break, end
end
% now output density values at grid
pdf = probfun(mesh,w,mu,Sig)/prod(scaling); % evaluate density
del=del*scaling; % adjust bandwidth for scaling
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pdf=probfun(x,w,mu,Sig)
[gam,d]=size(mu);
pdf=0;
for k=1:gam
    L=chol(Sig(:,:,k));s=diag(L);
    logpdf=-.5*sum(( bsxfun(@minus,x,mu(k,:))/L).^2,2)+log(w(k))...
        -sum(log(s))-d*log(2*pi)/2;
    pdf=pdf+exp(logpdf);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [w,mu,Sig,del,ent]=regEM(w,mu,Sig,del,X)
[gam,d]=size(mu);[n,d]=size(X);
log_lh=zeros(n,gam); log_sig=log_lh;
for i=1:gam
    L=chol(Sig(:,:,i));
    Xcentered = bsxfun(@minus, X, mu(i,:));
    xRinv = Xcentered /L; xSig = sum((xRinv /L').^2,2)+eps;
    log_lh(:,i)=-.5*sum(xRinv.^2, 2)-sum(log(diag(L)))...
        +log(w(i))-d*log(2*pi)/2-.5*del^2*trace((eye(d)/L)/L');
    log_sig(:,i)=log_lh(:,i)+log(xSig);
end
maxll = max (log_lh,[],2); maxlsig = max (log_sig,[],2);
p= exp(bsxfun(@minus, log_lh, maxll));
psig=exp(bsxfun(@minus, log_sig, maxlsig));
density = sum(p,2);  psigd=sum(psig,2);
logpdf=log(density)+maxll; logpsigd=log(psigd)+maxlsig;
p = bsxfun(@rdivide, p, density);
ent=sum(logpdf); w=sum(p,1);
for i=find(w>0)
    mu(i,:)=p(:,i)'*X/w(i);  %compute mu's
    Xcentered = bsxfun(@minus, X,mu(i,:));
    Xcentered = bsxfun(@times,sqrt(p(:,i)),Xcentered);
    Sig(:,:,i)=Xcentered'*Xcentered/w(i)+del^2*eye(d); % compute sigmas;
end
w=w/sum(w);curv=mean(exp(logpsigd-logpdf)); % estimate curvature
del=1/(4*n*(4*pi)^(d/2)*curv)^(1/(d+2));
end

repo: cluster


text:
function varargout = circleToPolygon(circle, varargin)
%CIRCLETOPOLYGON Convert a circle into a series of points.
%
%   PTS = circleToPolygon(CIRC, N);
%   Converts the circle CIRC into an array of  N-by-2 of double, containing
%   x and y positions of vertices.
%   CIRC is given as [x0 y0 r], where x0 and y0 are coordinate of center,
%   and r is the radius.
%
%   P = circleToPolygon(CIRCLE);
%   uses a default value of N=64 vertices.
%
%   Example
%     poly = circleToPolygon([30 20 15], 16);
%     figure; hold on;
%     axis equal;axis([0 50 0 50]);
%     drawPolygon(poly, 'b');
%     drawPoint(poly, 'bo');
%
%   See also:
%   circles2d, polygons2d, circleArcToPolyline, ellipseToPolygon
%

% ---------
% author : David Legland
% created the 06/04/2005.
% Copyright 2010 INRA - Cepia Software Platform.
%

% HISTORY
% 2007-04-20 return a closed polygon with N+1 vertices, use default N=64
% 2011-12-09 rename to 'circleToPolygon'
% 2017-08-31 returns N vertices instead of N+1

% determines number of points
N = 64;
if ~isempty(varargin)
    N = varargin{1};
end

% create circle
t = linspace(0, 2*pi, N+1)';
t(end) = [];

% coordinates of circle points
x = circle(1) + circle(3) * cos(t);
y = circle(2) + circle(3) * sin(t);

% foramt output
if nargout == 1
    varargout{1} = [x y];
elseif nargout == 2
    varargout{1} = x;
    varargout{2} = y;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = mapn(array, start, varargin)
% [OUTPUT] = MAPN(ARRAY, START, FUNC, OPTIONS)
%   This is MAP(ARRAY, FUNC', OPTIONS) where FUNC' takes the ordinal
%   position of the array element we are mapping across instead of the
%   actual value.  Positions start counting at start.
%

  [funclist, varargin] = select(varargin, @islambda);

  if nargin < 2, start = 1; end
  if nargin == 0
    error('MapN:Arguments', 'Array argument to MapN is required.');
  end

  if isempty(funclist)
    varargout{1} = start:(numel(array) - 1 + start);
    return;
  end

  if numel(funclist) > 1
    error('MapN:Arguments', 'Too many functions to mapn.');
  end

  func = funclist{1};
  function varargout = mapNFunc(index, varargin)
    [varargout{1:(max([nargout, 1]))}] = func(index);
  end

  [varargout{1:(max([nargout, 1]))}] = map(array, indexer(@mapNFunc, start, 1, 1), varargin{:});
end


text:
function[f,name]=tidefreq(str)
%TIDEFREQ  Frequencies of the eight major tidal components.
%
%   [F,NAME]=TIDEFREQ returns the frequencies of the eight major tidal
%   components together with a string matrix containing their names.  In
%   order of increasing frequency, these are
%
%        Mf O1 P1 K1 N2 M2 S2 K2.
%
%   The units are given in *radians* per hour. See Gill (1982) page 335.
%
%   F=TIDEFREQ(STR) returns the frequency for the component named STR,
%   if STR is a string.  F=TIDEFREQ(N) where N is a number also works.
%
%   Usage:   [f,name]=tidefreq;
%            f=tidefreq(str);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2005--2015 J.M. Lilly --- type 'help jlab_license' for details

name=('MFO1P1K1N2M2S2K2');

%Gill only gives first two decimal points
%I got the more exact numbers from
%http://oceanworld.tamu.edu/resources/ocng_textbook/chapter17/chapter17_04.htm

p(1)=327.85;
p(2)=25.8194;
p(3)=24.0659;
p(4)=23.9344;
p(5)=12.6584;
p(6)=12.4206;
p(7)=12.0000;
p(8)=11.9673;

f=1./p;
if nargin==1
  if ~ischar(str)
    f=f(str);
  else
    str=upper(str);
    ii=strfind(name,str);
    f=f((ii-1)/2+1);
  end
  clear name
else
  name=reshape(name,2,8)';
end

f=f*2*pi;

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [image_information_metric, rniirs, diagnostic] = RGIQE(filename, varargin)
% RGIQE Radar generalized image quality equation
%    [image_information_metric, rniirs] = RGIQE(bandwidth, nesz_db, graze_degrees)
%    [image_information_metric, rniirs] = RGIQE(filename, 'PropertyName', PropertyValue, ...)
%
% An implementation of a generalized image quality equation to provide a
% quantitative image quality metric and an RNIIRS prediction.
%
% RNIIRS (Radar National Imagery Interpretation Rating Scale) is an image
% quality scale from 1 to 9.  This "quality" refers to local image fidelity
% and a users ability to perform certain detection and recognition tasks
% with an image.  The "quality" that RNIIRS attempts to capture does NOT
% refer to image size, geopositioning accuracy, radiometric accuracy, or
% any other host of characterists that may improve an image's utility.
% Although RNIIRS is subjective, generally assessed by a trained analyst,
% an image quality equation (IQE) may be used to predict RNIIRS for a given
% SAR collect.
%
% The General Image Quality Equation (GIQE) provided here uses generic,
% fundamental, sensor-indepenedent parameters of radar collection to
% measure local image fidelity by quantifying the information density in
% the ground in terms of bits\m^2, using Shannon's channel capacity
% theorem.  Unlike RNIIRS, this value is quantitative and not subjective or
% based on analysts' assessments.  We can however compute an RNIIRS
% prediction based off of it by fitting the information metric to analysts
% ratings.  Both the information theory metric and the RNIIRS prediction
% are computed in this function.
%
%    [image_information_metric, rniirs] = RGIQE(bandwidth, nesz_db, graze_degrees)
% INPUTS:
%
%    bandwidth               Transmitted (and processed) bandwidth in Hz.
%    nesz_db                 Noise equivalent sigma-0 in dB.  Should
%                               capture all noise, including multiplicative
%                               sources.
%    graze_degrees           Grazing angle in degrees.
%
%
%    [image_information_metric, rniirs] = RGIQE(filename, 'PropertyName', PropertyValue, ...)
% INPUTS:
%
%    filename                Must either be a complex image in a format
%                               recognized by the MATLAB SAR Toolbox or an
%                               XML file with SICD-formatted XML.  If empty
%                               or not passed, the user will be prompted to
%                               select a file with a dialog box.
%
%    Property name           Description
%    multi_noise             Multiplicative noise.  This value is equal to:
%
%                               MNR * sigma_0_background
%
%                               where MNR is the multiplicative noise ratio
%                               and sigma_0_background is the assumed power
%                               level of the background and ambiguous
%                               areas.
%                               MNR can further be broken down into:
%
%                               MNR = ISLR + QNR + AMBR
%
%                               where ISLR is the integrated sidelobe
%                               ratio, QNR is the quantization noise ratio,
%                               and AMBR is the ambiguity ratio.
%                               (Default = 0)
%    interactive             "AUTO", "NOISE", "SIGNAL", "BOTH", or "NONE".
%                               Enable interactive measurement of noise and
%                               signal levels.  "AUTO" only runs
%                               interactive measurment if metadata from
%                               source file does not provide it. (Default
%                               is "AUTO").
%    frames                  If a file contains more than one complex
%                               dataset, this value specifies for which of
%                               those datasets to compute metrics. (Default
%                               is all datasets in the file.)
%    signal_sigma_override   Allows caller to override signal level used
%                               for computation (and INTERACTIVE input
%                               parameter).  In general, this is not
%                               recommended, as RNIIRS estimation function
%                               was trained using the default values. If
%                               the dataset provided is radiometrically
%                               calibrated, this value will be treated as
%                               calibrated signal level in the ground plane
%                               (sigma-0).  If the dataset is not
%                               calibrated, this value will be assume to be
%                               average pixel power. (Default is 1 square
%                               meter of radar return per square meter on
%                               the ground for co-pol, 0.25 for cross-pol.)
%
% Note that the selection of signal level and noise level are fundamentally
% different in this model.  Noise is measureable, quantifiable, and likely
% fairy accuractely predictable (at least for additive noise).  There is a
% right answer for what the noise level is for any given dataset. On the
% other hand, signal level is somewhat arbitrarily chosen, selected to be
% the power level of the scene content that one may be most interested in
% examining.  It is for this reason that this function allows for a signal
% level override, but not a direct noise level override (although noise can
% be adjusted with the MULTI_NOISE input parameter.)
%
% Author: Wade Schwartzkopf, NGA/Research

% Of the lines in this file, >95% is file and input parameter parsing and
% documention, and <5% actually computing the image quality metrics...

%% Bandwidth formulation
% Assumes azimuth bandwidth to match range in ground plane
% NESZ provided should capture all noise, including multiplicative sources
if nargin == 3 && isnumeric(filename) && isnumeric(varargin{1}) && isnumeric(varargin{2})
    bandwidth = filename;  % Hz
    nesz = varargin{1};  % in dB
    graze = varargin{2};  % in degrees
    bandwidth = bandwidth*2/SPEED_OF_LIGHT;  % Convert units from Hz to cycles/m
    bandwidth = bandwidth*cosd(graze);  % Convert from slant plane to ground plane
    bandwidth = bandwidth.^2;  % Convert from 1D bandwidth to 2D bandwidth area (assumes azimuth bandwidth same as range)
    image_information_metric = bandwidth.*log2(1+1/(10^(nesz/10)));
    rniirs = estimate_rniirs(image_information_metric);
    diagnostic = struct();  % No diagnostics for this version
    return;
end

%% SICD/SLC version.
%% Open input file
if ((nargin<1)||isempty(filename)) % If no filename was give, use dialog box to ask for one
    % Recall last interactively selected path used
    if ispref('matlab_sar_toolbox','last_used_directory')
        pathname = getpref('matlab_sar_toolbox','last_used_directory');
        if ~ischar(pathname)||~exist(pathname,'dir')
            pathname = pwd;
        end
    else
        pathname = pwd;
    end
    [filename,pathname]=uigetfile(sar_file_extensions(),...
        'Open SAR Data File',pathname);
    setpref('matlab_sar_toolbox','last_used_directory',pathname);
else % Path was already passed in with filenames
    pathname='';
end
if(filename)
    fullfilename=[pathname filename];
else % filename=0.  Cancel was pressed, instead of a file being chosen.
    return;
end
if ~isempty(pathname)
    setpref('matlab_sar_toolbox','last_used_directory',pathname); %store path
end
try
    ro = open_reader(fullfilename);
    % Check for multiple images
    if ~iscell(ro)
        ro = {ro};
    end
    meta = cell(numel(ro),1);
    for i = 1:numel(ro)
        meta{i} = ro{i}.get_meta();
    end
catch
    try
        fid = fopen(fullfilename);
        meta = {sicdxml2struct(xmlread( java.io.StringBufferInputStream( char(fread(fid)'))))};
        fclose(fid);
    catch
        error('RGIQE:InvalidFile','Invalid File.');
    end
end

%% Parse other input parameters
p = inputParser;
p.addParameter('multi_noise', 0, @(x) isnumeric(x) && isscalar(x));
p.addParameter('interactive', 'AUTO', @(x) ismember(upper(x), ...
    {'AUTO', 'NOISE', 'SIGNAL', 'BOTH', 'NONE'}));
p.addParameter('frames', 1:numel(meta), @(x) all(ismember(x, 1:numel(meta))));
p.addParameter('signal_sigma_override', [], @(x) isnumeric(x) && isscalar(x));
parse(p,varargin{:})

%% Iterate through datasets provided
meta = meta(p.Results.frames);
image_information_metric = zeros(numel(meta),1);
rniirs = zeros(numel(meta),1);
for i = 1:numel(meta)
    meta{i} = derived_sicd_fields(meta{i});  % Assure all radiometric fields are derived if not there

    %% Get noise level
    if isfield(meta{i}, 'Radiometric') && ...
            isfield(meta{i}.Radiometric,'NoiseLevel') && ...
            isfield(meta{i}.Radiometric.NoiseLevel,'NoiseLevelType') && ...
            strcmpi(meta{i}.Radiometric.NoiseLevel.NoiseLevelType,'ABSOLUTE') && ...
            isfield(meta{i}.Radiometric.NoiseLevel,'NoisePoly') && ...
            ~ismember(upper(p.Results.interactive), {'NOISE', 'BOTH'})
        % SICD metadata can only tell us additive noise.  Estimation of
        % multiplicative noise must come from elsewhere.
        noise_sigma = 10^(meta{i}.Radiometric.NoiseLevel.NoisePoly(1)/10) + p.Results.multi_noise;
    elseif ismember(upper(p.Results.interactive), {'AUTO', 'NOISE', 'BOTH'})
        % Measured noise contains both additive and multiplicative noise
        uiwait(msgbox('Please select a noise region for measurement.','Noise measurement'));
        aoi = mitm_viewer(fullfilename,'mode','aoi','closeAfterSelect',true,'initialFrame',i);
        if isempty(aoi), return; end  % Cancel was selected
        cdata = double(ro{i}.read_chip(...
            [aoi(1) aoi(1)+aoi(3)-1],[aoi(2) aoi(2)+aoi(4)-1]));
        noise_sigma = mean(abs(cdata(:)).^2);
    else
        noise_sigma = NaN;
    end
    if isfield(meta{i}, 'Radiometric') && ...
            isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly')
        % Scale to be equal to noise equivalent sigma zero
        noise_sigma = meta{i}.Radiometric.SigmaZeroSFPoly(1) * noise_sigma;
    end

    %% Get signal level
    % One issue with using the Shannon-Hartley channel capacity formulation
    % is that it requires a signal level.  Certainly the amount of
    % information a SAR image conveys is related to its signal level, so
    % this is how we hope it would behave.  A target with brighter returns
    % is always easier to interpret than a target that is in the noise.
    % However, it presents a problem to image quality prediction, since the
    % signal level can't be known prior to measurement.  This constraint is
    % implicitly true for traditional RNIIRS assessments as well, which
    % assumes specific types of targets.
    % One approach for setting the signal level used in the channel
    % capacity equation would be to assume a fixed radiometric signal level
    % that is representative of the sorts of scene content we are generally
    % interested in and use that same signal level across all images for
    % which we compute this metric.  This allows for an even
    % apples-to-apples comparison across images of varying content and is
    % consistent with the traditional RNIIRS viewpoint, which assumes
    % specific types of tasks and scene content.
    % Another approach would be to use the signal level in the actual scene
    % by measuring the average pixel power across a full image or AOI. This
    % approach may be more accurate at measuring the true information in
    % any given image area.  However, it may be less consistent with the
    % traditional RNIIRS viewpoint, and this approach can only be used
    % after the collection has been made, not prior to it for planning
    % collection and predicting image quality.
    % We also allow the caller of the function to pass their own signal
    % level. Passing a non-default signal level is not recommended for
    % RNIIRS estimation, since the RNIIRS function was trained on the
    % default values.  However, we provide this flexibility in case it
    % might be educational to see how it affects the information theory
    % metric.
    if ~ismember('signal_sigma_override', p.UsingDefaults)
        signal_sigma = p.Results.signal_sigma_override;
    elseif isfield(meta{i}, 'Radiometric') && ...
            isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly') && ...
            ~ismember(upper(p.Results.interactive), {'SIGNAL', 'BOTH'})
        % Default value: 1 square meter of return per square meter on
        % the ground (roughly similar to a large vehicle over its area).
        % Ulaby plots show cross-pol generally 5-10dB lower than
        % co-pol.  We pick 1/4 as simple round number.
        if isfield(meta{i},'ImageFormation') && ...
                isfield(meta{i}.ImageFormation,'TxRcvPolarizationProc') && ...
                numel(unique(split(meta{i}.ImageFormation.TxRcvPolarizationProc,':')))>1
            signal_sigma = 0.25;
        else
            signal_sigma = 1;
        end
    elseif ismember(upper(p.Results.interactive), {'AUTO', 'SIGNAL', 'BOTH'})
        ButtonName = questdlg(['Do you want to use average signal '...
            'power across entire image or select an AOI?'], ...
            'Signal level estimation', ...
            'Image', 'AOI', 'Image');
        if isempty(ButtonName), return; end %window was closed CANCEL
        if strcmpi(ButtonName,'AOI')
            aoi = mitm_viewer(fullfilename,'mode','aoi','closeAfterSelect',true,'initialFrame',i);
            if isempty(aoi), return; end  % Cancel was selected
            cdata = double(ro{i}.read_chip(...
                [aoi(1) aoi(1)+aoi(3)-1],[aoi(2) aoi(2)+aoi(4)-1]));
            signal_sigma = mean(abs(cdata(:)).^2) - noise_sigma;
        else
            samplesize=[1000 1000]; % Exract 1000x1000 array of samples to estimate mean
            subsample=ceil(double([meta{i}.ImageData.NumCols meta{i}.ImageData.NumRows])./samplesize);
            data = abs(single(ro{i}.read_chip([1 meta{i}.ImageData.NumCols],...
                    [1 meta{i}.ImageData.NumRows],subsample)));
            signal_sigma = mean(abs(data(:)).^2) - noise_sigma;
        end
        if isfield(meta{i}, 'Radiometric') && ...
                isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly')
            % Scale to be equal to noise equivalent sigma zero
            signal_sigma = meta{i}.Radiometric.SigmaZeroSFPoly * signal_sigma;
        end
        warning('RGIQE:SignalLevelMeasured',['More uncertainty with ' ...
            'RNIIRS estimates exists when signal level is sampled, ' ...
            'rather than when using a radiometric constant.']);
    else
        signal_sigma = NaN;
    end

    %% Save some diagnostic statistics
    if isfield(meta{i}, 'Radiometric') && ...
            isfield(meta{i}.Radiometric, 'SigmaZeroSFPoly')
        diagnostic.nesz(i) = 10*log10(noise_sigma);
    end
    diagnostic.snr(i) = signal_sigma/noise_sigma;
    diagnostic.sp_resolution(i) = sqrt(meta{i}.Grid.Row.ImpRespWid * ...
        meta{i}.Grid.Col.ImpRespWid);
    diagnostic.ellipicity(i) = max(meta{i}.Grid.Row.ImpRespWid, ...
        meta{i}.Grid.Col.ImpRespWid)/min(meta{i}.Grid.Row.ImpRespWid, ...
        meta{i}.Grid.Col.ImpRespWid);
    if numel(meta)>1
        diagnostic.meta{i} = meta{i};
    else
        diagnostic.meta = meta{i};
    end

    %% Compute bandwidth area
    % The cosine of the slope angle is proper scale factor to project the
    % bandwidth area into the ground plane.
    diagnostic.bandwidth_area(i) = ...
        meta{i}.Grid.Col.ImpRespBW * meta{i}.Grid.Row.ImpRespBW * ...
        cosd(meta{i}.SCPCOA.SlopeAng);
    % Computing area from vertices is much more work, but it may help
    % visualize.
    % Project slant plane spatial frequency bounds to ground plane
    % coords_slant_x = meta{i}.Grid.Col.KCtr + ...
    %     meta{i}.Grid.Col.ImpRespBW * [-1 1 1 -1 -1]/2;
    % coords_slant_y = meta{i}.Grid.Row.KCtr + ...
    %     meta{i}.Grid.Row.ImpRespBW * [1 1 -1 -1 1]/2;
    % ruvect = [meta{i}.Grid.Row.UVectECF.X meta{i}.Grid.Row.UVectECF.Y meta{i}.Grid.Row.UVectECF.Z];
    % cuvect = [meta{i}.Grid.Col.UVectECF.X meta{i}.Grid.Col.UVectECF.Y meta{i}.Grid.Col.UVectECF.Z];
    % coords_slant_3d = ruvect'*coords_slant_y + cuvect'*coords_slant_x;
    % Compute basis vectors for ground plane
    % gpn = wgs_84_norm([meta{i}.GeoData.SCP.ECF.X meta{i}.GeoData.SCP.ECF.Y meta{i}.GeoData.SCP.ECF.Z])';
    % gruvect = project_ground(ruvect, gpn); % Project range vector to ground plane
    % gruvect = gruvect/norm(gruvect);
    % gcuvect = cross(gruvect,gpn);
    % gcuvect = gcuvect/norm(gcuvect);
    % Project spatial frequencies to ground plane
    % coords_ground_3d = coords_slant_3d;
    % coords_ground_x = zeros(size(coords_slant_x));
    % coords_ground_y = zeros(size(coords_slant_y));
    % for j=1:numel(coords_slant_x)
    %     % Planes in ECF
    %     coords_ground_3d(:,j) = project_ground(coords_slant_3d(:,j)',gpn)';
    %     % Convert to 2D basis vectors in ground plane
    %     coords_ground_x(j) = dot(coords_ground_3d(:,j),gcuvect);
    %     coords_ground_y(j) = dot(coords_ground_3d(:,j),gruvect);
    % end
    % This is a way to compute the area of a any polygon from its vertices.
    % diagnostic.bandwidth_area = sum((coords_ground_x(1:(end-1)).*coords_ground_y(2:end)) - ...
    %     (coords_ground_y(1:(end-1)).*coords_ground_x(2:end)))/2;
    % Visualize spatial frequency projections in slant and ground planes
    % figure;
    % title('Spatial frequency bounds 3D');
    % g_plane = gruvect'*[100 100 -100 -100 100] + gcuvect'*[-100 100 100 -100 -100];
    % plot3(coords_slant_3d(1,:),coords_slant_3d(2,:),coords_slant_3d(3,:),...
    %     coords_ground_3d(1,:),coords_ground_3d(2,:),coords_ground_3d(3,:),...
    %     g_plane(1,:),g_plane(2,:),g_plane(3,:));
    % legend('Slant plane bounds', 'Ground plane bounds', 'Ground plane')
    % figure;
    % plot(coords_slant_x, coords_slant_y, coords_ground_x, coords_ground_y);
    % title('Spatial frequency bounds');
    % legend('Slant plane bounds', 'Ground plane bounds')

    %% After 300+ line of setup code and comments, we finally compute quality metric...
    % Shannon-Hartley theorem of channel capacity
    % C = B * log2(1 + (S/N))
    % We propose a 2-dimension information rate based on a 2-dimensional
    % bandwidth.  This 2D bandwidth is an area in a space where each axis
    % is cycles/m, as opposed to 1D distance in Hz (cycles per second).
    % Therefore the units of this metric are bits/m^2.
    image_information_metric(i) = diagnostic.bandwidth_area(i) * log2(1 + diagnostic.snr(i));

    %% Map bits/m^2 into RNIIRS
    if ismember('signal_sigma_override', p.UsingDefaults)
        rniirs(i) = estimate_rniirs(image_information_metric(i));
    else  % If assumed signal is not default, this RNIIRS mapping is not valid
        rniirs(i) = NaN;
    end
end
try
    for i = 1:numel(meta)
        ro{i}.close();
    end
end

end

function rniirs = estimate_rniirs(image_information_metric)
    % The coefficients here were derived by fitting hundreds of SAR
    % datasets with varying parameters to analysts ratings.  This
    % author prefers the information metric but also provides this as
    % well for comparison to this legacy scale.
    % coeffs = [0.3960 3.7555];  % Probably more digits of precision than we really have
    coeffs = [0.4357 3.4761];  % Coefficients were updated by NGA NIQU (Feb 1, 2022) to more accruately represent RNIIRS tables
    rniirs = polyval(coeffs,log2(image_information_metric));
    % For completeness, we handle the very low RNIIRS case where the
    % logarithmic equation would have been negative.  Negative is not
    % allowed in the RNIIRS scale. We compute where a line tangent to
    % the bits/m^2-to-RNIIRS function intersects with the origin.  We
    % will use this tangent line as the mapping for very low RNIIRS,
    % since it approaches zero as information approaches zero and
    % remains non-negative. This way the entire mapping function
    % between bits/m^2 and RNIIRS is continuous and has a continuous
    % derivative.  The function will be logarithmic with respect to the
    % information metric above this point (presumably nearly all data
    % will fit in here) and linear with respect to the metric below
    % this point (extremely low RNIIRS).  Note that since NGA has no
    % precedent or analysts ratings for RNIIRS<1, we are really free to
    % define this however we want.
    linlog_transition = exp(1 - coeffs(2)*log(2)/coeffs(1));
    % Derivative of linear portion (as well as intersection with log portion)
    d_rniirs = coeffs(1)/(log(2)*linlog_transition);
    % This is what the two lines would look line at very low RNIIRS:
    % x=linspace(0,.02,1000);
    % figure; plot(x, polyval(coeffs,log2(x)), x, x*d_rniirs);
    % That was a lot of explanation for a case that will likely rarely
    % be used, wasn't it?
    if image_information_metric<linlog_transition
        rniirs = d_rniirs * image_information_metric;
    end
end

function out = project_ground(point, gpn)
    % Projection of a point along a given direction to a plane is
    % just the intersection of the line defined by that point (l0)
    % and direction (l) and the plane defined by a point in the
    % plane (p0) and the normal (p):
    % l0 - ((l0 - p0).p/(l.p))*l
    % where . represents the dot product.
    % In this specific case, the origin is in the plane, and the projection
    % direction is the same as the plane normal, so this equation is much
    % simplified:
    out = point - dot(point,gpn)*gpn;
end

author: ngageoint


text:
function [fxs,fys,fzs,nxs,nys,nzs] = fdct3d_param(C)
  [m,n,p] = size(C{end}{1});
  nbscales = floor(log2(min([m,n,p])))-2;
  nbdstz_coarse = 8;
  allcurvelets = 0;
  [fxs,fys,fzs,nxs,nys,nzs] = fdct3d_param_mex(m,n,p,nbscales,nbdstz_coarse,allcurvelets);

author: thomaskuestner