text:
function bernstein_poly_01_values_test ( )

%*****************************************************************************80
%
%% BERNSTEIN_POLY_01_VALUES_TEST demonstrates the use of BERNSTEIN_POLY_01_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BERNSTEIN_POLY_01_VALUES_TEST:\n' );
  fprintf ( 1, '  BERNSTEIN_POLY_01_VALUES returns values of \n' );
  fprintf ( 1, '  the Bernstein polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N     K       X      BERNSTEIN(N,K)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, k, x, b ] = bernstein_poly_01_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %6d  %4d  %12f  %24.16f\n', n, k, x, b );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function www=trajt3(para,chrom)

x0=para(1);
y0=para(2);
z0=para(3);
vx0=0;
vy0=0;
vz0=0;
ax0=0;
ay0=0;
az0=0;
x2=para(4);
y2=para(5);
z2=para(6);
vx2=0;
vy2=0;
vz2=0;
ax2=0;
ay2=0;
az2=0;
x1=chrom(1);
y1=chrom(2);
z1=chrom(3);
vx1=chrom(4);
vy1=chrom(5);
vz1=chrom(6);
t1=chrom(7);
t2=chrom(8);

a00=x0;
a01=vx0;
a02=ax0/2;
a03=(4*x1-vx1*t1-4*x0-3*vx0*t1-ax0*t1^2)/t1^3;
a04=(vx1*t1-3*x1+3*x0+2*vx0*t1+ax0*t1^2/2)/t1^4;
ax1=2*a02+6*a03*t1+12*a04*t1^2;

c00=y0;
c01=vy0;
c02=ay0/2;
c03=(4*y1-vy1*t1-4*y0-3*vy0*t1-ay0*t1^2)/t1^3;
c04=(vy1*t1-3*y1+3*y0+2*vy0*t1+ay0*t1^2/2)/t1^4;
ay1=2*c02+6*c03*t1+12*c04*t1^2;

e00=z0;
e01=vz0;
e02=az0/2;
e03=(4*z1-vz1*t1-4*z0-3*vz0*t1-az0*t1^2)/t1^3;
e04=(vz1*t1-3*z1+3*z0+2*vz0*t1+az0*t1^2/2)/t1^4;
az1=2*e02+6*e03*t1+12*e04*t1^2;
%--------------------------------------------------------------------------
b10=x1;
b11=vx1;
b12=ax1/2;
b13=(20*x2-20*x1-(8*vx2+12*vx1)*t2-(3*ax1-ax2)*t2^2)/(2*t2^3);
b14=(30*x1-30*x2+(14*vx2+16*vx1)*t2+(3*ax1-2*ax2)*t2^2)/(2*t2^4);
b15=(12*x2-12*x1-(6*vx2+6*vx1)*t2-(ax1-ax2)*t2^2)/(2*t2^5);

d10=y1;
d11=vy1;
d12=ay1/2;
d13=(20*y2-20*y1-(8*vy2+12*vy1)*t2-(3*ay1-ay2)*t2^2)/(2*t2^3);
d14=(30*y1-30*y2+(14*vy2+16*vy1)*t2+(3*ay1-2*ay2)*t2^2)/(2*t2^4);
d15=(12*y2-12*y1-(6*vy2+6*vy1)*t2-(ay1-ay2)*t2^2)/(2*t2^5);

f10=z1;
f11=vz1;
f12=az1/2;
f13=(20*z2-20*z1-(8*vz2+12*vz1)*t2-(3*az1-az2)*t2^2)/(2*t2^3);
f14=(30*z1-30*z2+(14*vz2+16*vz1)*t2+(3*az1-2*az2)*t2^2)/(2*t2^4);
f15=(12*z2-12*z1-(6*vz2+6*vz1)*t2-(az1-az2)*t2^2)/(2*t2^5);
%--------------------------------------------------------------------------

t=linspace(0,t1,20);

x01=a00+a01*t+a02*t.^2+a03*t.^3+a04*t.^4;
vx01=a01+2*a02*t+3*a03*t.^2+4*a04*t.^3;
ax01=2*a02+6*a03*t+12*a04*t.^2;

y01=c00+c01*t+c02*t.^2+c03*t.^3+c04*t.^4;
vy01=c01+2*c02*t+3*c03*t.^2+4*c04*t.^3;
ay01=2*c02+6*c03*t+12*c04*t.^2;

z01=e00+e01*t+e02*t.^2+e03*t.^3+e04*t.^4;
vz01=e01+2*e02*t+3*e03*t.^2+4*e04*t.^3;
az01=2*e02+6*e03*t+12*e04*t.^2;
%--------------------------------------------------------------------------
t=linspace(0,t2,20);

x12=b10+b11*t+b12*t.^2+b13*t.^3+b14*t.^4+b15*t.^5;
vx12=b11+2*b12*t+3*b13*t.^2+4*b14*t.^3+5*b15*t.^4;
ax12=2*b12+6*b13*t+12*b14*t.^2+20*b15*t.^3;

y12=d10+d11*t+d12*t.^2+d13*t.^3+d14*t.^4+d15*t.^5;
vy12=d11+2*d12*t+3*d13*t.^2+4*d14*t.^3+5*d15*t.^4;
ay12=2*d12+6*d13*t+12*d14*t.^2+20*d15*t.^3;

z12=f10+f11*t+f12*t.^2+f13*t.^3+f14*t.^4+f15*t.^5;
vz12=f11+2*f12*t+3*f13*t.^2+4*f14*t.^3+5*f15*t.^4;
az12=2*f12+6*f13*t+12*f14*t.^2+20*f15*t.^3;
%--------------------------------------------------------------------------
ppx=[x01 x12];ppy=[y01 y12];ppz=[z01 z12];
q1=ppx;q2=ppy;q3=ppz;

vvx=[vx01 vx12];vvy=[vy01 vy12];vvz=[vz01 vz12];
vq1=vvx;vq2=vvy;vq3=vvz;

aax=[ax01 ax12];aay=[ay01 ay12];aaz=[az01 az12];
aq1=aax;aq2=aay;aq3=aaz;

www=[q1;q2;q3;vq1;vq2;vq3;aq1;aq2;aq3];

author: Sable

repo: mcbench-benchmarks


text:
function c = vertcat(varargin)
% function C=vertcat(A,B);
%
% DESCRIPTION
%   Vertical concatenation of polynomial objects.
%
% INPUTS
%   A,B: polynomials
%
% OUTPUTS
%   C:  vertical concatenation of input matrices.
%
% SYNTAX
%   [A; B]
%      Vertical concatenation of polynomial matrices A and B.
%      A and B must have the same number of columns.
%   [A1; A2; A3; ...]
%      Vertical concatenation of several polynomial matrices.
%   C = vertcat(A1,A2,A3,..);
%      Function-call form of vertical concatenation.
%
% See also horzcat

% 6/8/2002: PJS  Initial Coding

if nargin==1
    c = varargin{1};
