text:
function generateEvasionTest(lat_dev_m, speed_mps)

% Author: Alexander Wischnewski
%
% Description: creates a double lane change evasion scenario with the easiest trajectory planning
% possible to create challenging scenarios for the controller.
%
% Input:
%   lat_dev_m: Lateral deviation in m
%   speed_mps: Speed in mps

length = 3000;
idx_offset_start = 60;
idx_offset_end = 100;
n_blend = 10;
n_foresight = 4;
step_size = 5;
num_traj = length/step_size + 1;

% allocate memory
s.s_global_m = zeros(num_traj, 1);
s.debug_slow_tartraj_x_m = zeros(num_traj, 50);
s.debug_slow_tartraj_y_m = zeros(num_traj, 50);
s.debug_slow_tartraj_psi_rad = zeros(num_traj, 50) - pi/2;
s.debug_slow_tartraj_kappa_radpm = zeros(num_traj, 50);
s.debug_slow_tartraj_v_mps = speed_mps*ones(num_traj, 50);
s.debug_slow_tartraj_ax_mps2 = zeros(num_traj, 50);
s.debug_slow_tartraj_s_loc_m = zeros(num_traj, 50);
s.debug_slow_tartraj_LapCnt = ones(num_traj, 1);
s.debug_slow_tartraj_ax_lim_mps2 = 18.*ones(num_traj, 50);
s.debug_slow_tartraj_ay_lim_mps2 = 18.*ones(num_traj, 50);
s.debug_slow_tartraj_banking_rad = zeros(num_traj, 50);

% generate arrays with target trajectories in a step size of two meters
% the path is generated along the x axis
% after some time it is offset by lat_dev_m to the left side and back again
temp_traj = linspace(0, step_size*50, 50);
for idx = 1:1:num_traj
    s.s_global_m(idx) = (idx-1)*step_size;
    s.debug_slow_tartraj_x_m(idx, :) = temp_traj + (idx-1)*step_size;
    s.debug_slow_tartraj_s_loc_m(idx, :) = temp_traj;
    % if the current trajetory is one of the offset ones then adjust the y coordinate
    if idx >= idx_offset_start && idx < idx_offset_end
        % implement a sine wave over ten indices to move over to the new line
        steps_to_go = n_foresight+n_blend-(idx-idx_offset_start);
        current_step = n_foresight+n_blend-steps_to_go;
        if steps_to_go > 0
            if current_step < n_foresight
                s.debug_slow_tartraj_y_m(idx, 1:(n_foresight-current_step)) = 0;
                s.debug_slow_tartraj_psi_rad(idx, 1:(n_foresight-current_step)) = -pi/2;
                s.debug_slow_tartraj_kappa_radpm(idx, 1:(n_foresight-current_step)) = 0;
                s.debug_slow_tartraj_y_m(idx, (n_foresight-current_step+1):steps_to_go+1) = 0.5*lat_dev_m*(1-cos((0:n_blend)/n_blend*pi));
                s.debug_slow_tartraj_y_m(idx, (steps_to_go+2):end) = lat_dev_m;
                cos_dev_1 = 0.5*lat_dev_m*sin((0:n_blend)/n_blend*pi)/(n_blend*step_size/pi);
                cos_dev_2 = 0.5*lat_dev_m*cos((0:n_blend)/n_blend*pi)/(n_blend*step_size/pi)^2;
                s.debug_slow_tartraj_psi_rad(idx, (n_foresight-current_step+1):steps_to_go+1) = atan2(cos_dev_1, 1) - pi/2;
                s.debug_slow_tartraj_kappa_radpm(idx, (n_foresight-current_step+1):steps_to_go+1) = cos_dev_2./sqrt((1+cos_dev_1.^2).^3);
            else
                s.debug_slow_tartraj_y_m(idx, 1:steps_to_go+1) = 0.5*lat_dev_m*(1-cos(((current_step-n_foresight):n_blend)/n_blend*pi));
                s.debug_slow_tartraj_y_m(idx, (steps_to_go+2):end) = lat_dev_m;
                cos_dev_1 = 0.5*lat_dev_m*sin(((current_step-n_foresight):n_blend)/n_blend*pi)/(n_blend*step_size/pi);
                cos_dev_2 = 0.5*lat_dev_m*cos(((current_step-n_foresight):n_blend)/n_blend*pi)/(n_blend*step_size/pi)^2;
                s.debug_slow_tartraj_psi_rad(idx, 1:steps_to_go+1) = atan2(cos_dev_1, 1) - pi/2;
                s.debug_slow_tartraj_kappa_radpm(idx, 1:steps_to_go+1) = cos_dev_2./sqrt((1+cos_dev_1.^2).^3);
            end
        else
            s.debug_slow_tartraj_y_m(idx, :) = lat_dev_m;
        end
    elseif idx >= idx_offset_end
        % implement a sine wave over ten indices to move over to the new line
        steps_to_go = n_foresight+n_blend-(idx-idx_offset_end);
        current_step = n_foresight+n_blend-steps_to_go;
        if steps_to_go > 0
            if current_step < n_foresight
                s.debug_slow_tartraj_y_m(idx, 1:(n_foresight-current_step)) = lat_dev_m;
                s.debug_slow_tartraj_psi_rad(idx, 1:(n_foresight-current_step)) = -pi/2;
                s.debug_slow_tartraj_kappa_radpm(idx, 1:(n_foresight-current_step)) = 0;
                s.debug_slow_tartraj_y_m(idx, (n_foresight-current_step+1):steps_to_go+1) = lat_dev_m - 0.5*lat_dev_m*(1-cos((0:n_blend)/n_blend*pi));
                s.debug_slow_tartraj_y_m(idx, (steps_to_go+2):end) = 0;
                cos_dev_1 = -0.5*lat_dev_m*sin((0:n_blend)/n_blend*pi)/(n_blend*step_size/pi);
                cos_dev_2 = -0.5*lat_dev_m*cos((0:n_blend)/n_blend*pi)/(n_blend*step_size/pi)^2;
                s.debug_slow_tartraj_psi_rad(idx, (n_foresight-current_step+1):steps_to_go+1) = atan2(cos_dev_1, 1) - pi/2;
                s.debug_slow_tartraj_kappa_radpm(idx, (n_foresight-current_step+1):steps_to_go+1) = cos_dev_2./sqrt((1+cos_dev_1.^2).^3);
            else
                s.debug_slow_tartraj_y_m(idx, 1:steps_to_go+1) = lat_dev_m - 0.5*lat_dev_m*(1-cos(((current_step-n_foresight):n_blend)/n_blend*pi));
                cos_dev_1 = -0.5*lat_dev_m*sin(((current_step-n_foresight):n_blend)/n_blend*pi)/(n_blend*step_size/pi);
                cos_dev_2 = -0.5*lat_dev_m*cos(((current_step-n_foresight):n_blend)/n_blend*pi)/(n_blend*step_size/pi)^2;
                s.debug_slow_tartraj_psi_rad(idx, 1:steps_to_go+1) = atan2(cos_dev_1, 1) - pi/2;
                s.debug_slow_tartraj_kappa_radpm(idx, 1:steps_to_go+1) = cos_dev_2./sqrt((1+cos_dev_1.^2).^3);
            end
        else
            s.debug_slow_tartraj_y_m(idx, :) = 0;
        end
    end
end

%% Write to data dictionary
DDObj = ...
    Simulink.data.dictionary.open('TrajectoryPlanningEmulation.sldd');
dataSectObj = getSection(DDObj, 'Design Data');
try addEntry(dataSectObj, 'ltpl_log', s);
% if entry was there already modify it
catch e
    if isa(e, 'MSLException')
        ltplLogObj = getEntry(dataSectObj, 'ltpl_log');
        setValue(ltplLogObj, s);
    else
        print('Something went wrong.')
    end
end

% set switch to enable local log replay
switchObj = getEntry(dataSectObj, 'P_VDC_TrajEmulation_mode');
setValue(switchObj, 1);
% save & close
saveChanges(DDObj);
close(DDObj);

% set start position
Scenario_DD = Simulink.data.dictionary.open('raceline.sldd');
dataSectObj = getSection(Scenario_DD, 'Design Data');
x0_pose = getEntry(dataSectObj, 'x0_vehiclepose_stm');
setValue(x0_pose, [0, 0, -pi/2]);
saveChanges(Scenario_DD);
close(Scenario_DD);

author: TUMFTM


text:
%MDL_STANFORD_MDH Create model of Stanford arm using MDH conventions
%
%      mdl_stanford_mdh
%
% Script creates the workspace variable stanf which describes the
% kinematic and dynamic characteristics of the Stanford (Scheinman) arm
% using modified Denavit-Hartenberg parameters.
%
% Also defines the vectors:
%   qz   zero joint angle configuration.
%
% Notes::
% - SI units are used.
%
% References::
% - Kinematic data from "Modelling, Trajectory calculation and Servoing of
%   a computer controlled arm".  Stanford AIM-177.  Figure 2.3
% - Dynamic data from "Robot manipulators: mathematics, programming and control"
%   Paul 1981, Tables 6.5, 6.6
%
% See also SerialLink, mdl_puma560, mdl_puma560akb.

% MODEL: Stanford, Stanford arm, prismatic, 6DOF, modified_DH

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

clear L

qz = [0 0 0 0 0 0];

stanf = SerialLink([
            RevoluteMDH('d', 0.412)
            RevoluteMDH('d', 0.154, 'alpha', -pi/2)
            PrismaticMDH('alpha', pi/2, 'qlim', [0.2032 0.9144])
            RevoluteMDH()
            RevoluteMDH('alpha', -pi/2)
            RevoluteMDH('d', 0.263, 'alpha', pi/2)
        ], ...
    'name', 'Stanford arm MDH', ...
    'plotopt', {'workspace', [-2 2 -2 2 -2 2]} ...
    );

author: Allopart

repo: rbpf-gmapping


