text:
function sqi = xsqi(signal,qrs,fs,win)
%xSQI Extravagance SQI
%
% Despite the fancy name, this function does some pretty boring
% calculations, while attempting to describe how different a QRS complex is
% from the rest of the signal. This is particularly important for FECG
% signals, which are often buried into noise.
%
% Input:
%   signal:         single channel (F)ECG [1xN double]
%   qrs:            list with (F)QRS locations [1xNp double]
%   fs:             signal sampling frequency [Hz]
%   win:            half of the window length around (F)QRS complex [ms],
%                   same window is used to noise area
%
% Output:
%   sqi:            resulting xSQI for segment
%
% Fetal Extraction Toolbox, version 1.0, February 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014 Fernando Andreotti4
% Dresden University of Technology, Institute of Biomedical Engineering
% fernando.andreotti@mailbox.tu-dresden.de
%
% Last updated : 30-06-2016
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

%% Generate a template
% removing extremities detections
if size(qrs,2) > size(qrs,1); qrs = qrs';end
if size(signal,2) > size(signal,1); signal = signal';end

win =  2.*round(win*fs/2); % rounding to nearest even number
hwin = 0.5*win;
extremities = (qrs <= round(1.5*win) | qrs >= length(signal)-round(1.5*win));        % test if there are peaks on the border that may lead to error
qrs = round(qrs(~extremities));                                    % remove extremity peaks
if length(qrs) <3
    disp('xsqi: skipping due to low number of beats available')
    sqi = 0;
    return