else
    % Promote a to polynomial
    a = polynomial(varargin{1});
    [nra,nca] = size(a);

    % Promote b to polynomial
    b = polynomial(varargin{2});
    [nrb,ncb] = size(b);

    if isempty(b);
        c = a;
    elseif isempty(a);
        c = b;
    elseif nca == ncb
        % Get Dimensions
        nta = size(a.degmat,1);
        nva = length(a.varname);
        ntb = size(b.degmat,1);
        nvb = length(b.varname);

        if nva==0 && nvb==0
            % Combine constant terms
            ar = combine(a);
            coef1 = reshape(ar.coefficient,[nra nca]);
            br = combine(b);
            coef2 = reshape(br.coefficient,[nrb ncb]);

            % Stack up Coefficients and Form Polynomial
            coefficient = [coef1; coef2];
            c = polynomial(coefficient);
        else
            % Form Degmat, Varname, and Matdim
            adeg = a.degmat;
            bdeg = b.degmat;
            degmat = blkdiag(adeg,bdeg);
            varname = [a.varname(:); b.varname(:)];
            matdim = [nra+nrb nca];

            % Stack up Coefficients
            idx1 = [];
            idx2 = [];
            for i1 = 0:(nca-1);
                idx1 = [idx1 (1:nra)+i1*(nra+nrb)];
                idx2 = [idx2 nra+(1:nrb)+i1*(nra+nrb)];
            end;
            coef1 = spalloc(nta,(nra+nrb)*nca,nnz(a.coefficient));
            coef1(:,idx1) = a.coefficient;
            coef2 = spalloc(ntb,(nra+nrb)*ncb,nnz(b.coefficient));
            coef2(:,idx2) = b.coefficient;
            coefficient = [coef1; coef2];

            % Form Polynomial and combine terms
            chkval = 0; % skip validity check
            c = polynomial(coefficient,degmat,varname,matdim,chkval);
            c = combine(c);
        end

    else
        error('All columns must have the same row dimension')
    end

    if nargin>2
        c = vertcat(c,varargin{3:end});
    end
end

author: yu-jiang

repo: radpbook

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


text:
% Fit a piece-wise linear regression model.
% Here is the model
%
%  X \
%  | |
%  Q |
%  | /
%  Y
%
% where all arcs point down.
% We condition everything on X, so X is a root node. Q is a softmax, and Y is a linear Gaussian.
% Q is hidden, X and Y are observed.

X = 1;
Q = 2;
Y = 3;
dag = zeros(3,3);
dag(X,[Q Y]) = 1;
dag(Q,Y) = 1;
ns = [1 2 1]; % make X and Y scalars, and have 2 experts
dnodes = [2];
onodes = [1 3];
bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'observed', onodes);

IRLS_iter = 10;
clamped = 0;

bnet.CPD{1} = root_CPD(bnet, 1);

if 0
  % start with good initial params
  w = [-5 5];  % w(:,i) is the normal vector to the i'th decisions boundary
  b = [0 0];  % b(i) is the offset (bias) to the i'th decisions boundary

  mu = [0 0];
  sigma = 1;
  Sigma = repmat(sigma*eye(ns(Y)), [ns(Y) ns(Y) ns(Q)]);
  W = [-1 1];
  W2 = reshape(W, [ns(Y) ns(X) ns(Q)]);

  bnet.CPD{2} = softmax_CPD(bnet, 2, w, b,  clamped, IRLS_iter);
  bnet.CPD{3} = gaussian_CPD(bnet, 3, mu, Sigma, W2);
else
  % start with rnd initial params
  rand('state', 0);
  randn('state', 0);
  bnet.CPD{2} = softmax_CPD(bnet, 2, 'clamped', clamped, 'max_iter', IRLS_iter);
  bnet.CPD{3} = gaussian_CPD(bnet, 3);
end

load('/examples/static/Misc/mixexp_data.txt', '-ascii');
% Just use 1/10th of the data, to speed things up
data = mixexp_data(1:10:end, :);
%data = mixexp_data;

%plot(data(:,1), data(:,2), '.')

s = struct(bnet.CPD{2}); % violate object privacy
%eta0 = [s.glim.b1; s.glim.w1]';
eta0 = [s.glim{1}.b1; s.glim{1}.w1]';
s = struct(bnet.CPD{3}); % violate object privacy
W = reshape(s.weights, [1 2]);
theta0 = [s.mean; W]';

%figure(1)
%mixexp_plot(theta0, eta0, data);
%suptitle('before learning')

ncases = size(data, 1);
cases = cell(3, ncases);
cases([1 3], :) = num2cell(data');

engine = jtree_inf_engine(bnet);

% log lik before learning
ll = 0;
for l=1:ncases
  ev = cases(:,l);
  [engine, loglik] = enter_evidence(engine, ev);
  ll = ll + loglik;
end

% do learning
max_iter = 5;
[bnet2, LL2] = learn_params_em(engine, cases, max_iter);

s = struct(bnet2.CPD{2});
%eta2 = [s.glim.b1; s.glim.w1]';
eta2 = [s.glim{1}.b1; s.glim{1}.w1]';
s = struct(bnet2.CPD{3});
W = reshape(s.weights, [1 2]);
theta2 = [s.mean; W]';

%figure(2)
%mixexp_plot(theta2, eta2, data);
%suptitle('after learning')

fprintf('mixexp2: loglik before learning %f, after %d iters %f\n', ll, length(LL2),  LL2(end));

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
%% patchDetach
% Below is a demonstration of the features of the |patchDetach| function

%% Syntax
% |[Fs,Vs]=patchDetach(F,V,scalefactor);|

%% Description
% This function seperates the nodes (if shared) for all faces. If a
% constant or spatially varying scalefactor is provided the faces are
% scaled (around their mean) as well.

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
fontSize=15;
plotColor1=0.25.*ones(1,3);
plotColor2=0.75.*ones(1,3);
edgeWidth=2;
markerSize=25;

%% Example: Seperate and scale faces homogeneously
%Defining geodesic dome triangulation
r=1; %sphere radius
n=2; %Refinements
[F,V,~]=geoSphere(n,r);

%Detach nodes and shrink faces
scaleFactor=0.5;
[Fs,Vs]=patchDetach(F,V,scaleFactor);

%%
%Plotting results

cFigure;
subplot(1,2,1); hold on;
gpatch(F,V,'rw','k',1,edgeWidth);
axisGeom(gca,fontSize);
camlight headlight;
view(2);

subplot(1,2,2); hold on;
gpatch(F,V,'rw','none',0.5);
gpatch(Fs,Vs,'bw','k',1,edgeWidth);
axisGeom(gca,fontSize);
camlight headlight;
view(2);

drawnow;

%% Example: Seperate and shrink faces in a spatially varying way

%%
% Create a spatially varying shrink factor between 0 and 1
scaleFactor=V(:,2);
scaleFactor=mean(scaleFactor(F),2);
scaleFactor=scaleFactor-min(scaleFactor(:));
scaleFactor=scaleFactor./max(scaleFactor(:));

%%
% Detach nodes and shrink faces
[Fs,Vs]=patchDetach(F,V,scaleFactor);

%%
% Plotting results

cFigure;
subplot(1,2,1); hold on;
gpatch(F,V,'rw','k',1,edgeWidth);
axisGeom(gca,fontSize);
camlight headlight;
view(2);

subplot(1,2,2); hold on;
gpatch(F,V,'rw','none',0.5);
gpatch(Fs,Vs,'bw','k',1,edgeWidth);
axisGeom(gca,fontSize);
camlight headlight;
view(2);

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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 out = ndims(fa)
% Number of dimensions
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

out = size(fa);
out = length(out);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function scr_lin = linear_calibrate_scores(scr,key,obj_func,niters,prior,addzeros)
% A function for doing 'cheating' calibration on scores.  It trains
% a linear calibration (scaling and offset) on the scores and then
% applies the calibration to the same scores.
% Inputs:
%   scr: The object (of class Scores) containing the scores to be
%     calibrated.
%   key: The Key indicating target and non-target trials.
%   obj_func: The objective function for the calibration training.  The
%     default is cllr objective.
%   niters: The maximum number of iterations in the calibration
%     training.
%   prior: The effective target prior.
%   addzeros: If true, missing scores (required by key but not
%     present in scr) are added (with value zero).
% Outputs:
%   scr_lin: The calibrated version of the input scores.

assert(nargin==6)
assert(isa(scr,'Scores'))
assert(isa(key,'Key'))
assert(scr.validate())
assert(key.validate())

[tar,non] = scr.get_tar_non(key);
logprint(Logger.Info,'training calibration\n')
[flin,nxe,w_lin] = train_linear_calibration(tar,non,prior,obj_func,niters,true);

logprint(Logger.Info,'calibrating scores\n')
scr_lin = scr.align_with_ndx(key);
scr_lin = scr_lin.transform(flin);

logprint(Logger.Info,'displaying calibration results\n')
logprint(Logger.Info,'weights: scaling %f, offset %f\n',w_lin(1),w_lin(2))
res = Results(scr_lin,key);
logprint(Logger.Info,'norm_act_dcf = %g, norm_min_dcf = %g, prbep = %g, nxe = %g\n',res.get_norm_act_dcf(prior),res.get_norm_min_dcf(prior),res.get_prbep(),nxe);

if addzeros
    scr_lin = scr_lin.set_missing_to_value(key,0.0);
end

assert(scr_lin.validate())

end

author: nesl


text:
 function ob = Gblur2(mask, varargin)
%function ob = Gblur(mask, options)
%|
%| Construct Gblur object for image restoration.
%|
%| See Gblur_test.m for example usage.
%|
%| in
%|	mask	size(image)	logical array of object support.
%|
%| options
%|	'chat'		verbose printing of debug messages
%|	'psf'		point spread function (aka impulse response)
%|	'type'		type of blur:
%|				'conv,same'	usual case (default)
%|				'conv,per'	(periodic end conditions)
%|				'fft,same'	(periodic end conditions)
%|				todo: allow replicated end conditions!
%|				'imfilter,same'	(requires image toolbox)
%|				'imfilter,circ'	(requires image toolbox)
%|				'imfilter,mirror' (requires image toolbox)
%|					(do not use - adjoint fails)
%|				'sparse'	todo: return sparse matrix
%|	'imfilter_options'	options to 'imfilter' (if used)
%|
%| out
%|	ob [nd np]	np = sum(mask(:)), so it is already "masked"
%|			nd = prod(size(mask)) for 'conv,same' type
%|
%| Copyright 2005-4-22, Jeff Fessler, University of Michigan

if nargin == 1 && streq(mask, 'test'), Gblur_test, return, end
if nargin < 1, help(mfilename), error(mfilename), end

arg.mask = mask;

% option defaults
arg.chat = 0;
arg.psf = 1; % identity matrix
arg.type = 'conv,same';
arg.class = 'fatrix2';
arg.imfilter_options = {};

% options specified by name/value pairs
arg = vararg_pair(arg, varargin);

if ndims(arg.psf) ~= ndims(mask), error 'psf dim mismatch', end
if any(size(arg.psf) > size(mask)), error 'psf too large', end

if streq(arg.type, 'imfilter,circ') && exist('imfilter') ~= 2
	warn 'no imfilter (no image processing toolbox) so using fft'
	arg.type = 'fft,same';
end

if streq(arg.type, 'imfilter,mirror') && exist('imfilter') ~= 2
	fail 'no imfilter (no image processing toolbox)'
end

if streq(arg.type, 'imfilter,same') && exist('imfilter') ~= 2
	warn 'no imfilter (no image processing toolbox) so using slower convn'
	arg.type = 'conv,same';
end

arg.psf = Gblur_mask_psf_odd(arg.psf);

psf_flip = conj(flipdims2(arg.psf, 'odd', 1));
does_many = false;

switch arg.type

case 'conv,same'
	arg.odim = size(mask);
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) convn(x, arg.psf, 'same');
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			convn(y, arg.psf_flip, 'same'));
	does_many = true;