text:
% Program by Krish Mangal, R Harsha Vardhan and R Umesh Reddy, this program
% will be available for download in MATLAB central forever, on the Harsha's
% Account( http://goo.gl/cd7c ) along with other interesting programs we have
% developed or will develop in future.
%
% This program identifies all of the moving objects in a recorded video and
% shows those objects in color while displaying the other stationary
% objects in grayscale. This is a basic program which may be used with
% other programs which use applications of Motion Detection. Details of
% this technique are expained in the project report which comes along with
% this Program. We have used the function mmread also available on matlab
% central in our program, you must download this program also and add its
% location to MATLAB's path.
%
% see also DSBAM SSBAM Kirsch
%

clc
disp('Program by Krish Mangal, R Harsha Vardhan and R Umesh Reddy, IIT Patna');
clear all
close all
tic
h = waitbar(0,'Running......');
mov = (mmread('Input.wmv'));clc
disp('Program by Krish Mangal, R Harsha Vardhan and R Umesh Reddy, IIT Patna');
movi = mov.frames;
clear mov;
for i = 2:length(movi)
    movc = movi(i).cdata;
    waitbar(i/length(movi),h);
    movp = movi(i-1).cdata;
    u = motion(movc,movp);
    diff_img(i-1,:,:) = u;

end
close(h)

disp('Step 1 is completed Going for step 2 now!!!');
h = waitbar(0,'Running......');

movo = movi;
[a b c] = size(diff_img);
for k = 1:a
for i = 1:int32(b/5)
    for j = 1:int32(c/5)
        s = diff_img(k,(5*(i-1)+1):(5*(i-1)+5),(5*(j-1)+1):(5*(j-1)+5));
        s1(:,:) = s;
        q(i,j) = sum(sum(s1));
        if q(i,j)<500
            p(k,i,j) = false;
        else
            p(k,i,j) = true;
        end
        if(~p(k,i,j))
           temp = rgb2gray( movo(k).cdata((5*(i-1)+1):(5*(i-1)+5),(5*(j-1)+1):(5*(j-1)+5),1:3));
           for w = 1:3
           movo(k).cdata((5*(i-1)+1):(5*(i-1)+5),(5*(j-1)+1):(5*(j-1)+5),w) = temp;
           end
        end
    end
end
waitbar(k/a,h);
end
close(h)
clear movc movp diff_img temp q w  b u
disp('Step 2 is completed Going for final step now!!! This step takes a lot of time ');

for b = 1:a+1
movds(:,:,:,b) = movo(b).cdata;
end

[d e f ] = size(p);
t = 1;
img_path = movi(1).cdata;
for k = 1:d
    for i = 1:e
       for j = 1:f
        if(p(k,i,j))
            x_coe(t) = i;
            y_coe(t) = j;
            t = t+1;
        end
       end
    end
    xavg = int32(round((min(x_coe)+max(x_coe))/2));
    yavg = int32(round((min(y_coe)+ max(y_coe))/2));
    xcor = (2*xavg-1)*5/2;
    ycor = (2*yavg-1)*5/2;
   if k>1
    if xcor<xprev
        xsm = xcor;
        xlr = xprev;
    else

        xlr = xcor;
        xsm = xprev;
    end
    if ycor<yprev
        ysm = ycor;
        ylr = yprev;
    else

        ylr = ycor;
        ysm = yprev;
    end
    img_path(xsm:xlr,ysm:ylr,1) = 255;
    img_path(xsm:xlr,ysm:ylr,2) = 0;
    img_path(xsm:xlr,ysm:ylr,3) = 0;
    t=1;
    clear x_coe y_coe
   end
    xprev = xcor;
    yprev = ycor;

end
imshow(img_path)
disp('Execution Complete');
toc

author: Sable

repo: mcbench-benchmarks


text:
% ir_deblur_gcv1.m
% Edge-preserving regularized image deblurring (restoration)
% using (nonlinear!) generalized cross validation GCV (NGCV)
% to select the regularization parameter automatically
% 2013-10-01 Jeff Fessler, University of Michigan

if ~isvar('xtrue')
	ig = image_geom('nx', 100, 'ny', 128, 'dx', 2, 'down', 2);
	xell = [0 0 85 115 0 100;
		0 -60 30 20 0 -80; % mouth
		30 20 25 35 30 20; % eyes
		-30 20 25 35 -30 20;
		35 25 7 7 0 -100; % pupils
		-15 25 7 7 0 -100;
		0 75 60 15 0 -50; % hat
		]; % it is not north park, it is ...
	xtrue = ellipse_im(ig, xell, 'oversample', 3);
	xtrue = single(xtrue);
end

if ~isvar('y')
	psf = ones(3)/9;
	A = Gblur(ig.mask, 'psf', psf);
	ytrue = A * xtrue;

	%% add noise
	rng(0)
	snr = 20; % specify desired SNR of data in dB
	sig = 10^(-snr/20) * norm(ytrue(:)) / sqrt(numel(ytrue));
	y = ytrue + sig * randn(size(ytrue));

	snr_fun = @(x, xtrue) 20*log10(norm(xtrue(:)) / norm(x(:) - xtrue(:)));
	pr 'snr_fun(y, ytrue)'
	pr 'snr_fun(y, xtrue)'

	[nx ny] = size(y);
	nd = nx * ny; % # of data points
%	C = Cdiffs([], 'mask', ig.mask, 'offsets', [1 nx], 'type_diff', 'circshift');

	im plc 2 3
	clim = [-10 130]; % display all images on same gray scale
	im(xtrue, clim)
	im(ytrue, clim)
	im(y, clim)
prompt
end

delta = 5; % small compared to max(y) - min(y)
%dpot = @(t) t ./ sqrt(1 + (t / delta).^2); % derivative of potential
niter = 200; % max # of iter
tol = 1e-5; % stop tolerance for change in x

%% try several values of the regularization parameter

sig = 1; % pretend we don't know the noise level, using sigma=1

reglist = 2 .^ [-3:0.5:2];
nr = numel(reglist);
gcv = nan(nr,1);
rss = nan(nr,1);
snr = nan(nr,1); % solver only
gcv_best = inf;
x_best = [];
ir_best = [];

% emperically it seems better to generate w once outside the loop
% instead of making a new random bernoulli for every beta
w = 2 * (rand(nx,ny) > 0.5) - 1; % bernoulli +/- 1
weps = 0.05 * norm(y(:)) / norm(w(:));

for ir=1:nr
	reg = reglist(ir);
	step = sum(psf(:)) / (1 + reg * 8);

	R = Reg1(ig.mask, 'offsets', [1 nx], ...
		'beta', reg, 'pot_arg', {'hyper3', delta});

	% first run of PCG
	xinit = y;
	x = pwls_pcg1(xinit(ig.mask), A, 1, y(:), R, ...
		'niter', niter, 'stop_diff_tol', tol);
	x = ig.embed(x);

	snr(ir) = snr_fun(x, xtrue);

	% second run of PCG, now with a random perturbation
	xe = pwls_pcg1(xinit(ig.mask), A, 1, col(y + weps * w), R, ...
		'niter', niter, 'stop_diff_tol', tol);
	xe = ig.embed(xe);

	% todo: find reference for thie NGCV formula
	tr_ngcv = (1/nd) * ir_dot_double(w, A * (xe - x) / weps); % NGCV

	rss(ir) = norm(col(A*x) - y(:))^2 / sig^2; % RSS
	gcv(ir) = rss(ir) / (1 - tr_ngcv)^2;

	if gcv(ir) < gcv_best
		gcv_ir_best = ir;
		gcv_best = gcv(ir);
		x_best = x;
	end

	im(x, clim)
	xlabelf('$\log_2(\beta)$ = %4.1f', log2(reg))

%	im subplot 5
%	plot(log2(reglist), rss, '-o')
%	titlef 'RSS', xlabelf '$\log_2(\beta)$'

	im subplot 5
	plot(log2(reglist), snr, '-o', ...
		log2(reglist(gcv_ir_best)), snr(gcv_ir_best), '*')
	titlef 'SNR', xlabelf '$\log_2(\beta)$'

	im subplot 6
	plot(log2(reglist), gcv, '-o', ...
		log2(reglist(gcv_ir_best)), gcv(gcv_ir_best), '*')
	titlef 'NGCV', xlabelf '$\log_2(\beta)$'

	drawnow
end

	im(4, x_best, clim)
	xlabelf('$\log_2(\beta)$ = %4.1f', log2(reglist(gcv_ir_best)))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef TestTextUtilities < CoreTest
    % TestImageTextUtilities. Tests for the PTKTextUtilities class.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods
        function obj = TestTextUtilities
            % Test for StripFileparts method
            fs = filesep;
            filename_1 = ['Banana' fs 'Orange' fs 'MyName.ext'];
            filename_2 = ['Scott' fs 'Kubrick' fs 'NoName.bla'];
            both_filenames = {filename_1, filename_2};
            obj.Assert(strcmp('MyName', CoreTextUtilities.StripFileparts(filename_1)), 'Expected filename');
            obj.Assert(strcmp('NoName', CoreTextUtilities.StripFileparts(filename_2)), 'Expected filename');
            obj.Assert(strcmp({'MyName', 'NoName'}, CoreTextUtilities.StripFileparts(both_filenames)), 'Expected filename');
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%The aim of this simple benchmark is to illustrate the interest of restarting
%Nelder-Mead locally, from the last solution found, until no improvement is
%reached (to a given accuracy).
%Moreover, it shows that fminsearch has great difficulties at minimizing
%the most simple, smooth quadractic, objective function used. But
%restarting it locally corrects that.
%On the other hand, Nick Higham implementation of Nelder-Mead works fine,
%and the accuracy reached is also further improved by restarting it
%locally. Note that it may still happen that fminsearch performs better on
%other problems.

%Anyhow in theory, amongst direct search methods, one should not use even
%the restarted NM but rather MADS (C. Audet, J.E. Dennis, S. Le Digabel), which has
%guaranteed convergence even on non-smooth Clarke subdifferentiable objective functions.
%The restarted NM will also lead in practice to locally optimal solutions,
%although this is not theoretically guaranteed. It may fail in very
%particular of difficult situations. The reason for its good convergence
%properties in practice is that restarting it regenerates its search
%simplex and in the end many search directions are covered, which is a crude
%alternative to the POLL step of MADS (which is the step ensuring
%convergence). So the restarted NM will perform well even on non-smooth or
%discontinuous objective functions (not illustrated with this benchmark,
%other benchmarks are available on http://arxiv.org/abs/1104.5369).
%We put forward the restarted NM since it is easily available, and simple
%to use and will already work well enough in practice. But again, in
%theory, MADS should be used instead (to get a convergence certificate).

%For related works on optimization in systems and control, see the Matlab
%files http://www.mathworks.com/matlabcentral/fileexchange/33022 and
%http://www.mathworks.com/matlabcentral/fileexchange/33219 (and hyperlinked papers).

%Emile Simon, 18th october 2011

clear all
close all
clc

%Choice of the accuracies
precf = 10^-4;%'Objective' accuracy
precx = 10^-4;%'Simplex' accuracy
precs = 10^-4;%'Restarting' accuracy
Nmax = 20;%Number of variables tried% To be changed.
ntests = 10;%Number of tests for each N
disp = 0; %Turn to 1 to display the detail of the iterations of Nelder-Mead
if(disp==1) Disp = 'iter'; else Disp = 'off'; end

Objectif = @(x)sum(x.^2);%The objective function to be minimized, smooth quadratic. Perhaps the simplest objective function possible.
%func='sq';

options = optimset('TolF',precf,'TolX',precx,'MaxFunEvals',inf,'MaxIter',inf,'Display','off');%,'GradObj','on');%,'TolCon',precx

for s=2:Nmax
    for i =1:ntests

        x0=randn(s,1)./rand(s,1);%Difficult starting point

        tic
        [Sol1,Obj1(s,i)] = fminsearch(Objectif,x0,optimset('TolF',precf,'TolX',precx,'MaxFunEvals',inf,'MaxIter',inf,'Display',Disp));
        T1(s,i) = toc;

        acc = 1; Obj2(s,i) = Obj1(s,i); Sol2=Sol1;
        while (acc>precs)
            Objp(s,i) = Obj2(s,i);
            [Sol2,Obj2(s,i)] = fminsearch(Objectif,Sol2,optimset('TolF',precf,'TolX',precx,'MaxFunEvals',inf,'MaxIter',inf,'Display',Disp));
            acc = abs(abs(Objp(s,i)/Obj2(s,i))-1);
        end
        T2(s,i) = toc;

        %Nick Higham's implementation of Nelder-Mead
        tic
        [Sol3,Obj3(s,i),N3(s,i)]=nmsmax(@(x)-sum(x.^2),x0,[precx inf inf 0 disp]);
        T3(s,i) = toc;
        Obj3(s,i) = -Obj3(s,i);

        acc = 1; Obj4(s,i) = Obj3(s,i); Sol4=Sol3;
        while (acc>precs)
            Objp = Obj4(s,i);
            [Sol4,Obj4(s,i)] = nmsmax(@(x)-sum(x.^2),Sol4,[precx inf inf 0 disp]);
            T4(s,i)=toc;
            Obj4(s,i) = -Obj4(s,i);
            acc = Objp/Obj4(s,i)-1;
        end

    end
    s
end

figure
subplot(2,1,1)
semilogy(mean(Obj1'))
hold on
errorbar(mean(Obj1'),std(Obj1'))
errorbar(mean(Obj2'),std(Obj2'),'g')
errorbar(mean(Obj3'),std(Obj3'),'r')
errorbar(mean(Obj4'),std(Obj4'),'c')
legend('fminsearch','fminsearch','restarted fminsearch','nmsmax','restared nmsmax','Location','NorthWest')
title('Performance comparison of different versions of Nelder-Mead on \bf{min \Sigma_{i=1}^nx_i^2}')
ylabel('Avg. objective value')
subplot(2,1,2)
semilogy(mean(T1'))
hold on
errorbar(mean(T1'),std(T1'))
errorbar(mean(T2'),std(T2'),'g')
errorbar(mean(T3'),std(T3'),'r')
errorbar(mean(T4'),std(T4'),'c')
ylabel('Avg. computational time required (s)')
xlabel('Number of variables (n)')


text:
function [X, info] = IRhybrid_flsqr(A, b, varargin)
%IRhybrid_flsqr Hybrid version of FLSQR algorithm [...]
%
% options  = IRhybrid_flsqr('defaults')
% [X,info] = IRhybrid_flsqr(A,b)
% [X,info] = IRhybrid_flsqr(A,b,K)
% [X,info] = IRhybrid_flsqr(A,b,options)
% [X,info] = IRhybrid_flsqr(A,b,K,options)
%
% IRhybrid_flsqr is a hybrid iterative regularization method used for
% solving large-scale, ill-posed inverse problems of the form:
%               b = A*x + noise .
% The method combines FLSQR iteration (iterative regularization method)
% with a Tikhonov regularization method to stabilize the semiconvergence
% behavior that is characteristic of many iterative solvers applied to
% ill-posed problems.
%
% With 'defaults' as input returns the default options.  Otherwise outputs
% the iterates specified in K, using max(K) as MaxIter, and using all other
% default options.  With options as input: uses the user-specified options
% and all the other default options.
%
%
% Inputs:
%  A : either (a) a full or sparse matrix
%             (b) a matrix object that performs the matrix*vector operation
%             (c) user-defined function handle
%  b : right-hand side vector
%  K : (optional) integer vector that specifies which iterates are returned
%      in X; the maximum number of iterations is assumed to be max(K)
%      [ positive integer | vector of positive components ]
%  options : structure with the following fields (optional)
%      x0         - initial guess for the iterations; default = zero vector
%                   [ array | {'none'} ]
%      x_true     - true solution; allows us to returns error norms with
%                   respect to x_true at each iteration
%                   [ array | {'none'} ]
%      RegParam   - a value or a method to find the regularization
%                   parameter for the projected problems:
%                   [  non-neg. scalar | {'wgcv'} | 'gcv' | 'modgcv' |...
%                     'discrep' | 'discrepit' | 'optimal']
%                   This also determines which stopping rule is used
%                   If 'gcv', 'wgcv' or 'modgcv' is chosen, the iteration is
%                     stopped when the GCV function minimum stabilizes or
%                     increases within a certain window of iterations (see
%                     'stopGCV', 'FlatTol' and 'MinTol').
%                   If 'discrep' is chosen, and NoiseLevel is provided,
%                     then the discrepancy principle is used as stopping
%                     (see 'NoiseLevel' and 'eta').
%                   If 'discrepit' is chosen, and NoiseLevel is provided,
%                     then the discrepancy principle is used to set the
%                     parameter at each iteration, and the stabilization of
%                     successive parameters is used as stopping criterion
%                     (see 'NoiseLevel', 'eta' and 'regPflatTol').
%                   If 'optimal' is chosen, and x_true is provided,
%                     no stopping criterion is considered
%      stopGCV    - stopping criterion for the iterations when GCV is used
%                   [ 'GCVvalues' | {'resflat'} ]
%      resflatTol - tolerace for the stabilization of the residual
%                   (to be used if stopGCV is 'resflat')
%                   [ {0.05} | non-negative scalar ]
%      regPflatTol- tolerance for the stabilization of
%                   successive regularization parameters
%                   (to be used as stopping criterion)
%                   [ {0.9} | non-negative scalar ]
%      GCVflatTol - tolerance for detecting flatness (stabilization)
%                   in the GCV function as a stopping criterion
%                   [ {10^-6} | non-negative scalar ]
%      GCVminTol  - window of iterations - if the GCV minimum continues
%                   to increase over this window, then the iterations are
%                   stopped:
%                   [ {3} | positive integer ]
%      GCVweight  - weight to be used if RegParam is 'wgcv'
%                   [ 'adapt' | non-negative scalar ]
%      NoiseLevel - norm of noise in rhs divided by norm of rhs (must be
%                   assigned if RegParam is 'discrep')
%                   [ {'none'} | nonnegative scalar ]
%      eta        - safety factor for the discrepancy principle
%                   [ {1.01} | scalar greater than (and close to) 1 ]
%      RegParam0  - regularization parameter used in the first  projected
%                   problem (needed if RegParam is 'discrep')
%                   [ {1} | positive scalar ]
%      MaxIter    - maximum number of iterations
%                   [ {'none'} | positive integer ]
%      DecompOut  - returns the Golub-Kahan decomposition to the user
%                   [ 'on' | {'off'} ]
%      IterBar    - shows the progress of the iterations
%                   [ {'on'} | 'off' ]
%      NoStop     - specifies whether the iterations should proceed
%                   after a stopping criterion is satisfied
%                   [ 'on' | {'off'} ]
%   SparsityTrans - sparsity transform for the solution
%                   [ {'none'} | 'dwt' ]
%      wname      - discrete wavelet transform name (meaningful if
%                   SpartistyTrans is 'dwt')
%                   [ {'db1'} ]
%      wlevels    - discrete wavelet transform level (meaningful if
%                   SpartistyTrans is 'dwt')
%                   [ {2} | positive integer]
%   hybridvariant - kind of hybrid method to be implemented
%                   [ {'I'} | 'R' ]
%            tolX - tolerance for the weights (the modulus of the weights
%                   cannot be below this threshold
%                   [ {10^-10} | non-negative scalar ]
% Note: the options structure can be created using the function IRset.
%
% Outputs:
%   X : computed solutions, stored column-wise (at the iterations listed in K)
%   info: structure with the following fields:
%      its      - number of the last computed iteration
%      saved_iterations - iteration numbers of iterates stored in X
%      StopFlag - string that describes the output/stopping condition:
%                   * Flat GCV curve
%                   * Minimum of GCV function (within window of MinTol its)
%                   * Performed max number of iterations
%                   * Discrepancy principle satisfied
%                   * Breakdown of the Golub-Kahan bidiagonalization algorithm
%      StopReg  - structure with the following fields:
%                   * X: solution satisfying the stopping criterion
%                   * It: iteration satisfying the stopping criterion
%                   * RegP: regularization parameter at the iteration satisfying
%                     the stopping crierion
%                   * Xnrm: norm of the solution satisfying satisfying the
%                     stopping criterion
%                   * Rnrm: relative residual norm at the iteration
%                     satisfying the stopping criterion
%                   * Enrm: relative error norm at the iteration
%                     satisfying the stopping criterion (requires x_true)
%      Xnrm     - solution norms at each iteration
%      Rnrm     - relative residual norms at each iteration
%      Enrm     - relative error norms at each iteration (requires x_true)
%      RegP     - sequence of the regularization parameters
%      GCValues - GCV function evaluated at minimum point at each
%                 iteration (if RegParam is 'gcv', 'wgcv', 'modgcv')
%      V        - Golub-Kahan bidiagonalization basis vectors for the solution
%      U        - Golub-Kahan bidiagonalization basis vectors
%      B        - lower bidiagonal matrix computed by Golub-Kahan bidiagonalization
%
% See also: IRcgls, IRhybrid_fgmres, IRhybrid_gmres, IRhybrid_flsqr, IRget, IRset

% Julianne Chung, Virginia Tech
% Silvia Gazzola, University of Bath
% June, 2018.

% Initialization
defaultopt = struct('x0', 'none', 'MaxIter', 100 ,...
    'x_true', 'none', 'NoStop','off', 'IterBar', 'on',...
    'RegParam','wgcv',...
    'SparsityTrans', 'none', 'wname', 'db1', 'wlevels', 2,...
    'qnorm', 1, 'weight0', 'none',...
    'hybridvariant', 'I', 'tolX', 10^-10,...
    'GCVweight', 'adapt',...
    'GCVflatTol', 10^-6, 'GCVminTol', 3,...
    'stopGCV', 'GCVvalues', 'resflatTol', 0.05, 'regPflatTol', 0.9,...
    'NoiseLevel', 'none', 'eta', 1.01, 'RegParam0', 1, 'DecompOut', 'off');

if nargin == 0
    error('Not enough input arguments')
elseif nargin == 1
    % If input is 'defaults,' return the default options in X
    if nargout <= 1 && isequal(A,'defaults')
        X = defaultopt;
        return;
    else
        error('Not enough input arguments')
    end
end

defaultopt.restart = 'off';
defaultopt.verbosity = 'on';

% Check for acceptable number of optional input arguments
switch length(varargin)
    case 0
        K = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = [];
        else
            % no matter the order of appearance
            K = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = varargin{2};
        else
            % again, no matter the order of appearance
            K = varargin{2}; options = varargin{1};
        end
        if isfield(options, 'MaxIter') && ~isempty(options.MaxIter) && (~isempty(K) && options.MaxIter ~= max(K))
            warning('The value of MaxIter is discarded; the maximum value in K is taken as MaxIter')
        end
    otherwise
        error('Too many input parameters')
end

if isempty(options)
    options = defaultopt;
end

options = IRset(defaultopt, options);

MaxIter    = IRget(options, 'MaxIter',    [], 'fast');
RegParam   = IRget(options, 'RegParam',   [], 'fast');
x_true     = IRget(options, 'x_true',     [], 'fast');
NoStop     = IRget(options, 'NoStop',     [], 'fast');
IterBar    = IRget(options, 'IterBar',    [], 'fast');
omega      = IRget(options, 'GCVweight',  [], 'fast');
stopGCV    = IRget(options, 'stopGCV',    [], 'fast');
resdegflat = IRget(options, 'resflatTol', [], 'fast');
degflat    = IRget(options, 'GCVflatTol', [], 'fast');
regPflat   = IRget(options, 'regPflatTol', [], 'fast');
mintol     = IRget(options, 'GCVminTol',  [], 'fast');
NoiseLevel = IRget(options, 'NoiseLevel', [], 'fast');
eta        = IRget(options, 'eta',        [], 'fast');
RegParamk  = IRget(options, 'RegParam0',  [], 'fast');
restart    = IRget(options, 'restart',    [], 'fast');
verbose    = IRget(options, 'verbosity',  [], 'fast');
DecompOut  = IRget(options, 'DecompOut', [], 'fast');
tolX       = IRget(options, 'tolX',      [], 'fast');
q          = IRget(options, 'qnorm', [], 'fast');
SparsityTrans = IRget(options, 'SparsityTrans', [], 'fast');
hybridvariant = IRget(options, 'hybridvariant', [], 'fast');

verbose = strcmp(verbose, 'on');

adaptWGCV = strcmp(RegParam, {'wgcv'}) && strcmp(omega, {'adapt'});

% setting K
if isempty(K)
    K = MaxIter;
end
% sorting the iterations (in case they are shuffled in input)
K = K(:); K = sort(K,'ascend'); K = unique(K);
if ~((isreal(K) && (all(K > 0)) && all(K == floor(K))))
    error('K must be a vector of positive real integers')
end
if K(end) ~= MaxIter
    MaxIter = K(end);
end
% note that there is no control on K, as it does not go through IRset

if (strcmp(RegParam,'discrep') || strcmp(RegParam,'discrepit')) && ischar(NoiseLevel)
    error('The noise level (NoiseLevel) must be assigned')
end

if strcmp(RegParam,'optimal') && ischar(x_true)
    error('The true solution (x_true) must be assigned')
end

StopIt = MaxIter;

restart = strcmp(restart, 'on');
if restart
    ktotcount  = IRget(options, 'ktotcount', [],'fast');
    TotIterMax = IRget(options, 'TotIterMax',[],'fast');
    if strcmp(TotIterMax, 'none') || TotIterMax < MaxIter
        TotIterMax = MaxIter;
    end
    if strcmp(ktotcount, 'none')
        error('the total iteration counter must be assigned')
    end
    Ktot = IRget(options, 'Ktot', [], 'fast');
    % no checks on Ktot, it should be given from IRrestart
end

d = Atransp_times_vec(A, b);
n = length(d);
m = length(b);

if strcmp(SparsityTrans, 'none')
    Trans = speye(n);
elseif strcmp(SparsityTrans, 'dwt')
    wname   = IRget(options, 'wname',   [], 'fast');
    wlevels = IRget(options, 'wlevels', [], 'fast');
    if wlevels > 1/2*(log2(n))
        error('The assigned wavelet levels are too high. Make sure that wlevels <= 1/2*(log2(n))')
    end
    Trans = FreqMatrix('dwt', [sqrt(n) sqrt(n)], wname, wlevels);
end

% setting x0

x0 = IRget(options, 'x0', [], 'fast');

if strcmp(x0, 'none')
    x0 = zeros(n,1);
    r = b(:);
    precX = ones(n,1);
else
    try
        x0 = Trans*x0;
    catch
        error('Check the length of x0')
    end
    if max(abs(x0(:))) == 0
        r = b;
        x0 = zeros(n,1);
        precX = ones(n,1);
    else
        x0 = Trans'*x0;
        Ax0 = A_times_vec(A, x0);
        r = b(:) - Ax0;
        precX = x0;
        precX(abs(precX) < tolX) = eps;
        precX = (precX).^((2-q)/2);
    end
end

x = x0; % useful in case we have an immediate breakdown of the algorithm
beta = norm(r(:));
nrmb = norm(b(:));

% means no true solution
notrue = strcmp(x_true,'none');
% means we do not want to stop when the stopping criterion is satisfied
NoStop = strcmp(NoStop,'on');

% assessing if we want inner Tikhonov regularization
if strcmp(RegParam,'off')
    RegParam = 0;
end

Rfactor = strcmpi(hybridvariant, 'R');

% Declare matrices.
X                = zeros(n,length(K));
Xnrm             = zeros(max(K),1);
Rnrm             = zeros(max(K),1);
RegParamVect     = zeros(max(K),1);
M = zeros(max(K)+1,max(K));
T = zeros(max(K)+1);
Z = zeros(n, max(K));
V = zeros(n, max(K)+1);
U = zeros(m, max(K)+1);
rhs              = zeros(max(K)+1,1); % projected right-hand side
if restart
    saved_iterations = zeros(1, length(Ktot));
else
    saved_iterations = zeros(1, length(K));
end
GCV   = zeros(max(K),1);
Omega = zeros(max(K), 1);
warningGCV          = 0;
if notrue
    errornorms = false;
else
    errornorms = true;
    Enrm       = zeros(max(K),1);
    nrmtrue = norm(x_true(:));
    BestReg.RegP = [];
    BestReg.It = [];
    BestReg.X =[];
    BestReg.Enrm = [];
    BestReg.Xnrm = [];
    BestReg.Rnrm = [];
    BestEnrm = 1e10;
end
% Main Code Begins Here
u = r;
U(:,1) = u/beta; % no matter the preconditioning (because it's on the right)
rhs(1) = beta;
% Iterate
noIterBar = strcmp(IterBar,{'off'});
if ~noIterBar
  h_wait = waitbar(0, 'Running iterations, please wait ...');
end
j = 0;
for k=1:MaxIter
    if ~noIterBar
        waitbar(k/MaxIter, h_wait)
    end
    if restart, ktotcount = ktotcount + 1; end
    v = Atransp_times_vec(A, U(:,k)); v = v(:);
    v = Trans*v;
    for i = 1:k-1
        T(i,k)=V(:,i)'*v;
        v = v - T(i,k)*V(:,i);
    end
    T(k,k) = norm(v);
    v = v / T(k,k);
    %
    z = precX.*v;
    u = Trans'*z;
    u = A_times_vec(A, u); u = u(:);  %%%%
    for i = 1:k
        M(i,k) = U(:,i)'*u;
        u = u - M(i,k)*U(:,i);
    end
    M(k+1,k) = norm(u);
    u = u / M(k+1,k);
    U(:,k+1) = u;
    V(:,k) = v;
    Z(:,k) = z;

    rhsk = rhs(1:k+1); % current projected rhs

    if abs(T(k,k)) <= eps || abs(M(k+1,k)) <= eps
        if verbose
            disp('Flexible Golub-Kahan algorithm breaks down')
        end
        M = M(1:k+1,1:k);
        T = T(1:k,1:k);
        Z = Z(:,1:k);
        V = V(:,1:k);
        U = U(:,1:k+1);
        X(:,j+1) = x;
        X = X(:,1:j+1);
        if restart
            saved_iterations(j+1) = ktotcount-1;
        else
            saved_iterations(j+1) = k-1;
        end
        saved_iterations = saved_iterations(1:j+1);
        if k>1
            Xnrm    = Xnrm(1:k-1);
            Rnrm    = Rnrm(1:k-1);
            RegParamVect    = RegParamVect(1:k-1);
            if errornorms, Enrm = Enrm(1:k-1); end
        end
        % stop because the bidiagonal matrix is (numerically) singular
        % No chioce: even if NoStop is 'on'...we simpy cannot compute the solution, anymore
        if StopIt == MaxIter
            StopFlag = 'Breakdown of the Golub-Kahan algorithm';
            StopReg.X = x;
            StopReg.It = k-1;
            StopReg.RegP = RegParamk;
            StopReg.Xnrm = Xnrm(k-1);
            StopReg.Rnrm = Rnrm(k-1);
            if errornorms, StopReg.Enrm = Enrm(k-1); end
        end
        break
    end
    Mk = M(1:k+1,1:k);
    [Uk, Sk, Vk] = svd(Mk);
    if k==1
        Sk = Sk(1,1);
    else
        Sk = diag(Sk);
    end
    rhskhat = Uk'*rhsk;
    flsqr_res = abs(rhskhat(k+1))/nrmb;

    if Rfactor
        % update the Householder-QR factorization of Lk
        if k == 1
            [ZUk, ZRk] = householderQR(Z(:,1:k));
        else
            [ZUk, ZRk] = upd_householderQR(Z(:,1:k-1),...
            Z(:,k), ZUk, ZRk);
        end
        ZRksq = ZRk(1:k,1:k);
        [Uk, Vk, ~, Ck, Sk] = gsvd(Mk, ZRksq);
        rhskhat = Uk'*rhsk;
        if k==1
            gammak = Ck(1)/Sk(1);
        else
            gammak = sqrt(diag(Ck'*Ck)./diag(Sk'*Sk));
        end
    else
        ZRksq = eye(k);
    end

    % if tik
        if isscalar(RegParam)
            RegParamk = RegParam;
            RegParamVect(k) = RegParamk;
        elseif strcmp(RegParam,'discrep')
            if k==1
                RegParamVect(k) = RegParamk;
            end
        elseif strcmp(RegParam, 'discrepit')
            if flsqr_res > eta*NoiseLevel
                RegParamk = 0;
                RegParamVect(k) = RegParamk;
            else
                RegParamk = fzero(@(l)discrfcn(l, Mk, ZRksq, rhsk, nrmb, eta*NoiseLevel), [0, 1e10]);
                RegParamVect(k) = RegParamk;
            end
        elseif strcmp(RegParam, 'optimal')
            optfun = @(l)TikOptParam(l, Mk, ZRksq, rhsk, x0, Z(:,1:k), Trans, x_true);
            RegParamk = fmincon(optfun,0,[],[],[],[],0,.1);
            RegParamVect(k) = RegParamk;
        elseif strcmp(RegParam,'wgcv')
            if k>1
            if ~Rfactor
                if adaptWGCV
                    %Use the adaptive, weighted GCV method
                    Omega(k) = min(1, findomega(rhskhat, Sk));
                    omega = mean(Omega(1:k));
                end
                RegParamk = fminbnd('TikGCV', 0, Sk(1), [], rhskhat, Sk, omega);
                GCValk = GCVstopfun(RegParamk, Uk(1,:)', Sk, nrmb, m, n);
            else
                if adaptWGCV
                    %Use the adaptive, weighted GCV method
                    Omega(k) = min(1, findomega(rhskhat, gammak));
                    omega = mean(Omega(1:k));
                end
                RegParamk = fminbnd('TikGCV', 0, gammak(k), [], rhskhat, gammak, omega);
                GCValk = GCVstopfun(RegParamk, Uk(1,:)', gammak, nrmb, m, n);
            end
            RegParamVect(k) = RegParamk;
            GCV(k) = GCValk;
            else
            RegParamk = 0; GCValk = 0;
            RegParamVect(k) = RegParamk; GCV(k) = GCValk;
            end
        elseif strcmp(RegParam,'gcv')
            if ~Rfactor
                RegParamk = fminbnd('TikGCV', 0, Sk(1), [], rhskhat, Sk);
                GCValk = GCVstopfun(RegParamk, Uk(1,:)', Sk, nrmb, m, n);
            else
                RegParamk = fminbnd('TikGCV', 0, gammak(k), [], rhskhat, gammak);
                GCValk = GCVstopfun(RegParamk, Uk(1,:)', gammak, nrmb, m, n);
            end
            RegParamVect(k) = RegParamk;
            GCV(k) = GCValk;
        elseif strcmp(RegParam,'modgcv')
            if ~Rfactor
                RegParamk = fminbnd('TikGCV', 0, Sk(1), [], rhskhat, Sk, m);
                GCValk = GCVstopfun(RegParamk, Uk(1,:)', Sk, nrmb, m, n);
            else
                RegParamk = fminbnd('TikGCV', 0, gammak(k), [], rhskhat, gammak, m);
                GCValk = GCVstopfun(RegParamk, Uk(1,:)', gammak, nrmb, m, n);
            end
            RegParamVect(k) = RegParamk;
            GCV(k) = GCValk;
        else
            error('Invalid parameter choice method')
        end
        if ~Rfactor
            Dk = Sk.^2 + RegParamk^2;
            % Dk = Sk.^2 + RegParamk;
            rhskhat = Sk .* rhskhat(1:k);
            yhat = rhskhat(1:k)./Dk;
            y = Vk * yhat;
%             MZk = [Mk; RegParamk*eye(k)];
%             rhsZk = [rhsk; zeros(k,1)];
%             y = MZk\rhsZk;
        else
            MZk = [Mk; RegParamk*ZRksq(1:k,:)];
            % MZk = [Mk; sqrt(RegParamk)*ZRksq(1:k,:)];
            rhsZk = [rhsk; zeros(k,1)];
            y = MZk\rhsZk;
        end
        Rnrm(k) = norm(rhsk - Mk*y)/nrmb;
        d = Z(:,1:k)*y;
        x = x0 + d;
        precX = abs(x);
        precX(precX < tolX) = eps;
        precX = precX.^((2-q)/2);
        x = Trans'*x;
        % Compute norms
        Xnrm(k) = norm(x(:));
        if errornorms
            Enrm(k) = norm(x_true(:) - x(:))/nrmtrue;
            if Enrm(k)<BestEnrm
                BestReg.RegP = RegParamk;
                BestReg.It = k;
                BestReg.X = x;
                BestEnrm = Enrm(k);
                BestReg.Enrm = BestEnrm;
                BestReg.Xnrm = Xnrm(k);
                BestReg.Rnrm = Rnrm(k);
            end
        end
        AlreadySaved = 0;
        if any(k==K)
            j = j+1;
            X(:,j) = x;
            saved_iterations(j) = k;
            % this is used to save the last iteration, in the case
            % K = MaxIterIn
            % (when performing restarts, and the inner stopping criterion is not satisfied)
            if restart, saved_iterations(j) = ktotcount; end
            AlreadySaved = 1;
        end
        if restart
            if any(ktotcount == Ktot) && ~ AlreadySaved
                j = j+1;
                X(:,j) = x;
                saved_iterations(j) = ktotcount;
                AlreadySaved = 1;
            end
            if ktotcount == TotIterMax
                if ~ AlreadySaved
                    j = j+1;
                    saved_iterations(j) = ktotcount;
                    X(:,j) = x;
                end
                StopIt = k;
                StopReg.X = x;
                StopReg.It = k;
                StopReg.RegP = RegParamk;
                StopReg.Xnrm = Xnrm(k);
                StopReg.Rnrm = Rnrm(k);
                if errornorms
                    Enrm = Enrm(1:k);
                    StopReg.Enrm = Enrm(k);
                end
                Xnrm    = Xnrm(1:k);
                Rnrm    = Rnrm(1:k);
                RegParamVect    = RegParamVect(1:k);
                M = M(1:k+1,1:k);
                T = T(1:k,1:k);
                Z = Z(:,1:k);
                V = V(:,1:k);
                U = U(:,1:k+1);
                X = X(:,1:j);
                saved_iterations = saved_iterations(1:j);
                if verbose
                    disp('reached maximum number of iterations')
                end
                StopFlag = 'reached maximum number of iterations';
                break
            end
        end
        % update parameters, check stopping criteria
        if isscalar(RegParam)
        % Purely iterative method case.
        if strcmp(NoiseLevel, 'none')
            if k>1
            if abs((Rnrm(k)-Rnrm(k-1)))/Rnrm(k-1) < resdegflat && ...
                Rnrm(k) == min(Rnrm(1:k)) && StopIt == MaxIter
                if verbose
                    disp('The stopping criterion for flsqr is satisfied')
                end
                % Stop because the residual stabilizes.
                StopFlag = 'The residual norm stabilizes';
                if ~AlreadySaved && ~NoStop
                    j = j+1;
                    X(:,j) = x;
                    if restart
                        saved_iterations(j) = ktotcount;
                    else
                        saved_iterations(j) = k;
                    end
                    AlreadySaved = 1;
                end
                StopIt = k;
                StopReg.RegP = RegParamk;
                StopReg.It = k;
                StopReg.X = x;
                if errornorms, StopReg.Enrm = Enrm(k); end
                if ~ NoStop
                    Xnrm    = Xnrm(1:k);
                    Rnrm    = Rnrm(1:k);
                    RegParamVect = RegParamVect(1:k);
                    M = M(1:k+1,1:k);
                    T = T(1:k,1:k);
                    Z = Z(:,1:k);
                    V = V(:,1:k);
                    U = U(:,1:k+1);
                    if errornorms, Enrm = Enrm(1:k); end
                    X = X(:,1:j);
                    saved_iterations = saved_iterations(1:j);
                    break
                end
            end
            end
        else
            if Rnrm(k) < eta*NoiseLevel
            % Stopping criterion.
            if StopIt == MaxIter
                if verbose
                    disp('The discrepancy principle is satisfied')
                end
                StopFlag = 'The discrepancy principle satisfied';
                if ~AlreadySaved && ~NoStop
                    j = j+1;
                    X(:,j) = x;
                    if restart
                        saved_iterations(j) = ktotcount;
                    else
                        saved_iterations(j) = k;
                    end
                    AlreadySaved = 1;
                end
                StopIt = k;
                StopReg.RegP = RegParamk;
                StopReg.It = k;
                StopReg.X = x;
                if errornorms, StopReg.Enrm = Enrm(k); end
                if ~ NoStop
                    Xnrm    = Xnrm(1:k);
                    Rnrm    = Rnrm(1:k);
                    RegParamVect    = RegParamVect(1:k);
                    M = M(1:k+1,1:k);
                    T = T(1:k,1:k);
                    Z = Z(:,1:k);
                    V = V(:,1:k);
                    U = U(:,1:k+1);
                    if errornorms, Enrm = Enrm(1:k); end
                    X = X(:,1:j);
                    saved_iterations = saved_iterations(1:j);
                    % Stop because the discrepancy principle is satisfied.
                    break
                end
            end
            end
        end
        elseif strcmp(RegParam,'discrep')
            if Rnrm(k) < eta*NoiseLevel
                % stopping criterion
                if StopIt == MaxIter % the method has not stopped, yet
                    if verbose
                        disp('The discrepancy principle is satisfied')
                    end
                    StopFlag = 'discrepancy principle (secant update method) satisfied';
                    if ~AlreadySaved && ~NoStop
                        j = j+1;
                        X(:,j) = x;
                        if restart
                            saved_iterations(j) = ktotcount;
                        else
                            saved_iterations(j) = k;
                        end
                        AlreadySaved = 1;
                    end
                    StopIt = k;
                    StopReg.X = x;
                    StopReg.It = k;
                    StopReg.RegP = RegParamk;
                    StopReg.Xnrm = Xnrm(k);
                    StopReg.Rnrm = Rnrm(k);
                    if errornorms, StopReg.Enrm = Enrm(k); end
                    if ~ NoStop
                        Xnrm    = Xnrm(1:k);
                        Rnrm    = Rnrm(1:k);
                        RegParamVect    = RegParamVect(1:k);
                        M = M(1:k+1,1:k);
                        T = T(1:k,1:k);
                        Z = Z(:,1:k);
                        V = V(:,1:k);
                        U = U(:,1:k+1);
                        if errornorms, Enrm = Enrm(1:k); end
                        X = X(:,1:j);
                        saved_iterations = saved_iterations(1:j);
                        % stop because the discrepancy principle is satisfied
                        break
                    else
                        RegParamk = abs((eta*NoiseLevel - flsqr_res)/(Rnrm(k) - flsqr_res))*(RegParamk^2);
                        % RegParamk = RegParamk^2;
                        % RegParamk = abs((eta*NoiseLevel - flsqr_res)/(Rnrm(k) - flsqr_res))^2*RegParamk;
                        % RegParamk = abs((eta*NoiseLevel - flsqr_res)/(Rnrm(k) - flsqr_res))*RegParamk;
                        RegParamk = sqrt(RegParamk);
                        if k~=MaxIter, RegParamVect(k+1) = RegParamk; end
                    end
                else
                    RegParamk = abs((eta*NoiseLevel - flsqr_res)/(Rnrm(k) - flsqr_res))*(RegParamk^2);
                    % RegParamk = RegParamk^2;
                    % RegParamk = abs((eta*NoiseLevel - flsqr_res)/(Rnrm(k) - flsqr_res))*RegParamk;
                    RegParamk = sqrt(RegParamk);
                    if k~=MaxIter, RegParamVect(k+1) = RegParamk; end
                end
            else
                RegParamk = abs((eta*NoiseLevel - flsqr_res)/(Rnrm(k) - flsqr_res))*(RegParamk^2);
                % RegParamk = RegParamk^2;
                % RegParamk = abs((eta*NoiseLevel - flsqr_res)/(Rnrm(k) - flsqr_res))*RegParamk;
                RegParamk = sqrt(RegParamk);
                if k~=MaxIter, RegParamVect(k+1) = RegParamk; end
            end
        elseif strcmp(RegParam,'discrepit')
            if k>2
                % stopping criterion
                if StopIt == MaxIter % the method has not stopped, yet
                    if abs(RegParamVect(k)-RegParamVect(k-1))/RegParamVect(k-1) < regPflat && abs(RegParamVect(k-1)-RegParamVect(k-2))/RegParamVect(k-2)<regPflat
                        if verbose
                            disp('The stopping criterion for the discrepancy principle is satisfied')
                        end
                        StopFlag = 'discrepancy principle (stopping criterion) satisfied';
                        if ~AlreadySaved && ~NoStop
                            j = j+1;
                            X(:,j) = x;
                            if restart
                                saved_iterations(j) = ktotcount;
                            else
                                saved_iterations(j) = k;
                            end
                            AlreadySaved = 1;
                        end
                        StopIt = k;
                        StopReg.X = x;
                        StopReg.It = k;
                        StopReg.RegP = RegParamk;
                        StopReg.Xnrm = Xnrm(k);
                        StopReg.Rnrm = Rnrm(k);
                        if errornorms, StopReg.Enrm = Enrm(k); end
                        if ~ NoStop
                            Xnrm    = Xnrm(1:k);
                            Rnrm    = Rnrm(1:k);
                            RegParamVect    = RegParamVect(1:k);
                            M = M(1:k+1,1:k);
                            T = T(1:k,1:k);
                            Z = Z(:,1:k);
                            V = V(:,1:k);
                            U = U(:,1:k+1);
                            if errornorms, Enrm = Enrm(1:k); end
                            X = X(:,1:j);
                            saved_iterations = saved_iterations(1:j);
                            % stop because the discrepancy principle is satisfied
                            break
                        end
                    end
                end
            end
        elseif strcmp(RegParam,'wgcv') || strcmp(RegParam,'gcv') || strcmp(RegParam,'modgcv')
            % check the stopping criterion (all the possibilities)
            if k > 1
            if strcmpi(stopGCV, 'GCVvalues')
                if StopIt == MaxIter % the method has not stopped, yet
                if abs((GCV(k)-GCV(k-1)))/GCV(2) < degflat && StopIt == MaxIter
                % the method has not stopped, yet
                        if verbose
                            disp('The stopping criterion for GCV principle is satisfied')
                        end
                        % stop because the GCV curve is too flat
                        StopFlag = 'GCV curve too flat';
                        if ~AlreadySaved && ~ NoStop
                            j = j+1;
                            X(:,j) = x;
                            if restart
                                saved_iterations(j) = ktotcount;
                            else
                                saved_iterations(j) = k;
                            end
                            AlreadySaved = 1;
                        end
                        StopIt = k;
                        StopReg.X = x;
                        StopReg.It = k;
                        StopReg.RegP = RegParamk;
                        StopReg.Xnrm = Xnrm(k);
                        StopReg.Rnrm = Rnrm(k);
                        if errornorms, StopReg.Enrm = Enrm(k); end
                        if ~ NoStop
                            Xnrm    = Xnrm(1:k);
                            Rnrm    = Rnrm(1:k);
                            RegParamVect    = RegParamVect(1:k);
                            M = M(1:k+1,1:k);
                            T = T(1:k,1:k);
                            Z = Z(:,1:k);
                            V = V(:,1:k);
                            U = U(:,1:k+1);
                            if errornorms, Enrm = Enrm(1:k); end
                            X = X(:,1:j);
                            saved_iterations = saved_iterations(1:j);
                            % stop because the GCV stopping criterion is satisfied
                            break
                        end
                elseif GCV(k-1) < GCV(k) && ~ warningGCV && StopIt == MaxIter % Potential minimum reached.
                    warningGCV = 1;
                    % Save data just in case.
                    x_save = x;
                    k_save = k; % for computing the GCV stopping criterion
                    j_save = j;
                    AlreadySaved_save = AlreadySaved;
                    RegParamk_save = RegParamk;
                    if restart, ktotcount_save = ktotcount; end
                elseif warningGCV && k > min(k_save + mintol, MaxIter) && StopIt == MaxIter % Passed window
                    if GCV(k_save) < GCV(k_save+1:min(k_save + mintol, MaxIter))
                        if verbose
                            disp('The stopping criterion for GCV principle is satisfied')
                        end
                        StopFlag = 'increasing GCV minima';
                        StopIt = k_save;
                        StopReg.It = k_save;
                        StopReg.X = x_save;
                        StopReg.RegP = RegParamk_save;
                        StopReg.Xnrm = Xnrm(k_save);
                        StopReg.Rnrm = Rnrm(k_save);
                        if errornorms
                            StopReg.Enrm = Enrm(k_save);
                        end
                        if ~ NoStop
                            j = j_save;
                            saved_iterations = saved_iterations(1:j);
                            X = X(:,1:j);
                            if ~AlreadySaved_save
                                j = j+1;
                                X(:,j) = x_save;
                                if restart
                                    saved_iterations(j) = ktotcount_save;
                                else
                                    saved_iterations(j) = k_save;
                                end
                            end
                            Xnrm    = Xnrm(1:k_save);
                            Rnrm    = Rnrm(1:k_save);
                            RegParamVect    = RegParamVect(1:k_save);
                            M = M(1:k_save+1,1:k_save);
                            T = T(1:k_save,1:k_save);
                            Z = Z(:,1:k_save);
                            V = V(:,1:k_save);
                            U = U(:,1:k_save+1);
                            if errornorms
                                Enrm = Enrm(1:k_save);
                                if BestReg.It > k_save
                                    [BestReg.Enrm, BestReg.It] = min(Enrm);
                                    BestReg.RegP = RegParamVect(BestReg.It);
                                    BestReg.Xnrm = Xnrm(BestReg.It);
                                    BestReg.Rnrm = Rnrm(BestReg.It);
                                    % recompute the best solution again
                                    ktemp = BestReg.It;
                                    Mtemp = M(1:ktemp+1,1:ktemp);
                                    Ztemp = Z(:,1:ktemp);
                                    rhsk = rhs(1:ktemp+1);
                                    [Uk, Sk, Vk] = svd(Mtemp);
                                    if ktemp==1
                                        Sk = Sk(1,1);
                                    else
                                        Sk = diag(Sk);
                                    end
                                    rhskhat = Uk'*rhsk;
                                    if Rfactor
                                        [~, ZRksq] = qr(Ztemp,0);
                                    end
                                    if ~Rfactor
                                        Dk = Sk.^2 + RegParamk^2;
                                        rhskhat = Sk .* rhskhat(1:ktemp);
                                        yhat = rhskhat(1:ktemp)./Dk;
                                        y = Vk * yhat;
                                    else
                                        MLk = [Mtemp; RegParamk*ZRksq];
                                        rhsLk = [rhsk; zeros(ktemp,1)];
                                        y = MLk\rhsLk;
                                    end
                                    dtemp = Ztemp*y;
%                                     if precond, dtemp = P_solve(L, dtemp); end
                                    xtemp = x0 + dtemp;
                                    BestReg.X = xtemp;
                                end
                            end
                            X = X(:,1:j);
                            saved_iterations = saved_iterations(1:j);
                            if restart, ktotcount = ktotcount_save; end
                            % stop because the GCV stopping criterion is satisfied
                            break
                        end
                    else
                        warningGCV = 0;
                    end
                end
            elseif strcmpi(stopGCV, 'resflat')
                if abs((Rnrm(k)-Rnrm(k-1)))/Rnrm(k-1) < resdegflat && ...
                    Rnrm(k) == min(Rnrm(1:k)) && StopIt == MaxIter
                    if verbose
                        disp('The stopping criterion for GCV principle is satisfied')
                    end
                    % stop because discrepancy (i.e., residual for the
                    % regularized problem) stabilizes
                    StopFlag = 'the residual norm stabilizes';
                    if ~AlreadySaved && ~NoStop
                        j = j+1;
                        X(:,j) = x;
                        if restart
                            saved_iterations(j) = ktotcount;
                        else
                            saved_iterations(j) = k;
                        end
                        AlreadySaved = 1;
                    end
                    StopIt = k;
                    StopReg.X = x;
                    StopReg.It = k;
                    StopReg.RegP = RegParamk;
                    StopReg.Xnrm = Xnrm(k);
                    StopReg.Rnrm = Rnrm(k);
                    if errornorms, StopReg.Enrm = Enrm(k); end
                    if ~ NoStop
                        Xnrm    = Xnrm(1:k);
                        Rnrm    = Rnrm(1:k);
                        RegParamVect    = RegParamVect(1:k);
                        M = M(1:k+1,1:k);
                        T = T(1:k,1:k);
                        Z = Z(:,1:k);
                        V = V(:,1:k);
                        U = U(:,1:k+1);
                        if errornorms, Enrm = Enrm(1:k); end
                        X = X(:,1:j);
                        saved_iterations = saved_iterations(1:j);
                        break
                    end
                end
                end
            end
            end
        end
end
if k == MaxIter
    if StopIt == MaxIter
        % Stop because max number of iterations reached
        if verbose
            disp('Reached maximum number of iterations')
        end
        StopFlag = 'reached maximum number of iterations';
        if ~AlreadySaved
            j = j+1;
            X(:,j) = x;
            if restart
                saved_iterations(j) = ktotcount;
            else
                saved_iterations(j) = k;
            end
        end
        StopReg.X = x;
        StopReg.It = k;
        StopReg.RegP = RegParamk;
        StopReg.Xnrm = Xnrm(k);
        StopReg.Rnrm = Rnrm(k);
        if errornorms, StopReg.Enrm = Enrm(k); end
        Xnrm    = Xnrm(1:k);
        Rnrm    = Rnrm(1:k);
        RegParamVect    = RegParamVect(1:k);
        M = M(1:k+1,1:k);
        T = T(1:k,1:k);
        Z = Z(:,1:k);
        V = V(:,1:k);
        U = U(:,1:k+1);
        if errornorms, Enrm = Enrm(1:k); end
        X = X(:,1:j);
        saved_iterations = saved_iterations(1:j);
    end
end
if ~noIterBar, close(h_wait), end
if nargout==2
  if NoStop
      info.its = k;
  else
      info.its = StopIt;
  end
  info.saved_iterations = saved_iterations(1:j);
  info.StopReg = StopReg;
  info.StopFlag = StopFlag;
  if errornorms
    info.Enrm = Enrm;
    info.BestReg = BestReg;
  end
  info.Xnrm = Xnrm;
  info.Rnrm = Rnrm;
  info.RegP = RegParamVect;
  if strcmp(RegParam,'wgcv') || strcmp(RegParam,'gcv') || strcmp(RegParam,'modgcv')
    info.GCValues = GCV(1:k);
  end
  if strcmp(DecompOut,'on')
      info.V = V(:,1:k);
      info.U = U(:,1:k+1);
      info.Z = Z(:,1:k);
      info.T = T(1:k,1:k);
      info.M = M(1:k+1,1:k);
  end
  if restart
      info.ktotcount = ktotcount;
  end
end

%% ---------------SUBFUNCTIONS ---------------------------------------
function [U,R] = householderQR(L)
%
%  [U,R] = householderQR(L)
%  This function computes the Householder-QR factorization of L
%  (a "projected" regularization matrix), that will be used to define a
%  "projected" regularization matrix R to employ within the LSQR iterates.
%

[m, n] = size(L);
R = L;
U = zeros(m, n);
for k = 1:n
    x = L(k:m,k);
    e = zeros(length(x),1); e(1) = 1;
    u = sign(x(1))*norm(x(:))*e + x;
    u = u./norm(u(:));
    R(k:m, k:n) = R(k:m, k:n) -2*u*(u'*R(k:m, k:n));
    U(k:m,k) = u;
end

function [U,R] = upd_householderQR(L,ll,U,R)
%
% [U,R] = upd_householderQR(L, ll, U, R)
% This function updates the Householder-QR factorization of [L, ll].
%
% Input:
%   L  - matrix whose QR factorization is defined by U and R
%   ll - column appended to L, i.e., [L, ll]
%    U - matrix defining the orthogonal matrix Q, such that L = QR
%    R - upper triangular factor of L = QR
%

[m,n] = size(L);
Unew = zeros(m, n+1);
Unew(:,1:n) = U;
w = ll;
for i = 1:n
    u = U(i:m,i);
    w(i:m) = w(i:m) - 2*u*(u'*w(i:m));
end
v = w(1:n); x = w(n+1:m);
e = zeros(length(x),1); e(1) = 1;
u = sign(x(1))*norm(x(:))*e + x;
u = u./norm(u(:));
x = x -2*u*(u'*x);
Unew(n+1:m,n+1)=u;
U = Unew;
rr = [v; x];
R = [R, rr];

function omega = findomega(bhat, s)
%
%   omega = findomega(bhat, s)
%
%  This function computes a value for the omega parameter used in wGCV.
%
%  The method: Assume the 'optimal' regularization parameter to be the
%  smallest singular value.  Then we take the derivative of the GCV
%  function with respect to alpha, evaluate it at alpha_opt, set the
%  derivative equal to zero and then solve for omega.
%
%  Input:   bhat -  vector U'*b, where U = left singular vectors
%              s -  vector containing the singular values
%
%  Output:     omega - computed value for the omega parameter.

%
%   First assume the 'optimal' regularization parameter to be the smallest
%   singular value.
%

%
% Compute the needed elements for the function.
%
m = length(bhat);
n = length(s);

t0 = sum(abs(bhat(n+1:m)).^2);
alpha = s(end);
s2 = abs(s) .^ 2;
alpha2 = alpha^2;

tt = 1 ./ (s2 + alpha2);

t1 = sum(s2 .* tt);
t2 = abs(bhat(1:n).*alpha.*s) .^2;
t3 = sum(t2 .* abs((tt.^3)));

t4 = sum((s.*tt) .^2);
t5 = sum((abs(alpha2*bhat(1:n).*tt)).^2);

v1 = abs(bhat(1:n).*s).^2;
v2 = sum(v1.* abs((tt.^3)));

%
% Now compute omega.
%
omega = (m*alpha2*v2)/(t1*t3 + t4*(t5 + t0));

repo: IRtools


text:
classdef ZmapXYZcatalog < ZmapCatalog

end

author: CelsoReyes


text:
function [rf, info] = readrfvnmr(filename,verbosity)

% RF = READRFVNMR(FILENAME) reads a Varian VNMR format RF pulse
% description in FILENAME and generates an RF matrix with 3 columns
% specifying angle, magnitude, and duration.
%
% FILENAME can be in working directory, ~/vnmrsys/shapelib, or
% $vnmrsystem/shapelib and can be specified with or with trailing .RF

% Author: L. Martyn Klassen
% Copyright 2003 Robarts Research Institute
% This program is copyrighted worldwide by the Robarts Research
% Institute.  Any distribution, copying or redistribution is expressly
% forbidden.
% Created: April 17, 2003 14:22
% Modified: N/A

if nargin < 2
   verbosity = 0;
   if nargin < 1
      error('READRFVNMR requires RF file name.');
   end
end

% info defaults to an empty structure
info = struct;

% Get the vnmrsystem path
[status, vnmrpath] = unix('echo $vnmrsystem');
if status == 0
   % Strip end-of-line character
   vnmrpath = vnmrpath(1:end-1);
else
   % Make a guess - good for RRI
   vnmrpath = '/vnmr';
end

fid = findfile(filename, vnmrpath);
if fid < 0
   % Check if user forgot .RF or .rf ending
   orgfilename = filename;
   filename = [orgfilename '.RF'];
   fid = findfile(filename, vnmrpath);
   if fid < 0
      filename = [orgfilename '.rf'];
      fid = findfile(filename, vnmrpath);
      if fid < 0
         error('READRFVNMR unable to find %s', orgfilename);
      end
   end
end

if verbosity > 0
   % Report what file was actually opened for reading
   disp(sprintf('READRFVNMR reading file %s', fopen(fid)));
end

% Get the first line of the file
fileline = fgets(fid);
linenumber = 1;
rf = zeros(3,0);
while ischar(fileline)
   % Strip out any comment indicated by #
   idxcomment = strfind(fileline,'#');
   if ~isempty(idxcomment)
      % Look for information normally stored in header
      % It must be commented so we don't need to do this every line
      idx = strfind(fileline, 'VERSION');
      if ~isempty(idx)
         info.version = sscanf(fileline(idx+7:end), '%s', 1);
      end
      idx = strfind(fileline, 'TYPE');
      if ~isempty(idx)
         info.type = sscanf(fileline(idx+4:end), '%s', 1);
      end
      idx = strfind(fileline, 'MODULATION');
      if ~isempty(idx)
         info.modulation = sscanf(fileline(idx+10:end), '%s', 1);
      end
      idx = strfind(fileline, 'EXCITEWIDTH');
      if ~isempty(idx)
         info.excitewidth = sscanf(fileline(idx+11:end), '%f', 1);
      end
      idx = strfind(fileline, 'INVERTWIDTH');
      if ~isempty(idx)
         info.invertwidth = sscanf(fileline(idx+11:end), '%f', 1);
      end
      idx = strfind(fileline, 'INTEGRAL');
      if ~isempty(idx)
         info.integral = sscanf(fileline(idx+8:end), '%f', 1);
      end
      idx = strfind(fileline, 'RF_FRACTION');
      if ~isempty(idx)
         info.aref = sscanf(fileline(idx+11:end), '%f', 1);
      end
      idx = strfind(fileline, 'AREF');
      if ~isempty(idx)
         info.aref = sscanf(fileline(idx+4:end), '%f', 1);
      end
      idx = strfind(fileline, 'b1excite');
      if ~isempty(idx)
         info.b1excite = sscanf(fileline(idx+8:end), '%f', 1);
      end
      idx = strfind(fileline, 'b1invert');
      if ~isempty(idx)
         info.b1invert = sscanf(fileline(idx+8:end), '%f', 1);
      end

      % Strip out everything after the comment indicator
      fileline = fileline(1:idxcomment(1));
   end

   % Kill the last character: end-of-line or #
   fileline(end) = [];

   % If the line is empty skip it
   if ~isempty(fileline)
      % if it blank then skip
      if any(fileline ~= ' ')
         [A,count, errmsg, nextindex] = sscanf(fileline, '%f', inf);

         % If read failed, output the error
         if ~isempty(errmsg)
            fclose(fid)
            error('READRFVNMR failed with read error: %s', errmsg);
         end

         % Each line is to contain 3 and only 3 numbers
         if (count > 0) & (count ~= 3)
            fclose(fid);
            error('READRFVNMR malformed line number %d', linenumber);
         end

         % Store the read values into rf array
         rf(:,linenumber) = A;
         linenumber = linenumber + 1;
      end
   end
   % Get the next line
   fileline = fgets(fid);
end
fclose(fid);

% Add a default aref value if not found
if ~isfield(info, 'aref')
   info.aref = 0.5;
end

% Switch from row to column orientation
rf = rf';
return

function fid = findfile(filename, vnmrpath)
% Search local, user vnmrsys, and global vnmrsys for RF file
fid = fopen(filename, 'r');
if fid < 0
   orgfilename = filename;
   filename = ['~/vnmrsys/shapelib/' orgfilename];
   fid = fopen(filename, 'r');
   if fid < 0
      filename = [vnmrpath '/shapelib/' orgfilename];
      fid = fopen(filename, 'r');
   end
end
return

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [tg,theta] = detTG(im,radius,norient)
% function [tg,theta] = detTG(im,radius,norient)
%
% Compute smoothed but not thinned TG fields.

if nargin<2, radius=0.02; end
if nargin<3, norient=8; end

[h,w,unused] = size(im);
idiag = norm([h w]);
if isrgb(im), im=rgb2gray(im); end

% compute texture gradient
no = 6;
ss = 1;
ns = 2;
sc = sqrt(2);
el = 2;
k = 64;
fname = sprintf( ...
    'unitex_%.2g_%.2g_%.2g_%.2g_%.2g_%d.mat',no,ss,ns,sc,el,k);
textonData = load(fname); % defines fb,tex,tsim
tmap = assignTextons(fbRun(textonData.fb,im),textonData.tex);
[tg,theta] = tgmo(tmap,k,idiag*radius,norient,...
                  'smooth','savgol','sigma',idiag*radius);

author: Cloud-CV

repo: object-proposals


text:
function codeKeeper = syncPatternDetectNEV(NEVStruct)

% syncPatternDetectNEV
%
% Finds the code for the SYNC pattern on the SYNC output of Blackrock
% Microsystems Neural Signal Processor (NSP) recorded by analog input
% channel 16 and thresholded as extracted spikes and saved in the NEV file.
%
%   INPUT
%
%   NEVStruct:  The NEV structure containing the SYNC pulse.
%
%   OUTPUT
%
%   codeKeeper: A cell structure containing the unique code and their
%               respective timestamps.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   USAGE EXAMPLE:
%
%   codeKeeper = syncPatternDetectNEV(NEVStruct)
%
%   In the example above, the NEV structure, containing the continueus
%   signal from the SYNC pulse is passed to the function. The output will
%   contain the unique codes parsed from the signal file and their
%   respective timestamps.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Kian Torab
%   kian@blackrockmicro.com
%   Blackrock Microsystems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0: July 7, 2014
%   - Initial release.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Setup
separationVar = 'D';

%% Detect the rising edges
edgeTS = NEVStruct.Data.Spikes.TimeStamp(NEVStruct.Data.Spikes.Electrode == 144);
differenceTS = double(diff(edgeTS));
pulseDifferenceLenght = mode(differenceTS);

convertedChar = [];
for idx = 1:length(differenceTS)
    if differenceTS(idx) < pulseDifferenceLenght*1.1;
        convertedChar = [convertedChar, '1'];
    elseif differenceTS(idx) < 9*pulseDifferenceLenght;
        convertedChar = [convertedChar, repmat('0', 1, round((differenceTS(idx) - pulseDifferenceLenght)/pulseDifferenceLenght)), '1'];
    else
        convertedChar = [convertedChar, separationVar];
    end
end

begTS = edgeTS(find(differenceTS > pulseDifferenceLenght * 12) + 1);
%%
begTS = [edgeTS(1) begTS];

codeKeeper{1} = bin2dec(regexp(convertedChar, separationVar, 'split'))';
codeKeeper{2} = begTS;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% polyvalm2 - Evaluate polynomial with matrix argument.
%*************************************************************************************
%
%  MATLAB (R) is a trademark of The Mathworks (R) Corporation
%
%  Function:    polyvalm2
%  Filename:    polyvalm2.m
%  Programmer:  James Tursa
%  Version:     1.1
%  Date:        November 11, 2009
%  Copyright:   (c) 2009 by James Tursa, All Rights Reserved
%
%  This code uses the BSD License:
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.
%
%  polyvalm2 evaluates a polynomial with a square matrix argument faster
%  than the MATLAB built-in functions polyvalm or mpower.
%
%   Y = polyvalm2(P,X), when P is a vector of length N+1 whose
%   elements are the coefficients of a polynomial, is the value
%   of the polynomial evaluated with matrix argument X.  X must
%   be a square matrix.
%
%       Y = P(1)*X^N + P(2)*X^(N-1) + ... + P(N)*X + P(N+1)*I
%
%   Class support for inputs P, X:
%      float: double, single
%
%  The polyvalm2 speed improvements come from the following:
%
%  1) The MATLAB built-in function polyvalm uses Horner's method.
%     polyvalm2 uses a binary decomposition of the matrix powers
%     to do the calculation more efficiently, reducing the total
%     number of matrix multiplies used to calculate the answer.
%
%  2) polyvalm calculates the product of a scalar times a matrix
%     as the product of diag(scalar*ones(etc))*matrix ... i.e. it
%     does a matrix multiply. polyvalm2 will calculate this more
%     efficiently as the simple product scalar*matrix.
%
%  3) polyvalm does all of the calculations shown above, even if
%     the coefficient P(i) is zero. polyvalm2 does not do
%     calculations for P(i) coefficients that are zero.
%
%  4) polyvalm converts sparse matrix inputs into full matrices to
%     do the calculations, whereas polyvalm2 keeps the intermediate
%     calculations and the answer sparse.
%
%  An extreme case of speed difference can be found with a sparse
%  matrix example:
%
%  >> A = sprand(2500,2500,.01);
%  >> p = [1 2 3 4];
%  >> tic;polyvalm(p,A);toc
%     Elapsed time is 43.669362 seconds.
%  >> tic;polyvalm2(p,A);toc
%     Elapsed time is 4.240375 seconds.
%
%  The trade-off is that polyvalm2 uses more memory for intermediate
%  variables than polyvalm, so for very large matrices polyvalm2 can
%  run out of memory. In these cases polyvalm2 will abandon the
%  efficient calculation method and just call the built-in polyvalm.
%  For sparse matrix inputs, however, polyvalm2 will typically be more
%  memory efficient than the MATLAB polyvalm function.
%
%  Caution: Since polyvalm2 uses different calculations to form the matrix
%  powers, the end result may not match polyvalm exactly. Also, for the
%  case where only the leading coefficient is non-zero, polyvalm2 may not
%  match mpower exactly. But the answer will be just as accurate. And if
%  there are inf's or NaN's involved, then the end result will not, in
%  general, match polyvalm or mpower results. This should not be a great
%  drawback to using polyvalm2, however, since even the MATLAB built-in
%  functions polyvalm and mpower will not match each other in these cases.
%  (By reordering the calculations, the NaN's propagate differently)
%
%  Change Log:
%  Nov 11, 2009: Updated for sparse matrix input --> sparse result
%
%**************************************************************************

function Y = polyvalm2(p,X)
%\
% Check the arguments
%/
if( nargin ~= 2 )
    error('MATLAB:polyvalm2:InvalidNumberOfArgs','Need two input arguments.');
end
if( nargout > 1 )
    error('MATLAB:polyvalm2:TooManyOutputArgs','Too many output arguments.');
end
classname = superiorfloat(p,X);
if( ~(isvector(p) || isempty(p)) )
    error('MATLAB:polyvalm2:InvalidP','P must be a vector.');
end
z = size(X);
if( length(z) > 2 || z(1) ~= z(2) )
    error('MATLAB:polyvalm2:NonSquareMatrix','Matrix must be square.');
end
if( isempty(X) )
    if( issparse(X) )
        Y = sparse(z(1),z(2));
    else
        Y = zeros(z,classname);
    end
    return
end
%\
% Clear out any leading zeros and reverse the coefficients
%/
try
    f = find(p,1,'first');
    if( isempty(f) )
        if( issparse(X) )
            Y = sparse(z(1),z(2));
        else
            Y = zeros(z,classname);
        end
        return
    end
    p2 = p(end:-1:f);
%\
% Initialize return value with the constant term, and then set the
% constant term coefficient to 0 so we don't process it anymore.
%/
    if( issparse(X) )
        Y = diag(p2(1) * sparse(ones(z(1),1)));
    else
        Y = diag(p2(1) * ones(z(1),1,classname));
    end
    p2(1) = 0;
%\
% Special small cases use custom code for quick return
%/
    np = length(p2);
    if( np == 1 )
        return
    elseif( np == 2 )
        if( p2(2) == 1 )
            Y = Y + X;
        else
            Y = Y + p2(2) * X;
        end
        return
    end
%\
% Initialize the cell array that will hold the X powers
%/
    cp{np} = [];
%\
% Get the binary decomposition of the powers as rows of binary characters,
% each row representing a different power, and arranged from 0 at the top
% to the highest power at the bottom.
%/
    bp = dec2bin(0:(np-1));
%\
% Loop through the bit positions from least significant to most significant
%/
    P = X;
    zz = size(bp,2);
    for n=zz:-1:1
%\
% Only process those rows with non-zero coefficients. If the bit position
% for this row is 1, then we need to apply the current power of X to the
% cell for this row. Each cell is building up the appropriate power of X.
% cp{1} is for X^0 (not used or needed actually because we have that piece
% already calculated from above), cp{2} is for X^1, cp{3} is for X^2, etc.
% P is the current power of X, i.e. P = X^(2^(zz-n))
%/
        check = (p2 ~= 0);
        for m=1:np
            if( check(m) && bp(m,n) == '1' )
                if( isempty(cp{m}) )
                    cp{m} = P;
                else
                    cp{m} = cp{m} * P;
                end
%\
% Look at all the downstream bit patterns. If any match the current one for
% the bits processed so far, then just copy the current cell into the
% downstream cells (no need to repeat the calculation downstream because we
% already know the answer). Then reset the check flag for that downstream
% row so we don't process it for this particular loop index.
%/
                for k=m+1:np
                    if( check(k) && isequal(bp(m,n:zz),bp(k,n:zz)) )
                        cp{k} = cp{m};
                        check(k) = 0;
                    end
                end
%\
% If the remaining leftmost bits of the current row are 0, then we are done
% with this power of X. Apply the coefficient and add it to the result,
% then free up the memory for this cell position. Also, reset the
% coefficient for this row to 0 so we know not to process this row anymore.
%/
                if( all(bp(m,1:(n-1))=='0') )
                    Y = Y + p2(m) * cp{m};
                    p2(m) = 0;
                    cp{m} = [];
                end
            end
        end
%\
% Square the current power of X, but not if it is the last index in the
% loop because in that case it won't be used or needed. For some reason,
% P * P is a lot faster than P^2.
%/
        if( n ~= 1 )
            P = P * P;
        end
    end
%\
% The binary decomposition scheme used too much memory, so clear all of the
% local large variables and just call the built in polyvalm function
% instead. This is slower and computationally less efficient, but it is
% more memory efficient so it might work.
%/
catch
    warning('MATLAB:polyvalm2:OutOfMemory','Out Of Memory ... Resorting to built-in polyvalm');
    clear cp P bp p2 check;
    Y = polyvalm(p,X);
end
return
end


text:
function im_warp = warpImageFast(im,XXdense, YYdense)

%{
Citation:
J. Xiao, K. A. Ehinger, A. Oliva and A. Torralba.
Recognizing Scene Viewpoint using Panoramic Place Representation.
Proceedings of 25th IEEE Conference on Computer Vision and Pattern Recognition, 2012.
http://sun360.mit.edu
%}

minX = max(1,floor(min(min(XXdense)))-1);
minY = max(1,floor(min(min(YYdense)))-1);

maxX = min(size(im,2),ceil(max(max(XXdense)))+1);
maxY = min(size(im,1),ceil(max(max(YYdense)))+1);

im = im(minY:maxY,minX:maxX,:);

for c=1:size(im,3)
    % im_warp(:,:,c) = uint8(interp2(double(im(:,:,c)), XXdense-minX+1, YYdense-minY+1,'*cubic'));
    %im_warp(:,:,c) = interp2(im(:,:,c), XXdense-minX+1, YYdense-minY+1,'*cubic');
    im_warp(:,:,c) = interp2(im(:,:,c), XXdense-minX+1, YYdense-minY+1,'*linear');
end

author: jianxiongxiao

repo: ProfXkit


text:
function [px, py, gx, gy] = position_3_link(z,P)
%[px,py,gx,gy] = POSITION_3_LINK(Z,P)
%
%FUNCTION:  This function computes the cartesian positions
%    of the CoM and tip of each link
%INPUTS:
%
%
%OUTPUTS:
%    px = [nLink X nTime] position of the end of each link
%    py = [nLink X nTime] position of the end of each link
%    gx = [nLink X nTime] position of the CoM of each link
%    gy = [nLink X nTime] position of the CoM of each link
%
%NOTES:
%    This file was automatically generated by writePosition.m

g  = P.g ; %gravity
m1 = P.m(1); % Link 1 mass
m2 = P.m(2); % Link 2 mass
m3 = P.m(3); % Link 3 mass
l1 = P.l(1); % Link 1 length
l2 = P.l(2); % Link 2 length
l3 = P.l(3); % Link 3 length
I1 = P.I(1); % Link 1 moment of inertia about its center of mass
I2 = P.I(2); % Link 2 moment of inertia about its center of mass
I3 = P.I(3); % Link 3 moment of inertia about its center of mass
d1 = P.d(1); % Link 1 distance between center of mass and parent joint
d2 = P.d(2); % Link 2 distance between center of mass and parent joint
d3 = P.d(3); % Link 3 distance between center of mass and parent joint

th1 = z(1,:);
th2 = z(2,:);
th3 = z(3,:);

nTime = length(th1);
px = zeros(3,nTime);
py = zeros(3,nTime);
gx = zeros(3,nTime);
gy = zeros(3,nTime);

px(1,:) = l1.*cos(th1);
py(1,:) = l1.*sin(th1);
gx(1,:) = d1.*cos(th1);
gy(1,:) = d1.*sin(th1);

px(2,:) = l1.*cos(th1) + l2.*cos(th2);
py(2,:) = l1.*sin(th1) + l2.*sin(th2);
gx(2,:) = d2.*cos(th2) + l1.*cos(th1);
gy(2,:) = d2.*sin(th2) + l1.*sin(th1);

px(3,:) = l1.*cos(th1) + l2.*cos(th2) + l3.*cos(th3);
py(3,:) = l1.*sin(th1) + l2.*sin(th2) + l3.*sin(th3);
gx(3,:) = d3.*cos(th3) + l1.*cos(th1) + l2.*cos(th2);
gy(3,:) = d3.*sin(th3) + l1.*sin(th1) + l2.*sin(th2);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
i = 0;

i = i+1;
db(i).mouse_name    = 'M150329_MP009';
db(i).date          = '2015-04-29';
db(i).expts         = [4 5 6];
db(i).expred        = 3;
db(i).nchannels_red = 2;
db(i).nchannels     = 1;
db(i).gchannel      = 1;
db(i).nplanes       = 1;
db(i).planesToProcess = [1];
%
i = i+1;
db(i).mouse_name    = 'M150329_MP009';
db(i).date          = '2015-04-10';
db(i).expred        = 3;
db(i).expts         = [5 6 7 8 9 10 11];
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    = 'M150326_MP008';
db(i).date          = '2015-04-09';
db(i).expred        = 2;
db(i).expts         = [3 4 6 8];
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    = 'M150423_MP014';
db(i).date          = '2015-06-16';
db(i).expts         = [8];
db(i).expred        = 4;
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    =  'M150303_MP007';
db(i).date          = '2015-04-05';
db(i).expts         = [1 7 9 10 11 12 13 14];
db(i).expred        = 6;
db(i).nchannels_red = 2;
db(i).comments      = 'singleORIsat different pause lengths';
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    =  'M150303_MP007';
db(i).date          = '2015-04-04';
db(i).expts         = [3 5 6 10];
db(i).expred        = 4;
db(i).comments      = 'multiGRAT slow, fast and textures';
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    =  'M150326_MP008';
db(i).date          = '2015-04-11';
db(i).expts         = [4 6 7 9 10 11 12 14];
db(i).expred        = 5;
db(i).comments      = 'all sorts of crazy stims with SAT';
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    = 'M150326_MP008';
db(i).date          = '2015-04-16';
db(i).expts         = [1 3 5];
db(i).expred        = 2;
db(i).comments      = 'stimGRATfast not biased / biased';
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    = 'M150329_MP009';
db(i).date          = '2015-04-09';
db(i).expts         = [1 3];
db(i).expred        = 2;
db(i).comments      = 'multiGRAT, have not validated the channels';
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    = 'M150329_MP009';
db(i).date          = '2015-04-25';
db(i).expts         = [1 3 4];
db(i).expred        = 2;
db(i).comments      = 'fastgratings2';
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    = 'M150329_MP009';
db(i).date          = '2015-04-27';
db(i).expts         = [3 5 6];
db(i).expred        = 4;
db(i).comments      = 'fast gratings unbiased / biased';
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

i = i+1;
db(i).mouse_name    = 'M150422_MP012';
db(i).date          = '2015-05-20';
db(i).expts         = [4];
db(i).expred        = 4;
db(i).comments      = 'test registration types';
db(i).nchannels_red = 2;
db(i).nchannels = 1;
db(i).gchannel = 1;
db(i).nplanes = 1;

author: cortex-lab


text:
function	profile_end(profile_on, profile_dir)

if ~exist('profile_on','var'), profile_on = 0; end;

switch profile_on
case	0,
    ptime(toc)
case	1,
    ptime(toc)

    if str2num(strtok(version,'.')) > 6
%        profile off
		if ~exist('profile_dir','var'), profile_dir = 'profile'; end;
        profile viewer
        p = profile('info');
        profsave( p, profile_dir)
    else
		if ~exist('profile_dir','var'), profile_dir = './profile'; end;
	    pwd_dir = pwd;
		if ~exist(profile_dir,'dir'),
			flag=mkdir(profile_dir);
		end
	    cd(profile_dir)
        profile report profile_job
	    cd(pwd_dir)
    end

case	2,
    ptime(toc)
    profile viewer
end

author: KamitaniLab

repo: GenericObjectDecoding


text:
function [X]=tt_meshgrid_vert(varargin)
% Analogue of the meshgrid function for the TT-format,
% !!!!!!!!!!!! concatenating tensor trains with ones with tkron !!!!!!!!!!!
% !!!!!!!!!!!! contrarily to tt_meshgrid (which uses kron2)     !!!!!!!!!!!
%   X = TT_MESHGRID_VERT(A,B,C,...) Computes the meshgrid based on "1d"
%       representations
%   X = TT_MESHGRID_VERT(A) Computes the meshgrid based on the cell array A of
%       the representations
%   X = TT_MESHGRID_VERT(T,D) Computes the d-dimensional meshgrid, using T as a
%       one-dimensional grid
%   X = TT_MESHGRID_VERT(..., 'id', M) Expands with M instead of all-ones

vars = varargin;
M = [];
for i=1:numel(vars)
    if (isa(vars{i}, 'char'))&&(strcmp(vars{i},'id'))
        M = vars{i+1};
        vars(i:i+1) = [];
        break;
    end
end

if (numel(vars)==1)&&(isa(vars{1}, 'cell'))
    % We have a cell array of 1d points
    X = vars{1};
    d = numel(X);
    X = reshape(X, 1, d);
elseif (numel(vars)==2)&&(isa(vars{1}, 'tt_tensor')||isa(vars{1}, 'tt_matrix'))&&(isscalar(vars{2}))
    % First argument is a single x, second argument is the dimension
    d = vars{2};
    X = cell(1, d);
    for i=1:d
        X{i} = vars{1}; % make the same format as in the first case
    end
else
    % A set of independent x
    d = numel(vars);
    X = cell(1,d);
    for i=1:d
        if (isa(vars{i}, 'tt_tensor')||isa(vars{i}, 'tt_matrix'))
            X{i} = vars{i};
        else
            error('wrong input %d to tt_meshgrid_vert', i);
        end
    end
end

% Now X contains 1D tt_tensors. Expand them by ones

% Concat all n to the common storage
n = zeros(0,1);
% positions in n where each component starts
pos = ones(d+1,1);
for i=1:d
    n = [n; X{i}.n];
    pos(i+1) = pos(i)+X{i}.d;
end
% Expand
for i=1:d
    if (i>1)
        if (isempty(M))
            expand = tt_ones(n(1:pos(i)-1));
        else
            expand = M;
            if (i>2)
                expand = mtkron(repmat({M},1,i-1));
            end
        end
        X{i} = tkron(expand, X{i});
    end
    if (i<d)
        if (isempty(M))
            expand = tt_ones(n(pos(i+1):pos(d+1)-1));
        else
            expand = M;
            if (i<d-1)
                expand = mtkron(repmat({M},1,d-i));
            end
        end
        X{i} = tkron(X{i}, expand);
    end
end

end

author: oseledets

repo: TT-Toolbox


text:
function diagnostic = callkypd(F,h,options)

%F = constraint2kyp(F);
kyps = is(F,'kyp');
if all(kyps)
    kypConstraints = F;
    otherConstraints = lmi;
else
    kypConstraints   = F(find(kyps));
    otherConstraints = F(find(~kyps));
end

% Trivial case
if length(kypConstraints) == 0
    if ~isempty(options)
        options.solver = '';
    else
        options = sdpsettings;
    end
    diagnostic = solvesdp(F,h,options)
    return;
end

p_variables = [];
x_variables = [];
for i = 1:length(kypConstraints)
    [A{i},B{i},P{i},M{i},negated{i}] = extractkyp(sdpvar(kypConstraints(i)));
    if ~isa(P{i},'sdpvar') | ~isempty(intersect(p_variables,getvariables(P{i})))
        diagnostic.solvertime = 0;
        diagnostic.problem = -4;
        diagnostic.info = yalmiperror(-4);
        return;
    end
    P_base = getbase(P{i});
    [n,m] = size(P{i});
    if (n~=m) | (nnz(P_base)~=n*n) | (sum(sum(P_base))~=n*n)
        diagnostic.solvertime = 0;
        diagnostic.problem = -4;
        diagnostic.info = yalmiperror(-4);
        return;
    end
    p_variables = [p_variables getvariables(P{i})];
    x_variables = [x_variables getvariables(M{i})];
end

if intersect(p_variables,getvariables(h))
    diagnostic.solvertime = 0;
    diagnostic.problem = -4;
    diagnostic.info = yalmiperror(-4);
    return;
end

start = length(A);k = 1;
for i = 1:length(otherConstraints)
    if is(otherConstraints(i),'lmi')
        A{k+start}=[];
        B{k+start}=[];
        P{k+start}=[];
        M{k+start} = sdpvar(otherConstraints(i));
        x_variables = [x_variables getvariables(M{k+start})];
        k = k + 1;
    elseif is(otherConstraints(i),'elementwise')
        X = sdpvar(otherConstraints(i));
        for ii = 1:size(X,1)
            for jj = 1:size(X,2)
                A{k+start}=[];
                B{k+start}=[];
                P{k+start}=[];
                M{k+start} = X(ii,jj);
                x_variables = [x_variables getvariables(X(ii,jj))];
                k = k + 1;
            end
        end
    end
end

% Objective function variables
x_variables = unique([x_variables getvariables(h)]);

% Are P and x independent
if ~isempty(intersect(p_variables,x_variables))
    error
end

% Generate the M-matrices
for i = 1:length(M)
    m_variables = getvariables(M{i});
    n = length(M{i});
    M0{i} = getbasematrix(M{i},0);
    for j = 1:length(x_variables)
        Mi{i,j} = getbasematrix(M{i},x_variables(j));
    end
end

c = zeros(length(x_variables),1);
if ~isempty(h)
    for i = 1:length(x_variables)
        c(i) = getbasematrix(h,x_variables(i));
    end
end

for i = 1:length(M)
    matrixinfo.n(i) = length(A{i});
    matrixinfo.nm(i) = length(M{i});
end

matrixinfo.K  = length(x_variables);
matrixinfo.c  = c;
matrixinfo.N  = length(A);
matrixinfo.A  = A;
matrixinfo.B  = B;
matrixinfo.M0 = M0;
matrixinfo.M  = Mi;
matrixinfo.A  = A;

try
    solvertime = tic;
    [u,Popt,xopt,Zopt] = kypd_solver(matrixinfo,options);
    solvertime = toc(solvertime);

    % SAVE PRIMALS
    setsdpvar(recover(x_variables),xopt);
    for i = 1:length(kypConstraints)
        if negated{i}
            setsdpvar(P{i},-Popt{i});
        else
            setsdpvar(P{i},Popt{i});
        end
    end

    % SAVE DUALS
    lmi_index = [];
    j = 1;
    for i = 1:length(kypConstraints)
        %W = struct(kypConstraints(i));lmiid = W.LMIid;
        lmiid = getlmiid(kypConstraints(i));
        lmi_index = [lmi_index;lmiid];
        duals{j}=Zopt{j};j = j+1;
    end
    for i = 1:length(otherConstraints)
        %W = struct(otherConstraints(i));lmiid = W.LMIid;
        lmiid = getlmiid(otherConstraints(i));
        lmi_index = [lmi_index;lmiid];
        duals{j}=Zopt{j};j = j+1;
    end
    yalmip('setdual',lmi_index,duals);

    diagnostic.solvertime = solvertime;
    diagnostic.info = yalmiperror(0,'KYPD');
    diagnostic.problem = 0;
catch
    solvertime = etime(clock,solvertime);
    diagnostic.solvertime = solvertime;
    diagnostic.info = yalmiperror(9,'KYPD');
    diagnostic.problem = 9;
end

function Fout = constraint2kyp(F)
Fout = [];
for i = 1:length(F)
    if is(F(i),'lmi')
        [l,m,r] = factors(sdpvar(F(i)));
        [constants,general,singles,pairs] = classifyfactors(l,m,r);
        [constants,general,singles,pairs] = compressfactors2(constants,general,singles,pairs);
    else
        Fout = Foput + F(i);
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [ TP, FP, FN ] = measuresByConfidence( objects, min_confidence, all_TPs, all_FPs, TP_rep, FP )
%MEASURESBYCONFIDENCE Summary of this function goes here
%   Detailed explanation goes here

    countGT = 0;
    countFound = 0;

    % TP_rep: number of TP windows with repeated GT instances
    % TP: number of TP objects, without repeating GT instances

    nImgs = length(objects);
    for i = 1:nImgs

        %% Get TP valid samples and remove non-valid ones (by confidence)
        if(TP_rep > 0)
            this_pos = find(all_TPs(:,1) == i);
            more_conf = find(all_TPs(this_pos,4) >= min_confidence);

            nMoreConf = length(more_conf);
            nLess = length(this_pos) - nMoreConf;
            TP_rep = TP_rep - nLess;

            this_found = all_TPs(this_pos(more_conf),3); % label ids of detections with min_confidence
        else
            this_found = [];
        end

        %% Find unique TPs labels
        un_this_found = unique(this_found);
        nUnique = length(un_this_found);
        countFound = countFound + nUnique;
        lenGT = 0;

        for j = 1:length(objects(i).ground_truth)
            if(~isempty(objects(i).ground_truth(j).name))
                lenGT = lenGT+1;
            end
        end
        countGT = countGT + lenGT;
    end

    %% Remove non-valid FPs (by confidence)
    if(FP > 0)
        nLess = sum(all_FPs(:,4) < min_confidence);
        FP = FP - nLess;
    end

%     frac_NoObj = 1-countObjs/countTot;
%     detection_rate = countFound/countGT;

    FN = countGT - countFound;
    TP = countFound;

    % Add all repeated TPs to FPs
    FP =  FP + (TP_rep-TP);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function RHSdiv = divergenceTerm1D(F)
% This function calculates the divergence of a field using its face
% average value
%
% SYNOPSIS:
%   RHSdiv = divergenceTerm1D(F)
%
% PARAMETERS:
%   F: Face Variable
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nx = F.domain.dims(1);
G = 1:Nx+2;
DX = F.domain.cellsize.x(2:end-1);

% define the vector of cell index
row_index = reshape(G(2:Nx+1),Nx,1); % main diagonal (only internal cells)

% compute the divergence
div_x = (F.xvalue(2:Nx+1)-F.xvalue(1:Nx))./DX;

% define the RHS Vector
RHSdiv = zeros(Nx+2,1);

% assign the values of the RHS vector
RHSdiv(row_index) = reshape(div_x,Nx,1);

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%% Spherical Functions
%
% The |S2Fun| toolbox is a collection of Matlab classes for numerical computations with functions on the two-dimensional sphere.
% It overloads the default commands for vectors and matrices to compute the analogous operations for functions of the given type.
% The underlying mathematical approach is accomplished via spherical harmonics which form an orthonormal basis of the square-integrable functions on the two-dimensional sphere.
%
%%
%

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%% Twinning Tutorial
%
%%
%
%% Open in Editor
%
%% Contents
%


text:
function x = project_sample(x, P)

if ~isempty(P)
    x = cellfun(@(x, P) permute(mtimesx(permute(x, [4 3 1 2]), P, 'speed'), [3 4 2 1]), x, P, 'uniformoutput', false);
end

repo: CFWCR