end
% Stacking cycles
M = arrayfun(@(x) signal(x-hwin:x+hwin)'.^2,qrs,'UniformOutput',false);    % creates beat matrix
M = cell2mat(M);                                                        % converting cell output to array form (matrix is 2*width+1 x
N = arrayfun(@(x) signal([(x-hwin-win):(x-hwin) (x+hwin):(x+hwin+win)])'.^2,qrs,'UniformOutput',false);    % creates a surrounding matrix
N = cell2mat(N);                                                        % converting cell output to array form (matrix is 2*width+1 x

%% Check power
Psrd = median(median(N.^2)); % Power of sorroundings
Ppeak=median(median(M.^2)); % Power of peaks
sqi=Ppeak/(Psrd+Ppeak); % percentage of power that the peaks represent

end

author: fernandoandreotti


text:
%SPH   Implement spherical IBVS for point features
%
%  results = sph(T)
%  results = sph(T, params)
%
%  Simulate IBVS with for a square target comprising 4 points is placed
%  in the world XY plane. The camera/robot is initially at pose T and is
%  driven to the orgin.
%
%  Two windows are shown and animated:
%   1. The camera view, showing the desired view (*) and the
%      current view (o)
%   2. The external view, showing the target points and the camera
%
% The results structure contains time-history information about the image
% plane, camera pose, error, Jacobian condition number, error norm, image
% plane size and desired feature locations.
%
% The params structure can be used to override simulation defaults by
% providing elements, defaults in parentheses:
%
%   target_size    - the side length of the target in world units (0.5)
%   target_center  - center of the target in world coords (0,0,2)
%   niter          - the number of iterations to run the simulation (500)
%   eterm          - a stopping criteria on feature error norm (0)
%   lambda         - gain, can be scalar or diagonal 6x6 matrix (0.01)
%   ci             - camera intrinsic structure (camparam)
%   depth          - depth of points to use for Jacobian, scalar for
%                    all points, of 4-vector.  If null take actual value
%                    from simulation      ([])
%
% SEE ALSO: ibvsplot

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% IMPLEMENTATION NOTE
%
% 1.  As per task function notation (Chaumette papers) the error is
%     defined as actual-demand, the reverse of normal control system
%     notation.
% 2.  The gain, lambda, is always positive
% 3.  The negative sign is written into the control law

function results = sph(T0, params)

    %% default parameters

    % define default target geometry
    vw = 2; vh = 2;       % dimensions of target
    Tct_star = transl(0, 0, -2);    % desired pose of target wrt camera

    % define default IBVS parameters
    niter = 500;            % number of iterations
    eterm = 0;
    lambda = 0.04;         % control gain
    cam = scamera;
    depth = [];
    estim = [];

    if nargin == 2,
        % override default target size
        if isfield(params, 'target_size'),
            vw = params.target_size(1);
            vh = params.target_height(1);
        end
        % override default camera desired pose
        if isfield(params, 'target_center'),
            Tct_star = params.Tct_star;
        end
        % override default number of iterations
        if isfield(params, 'niter'),
            niter = params.niter;
        end
        % override default gain (can be scalar or diag 6x6)
        if isfield(params, 'lambda'),
            lambda = params.lambda;
        end
        % override default camera intrinsics
        if isfield(params, 'ci'),
            ci = params.ci;
        end
        % override default point depth
        if isfield(params, 'depth'),
            if length(params.depth) == 1,
                depth = params.depth * ones(4,1);
            else
                depth = params.depth;
            end
        end
        if isfield(params, 'estim'),
            estim = params.estim;
        end
    end

    % define feature points in XY plane
    %  make vertices of a square
    P = mkgrid(2, [vw vh]);

    % show the reference location, this is the view we wish to achieve
    % when Tc = Tct_star
    figure(1)
    clf
    P
    Tct_star
    % HACKuv_star = cam.plot(P, inv(Tct_star));    % create the camera view
    uv_star = cam.plot(P, (Tct_star));    % create the camera view
    k=1
    plot(uv_star(2,k), uv_star(1,k), '*');
    hold on
    for k=2:4
        plot(uv_star(2,k), uv_star(1,k), '*');
    end
    %plot(uv_star(:,1), uv_star(:,2), '*');  % show desired view
    %hold off
    %cam.plot(P, T0);    % show initial view
    pause(1)

    % this is the 'external' view of the points and the camera
    figure(2)
    plot3(P(1,:), P(2,:), P(3,:), '*')
    showcam(T0, P)
    grid
    xlabel('x');
    ylabel('y');
    zlabel('z');
    pause
    figure(1)
    %cam2 = showcamera(T0);
    %camup([0,-1,0]);

    %% initialize the vservo variables

    Tcam = T0;                % initial camera/robot pose

    % initialize some history variables, suffix _h, to hold time series data
    % about the IBVS run
    v_h = [];       % velocity demand to robot
    e_h = [];       % feature error
    en_h = [];      % scalar norm of the feature error
    c_h = [];       % Jacobian condition number
    uv_h = [];      % image plane feature coordinates
    Tcam_h = [];
    Z_h = [];
    PP = eye(1,1);
    theta = 0;
    smoothing = 0.95;

    for k=1:niter,
         % set and show the camera pose
        cam.Tcam = Tcam;
        %showcamera(cam2, Tcam);
        figure(2)
        plot3(P(1,:), P(2,:), P(3,:), '*')
        grid on
        xlabel('x');
        ylabel('y');
        zlabel('z');
        showcam(Tcam, P);
        figure(1)

        % compute the view
        uv = cam.plot(P);

        e = uv - uv_star;   % as per task function notation
        uv
        uv_star
        e
        if 1
            for i=1:numcols(e)
                % vectorize this
                if e(2,i) > pi
                    e(2,i) = e(2,i) - 2*pi;
                elseif e(2,i) < -pi
                    e(2,i) = 2*pi + e(2,i);
                end
            end
        end
        e
        e = reshape(e, numel(e), 1);

        if estim,
            % run the depth estimator
            if k > 1,

                % compute Jacobian for unit depth, r=1
                J = jac(uv, ones(numcols(uv),1));
                Jv = J(:,1:3);  % velocity part, depends on 1/z
                Jw = J(:,4:6);  % rotational part, indepedent of 1/z

                % estimate image plane velocity
                uv_d =  reshape(uv, 1, [])' - reshape(uv_p, 1, [])';

                % estimate coefficients for A (1/z) = B
                B = uv_d - Jw*v(4:6);
                A = Jv * v(1:3);

                eta = A\B;          % least squares solution

                % first order smoothing
                theta = (1-smoothing) * 1/eta + smoothing * theta;

                pt = transformp(inv(Tcam), P);

                fprintf('depth %.4g, est depth %.4g, rls depth %.4g\n', pt(3,1), 1/eta, theta);
                Z_h(k,:) = [theta pt(3,1)];   % estimated depth, true depth
            end
            uv_p = uv;
        end

        % compute the Jacobian
        if isempty(depth),
            % exact depth from simulation (not possible in practice)
            pt = transformp(inv(Tcam), P);
            J = jac(uv, colnorm(pt) );
        else
            J = jac(uv, depth );
        end
        J

        % compute the velocity of camera in camera frame

        v = -lambda * pinv(J) * e;
        fprintf('v: %.3f %.3f %.3f %.3f %.3f %.3f\n', v);
        norm(e)

        % update the camera pose
        Td = (eye(4,4) + delta2tr( v ) );    % differential motion
        Td = trnorm(Td);

        Tcam = Tcam * Td;       % apply it to current pose

        % update the history variables
        uv_h(k,:) = reshape(uv([2 1],:), 1, []);
        e_h(k,:) = e';
        en_h(k,:) = norm(e);
        c_h(k,:) = cond(J);
        v_h(k,:) = v';
        Tcam_h(:,:,k) = Tcam;

        if norm(e) < eterm,
            fprintf('completed on error tolerance\n');
            break;
        end
        drawnow
        pause(0.05)
        pause
    end
    fprintf('completed on iteration count\n');

    results.uv = uv_h;
    results.uv_labels = {'\phi (rad)', '\theta (rad)'};
    results.c = c_h;
    results.e = e_h;
    results.en = en_h;
    results.v = v_h;
    results.Tcam = Tcam_h;
    %results.limits = [0 cam.nu 0 cam.nv];
    results.limits = [-pi pi 0 pi];
    results.uv_star = uv_star([2 1],:);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function au_qmr_ilu_m_i(A,mc,max_it_qmr,tol_qmr,tol_ilu,info_qmr)
%
%  Generates the data used in 'au_qmr_ilu_m'. Data are stored in global
%  variables.
%
%  Moreover, it defines the parameters mc, max_it_qmr, tol_qmr,
%  tol_ilu, and info_qmr, which are needed later in 'au_qmr_ilu_l_i' and
%  'au_qmr_ilu_s_i'.
%
%  Calling sequence:
%
%    au_qmr_ilu_m_i(A,mc,max_it_qmr,tol_qmr,tol_ilu,info_qmr)
%
%  Input:
%
%    A          real matrix;
%    mc         (= 'M' or 'C') Optimize for memory ('M') or computation ('C')
%               when shifted systems are solved by 'au_qmr_ilu_l/s' later.
%               Optional (default value: 'M');
%    max_it_qmr maximal number of QMR steps allowed in 'au_qmr_ilu_l/s'.
%               Corresponds to parameter 'MAXIT' in MATLAB function 'qmr'.
%               Optional (default value: min([n-1,100]));
%    tol_qmr    stopping tolerance (w.r.t. normalized residual) for QMR
%               iteration in 'au_qmr_ilu_l/s'. Corresponds to parameter 'TOL'
%               in MATLAB function 'qmr'. Optional (default value: 1e-14);
%    tol_ilu    dropping tolerance for ILU preconditioner. Corresponds to
%               parameter 'DROPTOL' in MATLAB function 'luinc'. Optional
%               (default value: 1e-2);
%    info_qmr   (= 0, 1, ...) If 1, warnings are printed when parameter FLAG
%               of MATLAB function 'qmr' is not equal to zero in QMR
%               iteration in 'au_qmr_ilu_l/s'. If 0, no warnings are
%               printed. Optional (default value: 1).
%
%  Remark:
%
%    If mc = 'M', the incomplete LU factors for A+p(i)*I are computed
%    any time 'au_qmr_ilu_s' is called, but they are not stored as global
%    variables to save memory. That means they are in general computed
%    several times, which results in additional computational cost. If
%    'mc = C', the incomplete LU factors are computed once by calling
%    'au_qmr_ilu_s_i' and stored as global variables, which results is
%    an increased memory demand.
%
%
%  LYAPACK 1.0 (Thilo Penzl, August 1999)

na = nargin;

if na<1 | na>6
  error('Wrong number of input arguments.');
end

if nargin < 2, mc = []; end
if nargin < 3, max_it_qmr = []; end
if nargin < 4, tol_qmr = []; end
if nargin < 5, tol_ilu = []; end
if nargin < 6, info_qmr = []; end

if ~length(mc), mc = 'M'; end
if ~length(max_it_qmr), max_it_qmr = min([size(A,1)-1,100]); end
if ~length(tol_qmr), tol_qmr = 1e-14; end
if ~length(tol_ilu), tol_ilu = 1e-2; end
if ~length(info_qmr), info_qmr = 1; end

global LP_A LP_MC LP_MAX_IT_QMR LP_TOL_QMR LP_TOL_ILU LP_INFO_QMR

LP_A = A;
LP_MC = mc;
LP_MAX_IT_QMR = max_it_qmr;
LP_TOL_QMR = tol_qmr;
LP_TOL_ILU = tol_ilu;
LP_INFO_QMR = info_qmr;

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_ballfun( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = pref.techPrefs.chebfuneps;

% Example 1 :
f = ballfun(@(x,y,z)x);
exact = ballfun(@(r,lam,th)r.*sin(th).*cos(lam),'spherical');
pass(1) = norm( f - exact ) < tol;

% Example 2 :
f = ballfun(@(x,y,z)y);
exact = ballfun(@(r,lam,th)r.*sin(th).*sin(lam),'spherical');
pass(2) = norm( f - exact ) < tol;

% Example 3 :
f = ballfun(@(x,y,z)z);
exact = ballfun(@(r,lam,th)r.*cos(th),'spherical');
pass(3) = norm( f - exact ) < tol;

% Example 4 :
f = ballfun(@(x,y,z)x.*z);
exact = ballfun(@(r,lam,th)r.*sin(th).*cos(lam).*r.*cos(th),'spherical');
pass(4) = norm( f - exact ) < tol;

% Example 5 :
f = ballfun(@(x,y,z)sin(x.*y.*z));
exact = ballfun(@(r,lam,th)sin(r.*sin(th).*cos(lam).*r.*sin(th).*sin(lam).*r.*cos(th)),'spherical');
pass(5) = norm( f - exact ) < tol;

% Example 6 : (Test 'coeffs' flag.) with rcos(th)
m = 10; n = 11; p = 12;
F = zeros(m,n,p);
F(2,floor(n/2)+1,floor(p/2))=1/2;F(2,floor(n/2)+1,floor(p/2)+2)=1/2;
f = ballfun(F, 'coeffs');
g = ballfun(@(r,lam,th)r.*cos(th),'spherical');
pass(6) = norm( f - g ) < tol;

% Example 7 :
f = ballfun(@(x,y,z)cos(x.*y.*z));
g = ballfun(@(r,lam,th)cos(r.*sin(th).*cos(lam).*r.*sin(th).*sin(lam).*r.*cos(th)),'spherical');
pass(7) = norm( f - g ) < tol;

% Example 8
f = ballfun(@(x,y,z)cos(x.*y.*z));
g = ballfun(@(x,y,z)cos(x.*y.*z),'vectorize');
pass(8) = norm( f - g ) < tol;

% Example 9 :
f = ballfun(@(x,y,z)x.*z);
exact = ballfun(@(r,lam,th)r.*sin(th).*cos(lam).*r.*cos(th),'spherical');
pass(9) = norm( f - exact ) < tol;

% Example 10 :
f = ballfun(@(x,y,z)sin(x.*y.*z));
exact = ballfun(@(r,lam,th)sin(r.*sin(th).*cos(lam).*r.*sin(th).*sin(lam).*r.*cos(th)),'spherical');
pass(10) = norm( f - exact ) < tol;

% Example 11:
f = ballfun(@(x,y,z)x*y,'vectorize');
exact = ballfun(@(x,y,z)x.*y);
pass(11) = norm( f - exact ) < tol;

% Example 12:
f = ballfun(@(x,y,z)x*z*y,'vectorize');
exact = ballfun(@(x,y,z)x.*z.*y);
pass(12) = norm( f - exact ) < tol;

% Example 13
f = ballfun(@(x,y,z)cos(x.*y),[50,51,52]);
exact = ballfun(@(x,y,z)cos(x.*y));
pass(13) = norm( f - exact ) < tol;

% Example 14
exact = ballfun(@(x,y,z)cos(x.*y));
f = ballfun(exact, [51,50,49], 'vectorize');
pass(14) = norm( f - exact ) < tol;

% Example 15
exact = ballfun(@(x,y,z)sin(x.*y));
E = coeffs3(exact,49,50,51);
f = ballfun(E, 'coeffs', [51,53,51]);
pass(15) = norm( f - exact ) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [IMAGE_3D_DATA] = image3Ddata(M)
% function [IMAGE_3D_DATA]=image3Ddata(M)
% ------------------------------------------------------------------------
%
% This simple function creates a structure array containing coordinate and
% colour data for 3D images. It allows one to use the patch function to
% plot the whole image or a selection of voxels in 3D.
%
% N.B. The function has not been optimised for large images. Large images
% (function has been tested for images under 100x100x100) may produce
% memory problems.
%
%
% EXAMPLE
%
% M = rand(15,15,15);
% [IMAGE_3D_DATA] = image3D(M);
%
% Getting faces and vertices for full image:
% voxel_no=1:1:numel(M);
% voxel_face_no=IMAGE_3D_DATA.voxel_patch_face_numbers(voxel_no,:);
% M_faces=IMAGE_3D_DATA.voxel_patch_faces(voxel_face_no,:);
% M_vertices=IMAGE_3D_DATA.corner_coordinates_columns_XYZ;
%
% Getting faces and vertices for selection of voxels:
% voxel_no2=M>0.95;
% voxel_face_no2=IMAGE_3D_DATA.voxel_patch_face_numbers(voxel_no2,:);
% M_faces2=IMAGE_3D_DATA.voxel_patch_faces(voxel_face_no2,:);
% M_vertices2=IMAGE_3D_DATA.corner_coordinates_columns_XYZ;
%
% figure;fig=gcf; clf(fig); colordef (fig, 'white'); units=get(fig,'units'); set(fig,'units','normalized','outerposition',[0 0 1 1]); set(fig,'units',units);
% set(fig,'Color',[1 1 1]);
%
% subplot(1,2,1);
% hp=patch('Faces',M_faces,'Vertices',M_vertices,'EdgeColor','black', 'CData',IMAGE_3D_DATA.voxel_patch_CData(voxel_face_no,:),'FaceColor','flat');
% hold on; view(45,30); axis equal; axis tight; colormap jet; colorbar; caxis([0 1]);
% xlabel('J'); ylabel('I'); zlabel('K');
% title('Full image');
%
% subplot(1,2,2);
% hp2=patch('Faces',M_faces2,'Vertices',M_vertices2,'EdgeColor','black', 'CData',IMAGE_3D_DATA.voxel_patch_CData(voxel_face_no2,:),'FaceColor','flat');
% hold on; view(45,30); axis equal; axis tight; colormap jet; colorbar; caxis([0 1]); grid on;
% set(hp2,'FaceAlpha',0.8);
% xlabel('J'); ylabel('I'); zlabel('K');
% title('Selection of voxels');
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 11/05/2009
% ------------------------------------------------------------------------

% Copyright (c) 2009, Kevin Mattheus Moerman
% All rights reserved.
%
% 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.

% Setting up meshgrid of voxel centre coordinates
[X,Y,Z] = meshgrid(1:1:size(M,2),1:1:size(M,1),1:1:size(M,3));
IMAGE_3D_DATA.center_coordinates_meshgrid_X=X;
IMAGE_3D_DATA.center_coordinates_meshgrid_Y=Y;
IMAGE_3D_DATA.center_coordinates_meshgrid_Z=Z;

X=X(:); Y=Y(:); Z=Z(:);
IMAGE_3D_DATA.center_coordinates_columns_XYZ=[X Y Z];

% Creating coordinates for voxel corners
[X,Y,Z] = meshgrid(0.5:1:(size(M,2)+0.5),0.5:1:(size(M,1)+0.5),0.5:1:(size(M,3)+0.5));

IMAGE_3D_DATA.corner_coordinates_meshgrid_X=X;
IMAGE_3D_DATA.corner_coordinates_meshgrid_Y=Y;
IMAGE_3D_DATA.corner_coordinates_meshgrid_Z=Z;

X=X(:); Y=Y(:); Z=Z(:);
IMAGE_3D_DATA.corner_coordinates_columns_XYZ=[X Y Z];

clear X Y Z;

% Creating path face and color data
nodes_first_voxel = [    1 ...
                         2 ...
                        ( ( (size(M,1)+1)*(size(M,2)+1) ) +2) ...
                        ( ( (size(M,1)+1)*(size(M,2)+1) ) +1) ...
                        ( 1 + (size(M,1)+1) ) ...
                        ( 2 + (size(M,1)+1) ) ...
                        ( 2 + (size(M,1)+1) + ( (size(M,1)+1)*(size(M,2)+1) ) ) ...
                        ( 1 + (size(M,1)+1) + ( (size(M,1)+1)*(size(M,2)+1) ) )   ];

nodes_first_row_voxels=((0:1:(size(M,1)-1))' * ones(1,8)) + (ones(1,size(M,1))') * nodes_first_voxel;
A = repmat(nodes_first_row_voxels,(size(M,2)),1);
B = repmat(((size(M,1)+1)*(0:(size(M,2)-1))),(size(M,1)),1);
B = reshape(B,1, numel(B))' *ones(1,8);
nodes_first_slice_voxels=A+B;
A = repmat(nodes_first_slice_voxels,(size(M,3)),1);
B = repmat((((size(M,1)+1)*(size(M,2)+1))*(0:(size(M,3)-1))),((size(M,1))*(size(M,2))),1);
B = reshape(B,1,numel(B))'*ones(1,8);

IMAGE_3D_DATA.corner_numbers=A+B;
IMAGE_3D_DATA.voxel_patch_face_numbers=reshape(1:1:(6*numel(M)),6,numel(M))';
IMAGE_3D_DATA.voxel_patch_CData=reshape(((M(:)*ones(1,6)))',(6*numel(M)),1);

face_no=[1 2 3 4;1 2 6 5;2 3 7 6;3 4 8 7;1 4 8 5;5 6 7 8]';
faces=reshape((IMAGE_3D_DATA.corner_numbers(:,face_no))',4,[])';
IMAGE_3D_DATA.voxel_patch_faces=faces;

repo: Aero-matlab


text:
function bessel_i0_values_test ( )

%*****************************************************************************80
%
%% BESSEL_I0_VALUES_TEST demonstrates the use of BESSEL_I0_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_I0_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_I0_VALUES stores values of \n' );
  fprintf ( 1, '  the Bessel I0 function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            I0(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_i0_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = spwelch(s,w, varargin)
% PWELCH   (* DEPRECATED*) overloaded pwelch for spectralobjects
%   spwelch(spectralobject, waveform) - plots the spectral density
%       Pxx = pwelch(spectralobject, waveform) - returns the Power Spectral
%           Density (PSD) estimate, Pxx, of a discrete-time signal
%           vector X using Welch's averaged,  modified periodogram method.
%       [Pxx, Freqs] = pwelch(spectralobject, waveform) - returns spectral
%       density and associated frequency bins.
%
%       Options, spwelch(s,w, 'DEFAULT') - plots the spectral density using
%       pwelch's defaults (8 averaged windows, 50% overlap)
%   window is length of entire waveform..
%
% THIS FUNCTION HAS BEEN REPLACED WITH PWELCH, and will be removed from
% future versions.

% AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

error('Spectralobject:spwelch:Depricated',...
    ['spectralobject/spwelch has been deprecated. ',...
    'Please use spectralobject/pwelch instead.  the syntax is the same']);

author: geoscience-community-codes

repo: GISMO


text:
function bessel_i0_int_values_test ( )

%*****************************************************************************80
%
%% BESSEL_I0_INT_VALUES_TEST demonstrates the use of BESSEL_IO_INT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_I0_INT_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_I0_INT_VALUES stores values of \n' );
  fprintf ( 1, '  the integral of the Bessel function I0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_i0_int_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
function [infos, auxillary_toolboxes] = tapas_get_toolbox_infos()
	% Function storing information regarding toolboxes.
	%
	% Each toolbox is represented by a field with its lower-case-name.
	% The field has the following subfields:
	%   init_function   [character vector]
	%           If specified, use that function to initialze the toolbox (and do
	%           not anything else apart from adding init_dir to the path).
	%           If not specified, add init_dir and all subfolders
	%   init_dir        [character vector]
	%           If init function is specified; the directory, where it can be
	%           found (has to be added to MATLABPATH to call the function)
	%           If init function is not specified, the root dir of the toolbox,
	%           which (and all subdirs) will be added.
	%   dependencies   [character vector]
	%           Names of dependent toolboxes. Needs to be lower case (the field-
	%           name in the struct). At the moment, we cannot check for depen-
	%           dencies of dependencies, so they should be specified as depen-
	%           dencies as well.
	%   diagnose_files [character vector/cell array of character vectors]
        %           Name(s) of files to be checked by TAPAS_diagnose('toolboxName').
	%
	% muellmat@ethz.ch
	% copyright (C) 2022

	infos = struct();
    %% Order determines init-order: do that by alphabet for main toolboxes
	infos = tapas_default_toolbox_info(infos, 'ceode');

    infos.hgf.init_function = '';
	infos.hgf.init_dir = 'HGF';
	infos.hgf.dependencies = [];
	infos.hgf.diagnose_files = '';

	infos.huge.init_function = 'tapas_huge_compile';
	infos.huge.init_dir = 'huge';
	infos.huge.dependencies = [];
	infos.huge.diagnose_files = '';

	infos.physio.init_dir = strcat('PhysIO',filesep,'code');
	infos.physio.init_function = 'tapas_physio_init';
	infos.physio.dependencies = [];
	infos.physio.diagnose_files = 'tapas_physio_main_create_regressors';

	infos = tapas_default_toolbox_info(infos,'rDCM');

	infos = tapas_default_toolbox_info(infos,'sem');

	infos.uniqc.init_dir = 'UniQC';
	infos.uniqc.init_function = '';
	infos.uniqc.dependencies = [];
	infos.uniqc.diagnose_files = 'tapas_uniqc_demo_fmri_qa'; % in subfolder demo/MrSeries

	%% Auxillary toolboxes:
    infos.external.init_function = '';
	infos.external.init_dir = 'external';
	infos.external.dependencies = '';
	infos.external.diagnose_files = '';

    infos = tapas_default_toolbox_info(infos,'tools'); % Want to have that?

    auxillary_toolboxes = {'external','tools'}; % Just inform - let init decide what to do
end

function infos = tapas_default_toolbox_info(infos,folderName)
    % Setting default toolbox info (no init function, no dependencies).
    fld_name = lower(folderName);
    if ismember(fld_name,fieldnames(infos))
        warning('Infos for toolbox %s already there!\n',fld_name);
    end
    infos.(fld_name).init_function = '';
    infos.(fld_name).init_dir = folderName;
    infos.(fld_name).dependencies = '';
    infos.(fld_name).diagnose_files = '';

end

author: translationalneuromodeling

repo: tapas


text:
clear all
n = 0:40;
h = exp(-0.1*n);
x = exp(-0.2*n);
y = conv(h,x);
subplot(3,1,1);stem(h);title("h")
subplot(3,1,2);stem(x);title("x")
subplot(3,1,3);stem(y);title("y")

author: liyanluminary


text:
%%***************************************************************
%% linsysolve: solve linear system to get dy, and direction
%%             corresponding to unrestricted variables.
%%
%% [xx,coeff,L,resnrm] = linsysolve(schur,UU,Afree,EE,rhs);
%%
%% child functions: symqmr.m, mybicgstable.m, linsysolvefun.m
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***************************************************************

   function [xx,coeff,L,resnrm] = linsysolve(par,schur,UU,Afree,EE,rhs);

    global solve_ok existspcholsymb exist_analytic_term
    global nnzmat nnzmatold matfct_options matfct_options_old use_LU
    global Lsymb msg diagR diagRold perturb_opt

    spdensity  = par.spdensity;
    printlevel = par.printlevel;
    iter       = par.iter;

    m = length(schur);
    if (iter==1);
       use_LU = 0; matfct_options_old = '';
       diagR = ones(m,1); perturb_opt = 1;
    end
    if isempty(nnzmatold); nnzmatold = 0; end
    diagRold = diagR;
%%
%% schur = schur + rho*diagschur + lam*AAt
%%
    diagschur = full(abs(diag(schur)));
    if (min(diagR) < 1e-6) & (max(diagR) > 1e-2);
       perturb_opt = 2;
    end
    if (perturb_opt == 2) & (iter > 10)
       pertdiag = min(1e-2,1e-15*max(1,diagschur)./max(1e-12,diagR));
       idx = find(diagschur < 1e-6);
       if (norm(rhs(idx)) < 1e-6) & (max(diagschur) > 1e-2)
          pertdiag(idx) = 1e6*ones(length(idx),1);
       end
       mexschurfun(schur,pertdiag);
       if (printlevel); fprintf('*'); end
    else
       diagschur = max(0,full(diag(schur)));
       minrho(1) = max(1e-15, 1e-6/3.0^iter);
       minlam(1) = max(1e-10, 1e-4/2.0^iter);
       minrho(2) = max(1e-04, 1/1.5^iter);
       rho = min(minrho(1),minrho(2)*(1+norm(rhs))/(1+norm(diagschur.*par.y)));
       lam = min(minlam(1),0.1*rho*norm(diagschur)/par.normAAt);
       if (exist_analytic_term); rho = 0; end
       mexschurfun(schur,rho*diagschur);
       mexschurfun(schur,lam*par.AAt);
       %%if (printlevel); fprintf(' %2.1e %2.1e ',rho,lam); end
    end
%%
%% assemble coefficient matrix
%%
    len = size(Afree,2);
    if ~isempty(EE)
       EE(:,[1 2]) = len + EE(:,[1 2]); %% adjust for ublk
    end
    EE = [(1:len)' (1:len)' zeros(len,1); EE];
    if isempty(EE)
       coeff.mat22 = [];
    else
       coeff.mat22 = spconvert(EE);
    end
    coeff.mat12 = [Afree, UU];
    coeff.mat11 = schur; %% important to use perturbed schur matrix
    ncolU = size(coeff.mat12,2);
%%
%% pad rhs with zero vector
%% decide which solution methods to use
%%
    rhs = [rhs; zeros(m+ncolU-length(rhs),1)];
    if (ncolU > 300); use_LU = 1; end
%%
%% Cholesky factorization
%%
    L = []; resnrm = []; xx = inf*ones(m,1);
    if (~use_LU)
       solve_ok = 1;  solvesys = 1;
       nnzmat = mexnnz(coeff.mat11);
       nnzmatdiff = (nnzmat ~= nnzmatold);
       if (nnzmat > spdensity*m^2) | (m < 500)
          matfct_options = 'chol';
       else
          if (par.matlabversion < 7.3)
             matfct_options = 'spchol';
	  else
             matfct_options = 'spcholmatlab';
          end
       end
       if (printlevel>2); fprintf(' %s',matfct_options); end
       if strcmp(matfct_options,'chol')
          if issparse(schur); schur = full(schur); end;
          if (iter<=5); %% to fix strange anonmaly in Matlab
             mexschurfun(schur,1e-20,2);
          end
          L.matfct_options = 'chol';
          L.perm = [1:m];
          [L.R,indef] = chol(schur);
          diagR = diag(L.R).^2;
          %%semilogy(diagschur,'*'); hold on;
          %%semilogy(diagR,'or'); hold off; pause(0.1);
          if (indef)
             diagR = diagRold;
 	     solve_ok = -2; solvesys = 0;
             msg = 'linsysolve: Schur complement matrix not pos. def.';
             if (printlevel); fprintf('\n  %s',msg); end
          end
       elseif strcmp(matfct_options,'spcholmatlab')
          if ~issparse(schur); schur = sparse(schur); end;
          L.matfct_options = 'spcholmatlab';
          [L.R,indef,L.perm] = chol(schur,'vector');
          L.Rt = L.R';
          %%diagR(L.perm(1:length(L.R)),1) = diag(L.R).^2;
          if (indef)
             diagR = diagRold;
 	     solve_ok = -2; solvesys = 0;
             msg = 'linsysolve: Schur complement matrix not pos. def.';
             if (printlevel); fprintf('\n  %s',msg); end
          end
       elseif strcmp(matfct_options,'spchol')
          if ~issparse(schur), schur = sparse(schur); end;
          if (nnzmatdiff | ~strcmp(matfct_options,matfct_options_old))
             [Lsymb,flag] = symbcholfun(schur,par.cachesize);
             if (flag)
                solve_ok = -2;  solvesys = 0;
                existspcholsymb = 0;
                msg = 'linsysolve: symbolic factorization fails';
                if (printlevel); fprintf('\n  %s',msg); end
                use_LU = 1;
             else
                existspcholsymb = 1;
             end
          end
          if (existspcholsymb)
             L = sparcholfun(Lsymb,schur);
             L.matfct_options  = 'spchol';
             L.d(find(L.skip)) = 1e20;
             if any(L.skip) & (ncolU)
                solve_ok = -3; solvesys = 0;
                existspcholsymb = 0;
                use_LU = 1;
                if (printlevel)
                   fprintf('\n  L.skip exists but ncolU > 0.');
                   fprintf('\n  switch to LU factor.');
                end
             end
          end
       end
       if (solvesys)
          if (ncolU)
             tmp = coeff.mat12'*linsysolvefun(L,coeff.mat12)-coeff.mat22;
	     if issparse(tmp); tmp = full(tmp); end
             tmp = 0.5*(tmp + tmp');
             [L.Ml,L.Mu,L.Mp] = lu(tmp);
             tol = 1e-16;
             idx = find(abs(diag(L.Mu)) < tol);
             if ~isempty(idx) & (printlevel); fprintf('**'); end
          end
          [xx,resnrm,solve_ok] = symqmr(coeff,rhs,L);
          if (solve_ok<=0) & (printlevel)
             fprintf('\n  warning: symqmr fails: %3.1f.',solve_ok);
          end
       end
       if (solve_ok < 0)
          if (m < 5000 & strcmp(matfct_options,'chol')) | ...
             (m < 1e5 & strcmp(matfct_options,'spchol')) | ...
             (m < 1e5 & strcmp(matfct_options,'spcholmatlab'))
             use_LU = 1;
             if (printlevel); fprintf('\n  switch to LU factor.'); end
          end
       end
    end
%%
%% LU factorization
%%
    if (use_LU)
       nnzmat = mexnnz(coeff.mat11)+mexnnz(coeff.mat12);
       nnzmatdiff = (nnzmat ~= nnzmatold);
       solve_ok = 1; solvesys = 1;
       if ~isempty(coeff.mat22)
          raugmat = [coeff.mat11, coeff.mat12; coeff.mat12', coeff.mat22];
       else
          raugmat = coeff.mat11;
       end
       if (nnzmat > spdensity*m^2) | (m+ncolU < 500)
          matfct_options = 'lu';
       else
          matfct_options = 'splu';
       end
       if (printlevel > 2); fprintf(' %s ',matfct_options); end
       if strcmp(matfct_options,'lu')
          if issparse(raugmat); raugmat = full(raugmat); end
          [L.l,L.u,L.p] = lu(raugmat);
          L.matfct_options = 'lu';
          L.p = sparse(L.p);
          idx = find(abs(diag(L.u)) < 1e-20);
          if ~isempty(idx)
             msg = 'linsysolve: matrix is singular';
             if (printlevel); fprintf('\n  %s',msg); end
             solvesys = 0;
          end
          [ii,jj] = find(L.p); [dummy,idx] = sort(ii); L.perm = jj(idx);
       end
       if strcmp(matfct_options,'splu')
          if ~issparse(raugmat); raugmat = sparse(raugmat); end
          if (nnzmatdiff | ~strcmp(matfct_options,matfct_options_old))
             Lsymb.perm = symamd(raugmat);
          end
          L.perm = Lsymb.perm;
          L.matfct_options = 'splu';
          if (par.matlabversion >= 6.5)
             [L.l,L.u,L.p,L.q] = lu(raugmat(L.perm,L.perm));
          else
             [L.l,L.u,L.p] = lu(raugmat(L.perm,L.perm));
             L.q = speye(length(raugmat));
          end
       end
       if (solvesys)
          [xx,resnrm,solve_ok] = mybicgstab(coeff,rhs,L,[],[],printlevel);
          if (solve_ok<=0) & (printlevel)
             fprintf('\n  warning: bicgstab fails: %3.1f,',solve_ok);
          end
       end
    end
    if (printlevel>=3); fprintf(' %2.0d',length(resnrm)-1); end
%%
    nnzmatold = nnzmat; matfct_options_old = matfct_options;
%%***************************************************************

author: yu-jiang

repo: radpbook

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


text:
clear all; close all; clc

n=100; L=4;
x=linspace(0,L,n);
f=(x.^2).';  % parabola with 100 data points

M=20;  % polynomical degree
for j=1:M
  phi(:,j)=(x.').^(j-1);  % build matrix A
end

for j=1:4
    fn=(x.^2+0.1*randn(1,n)).';
    an=pinv(phi)*fn; fna=phi*an;  % least-square fit
    En=norm(f-fna)/norm(f);
    subplot(4,2,4+j),bar(an,'FaceColor',[.6 .6 .6],'EdgeColor','k')
end

subplot(2,1,1), plot(x,f,'k'), hold on

%% different regressions
subplot(2,1,1)
lambda=0.1; phi2=phi(:,2:end);
for jj=1:100
  f=(x.^2+0.2*randn(1,n)).';
  a1=pinv(phi)*f; f1=phi*a1; E1(jj)=norm(f-f1)/norm(f);
  a2=phi\f; f2=phi*a2; E2(jj)=norm(f-f2)/norm(f);
  [a3,stats]=lasso(phi,f,'Lambda',lambda); f3=phi*a3; E3(jj)=norm(f-f3)/norm(f);
  [a4,stats]=lasso(phi,f,'Lambda',lambda,'Alpha',0.8); f4=phi*a4; E4(jj)=norm(f-f4)/norm(f);
a5=robustfit(phi2,f);f5=phi*a5;E5(jj)=norm(f-f5)/norm(f);
a6=ridge(f,phi2,0.5,0);f6=phi*a6;E6(jj)=norm(f-f6)/norm(f);

  A1(:,jj)=a1;A2(:,jj)=a2;A3(:,jj)=a3;A4(:,jj)=a4;A5(:,jj)=a5;A6(:,jj)=a6;
  plot(x,f), hold on
end
Err=[E1; E2; E3; E4; E5; E6];
Err2=[E1; E2; E3; E4; E5];

figure(2)
subplot(3,2,1), boxplot(A1.'), set(gca,'Xtick',1:20,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20'})
subplot(3,2,2), boxplot(A2.'), set(gca,'Xtick',1:20,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20'})
subplot(3,2,3), boxplot(A3.'), set(gca,'Xtick',1:20,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20'})
subplot(3,2,4), boxplot(A4.'), set(gca,'Xtick',1:20,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20'})
subplot(3,2,5), boxplot(A5.'), set(gca,'Xtick',1:20,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20'})
subplot(3,2,6), boxplot(A6.'), set(gca,'Xtick',1:20,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20'})

figure(3)
subplot(3,3,1), boxplot(Err.'), axis([0 7 0.02 0.07])

%%
clear phi
M=10;
En=zeros(100,M);
for jj=1:M
  for j=1:jj
    phi(:,j)=(x.').^(j-1);
  end
  f=(x.^2).';
  for j=1:100
    fn=(x.^2+0.1*randn(1,n)).';
    an=pinv(phi)*fn; fna=phi*an;
    En(j,jj)=norm(f-fna)/norm(f);
  end
end

figure(3)
subplot(3,3,2)
boxplot(En)

subplot(3,3,3)
boxplot(En), axis([0.5 10.5 0 0.008])

author: dynamicslab


text:
function pc = ckf_packed_pc(x,fmmparam)
% CKF_PACKED_PC - Pack P and C for the Cubature Kalman filter transform
%
% Syntax:
%   pc = CKF_PACKED_PC(x,fmmparam)
%
% In:
%   x - Evaluation point
%   fmmparam - Array of handles and parameters to form the functions.
%
% Out:
%   pc - Output values
%
% Description:
%   Packs the integrals that need to be evaluated in nice function form to
%   ease the evaluation. Evaluates P = (f-fm)(f-fm)' and C = (x-m)(f-fm)'.

% Copyright (c) 2010 Hartikainen, Särkkä, Solin
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.
%%

  f  = fmmparam{1};
  m  = fmmparam{2};
  fm = fmmparam{3};
  if length(fmmparam) >= 4
      param = fmmparam{4};
  end

  if ischar(f) || strcmp(class(f),'function_handle')
      if ~exist('param','var')
         F = feval(f,x);
      else
         F = feval(f,x,param);
      end
  elseif isnumeric(f)
         F = f*x;
  else
      if ~exist('param','var')
         F = f(x);
      else
         F = f(x,param);
      end
  end
  d = size(x,1);
  s = size(F,1);

  % Compute P = (f-fm)(f-fm)' and C = (x-m)(f-fm)'
  % and form array of [vec(P):vec(C)]
  pc = zeros(s^2+d*s,size(F,2));
  P = zeros(s,s);
  C = zeros(d,s);
  for k=1:size(F,2)
    for j=1:s
      for i=1:s
          P(i,j) = (F(i,k)-fm(i)) * (F(j,k) - fm(j));
      end
      for i=1:d
          C(i,j) = (x(i,k)-m(i)) * (F(j,k) - fm(j));
      end
    end
    pc(:,k) = [reshape(P,s*s,1);reshape(C,s*d,1)];
  end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function p = problem( x )
p = x.problem_;

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

repo: REST


text:
%% Load ground truth image and SR results from real and simulated data.
load('G:\results\baselineDatasets\coffee\Uncoded\mat\global_pan_translation_xyz_inlier_bin4_sr0_f17_win01.mat');
I_gt = srImages.groundTruth;

load('G:\results\baselineDatasets\coffee\Uncoded\mat\global_pan_translation_xyz_inlier_bin4_sr19_f17_win01.mat');
I_real = srImages.vdsr;

load('G:\results\simulatedDatasets\coffee\Uncoded\mat\global_pan_translation_xyz_inlier_bin4_sr19_f17_win01.mat');
I_simulated = srImages.vdsr;

%% Error map for real data.
ROI = [570, 420, 290, 400];
rangeIn = [1e-4 2e-1];
rangeOut = [0.15 0.85];
gamma = 0.9;

errorMap = abs(I_gt - I_real);
errorMap = errorMap(ROI(1):(ROI(1) + ROI(3) - 1), ROI(2):(ROI(2) + ROI(4) - 1));

figure;
I_real = I_real(ROI(1):(ROI(1) + ROI(3) - 1), ROI(2):(ROI(2) + ROI(4) - 1));
I_real = imadjust(I_real, [], [], gamma);
imshow(I_real);
imwrite(I_real, 'D:\Paper\cvpr2017\latex\images\teaser_real_SR.png');
figure;
errorMap = imadjust(errorMap, rangeIn, rangeOut);
imwrite(gray2ind(errorMap), jet, 'D:\Paper\cvpr2017\latex\images\teaser_real_diff.png');
imshow(errorMap);
colormap(gca, jet);

%% Error map for simulated data.
errorMap = abs(I_gt - I_simulated);
errorMap = errorMap(ROI(1):(ROI(1) + ROI(3) - 1), ROI(2):(ROI(2) + ROI(4) - 1));

figure;
I_simulated = I_simulated(ROI(1):(ROI(1) + ROI(3) - 1), ROI(2):(ROI(2) + ROI(4) - 1));
I_simulated = imadjust(I_simulated, [], [], gamma);
imshow(I_simulated);
imwrite(I_simulated, 'D:\Paper\cvpr2017\latex\images\teaser_simulated_SR.png');
figure;
errorMap = imadjust(errorMap, rangeIn, rangeOut);
imwrite(gray2ind(errorMap), jet, 'D:\Paper\cvpr2017\latex\images\teaser_simulated_diff.png');
imshow(errorMap);
colormap(gca, jet)

author: thomas-koehler

repo: SupER


text:
function external(action)

% call external (not within matlab) application,
% e.g. netscape

if isunix
	[s,w] = unix([action ' & ']);
else

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function geometry_test199 ( )

%*****************************************************************************80
%
%% TEST199 tests SHAPE_RAY_INT_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  nside = 6;
  ntest = 4;

  center(1:2,1) = [ 3.0; 0.0 ];
  p1(1:2,1) = [ 5.0; 0.0 ];
  pa_test = [ ...
    3.0,  0.0; ...
    3.0,  0.0; ...
    3.0, -1.0; ...
    3.0, -1.0 ]';
  pb_test = [ ...
    4.0,  0.0; ...
    3.0,  1.0; ...
    3.0,  1.0; ...
    7.0,  5.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST199\n' );
  fprintf ( 1, '  For a shape in 2D,\n' );
  fprintf ( 1, '  SHAPE_RAY_INT_2D computes the intersection of\n' );
  fprintf ( 1, '    a shape and a ray whose origin is within\n' );
  fprintf ( 1, '    the shape.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sides:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %d\n', nside );

  r8vec_print ( 2, center, '  Hexagon center:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Hexagon vertex #1:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %12f  %12f\n', p1(1:2,1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I       XA          YA          ' );
  fprintf ( 1, 'XB          YB          XI          YI\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    pa(1:2,1) = pa_test(1:2,i);
    pb(1:2,1) = pb_test(1:2,i);

    pint = shape_ray_int_2d ( center, p1, nside, pa, pb );

    fprintf ( 1, '  %6d  %10f  %10f  %10f  %10f  %10f  %10f\n', ...
      i, pa(1:2,1), pb(1:2,1), pint(1:2,1) );

  end

  return
end


text:
function view = deleteAllROIs(view)
%
% view = deleteAllROIs(view)
%
% Deletes the ROIs in view.ROIs.
%
% djh, 1/26/98

view.ROIs = {};
view.selectedROI = 0;

% Set the ROI popup menu
setROIPopup(view);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
plotset.pfrequ.binArray   = binArray;
plotset.pfrequ.plotallbin = plotallbin;
plotset.pfrequ.chanArray  = chanArray;
plotset.pfrequ.plotallch  = plotallch;
plotset.pfrequ.chanArray_MGFP = chanArray_MGFP;
plotset.pfrequ.blcorr     = blcorr;
plotset.pfrequ.xscale     = xxscale;
plotset.pfrequ.yscale     = yyscale;
plotset.pfrequ.linewidth  = linewidth;
plotset.pfrequ.isiy       = isiy;
plotset.pfrequ.fschan     = fschan;
plotset.pfrequ.fslege     = fslege;
plotset.pfrequ.fsaxtick   = fsaxtick;
%plotset.pfrequ.meap       = meap;
plotset.pfrequ.pstyle     = pstyle;
plotset.pfrequ.errorstd   = errorstd;
plotset.pfrequ.stdalpha   = stdalpha;
plotset.pfrequ.pbox       = pbox;
plotset.pfrequ.counterwin = counterwin;
plotset.pfrequ.holdch     = holdch;
plotset.pfrequ.yauto      = yauto;
plotset.pfrequ.yautoticks = yautoticks;
plotset.pfrequ.xautoticks = xautoticks;
plotset.pfrequ.binleg     = binleg;
plotset.pfrequ.chanleg    = chanleg;
plotset.pfrequ.isMGFP     = isMGFP;
plotset.pfrequ.ismaxim    = ismaxim;
%plotset.pfrequ.istopo     = istopo;
plotset.pfrequ.axsize     = axsize;
plotset.pfrequ.minorticks = minorticks;
plotset.pfrequ.linespec   = linespeci;
plotset.pfrequ.legepos    = legepos;
plotset.pfrequ.posgui     = posgui;
if ~isfield(plotset.pfrequ, 'posminigui')
        plotset.pfrequ.posminigui = [];
end
plotset.pfrequ.posfig     = posfig;

% adjusting
plotset.pfrequ.binArray   = plotset.pfrequ.binArray(plotset.pfrequ.binArray<=ERP.nbin);
plotset.pfrequ.chanArray  = plotset.pfrequ.chanArray(plotset.pfrequ.chanArray<=ERP.nchan);
plotset.pfrequ.chanArray_MGFP = plotset.pfrequ.chanArray_MGFP(plotset.pfrequ.chanArray_MGFP<=ERP.nchan);

plotset.pfrequ.ibckground  = ibckground;

% if plotset.pfrequ.xscale(1) < ERP.xmin*1000
%       plotset.pfrequ.xscale(1) = ceil(ERP.xmin*1000);
% end
% if plotset.pfrequ.xscale(2) > ERP.xmax*1000
%       plotset.pfrequ.xscale(2) = floor(ERP.xmax*1000);
% end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [ source_directory, list_of_images ] = createListOfImagesFromVOCOpts( path_to_voc_devkit, set_name )

addpath([path_to_voc_devkit,'/VOCcode']);
VOCopts = VOCInitFrom( path_to_voc_devkit );

image_ids      = textread(sprintf(VOCopts.imgsetpath,set_name),'%s');
list_of_images = cell(numel(image_ids),1);

for i = 1:numel(image_ids)
    list_of_images{i} = sprintf(VOCopts.imgpath,image_ids{i});
end

end

author: gidariss

repo: mrcnn-object-detection


text:
function coef = sandia_sgmgg_coef_naive ( dim_num, point_num, sparse_index )

%*****************************************************************************80
%
%% SANDIA_SGMGG_COEF_NAIVE returns the combinatorial coefficients.
%
%  Discussion:
%
%    The coefficient of point I is calculated as follows:
%
%    *) point J is a "neighbor" of point I if every entry of the sparse
%       index for point J is either equal to, or 1 greater than, the
%       corresponding entry of the sparse index of point I.
%
%    *) If point J is a neighbor of point I, then it contributes
%       (-1)^D to the coefficient, where D is the sum of the differences
%       between the sparse indices of point I and point J.
%
%    This is a completely naive implementation of the calculation,
%    intended simply as a demonstration for small examples.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    An Anisotropic Sparse Grid Stochastic Collocation Method for Partial
%    Differential Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2411-2442.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the vector.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer SPARSE_INDEX(DIM_NUM,POINT_NUM),
%    the indices that define the points.
%
%    Output, integer COEF(POINT_NUM), the coefficients.
%
  coef = zeros ( point_num, 1 );

  for j1 = 1 : point_num

    for j2 = 1 : point_num

      neighbor = 1;
      term = + 1;

      for i = 1 : dim_num

        dif = sparse_index(i,j2) - sparse_index(i,j1);

        if ( dif == 0 )

        elseif ( dif == 1 )
          term = - term;
        else
          neighbor = 0;
          break;
        end

      end

      if ( neighbor )
        coef(j1) = coef(j1) + term;
      end

    end

  end

  return
end


text:
function is_kgraph = is_kuratowski_graph(A,varargin)
% IS_KURATOWSKI_GRAPH Test if a graph can be collapsed to K_3,3 or K_5
%
% is_k = is_kuratowski_graph(A) checks if A can be reduced to K_3,3 or K_5
% by repeated edge contractions.  If so, then is_k=1, otherwise is_k=0.
%
% ... = is_straight_line_drawing(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   No additional options for this function
%
% Example:
%   is_kuratowski_graph(clique_graph(4)) % false, K_4 is not kuratowski
%   is_kuratowski_graph(clique_graph(5)) % true, K_5 is kuratowski
%   is_kuratowski_graph(clique_graph([3,3])) % true, K_3,3 is kuratowski

% David Gleich
% Copyright, Stanford University, 2008

%% History
%  2008-10-05: Initial version
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

options = struct();
options = merge_options(options,varargin{:});

if check, check_matlab_bgl(A,struct('sym',1)); end

is_kgraph = planar_test_mex(A,1);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


function y = l2lossForward( x,r )

delta = x - r;

y = sum(delta(:).^2);

y = y / (size(x,1)*size(x,2));

end


text:
function v = WindowFeatureNameCompare(f1,f2)

v = [];

if isempty(f1) || isempty(f2),
  return;
end

iscell1 = iscell(f1{1});
iscell2 = iscell(f2{1});

if iscell1 && iscell2,
  if numel(f1) ~= numel(f2),
    error('Inputs f1 and f2 do not match');
  else
    v = false(size(f1));
    for i = 1:numel(f1),
      if numel(f1{i}) == numel(f2{i}),
        vcurr = true;
        for j = 1:numel(f1{i}),
          if strcmp(class(f1{i}{j}),class(f2{i}{j})),
            if isnumeric(f1{i}{j}),
              if ~all(f1{i}{j} == f2{i}{j}),
                vcurr = false;
                break;
              end
            else
              if ~strcmp(f1{i}{j},f2{i}{j}),
                vcurr = false;
                break;
              end
            end
          end
        end
        v(i) = vcurr;
      end
    end
  end
  return;
end

if iscell2,
  v = WindowFeatureNameCompare(f2,f1);
  return;
end

if iscell1,
  v = false(size(f1));
  for i = 1:numel(f1),
    if numel(f1{i}) == numel(f2),
      vcurr = true;
      for j = 1:numel(f1{i}),
        if strcmp(class(f1{i}{j}),class(f2{j})),
          if isnumeric(f1{i}{j}),
            if ~all(f1{i}{j} == f2{j}),
              vcurr = false;
              break;
            end
          else
            if ~strcmp(f1{i}{j},f2{j}),
              vcurr = false;
              break;
            end
          end
        end
      end
      v(i) = vcurr;
    end
  end
  return;
end

v = true;
for j = 1:numel(f1),
  if strcmp(class(f1{j}),class(f2{j})),
    if isnumeric(f1{j}),
      if ~all(f1{j} == f2{j}),
        v = false;
        break;
      end
    else
      if ~strcmp(f1{j},f2),
        v = false;
        break;
      end
    end
  end
end

author: kristinbranson

repo: JAABA


text:
function bool = VMCPresent
%"VMCPresent"
%   Returns true only if the .exe files necessary to run VMC are present
%   and in the proper location for the current copy of CERR.
%
%   JRA 10/15/04
%
%Usage:
%   function bool = VMCPresent
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

bool = 0;
IMRTPdir = fileparts(which('IMRTP'));
if exist(fullfile(IMRTPdir, 'vmc++', 'bin', 'vmc_Linux.exe'), 'file') == 2 && ...
   exist(fullfile(IMRTPdir, 'vmc++', 'bin', 'vmc_Windows.exe'), 'file') == 2
   bool = 1;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function mi = imMutualInformation(img1, img2)
%IMMUTUALINFORMATION  Mutual information between two images
%
%   MI = imMutualInformation(IMG1, IMG2)
%   Computes the mutual information between two images.
%   The mutual information can be related to entropy and joint entropy:
%     MI = H1 + H2 - H12
%   where H1 and H2 are the entropies computed on images IMG1 and IMG2, and
%   H12 is the joint entropy of images IMG1 and IMG2.
%
%   The mutual information between two independant images corresponds to
%   the sum of the entropies computed for each image. The mutual
%   information of an image with itself equals the entropy of this image.
%
%   Example
%   % compute mutual information between an image and a shifted copy of
%   % itself
%     img = imread('cameraman.tif');
%     img2 = circshift(img, [2 3]);
%     MI = imMutualInformation(img, img2)
%     img3 = circshift(img, [5, 7]);
%     MI3 = imMutualInformation(img, img2)
%
%   % Check that we get same results by computing entropies (could be wrong
%   % when computed on double images)
%     img1 = imread('rice.png');
%     img2 = circshift(img1, [3, 4]);
%     imMutualInformation(img1, img2)
%     ans =
%         1.0551
%     h1 = imEntropy(img1);
%     h2 = imEntropy(img2);
%     h12 = imJointEntropy(img1, img2);
%     h1 + h2 - h12
%     ans =
%         1.0551

%   See also
%   imJointEntropy, imJointHistogram, imEntropy
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-08-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% joint histogram, and reduced hisograms
hist12  = imJointHistogram(img1, img2);
hist1   = sum(hist12, 2);
hist2   = sum(hist12, 1);

% normalisation of histograms
hist12  = hist12(hist12 > 0) / sum(hist12(:));
hist1   = hist1(hist1 > 0) / sum(hist1(:));
hist2   = hist2(hist2 > 0) / sum(hist2(:));

% entropy of each image
h12 = -sum(hist12 .* log2(hist12));
h1  = -sum(hist1  .* log2(hist1));
h2  = -sum(hist2  .* log2(hist2));

% compute mutual information
mi = h1 + h2 - h12;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function varargout = ma_InsertAerobrakeGUI(varargin)
% MA_INSERTAEROBRAKE MATLAB code for ma_InsertAerobrake.fig
%      MA_INSERTAEROBRAKE, by itself, creates a new MA_INSERTAEROBRAKE or raises the existing
%      singleton*.
%
%      H = MA_INSERTAEROBRAKE returns the handle to a new MA_INSERTAEROBRAKE or the handle to
%      the existing singleton*.
%
%      MA_INSERTAEROBRAKE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MA_INSERTAEROBRAKE.M with the given input arguments.
%
%      MA_INSERTAEROBRAKE('Property','Value',...) creates a new MA_INSERTAEROBRAKE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ma_InsertAerobrake_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ma_InsertAerobrake_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ma_InsertAerobrake

% Last Modified by GUIDE v2.5 20-Jan-2019 16:22:40

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ma_InsertAerobrake_OpeningFcn, ...
                   'gui_OutputFcn',  @ma_InsertAerobrake_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before ma_InsertAerobrake is made visible.
function ma_InsertAerobrake_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ma_InsertAerobrake (see VARARGIN)

% Choose default command line output for ma_InsertAerobrake
handles.output = hObject;
handles.ma_MainGUI = varargin{1};

% Update handles structure
guidata(hObject, handles);

if(length(varargin)>1)
    event = varargin{2};
    populateGUIWithEvent(handles, event);
end

dragModelCombo_Callback(handles.dragModelCombo, [], handles);

% UIWAIT makes ma_InsertAerobrake wait for user response (see UIRESUME)
uiwait(handles.ma_InsertAerobrakeGUI);

function populateGUIWithEvent(handles, event)
    set(handles.titleLabel, 'String', 'Edit Aerobrake');
    set(handles.ma_InsertAerobrakeGUI, 'Name', 'Edit Aerobrake');
    set(handles.nameText, 'String', event.name);

    set(handles.dragCoeffText,'String', num2str(event.dragCoeff,10));

    colorStr = getStringFromLineSpecColor(event.lineColor);
    colorValue = findValueFromComboBox(colorStr, handles.lineColorCombo);
    set(handles.lineColorCombo,'value',colorValue);

    styleStr = getLineStyleFromString(event.lineStyle);
    styleValue = findValueFromComboBox(styleStr, handles.aerobrakeLineStyleCombo);
 	set(handles.aerobrakeLineStyleCombo,'Value',styleValue);

    contents = handles.lineWidthCombo.String;
    contentsDouble = str2double(contents);
    ind = find(contentsDouble == event.lineWidth);
    set(handles.lineWidthCombo,'Value',ind);

    if(~isfield(event,'dragModel'))
        event.dragModel = 'Stock';
    end
    value = findValueFromComboBox(event.dragModel, handles.dragModelCombo);
    set(handles.dragModelCombo,'Value',value);

% --- Outputs from this function are returned to the command line.
function varargout = ma_InsertAerobrake_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = [];
    else
        name = get(handles.nameText,'String');
        dragCoeff = str2double(get(handles.dragCoeffText,'String'));

        contents = cellstr(get(handles.dragModelCombo,'String'));
        dragModel = contents{get(handles.dragModelCombo,'Value')};

        contents = cellstr(get(handles.lineColorCombo,'String'));
        colorStr = contents{get(handles.lineColorCombo,'Value')};
        lineSpecColor = getLineSpecColorFromString(colorStr);

        contents = cellstr(get(handles.aerobrakeLineStyleCombo,'String'));
        lineStyleStr = contents{get(handles.aerobrakeLineStyleCombo,'Value')};
        lineStyle = getLineStyleStrFromText(lineStyleStr);

        contents = handles.lineWidthCombo.String;
        contentsDouble = str2double(contents);
        contensInd = get(handles.lineWidthCombo,'Value');
        lineWidth = contentsDouble(contensInd);

        varargout{1} = ma_createAerobrake(name, dragCoeff, dragModel, lineSpecColor, lineStyle, lineWidth);
        close(handles.ma_InsertAerobrakeGUI);
    end

% --- Executes on button press in saveAndCloseButton.
function saveAndCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveAndCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.ma_InsertAerobrakeGUI);
    else
        msgbox(errMsg,'Errors were found while inserting an aerobrake.','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    dragCoeff = str2double(get(handles.dragCoeffText,'String'));
    enteredStr = get(handles.dragCoeffText,'String');
    numberName = 'Drag Coefficient';
    lb = 0;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(dragCoeff, numberName, lb, ub, isInt, errMsg, enteredStr);

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.ma_InsertAerobrakeGUI);

function nameText_Callback(hObject, eventdata, handles)
% hObject    handle to nameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of nameText as text
%        str2double(get(hObject,'String')) returns contents of nameText as a double

% --- Executes during object creation, after setting all properties.
function nameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to nameText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function dragCoeffText_Callback(hObject, eventdata, handles)
% hObject    handle to dragCoeffText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of dragCoeffText as text
%        str2double(get(hObject,'String')) returns contents of dragCoeffText as a double
newInput = get(hObject,'String');
newInput = attemptStrEval(newInput);
set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function dragCoeffText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to dragCoeffText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on ma_InsertAerobrakeGUI or any of its controls.
function ma_InsertAerobrakeGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to ma_InsertAerobrakeGUI (see GCBO)
% eventdata  structure with the following fields (see FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.ma_InsertAerobrakeGUI);
    end

% --- Executes on button press in openCdCalcButton.
function openCdCalcButton_Callback(hObject, eventdata, handles)
% hObject    handle to openCdCalcButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    ma_AerobrakingCdCalcGUI(handles.ma_MainGUI);

% --- Executes on selection change in dragModelCombo.
function dragModelCombo_Callback(hObject, eventdata, handles)
% hObject    handle to dragModelCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns dragModelCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from dragModelCombo
    contents = cellstr(get(hObject,'String'));
    dragModel = contents{get(hObject,'Value')};

    switch(dragModel)
        case 'Stock'
%             set(handles.cdLabel,'String','');
            set(handles.cdLabel,'String','m^2');
        case 'FAR'
            set(handles.cdLabel,'String','m^2');
        case 'NEAR'
            set(handles.cdLabel,'String','m^2');
    end

% --- Executes during object creation, after setting all properties.
function dragModelCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to dragModelCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in lineColorCombo.
function lineColorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to lineColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns lineColorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from lineColorCombo

% --- Executes during object creation, after setting all properties.
function lineColorCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lineColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in aerobrakeLineStyleCombo.
function aerobrakeLineStyleCombo_Callback(hObject, eventdata, handles)
% hObject    handle to aerobrakeLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns aerobrakeLineStyleCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from aerobrakeLineStyleCombo

% --- Executes during object creation, after setting all properties.
function aerobrakeLineStyleCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to aerobrakeLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in lineWidthCombo.
function lineWidthCombo_Callback(hObject, eventdata, handles)
% hObject    handle to lineWidthCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns lineWidthCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from lineWidthCombo

% --- Executes during object creation, after setting all properties.
function lineWidthCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lineWidthCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function treepack_test11 ( )

%*****************************************************************************80
%
%% TREEPACK_TEST11 tests TREE_RB_ENUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TREEPACK_TEST11\n' );
  fprintf ( 1, '  TREE_RB_ENUM enumerates the rooted binary trees on a \n' );
  fprintf ( 1, '  given number of nodes.\n' );
  fprintf ( 1, '\n' );

  for nnode = 0 : 11

    num = tree_rb_enum ( nnode );

    fprintf ( 1, '  %8d  %8d\n', nnode, num );

  end

  return
end

author: johannesgerer


text:
function [xPred, SPred]=sqrtDiscEKFPred(xPrev,SPrev,f,FJacob,SQ)
%SQRTDISCEKFPRED Perform the discrete-time prediction step that comes with
%                the square-root implementation of the first-order Extended
%                Kalman Filter (EKF).
%
%INPUTS: xPrev The xDim X 1 state estimate at the previous time-step.
%        SPrev The xDim X xDim lower-triangular square root of the  state
%              covariance matrix at the previous time-step.
%            f A function handle for the state transition function that
%              takes the state as its parameter.
%       FJacob A function handle for calculating the xDim X xDim state
%              transition matrix. If an empty matrix is passed, then FJacob
%              will be found using numerical differentiation  via the
%              numDiff function with default parameters.
%           SQ The xDimX xDim lower-triangular square root of the  process
%              noise covariance matrix.
%
%OUTPUTS: xPred The xDim X 1 predicted state estimate.
%         SPred The xDim X xDim lower-triangular square root of the
%               predicted state covariance estimate.
%
%The first-order EKF is summarized in Figure 10.3.3-1 in Chapter 10.3.3 of
%[1].
%
%The partial derivatives in the Jacobian matrix returned by the function
%FJacob are ordered
%[dF/dx(1), dF/dx(2),...,dF/dx(xDim)]
%That is, column i consists of partial derivatives with respect to element
%i of the x vector.
%
%The mathematics behind the specific square root implementation used here
%are described in [2].
%
%REFERENCES:
%[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%[2] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%March 2015, David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(xPrev,1);
if(isempty(FJacob))
    FJacob=@(x)numDiff(x,f,xDim);
end

xPred=f(xPrev);
F=FJacob(xPrev);
SPred=tria([F*SPrev,SQ]);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% -------------------------------------------------------------------------
function fn = getBatchDagNNWrapper(opts, useGpu)
% -------------------------------------------------------------------------
fn = @(imdb,batch) getBatchDagNN(imdb,batch,opts,useGpu) ;

% -------------------------------------------------------------------------
function inputs = getBatchDagNN(imdb, batch, opts, useGpu)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
im = imdb_get_batch_bcnn(images, opts, ...
                            'prefetch', nargout == 0);
labels = imdb.images.label(batch) ;
numAugments = size(im{1},4)/numel(batch);

labels = reshape(repmat(labels, numAugments, 1), 1, size(im{1},4));

if nargout > 0
  if useGpu
    im1 = gpuArray(im{1}) ;
    im2 = gpuArray(im{2}) ;
  else
      im1 = im{1};
      im2 = im{2};
  end
  inputs = {'input', im1, 'netb_input', im2, 'label', labels} ;
end

repo: DBCNN


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

function pass = test_overlap(pref)

% Grab some preferences:
if ( nargin == 0 )
    pref = chebfunpref();
end

% TODO: Test for array-valued CHEBFUNS and quasimatrices.

% Test empty input.
f = chebfun();
g = chebfun();
[f2, g2] = overlap(f, g);
pass(1) = isempty(f2) && isempty(g2);

% Check behavior for input chebfuns with differing domains.
f = chebfun(@sin, [-1 -.5 0 1], pref);
g = chebfun(@sin, [-2 0 0.5 2], pref);

try
    [f2, g2] = overlap(f, g);
    pass(2) = false;
catch ME
    pass(2) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN:overlap:domains');
end

% Check behavior in the basic case.
f = chebfun(@sin, [-1 -.5 0 0.5 1], pref);
g = chebfun(@sin, [-1 0 0.5 1], pref);

[f2, g2] = overlap(f, g);
xx = linspace(-1, 1);
pass(3) = isequal(f2.domain, g2.domain) && ...
    norm(feval(f, xx) - feval(f2, xx), inf) < 10*eps*vscale(f) && ...
    norm(feval(g, xx) - feval(g2, xx), inf) < 10*eps*vscale(g);

% [TODO]: Are these the same tests as above?
% Check correct behavior for higher-order impulses.  [TODO]:  Use a function
% with real higher-order impulses instead of just faking them like this.
[f2, g2] = overlap(f, g);
xx = linspace(-1, 1);
pass(4) = isequal(f2.domain, g2.domain) && ...
    norm(feval(f, xx) - feval(f2, xx), inf) < 10*eps*vscale(f) && ...
    norm(feval(g, xx) - feval(g2, xx), inf) < 10*eps*vscale(g);

[g2, f2] = overlap(g, f);
xx = linspace(-1, 1);
pass(5) = isequal(f2.domain, g2.domain) && ...
    norm(feval(f, xx) - feval(f2, xx), inf) < 10*eps*vscale(f) && ...
    norm(feval(g, xx) - feval(g2, xx), inf) < 10*eps*vscale(g);

%% Test on singular function: piecewise smooth chebfun - splitting on.

% define the domain:
dom = [-2 7];
domCheck = [dom(1)+0.1 dom(2)-0.1];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

pow1 = -0.3;
pow2 = -0.5;
op1 = @(x) (x - dom(2)).^pow1.*sin(100*x);
op2 = @(x) (x - dom(2)).^pow2.*cos(300*x);
f = chebfun(op1, dom, 'exps', [0 pow1], 'splitting', 'on');
g = chebfun(op2, dom, 'exps', [0 pow2], 'splitting', 'on');
[fout, gout] = overlap(f,g);
vals_fout = feval(fout, x);
vals_gout = feval(gout, x);
vals_f = feval(op1, x);
vals_g = feval(op2, x);

check = zeros(1,4);
check(1) = all( fout.domain == gout.domain );
check(2) = all( fout.domain == unique([f.domain, g.domain]) );
check(3) = ( norm(vals_fout - vals_f, inf) < ...
    1e4*eps*norm(vals_fout, inf) );
check(4) = ( norm(vals_gout - vals_g, inf) < ...
    1e3*eps*norm(vals_gout, inf) );

pass(6) = all( check );

%% Test for function defined on unbounded domain:

% Function defined on [0 Inf]:

% Specify the domain:
dom = [0 Inf];
domCheck = [0 100];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

% A decaying function:
opf = @(x) 0.75+sin(10*x)./exp(x);
f = chebfun(opf, dom, 'splitting', 'on');

% Blow-up function:
opg = @(x) x.*(5+exp(-x.^3))./(x - dom(1));
g = chebfun(opg, dom, 'exps', [-1 0]);

[fout, gout] = overlap(f, g);

vals_fout = feval(fout, x);
vals_gout = feval(gout, x);
vals_f = feval(opf, x);
vals_g = feval(opg, x);

check = zeros(1,4);
check(1) = all( fout.domain == gout.domain );
check(2) = all( fout.domain == unique([f.domain, g.domain]) );
check(3) = ( norm(vals_fout - vals_f, inf) < ...
    1e6*eps*norm(vals_fout, inf) );
check(4) = ( norm(vals_gout - vals_g, inf) < ...
    1e3*eps*norm(vals_gout, inf) );

pass(7) = all( check );

end

author: chebfun

repo: chebfun


text:
clear all;

I=double(imread('data/lena.png'))/255;
% extract 8 x 8 patches
X=im2col(I,[8 8],'sliding');
X=X-repmat(mean(X),[size(X,1) 1]);
X=X ./ repmat(sqrt(sum(X.^2)),[size(X,1) 1]);

param.K=64;  % learns a dictionary with 64 elements
param.lambda=0.05;
param.numThreads=4; % number of threads
param.batchsize=400;
param.tol = 1e-3

param.iter=200;  %

if false
param.regul = 'l1';
fprintf('with Fista Regression %s\n',param.regul);
tic
D = mexStructTrainDL(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);
%
param.regul = 'l2';
fprintf('with Fista Regression %s\n',param.regul);
tic
D = mexStructTrainDL(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);
%
param.regul = 'elastic-net';
fprintf('with Fista %s\n',param.regul);
param.lambda2=0.1;
tic
D = mexStructTrainDL(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);

%%% GRAPH
param.lambda=0.1; % regularization parameter
param.tol=1e-5;
param.K = 10
graph.eta_g=[1 1 1 1 1];
graph.groups=sparse([0 0 0 1 0;
                     0 0 0 0 0;
                     0 0 0 0 0;
                     0 0 0 0 0;
                     0 0 1 0 0]);   % g5 is included in g3, and g2 is included in g4
graph.groups_var=sparse([1 0 0 0 0;
                         1 0 0 0 0;
                         1 0 0 0 0 ;
                         1 1 0 0 0;
                         0 1 0 1 0;
                         0 1 0 1 0;
                         0 1 0 0 1;
                         0 0 0 0 1;
                         0 0 0 0 1;
                         0 0 1 0 0]); % represents direct inclusion relations

param.graph = graph

param.regul = 'graph';
fprintf('with Fista %s\n',param.regul);
tic
D = mexStructTrainDL(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);

%%
%%% TREE
%?pause;
param = rmfield(param,'graph');

end

param.lambda=0.1; % regularization parameter
param.tol=1e-5;
param.K = 10

tree.own_variables=  int32([0 0 3 5 6 6 8 9]);   % pointer to the first variable of each group
tree.N_own_variables=int32([0 3 2 1 0 2 1 1]); % number of "root" variables in each group
tree.eta_g=[1 1 1 2 2 2 2.5 2.5];
tree.groups=sparse([0 0 0 0 0 0 0 0; ...
                    1 0 0 0 0 0 0 0; ...
                    0 1 0 0 0 0 0 0; ...
                    0 1 0 0 0 0 0 0; ...
                    1 0 0 0 0 0 0 0; ...
                    0 0 0 0 1 0 0 0; ...
                    0 0 0 0 1 0 0 0; ...
                    0 0 0 0 0 0 1 0]);  % first group should always be the root of the tree

param.tree = tree;

param.regul = 'tree-l0';
fprintf('with Fista %s\n',param.regul);

tic
D = mexStructTrainDL(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);

%
param.regul = 'tree-l2';
fprintf('with Fista %s\n',param.regul);
tic
D = mexStructTrainDL(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);

%
param.regul = 'tree-linf';
fprintf('with Fista %s\n',param.regul);

tic
D = mexStructTrainDL(X,param);
t=toc;
fprintf('time of computation for Dictionary Learning: %f\n',t);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function z = gt(x,y)
%GT Greater than for sptensors.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Observations for sparse matrix case.
% The result of a > 5 is sparse.
% The result of a > 0 is sparse.
% The result of a > full(a) is sparse.
% The result of a > zeros(a) is sparse.

%% Case 1: One argument is a scalar
if isscalar(y)
    subs1 = x.subs((x.vals > y),:);
    if y < 0
        subs2 = setdiff(allsubs(x),x.subs,'rows');
    else
        subs2 = [];
    end
    z = sptensor([subs1;subs2],true,size(x));
    return;
end

% Call back with the arguments reversed.
if isscalar(x)
    z = lt(y,x);
    return;
end

%% Case 2: Both x and y are tensors of some sort
% Check that the sizes match
if ~isequal(x.size,y.size)
    error('Size mismatch');
end

% Case 2a: Two sparse tensors
if isa(x,'sptensor') && isa(y,'sptensor')
    z = lt(y,x);
    return;
end

% Case 2b: One dense tensor
if isa(y,'tensor')

    % x zero and y < 0
    subs1 = find(y < 0);
    if ~isempty(subs1)
        subs1 = setdiff(subs1,x.subs,'rows');
    end

    % x and y nonzero
    subs2 = x.subs(x.vals > y(x.subs,'extract'),:);

    % assemble
    z = sptensor([subs1;subs2],true,size(x));

    return;

end

%% Otherwise
error('The arguments must be two sptensors or an sptensor and a scalar.');

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function area = sphere_cap_area_nd ( ndim, r, h )

%*****************************************************************************80
%
%% SPHERE_CAP_AREA_ND computes the area of a spherical cap in ND.
%
%  Discussion:
%
%    The spherical cap is a portion of the surface of the sphere:
%
%      sum ( X(1:N)^2 ) = R^2
%
%    which is no more than H units from the uppermost point on the sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Thomas Ericson and Victor Zinoviev,
%    Codes on Euclidean Spheres,
%    Elsevier, 2001, pages 439-441.
%    QA166.7 E75
%
%  Parameters:
%
%    Input, integer NDIM, the dimension of the space.
%
%    Input, real R, the radius of the sphere.
%
%    Input, real H, the "thickness" of the spherical cap,
%    which is normally between 0 and 2 * R.
%
%    Output, real AREA, the area of the spherical cap.
%
  if ( h <= 0.0 )
    area = 0.0;
    return
  end

  if ( 2.0 * r <= h )
    area = sphere_area_nd ( ndim, r );
    return
  end
%
%  For cases where R < H < 2 * R, work with the complementary region.
%
  haver_sine = sqrt ( ( 2.0 * r - h ) * h );

  theta = r8_asin ( haver_sine / r );

  if ( ndim < 1 )

    area = -1.0;

  elseif ( ndim == 1 )

    area = 0.0;

  elseif ( ndim == 2 )

    area = 2.0 * theta * r;

  else

    ti = theta;

    tj = ti;
    ti = 1.0 - cos ( theta );

    for i = 2 : ndim-2
      tk = tj;
      tj = ti;
      ti = ( ( i - 1 ) * tk - cos ( theta ) * sin ( theta )^( i - 1 ) ) / i;
    end

    area = sphere_k ( ndim-1 ) * ti * r^( ndim - 1 );

  end
%
%  Adjust for cases where R < H < 2R.
%
  if ( r < h )
    area2 = sphere_area_nd ( ndim, r );
    area = area2 - area;
  end

  return
end

repo: jburkardt-m


text:
function M = symmetricfactory(n, k)
% Returns a manifold struct to optimize over k symmetric matrices of size n
%
% function M = symmetricfactory(n)
% function M = symmetricfactory(n, k)
%
% Returns M, a structure describing the Euclidean space of n-by-n symmetric
% matrices equipped with the standard Frobenius distance and associated
% trace inner product, as a manifold for Manopt.
% By default, k = 1. If k > 1, points and vectors are stored in 3D matrices
% X of size nxnxk such that each slice X(:, :, i), for i = 1:k, is
% symmetric.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Jan. 22, 2014.
% Contributors:
% Change log:

    if ~exist('k', 'var') || isempty(k)
        k = 1;
    end

    M.name = @() sprintf('(Symmetric matrices of size %d)^%d', n, k);

    M.dim = @() k*n*(n+1)/2;

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d(:), 'fro');

    M.dist = @(x, y) norm(x(:)-y(:), 'fro');

    M.typicaldist = @() sqrt(k)*n;

    M.proj = @(x, d) multisym(d);

    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @(x, eg, eh, d) M.proj(x, eh);

    M.tangent = @(x, d) d;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

	M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.rand = @() multisym(randn(n, n, k));

    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = multisym(randn(n, n, k));
        u = u / norm(u(:), 'fro');
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(n, n, k);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2);

    % Elaborate list of indices of diagonal entries of an nxnxk matrix.
    single_diag_entries = (1:(n+1):n^2)';
    all_diag_entries = bsxfun(@plus, single_diag_entries, n^2*(0:(k-1)));
    all_diag_entries = all_diag_entries(:);

    % Likewise, elaborate list of indices of upper-triangular entries.
    single_upper_triangle = find(triu(ones(n), 1));
    all_upper_triangle = bsxfun(@plus, single_upper_triangle, n^2*(0:(k-1)));
    all_upper_triangle = all_upper_triangle(:);

    % To vectorize a matrix, we extract all diagonal entries, then all
    % upper-triangular entries, the latter being scaled by sqrt(2) to
    % ensure isometry, that is: given two tangent vectors U and V at a
    % point X, M.inner(X, U, V) is equal to u'*v, where u = M.vec(X, U) and
    % likewise for v. This construction has the advantage of providing a
    % vectorized representation of matrices that has the same length as the
    % intrinsic dimension of the space they live in.
    M.vec = @(x, u_mat) [u_mat(all_diag_entries) ; ...
                         sqrt(2)*u_mat(all_upper_triangle)];
    M.mat = @matricize;
    function u_mat = matricize(X, u_vec) %#ok<INUSL>
        u_mat = zeros(n, n, k);
        u_mat(all_upper_triangle) = u_vec((k*n+1):end) / sqrt(2);
        u_mat = u_mat + multitransp(u_mat);
        u_mat(all_diag_entries) = u_vec(1:(k*n));
    end
    M.vecmatareisometries = @() true;

end

% Former, easier versions for vec / mat. They had the disadvantage of
% giving vector representations of length k*n^2, instead of k*n*(n+1).
% M.vec = @(x, u_mat) u_mat(:);
% M.mat = @(x, u_vec) reshape(u_vec, [m, n]);
% M.vecmatareisometries = @() true;

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [ res, jac ] = opt07_rj ( x, flag )

%*****************************************************************************80
%
%% OPT07_RJ evaluates RES and JAC for test case #7.
%
%  Discussion:
%
%    This example is known as the helical valley function.
%
%    The optimizing value is
%
%      X* = (1,0,0)
%
%    for which
%
%      RES(X*) = (0,0,0).
%
%  Modified:
%
%    07 January 2008
%
%  Author:
%
%    Jeff Borggaard,
%    Gene Cliff,
%    Virginia Tech.
%
%  Reference:
%
%    John Dennis, Robert Schnabel,
%    Numerical Methods for Unconstrained Optimization
%    and Nonlinear Equations,
%    SIAM, 1996,
%    ISBN13: 978-0-898713-64-0,
%    LC: QA402.5.D44.
%
%  Parameters:
%
%    Input, real X(3), the evaluation point.
%
%    Input, string FLAG, indicates what must be computed.
%    'f' means only the value of RES is needed,
%    'g' means only the value of JAC is needed,
%    'all' means RES and JAC are needed.
%    It is acceptable to behave as though FLAG was 'all'
%    on every call.
%
%    Output, real RES(3,1), the function column vector.
%
%    Output, real JAC(3,3), the Jacobian matrix.
%
  n = length ( x );

  if ( n ~= 3 )
    fprintf ( '\n' );
    fprintf ( 'OPT07_RJ - Fatal error!\n' );
    fprintf ( '  The input vector X should have length 3.\n'),
    fprintf ( '  Instead, it has length = %d.\n', n );
    keyboard
  end

  if ( 0 < x(1) )
    theta = atan ( x(2) / x(1) ) / ( 2 * pi );
  elseif ( x(1) < 0 )
    theta = atan ( x(2) / x(1) ) / ( 2 * pi ) + 0.5;
  elseif ( 0 < x(2) )
    theta = 0.25;
  elseif ( x(2) < 0 )
    theta = - 0.25;
  else
    theta = 0;
  end

  res = zeros(3,1);

  res(1,1) = 10 * ( x(3) - 10 * theta );
  res(2,1) = 10 * ( sqrt ( x(1)^2 + x(2)^2 ) - 1 );
  res(3,1) = x(3);

  if ( x(1)^2 + x(2)^2 == 0 )
    dtdx1 = 0;
    dtdx2 = 0;
  else
    dtdx1 = - x(2) / ( 2 * pi * ( x(1)^2 + x(2)^2 ) );
    dtdx2 =   x(1) / ( 2 * pi * ( x(1)^2 + x(2)^2 ) );
  end

  jac = zeros(3,3);

  jac(1,1) = -100 * dtdx1;
  jac(1,2) = -100 * dtdx2;
  jac(1,3) = 10;

  if ( x(1)^2 + x(2)^2 == 0 )
    jac(2,1) = 0;
    jac(2,2) = 0;
  else
    jac(2,1) = 10 * x(1) / sqrt ( x(1)^2 + x(2)^2 );
    jac(2,2) = 10 * x(2) / sqrt ( x(1)^2 + x(2)^2 );
  end
  jac(2,3) = 0;

  jac(3,1) = 0;
  jac(3,2) = 0;
  jac(3,3) = 1;


text:
function [val,idx] = of_RMSE(obs,sim,idx)
% of_RMSE Calculates the Root Mean Squared Error of simulated streamflow.
% Ignores time steps with negative flow values.

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% In:
% obs       - time series of observations       [nx1]
% sim       - time series of simulations        [nx1]
% idx       - optional vector of indices to use for calculation, can be
%               logical vector [nx1] or numeric vector [mx1], with m <= n
%
% Out:
% val       - objective function value          [1x1]
% idx       - indices used for the calculation

%% Check inputs and select timesteps
if nargin < 2
    error('Not enugh input arguments')
end

if nargin < 3; idx = []; end
[sim, obs, idx] = check_and_select(sim, obs, idx);

%% Calculate metric
val = sqrt(mean((obs-sim).^2));
end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [x,err,k] = perform_conjugate_gradient(A,y,options)

% perform_conjugate_gradient - perform conjugate gradient
%
%     [x,err,k] = perform_conjugate_gradient(A,y,options);
%
%   Solves for A*x=y.
%   Works both for vector x,y and matrix x,y (parralel soving of many
%   linear equations).
%
%   Works only for semi-definite matrices.
%
%   A can be a matrix or a callback function y=A(x,options).
%   In this case, you have to set options.ncols as the number of columns of
%   A (if it is a callback).
%
%   err monitos the decreasing of |A*x-y|.
%   k is the total number of iterations.
%
%   You can set:
%       options.x is an initial guess
%       options.epsilon is maximum error
%       options.niter_max is maximum number of error
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
if isfield(options, 'niter_max')
    niter = options.niter_max;
else
    niter = 100;
end

use_callback = 0;
if not(isnumeric(A))
    use_callback = 1;
end
if isfield(options, 'x')
    x = options.x;
else
    if use_callback==0
        x = zeros(size(A,2),1);
    else
        if isfield(options, 'ncols')
            x = zeros(options.ncols, 1);
        else
            error('You have to specify options.ncols');
        end
    end
end
if isfield(options, 'epsilon')
    epsilon = options.epsilon;
else
    epsilon = 1e-5;
end

if norm(y, 'fro') == 0
    normb = epsilon;
else
    normb = epsilon * sum(y(:).^2);
end

if use_callback==0
    r  = y - A*x;
else
    r  = y - feval(A,x,options);
end
r0 = sum(r.^2);
err = [sum(r0)];
for it=1:niter
    if (it==1)
        p = r;
    else
        % search direction
        beta = r0./rprev;
        p = r + repmat(beta, [size(x,1) 1]).*p;
    end
    % auxiliary vector
    if use_callback==0
        w  = A*p;
    else
        w  = feval(A,p,options);
    end
    alpha = repmat( r0 ./ sum(p .* w), [size(x,1) 1] );              % found optimal alpha in line search
    x = x + alpha.*p;                       % new guess
    r = r - alpha.*w;                       % the residual is in fact r=b-A*x
    rprev = r0;                             % save norm of the old residual
    r0 = sum(r.^2);                         % compute norm of new residual
    err(end+1) = sum(r0);
    if err(end)<normb
        return;
    end
end

return;

err = [sqrt(r0)];
while ( sqrt(r0) > normb && k < kmax)
    if( k==1 )
        p = r;
    else
        beta = r0/rprev;
        % search direction
        p = r + beta*p;
    end
    if use_callback==1
        w = feval(A,p,options);
    else
        w = A * p;                          %% auxiliary vector
    end
    alpha = r0 / (p' * w);              %% found optimal alpha in line search
    x = x + alpha*p;                    %% new guess
    r = r - alpha*w;                    %% the residual is in fact r=b-A*x
    rprev = r0;                         %% ... save norm of the old residual
    r0 = ( norm(r, 'fro') )^2;                 %% ... compute norm of new residual
    k = k + 1;
    err(end+1) = sqrt(r0);
end

author: gpeyre

repo: matlab-toolboxes