case 'conv,per'
	arg.odim = size(mask);
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) ir_conv(x, arg.psf, 'per', true);
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			ir_conv(y, arg.psf_flip, 'per', true));

case 'fft,same'
	arg.odim = size(mask);
	[arg, handle_forw, handle_back] = Gblur_setup_fft_same(arg);

case 'imfilter,circ'
	arg.odim = size(mask);
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) imfilter(x, arg.psf, 'conv', 'same', ...
		'circular', arg.imfilter_options{:});
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			imfilter(y, arg.psf_flip, 'conv', 'same', ...
				'circular', arg.imfilter_options{:}));

case 'imfilter,mirror'
	arg.odim = size(mask);
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) imfilter(x, arg.psf, 'conv', 'same', ...
		'symmetric', arg.imfilter_options{:});
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			imfilter(y, arg.psf_flip, 'conv', 'same', ...
				'symmetric', arg.imfilter_options{:}));

case 'imfilter,same'
	arg.odim = size(mask);
	arg.psf_flip = psf_flip;
	handle_forw = @(arg,x) imfilter(x, arg.psf, 'conv', 'same', ...
		arg.imfilter_options{:});
	handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			imfilter(y, arg.psf_flip, 'conv', 'same', ...
				arg.imfilter_options{:}));

case 'sparse'
	ob = Gblur_sparse(arg.psf, arg.idim);
	ob = ob(:,arg.mask(:));
	return

otherwise
	error 'unknown blur type'
end

switch arg.class

case 'Fatrix'

	arg.nd = prod(arg.odim);
	arg.np = sum(mask(:));
	dim = [arg.nd arg.np]; % trick: make it masked by default!

	ob = Fatrix(dim, arg, 'caller', 'Gblur', ...
		'abs', @Gblur_abs, ...
		'forw', @Gblur_forw_Fatrix, 'back', @Gblur_back_Fatrix);

case 'fatrix2'

	ob = fatrix2('idim', arg.odim, 'arg', arg, 'does_many', does_many, ...
		'abs', @Gblur_abs, 'forw', handle_forw, 'back', handle_back);

otherwise
	fail('bug')
end

% Gblur_mask_psf_odd(psf)
% having an odd-sized psf facilites adjoint
function psf = Gblur_mask_psf_odd(psf)
nd = ndims(psf);
for id=1:nd
	sz = size(psf, id);
	if ~mod(sz, 2) % even
		warn('size(psf, %d) = %d is even; appending 0', id, sz)
		tmp = size(psf);
		tmp(id) = 1;
		z = zeros(tmp);
		%psf = cat(id, psf, z);
		psf = cat(id, z, psf);
	end
end

% Gblur_sparse()
% a_ij = b[ n(i) - n(j), m(i) - m(j) ]
% n(i) = (i-1) mod N
% m(i) = floor((i-1) / N)
function sp = Gblur_sparse(psf, idim)
if numel(idim) ~= 2
	fail 'not done'
end
fail 'todo: not done'
%sp = sparse(i, j, s, length(ii), size(ob.arg.G, 2));

% Gblur_abs(): |A| for abs(A)
function ob = Gblur_abs(ob)
ob.arg.psf = abs(ob.arg.psf);
if isfield(ob.arg, 'psf_flip')
	ob.arg.psf_flip = abs(ob.arg.psf_flip);
end
if isfield(ob.arg, 'psf_fft') % fixed 2012-04-08
	ob.arg.psf_fft = Gblur_setup_fft_same_fft(ob.arg.psf, ob.mask);
end

% Gblur_setup_fft_same()
function [arg, handle_forw, handle_back] = Gblur_setup_fft_same(arg)

arg.psf_fft = Gblur_setup_fft_same_fft(arg.psf, arg.mask);

handle_forw = @(arg,x) ifftn(fftn(x) .* arg.psf_fft);
handle_back = @(arg,y) fatrix2_maskit(arg.mask, ...
			ifftn(fftn(y) .* conj(arg.psf_fft)));

% Gblur_setup_fft_same_fft()
function psf_fft = Gblur_setup_fft_same_fft(psf, mask)

% % put psf in center of array
% tmp = zeros(size(mask));
% switch ndims(psf)
% case 2
% 	[n1 n2] = size(mask);
% 	[p1 p2] = size(psf);
% %	h1 = (size(psf,1)-1)/2;
% %	h2 = (size(psf,2)-1)/2;
% %	i1 = floor(n1/2) + 1 + [-h1:h1]
% %	i2 = floor(n2/2) + 1 + [-h2:h2]
% 	i1 = floor(n1/2) + 1 + (1:p1) - floor((p1+1)/2);
% 	i2 = floor(n2/2) + 1 + (1:p2) - floor((p2+1)/2);
% 	tmp(i1,i2) = psf;
% case 3
% 	if any(~rem(size(psf),2)), error 'psf size must be odd', end
% 	h1 = (size(psf,1)-1)/2;
% 	h2 = (size(psf,2)-1)/2;
% 	h3 = (size(psf,3)-1)/2;
% 	i1 = floor(size(mask,1)/2) + 1 + [-h1:h1];
% 	i2 = floor(size(mask,2)/2) + 1 + [-h2:h2];
% 	i3 = floor(size(mask,3)/2) + 1 + [-h3:h3];
% 	tmp(i1,i2,i3) = psf;
% otherwise
% 	error 'only 2d & 3d done'
% end
% psf_fft = fftn(fftshift(tmp));

% pad psf with 0, and make size(psf)==size(mask)
sz_mask = size(mask);
sz_psf = ones(1,ndims(mask));
for ii=1:length(sz_psf), sz_psf(ii) = size(psf,ii); end
if any(sz_mask>sz_psf)
    psf = padarray(psf, sz_mask-sz_psf, 'post');
end
center = ceil(sz_psf/2-1/2) + 1; % or center = floor(sz_psf/2) + 1;
psf_fft = fftn( circshift(psf, 1-center) ); % by yusheng

% Gblur_forw_Fatrix(): y = A * x
function y = Gblur_forw_Fatrix(arg, x)

[x ei] = embed_in(x, arg.mask, arg.np);

switch arg.type
case 'conv,same'
	y = convn(x, arg.psf, 'same');
case 'conv,per'
	y = ir_conv(x, arg.psf, 'per', true);
case 'fft,same'
	if ndims(x) > ndims(arg.mask);
		y = fft_conv_multi(x, arg.psf_fft);
	else
		y = ifftn(fftn(x) .* arg.psf_fft);
	end
case 'imfilter,circ'
	y = imfilter(x, arg.psf, 'conv', 'same', ...
		'circular', arg.imfilter_options{:});
case 'imfilter,same'
	y = imfilter(x, arg.psf, 'conv', 'same', arg.imfilter_options{:});
otherwise
	error 'bug'
end

y = ei.shape(y);

% Gblur_back_Fatrix(): x = A' * y
function x = Gblur_back_Fatrix(arg, y)

[y eo] = embed_out(y, arg.odim);

switch arg.type
case 'conv,same'
	x = convn(y, arg.psf_flip, 'same');
case 'conv,per'
	x = ir_conv(y, arg.psf_flip, 'per', true);
case 'fft,same'
	if ndims(y) > ndims(arg.mask);
		x = fft_conv_multi(y, conj(arg.psf_fft));
	else
		x = ifftn(fftn(y) .* conj(arg.psf_fft));
	end
case 'imfilter,circ'
	x = imfilter(y, arg.psf_flip, 'conv', 'same', ...
		'circular', arg.imfilter_options{:});
case 'imfilter,same'
	x = imfilter(y, arg.psf_flip, 'conv', 'same', arg.imfilter_options{:});
otherwise
	error 'bug'
end

x = eo.shape(x, arg.mask, arg.np);

% fft_conv_multi()
function y = fft_conv_multi(x, H)
dim = size(x);
y = zeros(size(x));
y = [];
for ii=1:dim(end)
	tmp = ifftn(fftn(stackpick(x, ii)) .* H);
	y = cat(length(dim), y, tmp); % slow: keeps enlarging y.
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function xs = searchOptim(x0,gpstruct,LB,UB,Scale,optimState,options)
%SEARCHOPTIM Search step by local maximization of expected improvement.

if nargin < 1
    xs = 'ei';
    return;
end

optoptions = optimset('Display','off','GradObj','on','DerivativeCheck','off',...
    'TolX',optimState.TolMesh,'TolFun',options.TolFun);

[ymui,ys2i,fmui,fs2i,~,post] = gp(gpstruct.hyp(1),gpstruct.inf,gpstruct.mean,gpstruct.cov, ...
    gpstruct.lik,gpstruct.x,gpstruct.y,x0);
gpstruct.post = post;

idx = 1;
funcCount = 0;
while funcCount < options.Nsearch && idx <= 10
    x0 = x0 + 0.1*Scale*randn(size(x0));
    try
        optoptions.MaxFunEval = options.Nsearch - funcCount;
        % [xs(idx,:),~,~,output] = fmincon(@(x_) LowestUpperBound(x_,gpstruct,Scale),x0,[],[],[],[],LB,UB,[],optoptions);
        [xs(idx,:),~,~,output] = fmincon(@(x_) NegExpectedImprovement(x_,optimState.ftarget,gpstruct),x0,[],[],[],[],LB,UB,[],optoptions);
        funcCount = funcCount + output.funcCount;
    catch
        warning('ah');
        xs(idx,:) = x0;
    end
    idx = idx + 1;
end

end

%--------------------------------------------------------------------------
function [y,dy] = NegExpectedImprovement(xi,target,gpstruct,scale)
%NEGEXPECTEDIMPROVEMENT Return NFEI

if ~isempty(gpstruct.x0)
    xi = bsxfun(@minus,xi,gpstruct.x0);
    gpstruct.x = bsxfun(@minus,gpstruct.x,gpstruct.x0);
end

Nhyp = length(gpstruct.hyp);
hypw = zeros(Nhyp,1);

for i = 1:Nhyp; hypw(i) = gpstruct.hypweight(i); end
[ymu,ys2,fmu,fs2,dymu,dys2,dfmu,dfs2] = gpgrad(xi,gpstruct,'central');
fs = sqrt(fs2);
ys = sqrt(ys2);

gammaz = (target - fmu)./fs;
fpi = 0.5*erfc(-gammaz/sqrt(2));    % Probability of improvement
y = -fs.*(gammaz.*fpi + exp(-0.5*gammaz.^2)/sqrt(2*pi));
y = sum(bsxfun(@times,hypw,y),1);

if any(isnan(y) | isinf(y))
    y = 0;
    dy = zeros(1,size(xi,2));
    return;
end

dfs = 0.5*dfs2./fs;
dgammaz = -(dfmu.*fs + (target - fmu).*dfs)./fs2;
dfpi = -0.5*dgammaz/sqrt(2)*(-2*exp(-gammaz.^2/2)/sqrt(pi));

dy = -(dfs.*gammaz.*fpi + dgammaz.*fs.*fpi + dfpi.*gammaz.*fs) ...
    + (fs.*gammaz.*dgammaz - dfs).*exp(-0.5*gammaz.^2)/sqrt(2*pi);
dy = sum(bsxfun(@times,hypw,dy),1);

end

%--------------------------------------------------------------------------
function [y,dy] = LowestUpperBound(xi,gpstruct,scale,kappa)
%LOWESTUPPERBOUND Return LUB

if nargin < 4 || isempty(kappa); kappa = 1; end

if ~isempty(gpstruct.x0)
    xi = bsxfun(@minus,xi,gpstruct.x0);
    gpstruct.x = bsxfun(@minus,gpstruct.x,gpstruct.x0);
end

Nhyp = length(gpstruct.hyp);
hypw = zeros(Nhyp,1);

for i = 1:Nhyp; hypw(i) = gpstruct.hypweight(i); end
if nargout > 1
    [ymu,ys2,fmu,fs2,dymu,dys2,dfmu,dfs2] = gpgrad(xi,gpstruct,'central');
else
    [ymu,ys2,fmu,fs2] = dgp(xi,gpstruct);
end
fs = sqrt(fs2);
% ys = sqrt(ys2);

y = fmu + kappa.*fs;
y = sum(bsxfun(@times,hypw,y),1);

if any(isnan(y) | isinf(y))
    y = 0;
    dy = zeros(1,size(xi,2));
    return;
end

if nargout > 1
    dy = dfmu + 0.5*kappa.*dfs2./fs;
    dy = sum(bsxfun(@times,hypw,dy),1);
    if any(isnan(dy) | isinf(dy)); dy = zeros(1,size(xi,2)); end
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [featureVector,calcTimes,calcQuality] = TS_CalculateFeatureVector(tsStruct,doParallel,Operations,MasterOperations,codeSpecial,howVocal)
% TS_CalculateFeatureVector	Compute a feature vector from an input time series
%
%---INPUTS:
% tsStruct, either (i) a vector of time-series data to compute
% 				or (ii) an hctsa-style table of time-series data
% doParallel, (binary) whether to compute the features using parallel processing
% Operations, an hctsa-style table of Operations
% MasterOperations, an hctsa-style table of MasterOperations
% codeSpecial, whether to code special values with quality labels (for mySQL database)
% 				this makes sure the featureVector is all real numbers, and any
%				special values (like NaNs, Infs, etc.) are coded with corresponding
% 				labels in calcQuality.
% 				codeSpecial = false [default]: special values are kept in the
% 							feature vector, and any errors are coded as NaN.
% 				codeSpecial = true: featureVector is all real numbers, and is set to
% 							     zero where any special-valued outputs occur.
% howVocal, {'fast', 'minimal', or 'full'}: how to give user feedback on the computation.
%
%---OUTPUTS:
% featureVector, the feature vector obtained by running MasterOperations and
% 					retrieving all features defined in the Operations structure
% 					array on the time series data given in tsStruct.
% calcTimes, corresponding calculation times for each feature in featureVector
% calcQuality, quality labels of each calculation (e.g., coding for NaNs, Infs, etc.)
%
%---USAGE:
% Quickly compute a feature vector of a 500-long random time-series using the
% default operation library using parallel processing:
% >> features = TS_CalculateFeatureVector(randn(500,1));

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check Inputs
%-------------------------------------------------------------------------------
if isnumeric(tsStruct)
	% Provided data only without metadata:
	tsData = tsStruct;
	tsStruct = struct('Name','Input Timeseries', ...
						'Data',tsData, ...
						'ID',1, ...
						'Length',length(tsData));
elseif istable(tsStruct)
	tsStruct = table2struct(tsStruct);
end

if nargin < 2
    fprintf(1,'Computing features in parallel by default.\n');
	doParallel = true;
end

if nargin < 3 || isempty(Operations)
    fprintf(1,'Importing the default hctsa set of time-series features!\n');
    theINPfile = 'INP_ops_hctsa.txt';
    Operations = SQL_Add('ops',theINPfile,false,false);
elseif ischar(Operations)
    theINPfile = Operations;
    Operations = SQL_Add('ops',theINPfile,false,false);
end
if isnumeric(Operations)
	error('Provide an input file or a structure array of Operations.');
end

if nargin < 4 || isempty(MasterOperations)
	% Use the default library:
	MasterOperations = SQL_Add('mops','INP_mops_hctsa.txt',false,false);
	% Need to link operations to masters if not already supplied:
	[Operations,MasterOperations] = TS_LinkOperationsWithMasters(Operations,MasterOperations);
end

% Whether to code up special-valued outputs
if nargin < 5
	codeSpecial = false;
end

% Whether to give information out to screen
if nargin < 6
	howVocal = 'minimal';
end

%-------------------------------------------------------------------------------
% Check Statistics toolbox is available (needed throughout hctsa, including for
% z-scoring)
%-------------------------------------------------------------------------------
BF_CheckToolbox('statistics_toolbox');

% ------------------------------------------------------------------------------
%% Open parallel processing worker pool
% ------------------------------------------------------------------------------
if doParallel
    % Check that a parallel worker pool is open (if not, attempt to initiate it):
	doParallel = TS_InitiateParallel(false);
end
% Tell the user about it:
if strcmp(howVocal,'full')
    if doParallel
        fprintf(1,['Computation will be performed across multiple cores' ...
                ' using Matlab''s Parallel Computing Toolbox.\n']);
    else % use single-threaded for loops
        fprintf(1,'Computations will be performed serially without parallelization.\n');
    end
end

% --------------------------------------------------------------------------
%% Basic checking on the data
% --------------------------------------------------------------------------
% (x a N x 1 column vector)
x = tsStruct.Data;
if size(x,2) ~= 1
	if size(x,1) == 1
		fprintf(1,['***** The time series %s is a row vector. Not sure how it snuck',...
								' through the processing cracks, but I need' ...
								' a column vector...\n'],tsStruct.Name);
		fprintf(1,'I''ll transpose it for you for now....\n');
		x = x';
	else
		fprintf(1,'******************************************************************************************\n');
		error('ERROR WITH ''%s'' -- is it multivariate or something weird? Skipping!\n',tsStruct.Name);
	end
end
% Basics checks on data type:
if ~isa(x,'numeric')
    error('The time series provided must be numerical data')
elseif isa(x,'integer')
    error('Methods in hctsa are generally not well-suited to integer-valued time-series data. Convert to double.')
elseif isa(x,'single')
    warning('Your data is provided in single precision. Converting to double for compatibility with hctsa methods.')
    x = double(x);
end
% (x contains no special values)
if ~all(isfinite(x))
	error('ERROR WITH ''%s'' -- contains non-finite values',tsStruct.Name);
end
if all(x==x(1))
	warning('Data are a constant -- there is no information to derive from the time series; many features will fail')
end

% --------------------------------------------------------------------------
%% Display information
% --------------------------------------------------------------------------
numCalc = height(Operations); % Number of features to calculate
if numCalc == 0
	error('Nothing to calculate :-/')
end

if strcmp(howVocal,'full')
    fprintf(1,'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=\n');
    fprintf(1,'Preparing to calculate %s\nts_id = %u, N = %u samples\nComputing %u features.\n', ...
                    tsStruct.Name,tsStruct.ID,tsStruct.Length,numCalc);
    fprintf(1,'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=\n\n');
end

% Initialize variables:
featureVector = zeros(numCalc,1); % Output of each operation
calcQuality = zeros(numCalc,1); % Quality of output from each operation
calcTimes = nan(numCalc,1); % Calculation time for each operation

% --------------------------------------------------------------------------
%% Pre-Processing
% --------------------------------------------------------------------------
% x_z is a z-scored transformation of the time series:
x_z = zscore(x);

% So we now have the raw time series x and the z-scored time series x_z.
% Operations take these as inputs.

fullTimer = tic;

% --------------------------------------------------------------------------
%% Evaluate all master operation functions (maybe in parallel)
% --------------------------------------------------------------------------
% Because of parallelization, we have to evaluate all the master functions *first*
% Check through the metrics to determine which master functions are relevant for this run

% Put the output from each Master operation in an element of MasterOutput
MasterOutput = cell(height(MasterOperations),1); % Ouput structures
MasterCalcTime = zeros(height(MasterOperations),1); % Calculation times for each master operation

Master_IDs_calc = unique(Operations.MasterID); % Master_IDs that need to be calculated
Master_ind_calc = arrayfun(@(x)find(MasterOperations.ID==x,1),Master_IDs_calc); % Indicies of MasterOperations that need to be calculated
numMopsToCalc = length(Master_IDs_calc); % Number of master operations to calculate

% Index sliced variables to minimize the communication overhead in the parallel processing
par_MasterOpCodeCalc = MasterOperations.Code(Master_ind_calc); % String array of strings of Code to evaluate
par_mop_ids = MasterOperations.ID(Master_ind_calc); % mop_id for each master operation

% Store in temporary variables for parfor loop then map back later
MasterOutput_tmp = cell(numMopsToCalc,1);
MasterCalcTime_tmp = zeros(numMopsToCalc,1);

% ----
% Evaluate all the master operations
% ----
if strcmp(howVocal,'full')
    fprintf(1,'Evaluating %u master operations...\n',length(Master_IDs_calc));
end
TimeSeries_i_ID = tsStruct.ID; % Make a PARFOR-friendly version of the ID
masterTimer = tic;
if doParallel
    % PARFOR Loop (parallel)
	parfor jj = 1:numMopsToCalc
		[MasterOutput_tmp{jj},MasterCalcTime_tmp(jj)] = ...
			TS_ComputeMasterLoop(x,x_z,par_MasterOpCodeCalc{jj}, ...
				par_mop_ids(jj),numMopsToCalc,howVocal,TimeSeries_i_ID,jj);
	end
else
    % Normal FOR Loop (serial)
    if strcmp(howVocal,'minimal')
        BF_ProgressBar('new');
    end
	for jj = 1:numMopsToCalc
		[MasterOutput_tmp{jj},MasterCalcTime_tmp(jj)] = ...
			TS_ComputeMasterLoop(x,x_z,par_MasterOpCodeCalc{jj}, ...
				par_mop_ids(jj),numMopsToCalc,howVocal,TimeSeries_i_ID,jj);

        if strcmp(howVocal,'minimal')
            BF_ProgressBar(jj/numMopsToCalc);
        end
	end
    if strcmp(howVocal,'minimal')
        BF_ProgressBar('close');
    end
end

% Map from temporary versions to the full versions:
MasterOutput(Master_ind_calc) = MasterOutput_tmp;
MasterCalcTime(Master_ind_calc) = MasterCalcTime_tmp;

if strcmp(howVocal,'full')
    fprintf(1,'%u time-series features evaluated in %s///\n\n',...
                        numMopsToCalc,BF_TheTime(toc(masterTimer)));
end
clear('masterTimer')

% --------------------------------------------------------------------------
%% Assign all the results to the corresponding operations
% --------------------------------------------------------------------------
% Set sliced version of matching indicies across the range toCalc
% Indices of MasterOperations corresponding to each Operation (i.e., each index of toCalc)
MasterOp_ind = arrayfun(@(x)find(MasterOperations.ID==x,1),Operations.MasterID);

for jj = 1:numCalc
	try
		[featureVector(jj),calcQuality(jj),calcTimes(jj)] = ...
			TS_ComputeOpLoop(MasterOutput{MasterOp_ind(jj)}, ...
				MasterCalcTime(MasterOp_ind(jj)), ...
				MasterOperations.Label{MasterOp_ind(jj)}, ... % Master label
				Operations.CodeString{jj}); % Code string for each operation to calculate (i.e., in toCalc)
	catch
		fprintf(1,'---Error with %s\n',Operations.CodeString{jj});
		if (MasterOperations.ID(MasterOp_ind(jj)) == 0)
			error(['The operations database is corrupt: there is no link ' ...
					'from ''%s'' to a master code'],Operations.CodeString{jj});
		else
			fprintf(1,'Error retrieving element %s from %s.\n', ...
				Operations.CodeString{jj},MasterOperations.Label{MasterOp_ind(jj)});
		end
	end
end

% --------------------------------------------------------------------------
%% Code special values:
% --------------------------------------------------------------------------
% (*) Errorless calculation: calcQuality = 0, featureVector = <real number>

% (*) Output = NaN: calcQuality = 2, set featureVector = 0
RR = isnan(featureVector); % NaN
if any(RR)
	calcQuality(RR) = 2;
	if codeSpecial
		featureVector(RR) = 0;
	end
end

% (*) Output = Inf: calcQuality = 3, set featureVector = 0
RR = (isinf(featureVector) & featureVector > 0); % Inf
if any(RR)
	calcQuality(RR) = 3;
	if codeSpecial
		featureVector(RR) = 0;
	end
end

% (*) Output = -Inf: calcQuality = 4, set featureVector = 0
RR = (isinf(featureVector) & featureVector < 0);
if any(RR)
	calcQuality(RR) = 4;
	if codeSpecial
		featureVector(RR) = 0;
	end
end

% (*) Output is a complex number: calcQuality = 5, set featureVector = 0
RR = (imag(featureVector)~=0);
if any(RR)
	calcQuality(RR) = 5;
	if codeSpecial
		featureVector(RR) = 0;
	end
end

% (*) Fatal error: calcQuality = 1, featureVector = 0
% 		(this is done already in the code above)
if ~codeSpecial
	% If you want special values in your feature vector (~codeSpecial), then
	% we need to put errors back into the feature vector as NaNs:
	RR = (calcQuality==1);
	featureVector(RR) = NaN;
end

% --------------------------------------------------------------------------
%% Calculation complete: display information about this time series calculation
% --------------------------------------------------------------------------

% Calculate statistics for writing to file/screen
% The number of calculated operations that returned real outputs without errors, numGood:
numGood = sum(calcQuality == 0);
% The number of errors encountered, numErrors:
numErrors = sum(calcQuality == 1);
% The number of other special outputs, numSpecial:
numSpecial = sum(calcQuality > 1);

switch howVocal
case 'full'
    fprintf(1,'********************************************************************\n');
    fprintf(1,'; ; ; : : : : ; ; ; ;   %s    ; ; ; ; : : : ; ; ;\n',datestr(now));
    fprintf(1,'Calculation complete for %s (ts_id = %u, N = %u)\n', ...
                        tsStruct.Name,tsStruct.ID,tsStruct.Length);
    fprintf(1,'%u real-valued outputs, %u errors, %u special-valued outputs stored.\n',...
                        numGood,numErrors,numSpecial);
    fprintf(1,'All %u calculations for this time series took %s.\n',numCalc,BF_TheTime(toc(fullTimer),1));
    fprintf(1,'********************************************************************\n');
case 'minimal'
    fprintf(1,'Calculation complete for %s (ts_id = %u, N = %u) in %s\n', ...
                        tsStruct.Name,tsStruct.ID,tsStruct.Length,BF_TheTime(toc(fullTimer),1));
    fprintf(1,'%u real-valued outputs, %u errors, %u special-valued outputs stored.\n\n',...
                        numGood,numErrors,numSpecial);
case 'fast'
    % (none)
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function compiler()

mex -g Seg2_Stim.c levelset.c Img_new.c gauss_seg.c gauss_vect_part_seg.c gauss_vect_seg.c gamma_seg.c ker_part_seg.c ker_seg.c levelutil.c mean_clust_seg.c mean_part_seg.c mean_seg.c mean_vect_seg.c;


text:
function z = mtimes( x, y, oper )

%   Disciplined convex programming information for MTIMES:
%      True matrix multiplications Z = X * Y---that is, where neither X
%      nor Y is a scalar---require both multiplications and additions.
%      For example, element (i,j) of Z is given by
%         X(i,1) * Y(1,j) + X(i,2) * Y(2,j) + ... + X(i,k) * Y(k,j)
%      Therefore, matrix multiplications must satisfy *both* the
%      "no-product" rule of multiplication and the "same curvature" rule
%      of addition. See the help for CVX/TIMES and CVX/PLUS,
%      respectively, for individual descriptions of these rules.
%
%      An exception is made to these general rules for quadratic forms.
%      That is, two affine expressions may be multiplied together if the
%      result can immediately be verified as a convex quadratic form.
%      For example, the construction
%         variable x(n)
%         x' * Q * x  <= 1;
%      would be permitted if Q is *constant* and positive semidefinite.
%
%   Disciplined geometric programming information for TIMES:
%      As mentioned above, true matrix multiplies Z = X * Y require both
%      multiplications and additions. Since only log-convex terms can be
%      summed, both X and Y must be elementwise log-convex/affine.

persistent remap
if nargin < 3, oper = 'times'; end

%
% Check sizes
%

sx = size( x );
sy = size( y );
if all( sx == 1 ) || all( sy == 1 ),
    z = feval( oper, x, y );
    return
elseif length( sx ) > 2 || length( sy ) > 2,
    error( 'Input arguments must be 2-D.' );
elseif sx( 2 ) ~= sy( 1 ),
    error( 'Inner matrix dimensions must agree.' );
else
    sz = [ sx( 1 ), sy( 2 ) ];
end
nz = prod( sz );

%
% Check expression types
%

if cvx_isconstant( x ),

    xC = cvx_constant( x );
    if nnz( isnan( xC ) ),
        error( 'Disciplined convex programming error:\n    Invalid numeric values (NaNs) may not be used in CVX expressions.', 1 ); %#ok
    elseif cvx_isconstant( y ),
        yC = cvx_constant( y );
        if nnz( isnan( yC ) ),
            error( 'Disciplined convex programming error:\n    Invalid numeric values (NaNs) may not be used in CVX expressions.', 1 ); %#ok
        end
        z = feval( [ 'm', oper ], xC, yC );
        if nnz( isnan( z ) ),
            error( 'Disciplined convex programming error:\n    This expression produced one or more invalid numeric values (NaNs).', 1 ); %#ok
        end
        z = cvx( z );
        return
    elseif isequal( oper, 'rdivide' ),
        error( 'Disciplined convex programming error:\n    Matrix divisor must be constant.', 1 ); %#ok
    end
    yA   = cvx_basis( y );
    laff = true;
    cnst = false;
    raff = false;
    quad = false;
    posy = false;
    vpos = false;

elseif cvx_isconstant( y ),

    yC = cvx_constant( y );
    if nnz( isnan( yC ) ),
        error( 'Disciplined convex programming error:\n    Invalid numeric values (NaNs) may not be used in CVX expressions.', 1 ); %#ok
    elseif isequal( oper, 'ldivide' ),
        error( 'Disciplined convex programming error:\n    Matrix divisor must be constant.', 1 ); %#ok
    end
    xA   = cvx_basis( x );
    raff = true;
    laff = false;
    quad = false;
    cnst = false;
    posy = false;
    vpos = false;

else

    if isempty( remap ),
        remap_0 = cvx_remap( 'zero' );
        remap_1 = cvx_remap( 'nonzero', 'complex' );
        temp    = ~( remap_0 | remap_1 );
        remap_2 = cvx_remap( 'affine' ) & temp;
        remap_4 = cvx_remap( 'log-convex' );
        remap_5 = cvx_remap( 'log-concave' ) & ~remap_4;
        remap_4 = remap_4 & temp;
        remap_3 = cvx_remap( 'valid' ) & ~( remap_0 | remap_1 | remap_2 | remap_4 | remap_5 );
        remap   = remap_1 + 2 * remap_2 + 3 * remap_3 + 4 * remap_4 + 5 * remap_5 - cvx_remap( 'invalid' );
    end
    vx = remap( cvx_classify( x ) );
    vy = remap( cvx_classify( y ) );
    xA = cvx_basis( x );
    yA = cvx_basis( y );
    xC = cvx_reshape( xA( 1, : ), sx );
    yC = cvx_reshape( yA( 1, : ), sy );
    vx = reshape( vx, sx );
    vy = reshape( vy, sy );
    cx = xC ~= 0;
    cy = yC ~= 0;
    ax = vx == 2;
    ay = vy == 2;
    px = vx == 4;
    py = vy == 4;
    gx = vx == 5;
    gy = vy == 5;
    quad = +ax * +ay;
    if nnz( quad ) ~= 0,
        if length( quad ) ~= 1,
            error( 'Disciplined convex programming error:\n    Only scalar quadratic forms can be specified in CVX\n.', 1 ); %#ok
        else
            cx = cx & ~ax;
            cy = cy & ~ay;
            xC( ax ) = 0;
            yC( ay ) = 0;
        end
    end
    cnst = +cx * +cy; %#ok
    laff = +cx * +( vy > 1 );
    raff = +( vx > 1 ) * cy;
    posy = +px * +py;
    vpos = +gx * +gy;
    if nnz( raff ) ~= 0,
        raff = true;
        cnst = false;
    elseif nnz( laff ) ~= 0,
        laff = true;
        cnst = false;
    else
        laff = false;
        raff = false;
        cnst = true;
    end
    othr = +( vx > 1 | vx < 0 ) * +( vy > 1 | vy < 0 ) - quad - posy - vpos;
    if nnz( othr ) ~= 0,
        error( 'Disciplined convex programming error:\n    Cannot perform the operation {%s}*{%s}', cvx_class( x ), cvx_class( y ) );
    end
    quad = nnz( quad ) ~= 0;
    posy = nnz( posy ) ~= 0;

end

first = true;

if cnst,
    switch oper,
    case 'ldivide', z2 = xC \ yC;
    case 'rdivide', z2 = xC / yC;
    otherwise,      z2 = xC * yC;
    end
    if first, z = z2; first = false; else z = z + z2; end %#ok
end

if raff,
    % everything * constant
    nA = size( xA, 1 );
    z2 = cvx_reshape( xA, [ nA * sx( 1 ), sx( 2 ) ] );
    if issparse( z2 ),
        tt = any( z2, 2 );
        if cvx_use_sparse( size( z2 ), nnz( tt ) * sx( 2 ), isreal( z2 ) & isreal( yC ) ),
            z2 = z2( tt, : );
        else
            tt = [];
        end
    else
        tt = [];
    end
    switch oper,
    case 'rdivide', z2 = z2 / yC;
    otherwise,      z2 = z2 * yC;
    end
    z2 = cvx_reshape( z2, [ nA, nz ], tt );
    z2 = cvx( sz, z2 );
    if first, z = z2; first = false; else z = z + z2; end
end

if laff,
    % constant * everything
    nA = size( yA, 1 );
    t1 = reshape( 1 : prod( sy ), sy )';
    t2 = reshape( 1 : prod( sz ), [ sz(2), sz(1) ] )';
    z2 = yA; if raff, z2( 1, : ) = 0; end
    z2 = cvx_reshape( z2, [ nA * sy( 2 ), sy( 1 ) ], [], t1 );
    if issparse( z2 ),
        tt = any( z2, 2 );
        if cvx_use_sparse( size( z2 ), nnz( tt ) * sy( 1 ), isreal( z2 ) & isreal( xC ) ),
            z2 = z2( tt, : );
        else
            tt = [];
        end
    else
        tt = [];
    end
    switch oper,
    case 'ldivide', z2 = z2 / xC.';
    otherwise,      z2 = z2 * xC.';
    end
    z2 = cvx_reshape( z2, [ nA, nz ], tt, [], t2 );
    z2 = cvx( sz, z2 );
    if first, z = z2; first = false; else z = z + z2; end
end

if quad,
    % affine * affine
    tt = ax( : ) & ay( : );
    xA = xA( :, tt ); xB = xA( 1, : ); xA( 1, : ) = 0;
    yA = yA( :, tt ); yB = yA( 1, : ); yA( 1, : ) = 0;
    xM = size( xA, 1 ); yM = size( yA, 1 );
    if xM < yM, xA( yM, end ) = 0;
    elseif yM < xM, yA( xM, end ) = 0; end
    %
    % Quadratic form test 1: See if x == a conj( y ) + b for some real a, b,
    % so that the quadratic form involves a simple squaring (or sum of squares)
    %
    cyA   = conj( yA );
    alpha = sum( sum( real( xA .* yA ) ) ) ./ max( sum( sum( cyA .* yA ) ), realmin );
    if sum( sum( abs( xA - alpha * cyA ) ) ) <= 2 * eps * sum( sum( abs( xA ) ) ),
        beta = xB - alpha * conj( yB );
        yt = cvx( [ 1, size( yA, 2 ) ], yA ) + yB;
        if isreal( yA ) && isreal( yB ) && isreal( beta ),
            beta = ( 0.5 / alpha ) * beta;
            z2 = alpha * ( sum_square( yt + beta ) - sum_square( beta ) );
        elseif all( abs( beta ) <= 2 * eps * abs( xB ) ),
            z2 = alpha * sum_square_abs( yt );
        else
            error( 'Disciplined convex programming error:\n    Invalid quadratic form: product is not real.\n', 1 ); %#ok
        end
    else
        %
        % Quadratic form test 2: Extract the quadratic coefficient matrix
        % and test it for semidefiniteness
        %
        dx = find( any( xA, 2 ) | any( yA, 2 ) );
        zb = length( dx );
        cxA = conj( xA( dx, : ) );
        cyA = cyA( dx, : );
        P  = cxA * cyA.';
        Q  = cxA * yB.' + cyA * xB.';
        R  = xB * yB.';
        P  = 0.5 * ( P + P.' );
        if ~isreal( R ) || ~isreal( Q ) || ~isreal( P ),
            error( 'Disciplined convex programming error:\n   Invalid quadratic form: product is complex.', 1 ); %#ok
        else
            xx = cvx( zb, sparse( dx, 1 : zb, 1 ) );
            [ z2, success ] = quad_form( xx, P, Q, R );
            if ~success,
                error( 'Disciplined convex programming error:\n   Invalid quadratic form: neither convex nor concave.', 1 ); %#ok
            end
        end
    end
    if first, z = z2; first = false; else z = z + z2; end
end

if posy,
    [ ix, jx ] = find( reshape( px, sx ) );
    vx = log( cvx_subsref( x, px ) );
    [ iy, jy ] = find( reshape( py, sy ) );
    vy = log( cvx_subsref( y, py ) );
    [ iz, jz ] = find( sparse( 1 : nnz(px), jx, 1 ) * sparse( iy, 1 : nnz(py), 1 ) );
    z2 = exp( vec( cvx_subsref( vx, iz ) ) + vec( cvx_subsref( vy, jz ) ) );
    z2 = sparse( ix(iz), jy(jz), z2, sz(1), sz(2) );
    if first, z = z2; first = false; else z = z + z2; end
end

if vpos,
    [ ix, jx ] = find( reshape( gx, sx ) );
    vx = log( cvx_subsref( x, gx ) );
    [ iy, jy ] = find( reshape( gy, sy ) );
    vy = log( cvx_subsref( y, gy ) );
    [ iz, jz ] = find( sparse( 1 : nnz(gx), jx, 1 ) * sparse( iy, 1 : nnz(gy), 1 ) );
    z2 = exp( vec( cvx_subsref( vx, iz ) ) + vec( cvx_subsref( vy, jz ) ) );
    z2 = sparse( ix(iz), jy(jz), z2, sz(1), sz(2) );
    if first, z = z2; first = false; else z = z + z2; end %#ok
end

%
% Check that the sums are legal
%

v = cvx_vexity( z );
if any( isnan( v( : ) ) ),
    temp = 'Disciplined convex programming error:';
    tt = isnan( cvx_constant( z ) );
    if any( tt ),
        temp = [ temp, '\n    This expression produced one or more invalid numeric values (NaNs).' ];
    end
    if any( isnan( v( ~tt ) ) ),
        temp = [ temp, '\n   Illegal affine combination of convex and/or concave terms detected.' ];
    end
    error( temp, 1 );
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function VDSR(datasetName, SF, model, outRoute)

if isempty(model)
    error('no model');
else
    modelPath = ['VDSR model/',model];
    %gpu = 1;
    gpu = 0; % MAGI-ADAPT
end

model = load(modelPath);
net = model.net;
net = vl_simplenn_tidy(net);

managableMax = 300000;

dataDir = fullfile('data', datasetName);
f_lst = dir(fullfile(dataDir, '*.*'));
evaltable = zeros(numel(f_lst),2);
for f_iter = 1:numel(f_lst)
    f_info = f_lst(f_iter);
    if f_info.isdir, continue; end
    [~,imgName,~] = fileparts(f_lst(f_iter).name);
    imGT = imread(fullfile(dataDir, f_info.name));
    if size(imGT,3) > 1
        im = rgb2ycbcr(imGT);
    else
        im = imGT;
    end
    imhigh = modcrop(im, SF);
    imhigh = single(imhigh)/255;
    imlow = imresize(imhigh, 1/SF, 'bicubic');
    imlow = imresize(imlow, SF, 'bicubic');
    if size(imlow,3)>1
        imlowy = imlow(:,:,1);
        imlowy = max(16.0/255, min(235.0/255, imlowy));
        imlowcb = imlow(:,:,2);
        imlowcr = imlow(:,:,3);
    else
        imlowy = imlow;
    end

    if size(imlowy,1)*size(imlowy,2) > managableMax
        impred = runPatchVDSR(net, imlowy, gpu, 20);
    else
        if gpu, imlowy = gpuArray(imlowy); end;
        impred = runVDSR(net, imlowy, gpu);
    end
    if size(imGT,3) > 1
        impredColor = cat(3,impred,imlowcb,imlowcr);
        imwrite(ycbcr2rgb(uint8(impredColor*255)), fullfile(outRoute, [imgName, '.png']));
    else
        impredColor = impred;
        imwrite(uint8(impredColor*255), fullfile(outRoute, [imgName, '.png']));
    end
    imtest = imread(fullfile(outRoute, [imgName, '.png']));
    [psnr, ssim] = compute_diff(imGT,imtest,SF);
    evaltable(f_iter,1) = psnr; evaltable(f_iter,2) = ssim;
end

save(fullfile(outRoute,'psnr_ssim.mat'),'evaltable');

author: thomas-koehler

repo: SupER


text:
function A = ambiguityfunction(f,g)
%AMBIGUITYFUNCTION Ambiguity function
%   Usage: A = ambiguityfunction(f);
%          A = ambiguityfunction(f,g);
%
%   Input parameters:
%         f,g    : Input vector(s).
%
%   Output parameters:
%         A      : ambiguity function
%
%   `ambiguityfunction(f)` computes the (symmetric) ambiguity function of f.
%   The ambiguity function is computed as the two-dimensional Fourier transform
%   of the Wigner-Ville distribution |wignervilledist|.
%
%   **WARNING**: The quadratic time-frequency distributions are highly
%   redundant. For an input vector of length L, the quadratic time-frequency
%   distribution will be a $L \times L$ matrix.

% AUTHOR: Jordy van Velthoven
% TESTING: TEST_AMBIGUITYFUNCTION
% REFERENCE: REF_AMBIGUITYFUNCTION

upfname = upper(mfilename);
complainif_notenoughargs(nargin, 1, upfname);
complainif_toomanyargs(nargin, 2, upfname);

[f,~,W]=comp_sigreshape_pre(f,upfname);

if W>1
   error('%s: Only one-dimensional vectors can be processed.',upfname);
end

if (nargin == 1)
  if isreal(f)
    z1 = comp_fftanalytic(f);
  else
    z1 = f;
  end
  z2 = z1;

elseif (nargin == 2)
  [g,~,W]=comp_sigreshape_pre(g,upfname);
  if W>1
     error('%s: Only one-dimensional vectors can be processed.',upfname);
  end

  if ~all(size(f)==size(g))
  	error('%s: f and g must have the same length.', upfname);
  end;

  if xor(isreal(f), isreal(g))
      error('%s: One input is real, the other one must be real too. ',...
            upfname);
  end

  if isreal(f) || isreal(g)
    z1 = comp_fftanalytic(f);
    z2 = comp_fftanalytic(g);
  else
    z1 = f;
    z2 = g;
  end;
end

R = comp_instcorrmat(z1, z2);

A = fftshift(fft2(fft(R)));

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat