text:
function [movie_frame] = project_s3_partition(N,varargin)
%PROJECT_S3_PARTITION Use projection to illustrate an EQ partition of S^3
%
%Syntax
% [movie_frame] = project_s3_partition(N,options);
%
%Description
% PROJECT_S3_PARTITION(N) uses projection to illustrate the partition of
% the unit sphere S^3 into N regions.
%
% MOVIE_FRAME = PROJECT_S3_PARTITION(N) sets MOVIE_FRAME to be an array of
% movie frames for use with MOVIE. The movie frames will contain the region by
% region build-up of the illustration.
%
% PROJECT_S3_PARTITION(N,'offset','extra') uses experimental extra offsets.
% For more detail on partition options, see HELP PARTITION_OPTIONS.
%
% PROJECT_S3_PARTITION(N,options) also recognizes a number of illustration
% options, which are specified as name, value pairs.
% Any number of pairs can be used, in any order.
%
% The following illustration options are used.
%
% PROJECT_S3_PARTITION(N,'fontsize',size)
% Font size used in titles (numeric, default 18).
%
% PROJECT_S3_PARTITION(N,'title','long')
% PROJECT_S3_PARTITION(N,'title','short')
% Use long or short titles (default 'long').
%
% PROJECT_S3_PARTITION(N,'proj','stereo')
% PROJECT_S3_PARTITION(N,'proj','eqarea')
% Use stereographic or equal area projection (default 'stereo').
%
% PROJECT_S3_PARTITION(N,'points','show')
% PROJECT_S3_PARTITION(N,'points','hide')
% Show or hide center points (default 'show').
%
% PROJECT_S3_PARTITION(N,'surf','show')
% PROJECT_S3_PARTITION(N,'surf','hide')
% Show or hide surfaces of regions (default 'show').
%
% For more detail on illustration options, see HELP ILLUSTRATION_OPTIONS.
%
%Examples
% > project_s3_partition(10)
% > frames=project_s3_partition(9,'offset','extra','proj','eqarea')
% frames =
% 1x18 struct array with fields:
%     cdata
%     colormap
% > project_s3_partition(99,'proj','eqarea','points','hide')
%
%See also
% MOVIE, PARTITION_OPTIONS, ILLUSTRATION_OPTIONS, PROJECT_S2_PARTITION

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Function changed name from s2x to polar2cart
% Function changed name from x2s2 to cart2polar2
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-13 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

pdefault.extra_offset =  false;

popt = partition_options(pdefault, varargin{:});

gdefault.fontsize = 18;
gdefault.stereo =        true;
gdefault.long_title =    true;
gdefault.show_points =   true;
gdefault.show_surfaces = true;

gopt = illustration_options(gdefault, varargin{:});

dim = 3;

[X,Y,Z] = sphere(90);
if gopt.stereo
    r = 0;
else
    r = (area_of_sphere(dim)/volume_of_ball(dim)).^(1/dim);
end

hold off

if gopt.show_surfaces
    surf(r*X,r*Y,r*Z,zeros(size(Z)),...
        'FaceAlpha',1/20,'FaceColor','interp','FaceLighting','phong','EdgeColor','none')
else
    plot3(0,0,0,'w.')
end
axis equal;hold on
camlight right
colormap jet
grid off
axis off

if gopt.long_title
    if gopt.stereo
        s = 'Stereographic';
    else
        s = 'Equal volume';
    end
    if gopt.show_points
        pointstr = ', showing the center point of each region';
    else
        pointstr = '';
    end

    title(sprintf(...
        '\n%s projection of recursive zonal equal area partition of {S^3} \n into %d regions%s.',...
        s,N,pointstr),'FontSize',gopt.fontsize);
else
    title(sprintf('\nEQ(3,%d)',N),'FontSize',gopt.fontsize);
end

axis equal
grid off
axis off

pause(0);
if nargout > 0
    movie_frame(1) = getframe(gcf);
end

if gopt.stereo && (N == 1)
    return;
end

if popt.extra_offset
    [R,dim_1_rot] = eq_regions(dim,N,popt.extra_offset);
else
    R = eq_regions(dim,N);
end

for i = N:-1:2
    if popt.extra_offset
        project_s3_region(R(:,:,i),N,gopt.stereo,gopt.show_surfaces,dim_1_rot{i});
    else
        project_s3_region(R(:,:,i),N,gopt.stereo,gopt.show_surfaces);
    end
    pause(0);
    if nargout > 0
        movie_frame(N-i+2) = getframe(gcf);
    end
end

if gopt.show_points
    project_s3_eq_point_set(N,popt.extra_offset,gopt.stereo);
    if nargout > 0
        for k=1:min(N,40)
            movie_frame(N+k) = getframe(gcf);
        end
    end
end

hold off
%
% end function

function project_s3_region(region, N, stereo, show_surfaces, rot_matrix)
%PROJECT_S3_REGION Use projection to illustrate an EQ region of S^3
%Syntax
% project_s3_region(region, stereo, show_surfaces, rot_matrix);
%
%Notes
% The region is given as a 3 x 2 matrix in spherical polar coordinates
%
% The default is to use stereographic projection
% If the optional second argument, stereo is false,
% then use a equal area projection.

if nargin < 3
    stereo = true;
end
if stereo
    projection = 'x2stereo';
else
    projection = 'x2eqarea';
end
if nargin < 4
    show_surfaces = true;
end

offset_regions = (nargin >= 5);

tol = eps*2^5;

dim = size(region,1);
t = region(:,1);
b = region(:,2);

if abs(b(1)) < tol
    b(1) = 2*pi;
end
pseudo = 0;
if abs(t(1)) < tol && abs(b(1)-2*pi) < tol
    pseudo = 1;
end
n = 33;
delta = 1/(n-1);
h = 0:delta:1;
[h1, h2]  = meshgrid(h,h);
t_to_b = zeros(dim,n,n);
b_to_t = t_to_b;
r = N^(-1/3)/32;
for k = 1:dim
    if ~pseudo || k < 3
        L = 1:dim;
        j(L) = mod(k+L,dim)+1;
        t_to_b(j(1),:,:) = t(j(1))+(b(j(1))-t(j(1)))*h1;
        t_to_b(j(2),:,:) = t(j(2))+(b(j(2))-t(j(2)))*h2;
        t_to_b(j(3),:,:) = t(j(3))*ones(n,n);
        t_to_b_v = reshape(t_to_b,dim,n*n);
        if offset_regions
            t_to_b_x = polar2cart([cart2polar2(rot_matrix*polar2cart(t_to_b_v(1:dim-1,:)));t_to_b_v(dim,:)]);
        else
            t_to_b_x = polar2cart(t_to_b_v);
        end
        s = reshape(feval(projection,t_to_b_x),dim,n,n);
        degenerate = (norm(s(:,1,1)-s(:,1,2)) < tol);
        if ~degenerate && (~pseudo || k > 1)
            [X,Y,Z] = fatcurve(squeeze(s(:,1,:)),r);
            surface(X,Y,Z,zeros(size(Z)),...
                'FaceAlpha',(t(dim)/pi)/2,'FaceColor','interp','FaceLighting','phong','EdgeColor','none')
            axis equal; hold on
        end
        if show_surfaces
            surf(squeeze(s(1,:,:)),squeeze(s(2,:,:)),squeeze(s(3,:,:)),t(3)*ones(n,n),...
                'FaceAlpha',(t(dim)/pi)/2,'FaceColor','interp','FaceLighting','phong','EdgeColor','none')
        end
        axis equal; hold on
        camlight right
        b_to_t(j(1),:,:) = b(j(1))-(b(j(1))-t(j(1)))*h1;
        b_to_t(j(2),:,:) = b(j(2))-(b(j(2))-t(j(2)))*h2;
        b_to_t(j(3),:,:) = b(j(3))*ones(n,n);
        b_to_t_v = reshape(b_to_t,dim,n*n);
        if offset_regions
            b_to_t_x = polar2cart([cart2polar2(rot_matrix*polar2cart(b_to_t_v(1:dim-1,:)));b_to_t_v(dim,:)]);
        else
            b_to_t_x = polar2cart(b_to_t_v);
        end
        s = reshape(feval(projection,b_to_t_x),dim,n,n);
        degenerate = (norm(s(:,1,1)-s(:,1,2)) < tol);
        if ~degenerate && (~pseudo || (k > 1 && abs(b(2)-pi) > tol))
            [X,Y,Z] = fatcurve(squeeze(s(:,1,:)),r);
            surface(X,Y,Z,zeros(size(Z)),...
                'FaceAlpha',(t(dim)/pi)/2,'FaceColor','interp','FaceLighting','phong','EdgeColor','none')
        end
        if show_surfaces && k < 2
            surf(squeeze(s(1,:,:)),squeeze(s(2,:,:)),squeeze(s(3,:,:)),t(3)*ones(n,n),...
                'FaceAlpha',(t(dim)/pi)/2,'FaceColor','interp','FaceLighting','phong','EdgeColor','none')
            camlight right
        end
    end
end
colormap jet
grid off
axis off
%
% end function

function project_s3_eq_point_set(N,extra_offset,stereo)
%PROJECT_S3_EQ_POINT_SET Use projection to illustrate an EQ point set of S^3
%
%Syntax
% project_s3_eq_point_set(N,min_energy,stereo);

if nargin < 2
    extra_offset = true;
end
if nargin < 3
    stereo = true;
end
if stereo
    projection = 'stereo';
else
    projection = 'eqarea';
end

x = eq_point_set(3,N,extra_offset);
project_point_set(x,'title','hide','proj',projection);
%
% end function

repo: gspbox


text:
function [y,fs,h]=v_readcnx(filename,mode)
%V_READCNX  Read a .CNX format sound file [Y,FS,H]=(FILENAME)
%
% Inputs:
%           filename : character string containing filename (default extension .cnx)
%           mode:    't' = trim leading and trailing silences
%                    'h' = read header only
%             'd'    Look in data directory: v_voicebox('dir_data')
%
% Outputs:
%           y  : column vector containing waveform
%           fs : sample frequency
%           h  : parameter array:
%                h(1) = number of samples in file
%                h(2) = status: 0=good, 1=bad
%                h(3) = start sample number
%                h(4) = ending sample number
%                h(5) = speaker identification number
%                h(6) = speaker age group
%                h(7) = speaker sex: 0=male, 1=female
%                h(8) = ascii character
%                h(9) = repetition number
%
% This is the format of the BT Connex-S1 alphabet database
% Note: the decoding is not particularly robust and assumes
%       that all headers contain the same sequence of fields

%	   Copyright (C) Mike Brookes 1998
%      Version: $Id: v_readcnx.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ix=[17 71; 18 0; 19 0; 10 0; 12 0; 13 77; 15 -1; 16 0; ];

if nargin<2  mode='0'; end
if any(mode=='d')
    filename=fullfile(v_voicebox('dir_data'),filename);
end
fid=fopen(filename,'rb','l');
if fid == -1
    fn=[filename,'.cnx'];
    fid=fopen(fn,'rb','l');
    if fid ~= -1 filename=fn; end
end
if fid == -1 error(sprintf('Can''t open %s for input',filename)); end
[hdr,n]=fread(fid,512,'uchar');
if n ~= 512 error(sprintf('Can''t read header from connex file %s',filename)); end
del=find(hdr(5:end)=='|')';
fs=sscanf(char(hdr(17:del(1)+3)),'%f');
h=zeros(size(ix,1),1);
for i=1:length(h)
    e=find(hdr(del(ix(i)-1)+5:del(ix(i))+3)=='=');
    if ix(i,2)
        h(i)=hdr(del(ix(i)-1)+e+5);
        if ix(i,2)>0
            h(i)=1-(h(i)==ix(i,2));
        end
    else
        h(i)=sscanf(char(hdr(del(ix(i)-1)+e+5:del(ix(i))+3)),'%d');
    end
end
if any(mode =='h')
    y=[];
elseif any(mode =='t')
    fseek(fid,2*h(2),0);
    [y,n]=fread(fid,h(3)-h(2)+1,'short');
    if n ~= h(3)-h(2)+1 error(sprintf('Error reading data from connex file %s',filename)); end
else
    y=fread(fid,inf,'short');
end
fseek(fid,0,1);
h=[ftell(fid)/2-256; h];
fclose(fid);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% Script for internal menpo data split, convert the original training data
% to 2/3rds training and 1/3rd validation

menpo_root = 'D:\Datasets\menpo/';

pts_files = dir([menpo_root, '/*.jpg']);
jpg_files = dir([menpo_root, '/*.pts']);

% Do the actual copying to respective folders
out_train = [menpo_root, '/train/'];
out_valid = [menpo_root, '/valid/'];
mkdir(out_train);
mkdir(out_valid);

load('train_valid_split');

for i=1:numel(train_imgs)

    copyfile([menpo_root, train_pts(i).name], [out_train, train_pts(i).name]);
    copyfile([menpo_root, train_imgs(i).name], [out_train, train_imgs(i).name]);

end

for i=1:numel(valid_imgs)

    copyfile([menpo_root, valid_pts(i).name], [out_valid, valid_pts(i).name]);
    copyfile([menpo_root, valid_imgs(i).name], [out_valid, valid_imgs(i).name]);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [interp_cmap]=interpolate_cbrewer(cbrew_init, interp_method, ncolors)
%
% INTERPOLATE_CBREWER - interpolate a colorbrewer map to ncolors levels
%
% INPUT:
%   - cbrew_init: the initial colormap with format N*3
%   - interp_method: interpolation method, which can be the following:
%               'nearest' - nearest neighbor interpolation
%               'linear'  - bilinear interpolation
%               'spline'  - spline interpolation
%               'cubic'   - bicubic interpolation as long as the data is
%                           uniformly spaced, otherwise the same as 'spline'
%   - ncolors=desired number of colors
%
% Author: Charles Robert
% email: tannoudji@hotmail.com
% Date: 14.10.2011

% just to make sure, in case someone puts in a decimal
ncolors=round(ncolors);

% How many data points of the colormap available
nmax=size(cbrew_init,1);

% create the associated X axis (using round to get rid of decimals)
a=(ncolors-1)./(nmax-1);
X=round([0 a:a:(ncolors-1)]);
X2=0:ncolors-1;

z=interp1(X,cbrew_init(:,1),X2,interp_method);
z2=interp1(X,cbrew_init(:,2),X2,interp_method);
z3=interp1(X,cbrew_init(:,3),X2, interp_method);
interp_cmap=round([z' z2' z3']);

end
% Copyright (c) 2011, Charles Robert
% 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.
%

author: alecjacobson

repo: gptoolbox


text:
function test_suite = test_grRemoveNodes
%TEST_GRREMOVENODES  Test case for the file grRemoveNodes
%
%   Test case for the file grRemoveNodes

%   Example
%   test_grRemoveNodes
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2014-03-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2014 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_RemoveOneVertex(testCase) %#ok<*DEFNU>
% Test call of function without argument
[nodes, edges] = makeSimpleGraph();
[nodes2, edges2] = grRemoveNodes(nodes, edges, 3);
testCase.assertEqual(5, size(nodes2, 1));
testCase.assertEqual(7, size(edges2, 1));

function test_RemoveTwoMiddleVertices(testCase)
% Test call of function without argument
[nodes, edges] = makeSimpleGraph();
[nodes2, edges2] = grRemoveNodes(nodes, edges, [3 4]);
testCase.assertEqual(4, size(nodes2, 1));
testCase.assertEqual(5, size(edges2, 1));

function test_RemoveTwoExtremeVertices(testCase)
% Test call of function without argument
[nodes, edges] = makeSimpleGraph();
[nodes2, edges2] = grRemoveNodes(nodes, edges, [1 6]);
testCase.assertEqual(4, size(nodes2, 1));
testCase.assertEqual(3, size(edges2, 1));

function [nodes, edges] = makeSimpleGraph()
nodes = [...
    10 10; 20 10; 30 10; ...
    10 20; 20 20; 30 20];
edges = [...
    1 2; 1 4; 1 5; ...
    2 3; 2 5; 2 6; ...
    3 6; 4 5; 5 6];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [p1,p2,g1,g2] = doublePendulumPosition(z,P)
%[p1,p2,g1,g2] = DOUBLEPENDULUMPOSITION(Z,P)
%
%FUNCTION:  This function computes the energy of a double
%    pendulum.
%INPUTS:
%    z = [4xn] matrix of states.
%    P = struct of parameters
%OUTPUTS:
%    p1 = [2xn] position of the end of the first link
%    p2 = [2xn] position of the end of the second link
%    g1 = [2xn] position of the CoM of link one
%    g2 = [2xn] position of the CoM of link two
%
%NOTES:
%    This file was automatically generated by writeDoublePendulumPosition

m1 = P.m1; %link one mass
m2 = P.m2; %link two mass
g  = P.g ; %gravity
l1 = P.l1; %link one length
l2 = P.l2; %link two length
I1 = P.I1; %link one moment of inertia about its center of mass
I2 = P.I2; %link two moment of inertia about its center of mass
d1 = P.d1; %distance between link one center of mass and parent joint
d2 = P.d2; %distance between link two center of mass and parent joint

th1 = z(1,:); %link one absolute angle
dth1 = z(2,:); %link one angular rate
th2 = z(3,:); %link two absolute angle
dth2 = z(4,:); %link two angular rate

n = length(th1);
p1 = zeros(2,n);
p2 = zeros(2,n);
g1 = zeros(2,n);
g2 = zeros(2,n);
p1(1,:) = l1.*cos(th1);
p1(2,:) = l1.*sin(th1);
p2(1,:) = l1.*cos(th1) + l2.*cos(th2);
p2(2,:) = l1.*sin(th1) + l2.*sin(th2);
g1(1,:) = d1.*cos(th1);
g1(2,:) = d1.*sin(th1);
g2(1,:) = d2.*cos(th2) + l1.*cos(th1);
g2(2,:) = d2.*sin(th2) + l1.*sin(th1);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function[cv,ca]=latlon2uv(varargin)
%LATLON2UV  Converts latitude and longitude to horizontal velocity.
%
%   CV=LATLON2UV(NUM,LAT,LON) where NUM is the data in DATENUM format
%   and LAT and LON are the latitude and longitude in degrees, outputs the
%   complex-valued velocity CV=U+SQRT(-1)*V in cm/s, computed using the
%   first central difference.
%
%   NUM is a column vector or a matrix of the same size as LAT and LON.
%   LAT and LON are matices having SIZE(NUM,1) rows.  CV is the same size
%   as LAT and LON. NANs in LAT or LON become NAN+SQRT(-1)*NAN.
%
%   LATLON2UV computes the velocity components from the distance travelled
%   across the surface of the sphere and the heading, taking account of
%   the sphericity of the earth.
%
%   The first and last points must be treated differently, as the central
%   difference is not defined there.  At the first and last point we use,
%   respectively, the first forward and first backward difference.
%
%   The radius of the earth is given by the function RADEARTH.
%
%   [CV,CA]=LATLON2UV(...,'acceleration') also returns CA, a measure of the
%   complex-valued acceleration based on the difference between the first
%   forward difference velocity and the first backward difference velocity.
%   The units of acceleration are centimeters per inverse second squared.
%   ___________________________________________________________________
%
%   Cell array input / output
%
%   LATLON2UV returns cell array output given cell array input.
%
%   That is, if NUM, LAT, and LON, are all cell arrays of length K,
%   containing K different numerical arrays, then the output will also be
%   cell arrays of length K.
%   ___________________________________________________________________
%
%   Difference algorithm
%
%   By default, LATLON2UV uses a first central difference algorithm which
%   is defined to be the average of a forward and a backward difference.
%
%   CV=LATLON2UV(...,'forward') instead uses a first forward difference.
%   Velocity is computed based on the great circle distance and bearing
%   from each point to the next point.
%
%   CV=LATLON2UV(...,'backward') instead uses a first backward difference.
%   Velocity is computed based on the great circle distance and bearing
%   from each point to the previous point.  This is equivalent to the first
%   forward difference computed backward in time.
%
%   The first forward and first backward difference can be thought of as
%   giving the Cartesian velocities of a departing and arriving particle,
%   respectively.
%
%   Note that acceleration is not output when the 'forward' or 'backward'
%   algorithms are specified, as acceleration is computed as a part of the
%   first central difference operation.
%   ___________________________________________________________________
%
%   Parallelization
%
%   LATLON2UV(...,'parallel') with cell array input parallelizes the
%   computation using a PARFOR loop.  This requires that Matlab's Parallel
%   Computing Toolbox be installed, and is useful for very large datasets.
%   __________________________________________________________________
%
%   See also XY2LATLON, LATLON2XY, UV2LATLON.
%
%   Usage:  cv=latlon2uv(num,lat,lon);
%           [cv,ca]=latlon2uv(num,lat,lon);
%
%   'latlon2uv --t' runs a test.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin, '--t')
    latlon2uv_test,return
end

str='centered';
cores='serial';

for i=1:2
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'par')||strcmpi(varargin{end}(1:3),'ser')
            cores=varargin{end};
        else
            str=varargin{end};
        end
    end
end

num=varargin{1};
lat=varargin{2};
lon=varargin{3};
varargin=varargin(4:end);

if strcmpi(cores(1:3),'par')
    if exist('parpool')~=2
        disp('Sorry, parallel algorithm requires the Parallel Computing Toolbox.')
        disp('Defaulting to the standard algorithm.')
        cores='serial';
    end
end

if nargout==2
    str='acc';  %Output acceleration
end

if ~iscell(num)
    [cv,ca]=latlon2uv_celloop(num,lat,lon,str);
else
    if strcmpi(cores(1:3),'ser')
        for i=1:length(num)
            [cv{i,1},ca{i,1}]=latlon2uv_celloop(num{i},lat{i},lon{i},str);
        end
    elseif strcmpi(cores(1:3),'par')
        parfor i=1:length(num)
            [cv{i,1},ca{i,1}]=latlon2uv_celloop(num{i},lat{i},lon{i},str);
        end
    end
end

if ~strcmpi(str(1:3),'acc')
    clear ca
end

function[cv,ca]=latlon2uv_celloop(num,lat,lon,str)

ca=[];

if length(lat)<=2
    cv=nan+0*lat;
    ca=nan+0*lat;
else
    if size(num,1)==1
        if size(num')==size(lat)
            num=num';
        end
    end
    if size(num,2)==1
        num=vrep(num,size(lat,2),2);
    end
    if strcmpi(str(1:3),'cen')||strcmpi(str(1:3),'acc')
        [u1,v1]=latlon2uv_celloop_one(num,lat,lon,'forward');
        [u2,v2]=latlon2uv_celloop_one(num,lat,lon,'backward');

        %Correction for first and last points
        u2(1,:)=u1(1,:);
        v2(1,:)=v1(1,:);
        u1(end,:)=u2(end,:);
        v1(end,:)=v2(end,:);

        cv=frac(1,2)*(u1+u2)+1i*frac(1,2)*(v1+v2);

        if strfind(str,'acc')
            dt=vdiff(num*24*3600,1,1);
            ca=((u1-u2)+1i*(v1-v2))./dt;
            %Correction for first and last points
            ca([1 end],:)=ca([2 end-1],:);
        end
    else
        [u,v]=latlon2uv_celloop_one(num,lat,lon,str);
        cv=u+1i*v;
    end
end

function[u,v]=latlon2uv_celloop_one(num,lat,lon,str)

if strfind(str,'old')
    [dr,dt,gamma]=latlon2uv_old(num,lat,lon,str);
elseif strfind(str,'for')
    [dr,dt,gamma]=latlon2uv_forward(num,lat,lon,str);
elseif strfind(str,'bac')
    [dr,dt,gamma]=latlon2uv_forward(flip(num,1),flip(lat,1),flip(lon,1),str);
    dr=flip(dr,1);
    dt=flip(dt,1);
    gamma=flip(gamma,1);
end

%Convert to centimeters
c=100*1000;
u=c.*dr./dt.*cos(gamma);
v=c.*dr./dt.*sin(gamma);
index=find(isnan(u)|isnan(v));
if ~isempty(index)
    u(index)=nan;
    v(index)=nan;
end

function[dr,dt,gamma]=latlon2uv_old(num,lat,lon,str)
%Former version of a first central difference algorithm

%[phi,th]=jdeg2rad(lat,lon);
[x,y,z]=latlon2xyz(lat,lon);

dt=vdiff(num*24*3600,1);
[dx,dy,dz]=vdiff(x,y,z,1);
xm=vshift(x,1,1)/2+vshift(x,-1,1)/2;
ym=vshift(y,1,1)/2+vshift(y,-1,1)/2;
zm=vshift(z,1,1)/2+vshift(z,-1,1)/2;

%Correction for first and last point
xm([1 end],:)=x([1 end],:);
ym([1 end],:)=y([1 end],:);
zm([1 end],:)=z([1 end],:);

[latm,lonm]=xyz2latlon(xm,ym,zm);
[phim,thm]=jdeg2rad(latm,lonm);

%Now perform a local rotation
dx2=dx.*cos(-thm)-dy.*sin(-thm);
dy2=dx.*sin(-thm)+dy.*cos(-thm);
dz2=dz;

%dx3=dx2.*cos(-phim)-dz2.*sin(-phim);    %Correct but not needed
dy3=dy2;
dz3=dx2.*sin(-phim)+dz2.*cos(-phim);

[lat1,lon1]=vshift(lat,lon,1,1);
[lat2,lon2]=vshift(lat,lon,-1,1);
dr=spheredist(lat1,lon1,lat2,lon2)/2;

%Treat first and last points differently
dr(1,:)=spheredist(lat(1,:),lon(1,:),lat(2,:),lon(2,:));
dr(end,:)=spheredist(lat(end-1,:),lon(end-1,:),lat(end,:),lon(end,:));

gamma=imlog(dy3+sqrt(-1)*dz3);

function[dr,dt,gamma]=latlon2uv_forward(num,lat,lon,str)

%lon=degunwrap(lon);
%First forward difference
dt=diff(num*24*3600,1,1);

[lat2,lon2]=vshift(lat,lon,1,1);
dr=spheredist(lat,lon,lat2,lon2); %figure,plot(dr)

xx=sind(lon2-lon).*cosd(lat2);
yy=cosd(lat).*sind(lat2)-sind(lat).*cosd(lat2).*cosd(lon2-lon);

gamma=atan2(yy,xx);
%figure,uvplot(diff(rot(gamma)))
%gamma2=imag(log(xx+1i*yy));
%aresame(gamma,gamma2,1e-10)

%From http://www.movable-type.co.uk/scripts/latlong.html
%ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1),SIN(lon2-lon1)*COS(lat2))
%But note, that source actually reverses the inputs to the inverse tangent function

dt(end+1,:)=dt(end,:);
dr(end,:)=dr(end-1,:);
gamma(end,:)=gamma(end-1,:);

function[]=latlon2uv_test
latlon2uv_dlat_test
latlon2uv_dlon_test
latlon2uv_small_test
latlon2uv_displacement_test
latlon2uv_npg2006_test
latlon2uv_npg2006_nonuniform_test
latlon2uv_testformer

function[]=latlon2uv_dlat_test

N=100;
tol=1e-3;

lono=2*pi*rand(1,N)-pi;
lato=pi*rand(1,N)-pi/2;
[lato,lono]=jrad2deg(lato,lono);

dlat=(rand(1,N)-1/2)/10;

lat=[lato-dlat;lato;lato+dlat];
lon=[lono;lono;lono];

num=[-1+0*lato;0*lato;1+0*lato];
index=find(max(lat)<90&min(lat)>-90);
vindex(num,lat,lon,dlat,index,2);

cv1=0*dlat+1i*100*1000*2*pi*radearth/360.*dlat/(3600*24);
cv2=latlon2uv(num,lat,lon);
b=aresame(cv2(2,:),cv1,tol);
reporttest('LATLON2UV small delta latitude',b);

function[]=latlon2uv_dlon_test

N=100;
tol=1e-3;

lono=2*pi*rand(1,N)-pi;
lato=pi*rand(1,N)-pi/2;
[lato,lono]=jrad2deg(lato,lono);

dlon=(rand(1,N)-1/2)/10;

lat=[lato;lato;lato];
lon=[lono-dlon;lono;lono+dlon];

num=[-1+0*lato;0*lato;1+0*lato];

cv1=100*1000*2*pi*radearth/360.*dlon/(3600*24).*cos(jdeg2rad(lato))+1i*0*lato;
cv2=latlon2uv(num,lat,lon);
b=aresame(cv2(2,:),cv1,tol);
reporttest('LATLON2UV small delta longitude',b);

function[]=latlon2uv_small_test

N=100;
tol=5e-3;

lono=2*pi*rand(1,N)-pi;
lato=pi*rand(1,N)-pi/2;
[lato,lono]=jrad2deg(lato,lono);

dlon=(rand(1,N)-1/2)/10;
dlat=(rand(1,N)-1/2)/10;

lat=[lato-dlat;lato;lato+dlat];
lon=[lono-dlon;lono;lono+dlon];

num=[-1+0*lato;0*lato;1+0*lato];
index=find(max(lat)<90&min(lat)>-90);
vindex(num,lat,lon,dlat,dlon,lato,index,2);

u1=100*1000*2*pi*radearth/360.*dlon/(3600*24).*cos(jdeg2rad(lato));
v1=100*1000*2*pi*radearth/360.*dlat/(3600*24);
cv1=u1+1i*v1;
cv=latlon2uv(num,lat,lon);

b=aresame(cv(2,:),cv1,tol);
reporttest('LATLON2UV small displacements',b);

function[]=latlon2uv_displacement_test

N=100;
tol=2e-3;

lon=2*pi*rand(N,1)-pi;
lat=pi*rand(N,1)-pi/2;

num=(1:N)';
cv=latlon2uv(num,lat,lon);

dr=abs(cv).*(3600.*24)./(100*1000);

[lat1,lon1]=vshift(lat,lon,1,1);
[lat2,lon2]=vshift(lat,lon,-1,1);

dr1=spheredist(lat1,lon1,lat2,lon2)/2;

b=max(frac(abs(dr(2:end-1)-dr1(2:end-1)),dr(2:end-1)))<tol;
reporttest('LATLON2UV displacement matches SPHEREDIST',b);

function[]=latlon2uv_npg2006_test
load npg2006
use npg2006
cv1=latlon2uv(num,lat,lon);
cv2=vdiff(cx,1).*100.*1000./(3600.*24.*dt);
tol=1/2;
b=aresame(cv1(2:end-1,:),cv2(2:end-1,:),tol);
reporttest('LATLON2UV matches VDIFF of Cartesian position for NPG2006 data',b);

function[]=latlon2uv_npg2006_nonuniform_test
load npg2006
use npg2006
cv1=latlon2uv([num 2*num],[lat lat],[lon lon]);
cv2=vdiff(cx,1).*100.*1000./(3600.*24.*dt);

tol=1/2;
b=aresame(cv1(2:end-1,1),cv2(2:end-1,1),tol)&&aresame(cv1(2:end-1,2),cv2(2:end-1,1)./2,tol);
reporttest('LATLON2UV matches NPG2006 data for DT different by columns, central diff',b);

cv1=latlon2uv([num 2*num],[lat lat],[lon lon],'forward');
cv2=diff(cx,1,1).*100.*1000./(3600.*24.*dt);

tol=1/2;
b=aresame(cv1(1:end-1,1),cv2(:,1),tol)&&aresame(cv1(1:end-1,2),cv2(:,1)./2,tol);
reporttest('LATLON2UV matches NPG2006 data for DT different by columns, forward diff',b);

cv1=latlon2uv(num,lat,lon,'forward');
cv2=latlon2uv(num,lat,lon,'backward');

tol=1/2;
b=aresame(cv1(1:end-1,:),cv2(2:end,:),tol);
reporttest('LATLON2UV forward and shifted backward difference agree approximately',b);

cv3=latlon2uv(num,lat,lon,'central');

b=aresame(frac(cv1(2:end-1,:)+cv2(2:end-1,:),2),cv3(2:end-1,:));
reporttest('LATLON2UV central difference is exactly average of forward and backward',b);

cv4=latlon2uv(num,lat,lon,'old');
b=aresame(cv3(2:end-1,:),cv4(2:end-1,:),0.02);
reporttest('LATLON2UV central difference matches old algorithm for NPG2006 data',b);

% cx=10*rot([0:0.001:4*pi]');
% [lat,lon]=xy2latlon(cx,65,0);
% num=[1:length(cx)]'/24/100;
% cv=latlon2uv(num,lat,lon);
% ca=latlon2uv(num,lat,lon,'acc');
% ca2=vdiff(cv,1)./(1/24/100*24*3600);
%
% uvplot(num,ca), hold on, uvplot(num,ca2)
%
%But there is a problem with eastward acceleration for very small amplitude
%circles.  Keep an eye out for noise coming from latlon2uv.

function[]=latlon2uv_testformer

lon1=[1:360]';
lat1=60 +0*lon1;

lon2=180+0*lon1;
lat2=0+0*lon1;

lon3=[360:-1:1]';
lat3=-60 +0*lon1;

lon=[lon1';lon2';lon3'];lon=lon(:);
lat=[lat1';lat2';lat3'];lat=lat(:);

num=[1:length(lat)]';

clear cv
cv(:,1)=latlon2uv(num,lat,lon,'forward');
cv(:,2)=latlon2uv(num,lat,lon,'backward');
cv(:,3)=latlon2uv(num,lat,lon,'central');
cv(:,4)=latlon2uv(num,lat,lon,'old');

cv=cv(2:3:end,:);
lat=lat(1:3:end,:);
lon=lon(1:3:end,:);

ii=[91:269]';
b=aresame(1+0*ii,cv(ii,4)./cv(ii,3),0.0001);

reporttest('LATLON2UV former algorithm and new algorthm agree to 1e-4 when Delta Lon<90',b);

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% Check correspondence between images in items_classification.mat
%   and data_evaluation.mat

addpath(genpath('..'));
fprintf(1,'Loading dataset...\n');
load('data_evaluation','drawings','images','names','timing');
load('items_classification','nrun','ntrain','ntest','cell_train','cell_test','cell_Y');
D = Dataset(drawings,images,names,timing);

fprintf(1,'Images are loaded...\n');
for r=1:nrun
     trainset = cell_train{r};
     for i=1:ntrain
         item = trainset{i};
         assert(isequal(item.img,D.get('image',item.a,item.c,item.r)));
     end

     testset = cell_test{r};
     for i=1:ntest
         item = testset{i};
         assert(isequal(item.img,D.get('image',item.a,item.c,item.r)));
     end

end
fprintf(1,'Images are verified as the same as in evaluation set...\n');

author: brendenlake

repo: omniglot


text:
function yp = p12_fun ( neqn, t, y )

%*****************************************************************************80
%
%% P12_FUN evaluates the function for problem P12.
%
%  Discussion:
%
%    10 equations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Wayne Enright, John Pryce,
%    Algorithm 648,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 1, pages 28-34.
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the derivative
%    function.
%
%    Output, real YP(NEQN), the value of the derivative function.
%
  yp = zeros ( neqn, 1 );

  yp(1) = - y(1);
  for i = 2 : neqn-1
    yp(i) = ( i - 1 ) * y(i-1) - i * y(i);
  end
  yp(neqn) = ( neqn - 1 ) * y(neqn-1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = std( varargin )

%STD    Internal cvx version.
%This actually implements VAR and STD, controlled by the square_it flag.

persistent P
if isempty( remap ),
    P.map      = cvx_remap( 'affine' );
    P.funcs    = { @std_1 };
    P.zero     = NaN;
    P.reduce   = true;
    P.reverse  = false;
    P.constant = [];
    P.dimarg   = [];
end
[ sx, x, w, dim, square_it ] = cvx_get_dimension( varargin, 3 );
if ~isempty( w ),
    w = false;
elseif ~( isnumeric(w) && numel(w) ~= length(w) && isreal(w) ),
    cvx_throw( 'Second argument must be a logical scalar or a nonnegative vector.' );
elseif numel(w) == 1,
    w = w ~= 0;
elseif numel(w) ~= sx(dim),
    cvx_throw( 'Weight vector expected to have %d elements.', nx );
elseif any( w < 0 ),
    cvx_throw( 'Weight vector must be nonnegative.' );
elseif ~any( w ),
    cvx_throw( 'Weight vector must not be all zeros.' );
else
    w = w(:) / sum( w );
end
if isempty( square_it ),
    P.name = 'var';
else
    P.name = 'std';
end
y = cvx_reduce_op( P, x, dim, w, square_it );

function y = std_1( x, w, square_it ) %#ok
[nx,nv] = size( x );
if nx == 1,
    y = cvx( [nx,nv], [] );
else
    cvx_begin
        variable xbar( 1, nv )
        nw = numel(w);
        if nw <= 1,
            nx * xbar == sum( x, 1 ); %#ok
            xmid = bsxfun( @minus, x, xbar ) / sqrt( nx - ~(nw&&w) );
        else
            xbar == sum( w' * x ); %#ok
            xmid = bsxfun( @times, sqrt( w ), bsxfun( @minus, x, xbar ) );
        end
        if square_it
            minimize( sum_square( xmid, 1 ) );
        else
            minimize( norms( xmid, 2, 1 ) );
        end
    cvx_end
    y = cvx_optval;
end

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

author: yu-jiang

repo: radpbook

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


text:
function [fx] = f_vgo(x,P,u,in)
fx = x + VBA_sigmoid(P(1))*(P(2)*u-x);

author: MBB-team

repo: VBA-toolbox


text:
function res=tesla(A, y, lambda1, lambda2, opts)
%
%%
% Function mtLogisticR:
%      Logistic Loss for ...... (to be added)
%
%% Problem
%
%   to be added
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size m x 1)
%  lambda1 -
%  lambda2 -
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x k
%  c-         The obtained intercept if size 1 x k
%  funVal-    Function value during iterations
%
%% Related papers
%
% [1]  Jun Liu, Lei Yuan, and Jieping Ye, An Efficient Algorithm for
%      a Class of Fused Lasso Problems, KDD, 2010.
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified 3 November 2009.
%
% Related functions:
%  ....
%%

%% Verify and initialize the parameters
%%
if (nargin <5)
    error('\n Inputs: A, y, lambda1,  lambda2, and opts.ind should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (lambda1<0 || lambda2 <0)
    error('\n lambda1 should be nonnegative!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%%

% Initialize ind and q
if ~isfield(opts,'ind')
    error('\n In tesla, .ind should be specified');
else
    ind=opts.ind;
    k=length(ind)-1;

    if ind(k+1)~=m
        error('\n Check opts.ind');
    end
end

%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Starting point initialization

p_flag=(y==1);                  % the indices of the postive samples

for i=1:k
    ind_i=(ind(i)+1):ind(i+1);  % indices for the i-th group

    m1(1,i)=sum(p_flag(ind_i));      % the total number of the positive samples
    m2(1,i)=length(ind_i)-m1(1,i);   % the total number of the negative samples
end

% process the regularization parameter
if (opts.rFlag==1) % lambda1 here is the scaling factor lying in [0,1]
    if (lambda1<0 || lambda1>1)
        error('\n opts.rFlag=1, and lambda1 should be in [0,1]');
    end

    % we compute ATb for computing lambda_max, when the input lambda1 is a ratio

    p_flag=(y==1);                  % the indices of the postive samples

    for i=1:k
        ind_i=(ind(i)+1):ind(i+1);  % indices for the i-th group

        b(ind_i,1)=p_flag(ind_i)*m1(i)/(m1(i)+m2(i))-...
            (~p_flag(ind_i))*m2(i)/(m1(i)+m2(i));
    end

    ATb=zeros(n, k);
    % compute AT b
    for i=1:k
        ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group

        if (opts.nFlag==0)
            tt =A(ind_i,:)'*b(ind_i,1);
        elseif (opts.nFlag==1)
            tt= A(ind_i,:)'*b(ind_i,1) - sum(b(ind_i,1)) * mu';
            tt=tt./nu(ind_i,1);
        else
            invNu=b(ind_i,1)./nu(ind_i,1);
            tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
        end

        ATb(:,i)= tt;
    end

    q_bar=Inf;

    lambda_max=0;
    for i=1:n
        lambda_max=max(lambda_max,...
            norm(  ATb(i,:), q_bar) );
    end

    lambda_max=lambda_max / m;
    lambda1=lambda1*lambda_max;

    if isfield(opts,'rFlag2')
        if (opts.rFlag2==1)
            lambda2=lambda1*lambda_max;
        end
    end
end

% initialize a starting point
if opts.init==2
    x=zeros(n,k); c=zeros(1,k);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if ( size(x,1)~=n && size(x,2)~=k )
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,k);
    end

    if isfield(opts,'c0')
        c=opts.c0;

        if ( length(c)~=k )
            error('\n Check the input .c0');
        end
    else
        c=log(m1./m2);
    end
end

Ax=zeros(m,1); % m x 1
% compute Ax: Ax_i= A_i * x_i
for i=1:k
    ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
    m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

    if (opts.nFlag==0)
        Ax(ind_i,1)=A(ind_i,:)* x(:,i);
    elseif (opts.nFlag==1)
        invNu=x(:,i)./nu; mu_invNu=mu * invNu;
        Ax(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
    else
        Ax(ind_i,1)=A(ind_i,:)*x(:,i)-repmat(mu*x(:,i), m, 1);
        Ax(ind_i,1)=Ax./nu(ind_i,1);
    end
end

%% The main program
% The Armijo Goldstein line search schemes + accelearted gradient descent

z0=zeros(n-1,k);

if (opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m; % the intial guess of the Lipschitz continuous gradient

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);
    cp=c;         ccp=zeros(1,k);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;   sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        vec_sc=zeros(m,1);
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
            vec_sc(ind_i,1)=sc(i);
        end
        aa=- y.*(As+ vec_sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y) * (1 - prob)
        b= -y.*(1-prob) / m;

        gc=zeros(1,k);                     % the gradient of c
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
            gc(1,k)=sum(b(ind_i));
        end

        % compute g= AT b, the gradient of x
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group

            if (opts.nFlag==0)
                tt =A(ind_i,:)'*b(ind_i,1);
            elseif (opts.nFlag==1)
                tt= A(ind_i,:)'*b(ind_i,1) - sum(b(ind_i,1)) * mu';
                tt=tt./nu(ind_i,1);
            else
                invNu=b(ind_i,1)./nu(ind_i,1);
                tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
            end

            g(:,i)= tt;
        end

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            v=s-g/L; c= sc- gc/L;

            [x, z, gap]=tesla_proj(v, z0, ...
                lambda1/L, lambda2/L, n, k,...
                1000, 1e-8, 1, 6);
            z0=z;

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Ax: Ax_i= A_i * x_i
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

                if (opts.nFlag==0)
                    Ax(ind_i,1)=A(ind_i,:)* x(:,i);
                elseif (opts.nFlag==1)
                    invNu=x(:,i)./nu; mu_invNu=mu * invNu;
                    Ax(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
                else
                    Ax(ind_i,1)=A(ind_i,:)*x(:,i)-repmat(mu*x(:,i), m, 1);
                    Ax(ind_i,1)=Ax./nu(ind_i,1);
                end
            end

            % aa= - diag(y) * (A * x + c)
            vec_sc=zeros(m,1);
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                vec_sc(ind_i,1)=c(i);
            end
            aa=- y.*(Ax+ vec_sc);

            % fun_s is the logistic loss at the search point
            bb=max(aa,0);
            fun_x= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

            r_sum=(norm(v,'fro')^2 + norm(c-sc,2)^2) / 2;
            l_sum=fun_x - fun_s - sum(sum(v.* g)) - (c-sc)* gc';

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + <v, g> + <c ,gc>
            %                           + L/2 * (<v,v> + <c-sc,c-sc> )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;

        xxp=x-xp;   ccp=c-cp;

        funVal(iterStep)=fun_x;

        funVal(iterStep)=funVal(iterStep)+ lambda1 * sum(abs(x(:))) + lambda2 * sum( sum( abs(x(:,1:(k-1)) - x(:,2:k) ) ) );

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
        res.gap(:,iterStep)=gap;
    end
    res.x=x;
    res.c=c;
    res.funVal=funVal;
    res.ValueL=ValueL;
    res.lambda_max=lambda_max;
end

author: thomaskuestner


text:
function [Y, X, Z, n, m, p, T, k1, k3, q1, q2, q3,X1,Y1,data_exo]=TVEmaprelim(data_endo,data_exo,const,lags,regimeperiods,names)

% function [Y X Z n m p T k1 k3 q1 q2 q3]=maprelim(data_endo,data_exo,const,lags)
% creates the preliminary values necessary for all the subsequent computations
% inputs:  - matrix 'data_endo': the matrix storing the endogenous time series data used to estimate the model
%          - matrix 'data_exo': the matrix storing the exogenous time series data used to estimate the model
%          - integer 'const': 0-1 value determining whether a constant term should be included in the model
%          - integer 'lags': the number of lags to include in the model
% outputs: - matrix 'Y': the matrix of endogenous variables, defined in (3.5.10)
%          - matrix 'X': the matrix of endogenous regressors, defined in (3.5.10)
%          - matrix 'Z': the matrix of exogenous regressors, defined in (3.5.10)
%          - integer 'n': the number of endogenous variables in the model
%          - integer 'm': the number of exogenous variables in the model
%          - integer 'p': the number of lags in the model
%          - integer 'T': the sample size, i.e. the number of time periods used to estimate the model
%          - integer 'k1': the number of coefficients related to the endogenous variables for each equation in the model
%          - integer 'k3': the number of coefficients related to the exogenous variables for each equation, in the reformulated model (3.5.5)
%          - integer 'q1': the total number of VAR coefficients related to the endogenous variables
%          - integer 'q2': the total number of VAR coefficients related to the exogenous variables
%          - integer 'q3': the total number of VAR coefficients related to the exogenous variables, in the reformulated model (3.5.5)

% first compute p, the number of lags in the model, defined p77
p=lags;

% then compute n, the number of endogenous variables in the model; it is simply the number of columns in the matrix 'data_endo'
n=size(data_endo,2);

% augment the matrix of exogenous with a column of ones to account for the constant
%data_exo=[ones(size(data_endo,1),1) data_exo];
%alternative data exo

if isempty(regimeperiods)
    data_exo=[ones(size(data_endo,1),1)];
else
data_exo=[ones(size(data_endo,1),1) zeros(size(data_endo,1),1)];
data_exo(find(strcmp(names(2:end,1),regimeperiods(1))):find(strcmp(names(2:end,1),regimeperiods(2))),2)=1;
data_exo(find(strcmp(names(2:end,1),regimeperiods(1))):find(strcmp(names(2:end,1),regimeperiods(2))),1)=0;
end

% then compute m, the number of exogenous variables in the model, defined p77
% if data_exo is empty, set m=0
m=size(data_exo,2);

%%%%% Also, trim a number initial rows equal to the number of lags, as they will be suppressed from the endogenous as well to create initial conditions
data_exo1=data_exo(p+1:end,:);

% estimate k1, the number of parameters related to endogenous variables in each equation, defined p77
k1=n*p;

% estimate q1, the total number of parameters related to endogenous variables, defined p77
q1=n*k1;

% estimate q2, the total number of parameters related to exogenous variables, defined p77
q2=n*m;

% estimate k3, the number of parameters related to exogenous variables in each equation of the modified system, defined p78
k3=m*(p+1);

% estimate q3, the total number of parameters related to endogenous variables in the modified system, defined p78
q3=n*k3;

% obtain the matrices Y and X, defined in (3.6.10)
% to do so, use the lagx function on the data matrix
temp=bear.lagx(data_endo,lags);

% to build X, take off the n initial columns of current data
X=temp(:,n+1:end);

%%%%% to build X, take off the n initial columns of current data, and concatenate the exogenous
X1=[temp(:,n+1:end) data_exo1];

% save the n first columns of temp as Y
Y=temp(:,1:n);

%%%%% save the n first columns of temp as Y
Y1=temp(:,1:n);

% obtain the matrix Z, defined in (3.6.10)
temp=bear.lagx(data_exo,lags);
temp(:,m+1:end)=-temp(:,m+1:end);
Z=temp;

% Define T, the number of periods of the model, as the number of rows of Y
T=size(Y,1);

data_exo=Z;

% obtain eventually y, Xbar and Zbar, as defined in (XXX)
% % % y=Y(:);
% % % Xbar=kron(eye(n),X);
% % % Zbar=kron(eye(n),Z);

author: european-central-bank

repo: BEAR-toolbox


text:
function calpak_test389 ( )

%*****************************************************************************80
%
%% CALPAK_TEST389 tests MONTH_TO_MONTH_NAME_BAHAI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST389\n' );
  fprintf ( 1, '  For the Bahai calendar,\n' );
  fprintf ( 1, '  MONTH_TO_MONTH_NAME_BAHAI names the months.\n' );
  fprintf ( 1, '\n' );

  y = 1;
  months = year_length_months_bahai ( y );

  for m = 1 : months
    month_name = month_to_month_name_bahai ( m );
    fprintf ( 1, '  %2d  %s\n', m, month_name );
  end

  return
end


text:
function jnii=loadjnifti(filename, varargin)
%
%    jnii=loadjnifti(inputfile)
%       or
%    jnii=loadjnifti(inputfile, 'Param1',value1, 'Param2',value2,...)
%
%    Load a standard NIFTI-1/2 file or text or binary JNIfTI file with
%    format defined in JNIfTI specification: https://github.com/fangq/jnifti
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        inputfile: the output file name to the JNIfTI or NIFTI-1/2 file
%                *.bnii for binary JNIfTI file
%                *.jnii for text JNIfTI file
%                *.nii  for NIFTI-1/2 files
%        options: (optional) if loading from a .bnii file, please see the options for
%               loadbj.m (part of JSONLab); if loading from a .jnii, please see the
%               supported options for loadjson.m (part of JSONLab).
%
%    output:
%        jnii: a structure (array) or cell (array). The data structure can
%            be completely generic or auxilary data without any JNIfTI
%            constructs. However, if a JNIfTI object is included, it shall
%            contain the below subfields (can appear within any depth of the
%            structure)
%                jnii.NIFTIHeader -  a structure containing the 1-to-1 mapped NIFTI-1/2 header
%                jnii.NIFTIData - the main image data array
%                jnii.NIFTIExtension - a cell array contaiing the extension data buffers
%
%    example:
%        jnii=jnifticreate(uint8(magic(10)),'Name','10x10 magic matrix');
%        savejnifti(jnii, 'magic10.jnii')
%        newjnii=loadjnifti('magic10.jnii');
%
%    this file is part of JNIfTI specification: https://github.com/fangq/jnifti
%
%    License: Apache 2.0, see https://github.com/fangq/jnifti for details
%

if(nargin<1)
    error('you must provide data and output file name');
end

if(~exist('savejson','file'))
    error('you must first install JSONLab from http://github.com/fangq/jsonlab/');
end

if(regexp(filename,'\.nii$'))
    jnii=nii2jnii(filename,'jnii');
elseif(regexp(filename,'\.jnii$'))
    jnii=loadjson(filename,varargin{:});
elseif(regexp(filename,'\.bnii$'))
    jnii=loadbj(filename,varargin{:});
else
    error('file suffix must be .jnii for text JNIfTI, .bnii for binary JNIfTI or .nii for NIFTI-1/2 files');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [rVect] = getAbsPositBetweenSpacecraftAndStation(time, rVectSC, bodySC, station, celBodyData)
%getAbsPositBetweenSpacecraftAndStation Summary of this function goes here
%   Detailed explanation goes here

    stnBodyInfo = getBodyInfoByNumber(station.parentID, celBodyData);
    dVect = getAbsPositBetweenSpacecraftAndBody(time, rVectSC, bodySC, stnBodyInfo, celBodyData);
    stnRVectECIRelToParent = getInertialVectFromLatLongAlt(time, station.lat, station.long, station.alt, stnBodyInfo, [NaN;NaN;NaN]);
    rVect = dVect + stnRVectECIRelToParent;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% EX_LAPLACE_EIG_SQUARE: Compute eigenvalues and eigenvectors for the Laplace operator in the square.
%
% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_square.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [];
problem_data.drchlt_sides = [1 2 3 4];

% Physical parameters
problem_data.c_diff = @(x, y) ones(size(x));
problem_data.c_mass = @(x, y) ones(size(x));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];     % Degree of the bsplines
method_data.regularity = [2 2];     % Regularity of the splines
method_data.nsub       = [8 8];     % Number of subdivisions
method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER

[geometry, msh, space, lambda, u] = ...
                 solve_laplace_eig (problem_data, method_data);

% 4) POST-PROCESSING

if (all (imag (lambda) == 0))
  [lambda, perm] = sort (lambda);
elseif (any (abs (imag (lambda)) > 1e-9))
  error ('Complex eigenvalues appeared. I skip the postprocess.')
else
  warning ('Complex eigenvalues appeared, with small imaginary part. Only the real part is used for postprocessing')
  [lambda, perm] = sort (real (lambda));
end

% Plot of the 11th eigenfunction
figure
subplot(1,2,1)
sp_plot_solution (u(:,perm(11)), space, geometry, [31 31]);
title ('Plot of the 11^{th} eigenfunction')

% Comparison with the exact eigenvalues
ndofs_1 = repmat ([1:space.ndof_dir(1)-2], space.ndof_dir(2)-2, 1);
ndofs_2 = repmat ([1:space.ndof_dir(2)-2]', 1, space.ndof_dir(1)-2);
exact = pi * sqrt (ndofs_1.^2 + ndofs_2.^2);
exact = sort (exact(:));
spectrum = sqrt (lambda) ./ exact;
subplot(1,2,2)
plot (linspace (0, 1, numel (spectrum)), spectrum, '*')
title ('Ratio of numerical to exact eigenvalues')

%!demo
%! ex_laplace_eig_square

%!test
%! problem_data.geo_name = 'geo_square.txt';
%! problem_data.nmnn_sides   = [];
%! problem_data.drchlt_sides = [1 2 3 4];
%! problem_data.c_diff = @(x, y) ones(size(x));
%! problem_data.c_mass = @(x, y) ones(size(x));
%! method_data.degree     = [3 3];     % Degree of the bsplines
%! method_data.regularity = [2 2];     % Regularity of the splines
%! method_data.nsub       = [8 8];     % Number of subdivisions
%! method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, lambda, u] = ...
%!                  solve_laplace_eig (problem_data, method_data);
%! [lambda, perm] = sort (lambda);
%! assert (msh.nel, 64);
%! assert (space.ndof, 121);
%! assert (lambda(1:6)/pi^2, [2.00000025982960, 5.00004032992316, 5.00004032992316, 8.00008040001672, 10.00134779087585, 10.00134779087586]', 5e-13)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function computeModelCoefficients_Aerts(pathExperiments,imbalance,seed)
% -------------------------------------------------------------------------
% function computeModelCoefficients_batchHN(pathExperiments,nExp,fSetNames,imbalance,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the final logistic regression coefficients and
% bootstrap confidence intervals of the final models obtained for all
% outcomes analyzed in the HN study. See ref.[1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to the directory containing all experiments.
%                     --> Ex: '/myProject/WORKSPACE/CV-BASED_RESULTS'
% 2. nExp: Numerical value specifying the number of experiments to analyze.
%          --> Ex: 10
% 3. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (see ref.[2]).
%               --> Ex: 'IALR'
% 4. nBatch: Number of parallel batch.
%            --> Ex: 8
% 5. matlabPATH: Full path to the MATLAB excutable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% -------------------------------------------------------------------------
% OUTPUTS: Final coefficients, model response and confidence intervals
%          saved in the corresponding ../experiment/outcome/fSetName
%          folder.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

% INITIALIZATIONS
cd(pathExperiments), load('training')
cd('FINAL_MODELS'), mkdir('DeathSign'), cd('DeathSign'), mkdir('CT'), cd('CT'), pathFinalModel = pwd;
finalModel.Data = training.Death.sign.CT;
finalModel.Name = {'CT_Energy';'CT_Compactness';'CT_GLN';'CT_GLN_HLH'};
finalModel.Order = 4;
finalModel.outcome = 'Death';
save('finalModel','finalModel')

% COMPUTATIONS
[coeff,response,modelCI] = computeModelCoefficients_HN(finalModel.Data,training.Death.outcome,imbalance,seed);
save('coeff','coeff'), save('response','response'), save('modelCI','modelCI')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [theta, alpha, minError, inliers] = ObtainPlanePoint(pc,point)

    %find 10 nearest points to the selected point
    %one is the point ifself
    dis = sqrt(sum(bsxfun(@minus,pc,point).^2,2));
%     [~,index] = sort(dis);
%     n_neighbour = 200;
%     nearest_points = pc(index(1:n_neighbour),:);
    nearest_points = pc(dis<0.05,:);
    n_neighbour = size(nearest_points,1);
    n_samples = 10;

    max_inliers = 0;
    inliers = [];

    for i=1:500
        temp = randperm(n_neighbour);
        index_samples = temp(1:n_samples);
        samples = nearest_points(index_samples,:);
        [n,p] = affine_fit(samples);

        index = find(abs((bsxfun(@minus,pc,p))*n)<0.002);
        if length(index)>max_inliers
            max_inliers = length(index);
            inliers = pc(index,:);
            theta = n;
            alpha = p*n;
        end

    end

    index = find(sqrt(sum(bsxfun(@minus,inliers,p).^2,2))<0.1);
    inliers = inliers(index,:);

    if(alpha<0)
        theta=-theta;
        alpha=-alpha;
    end

    minError = median(abs(bsxfun(@minus,inliers,mean(inliers,1))*n));

end

function [n,p] = affine_fit(X)
    %Computes the plane that fits best (lest square of the normal distance
    %to the plane) a set of sample points.
    %INPUTS:
    %
    %X: a N by 3 matrix where each line is a sample point
    %
    %OUTPUTS:
    %
    %n : a unit (column) vector normal to the plane
    %V : a 3 by 2 matrix. The columns of V form an orthonormal basis of the
    %plane
    %p : a point belonging to the plane
    %
    %NB: this code actually works in any dimension (2,3,4,...)
    %Author: Adrien Leygue
    %Date: August 30 2013

    %the mean of the samples belongs to the plane
    p = mean(X,1);

    %The samples are reduced:
    R = bsxfun(@minus,X,p);
    %Computation of the principal directions if the samples cloud
    [V,~] = eig(R'*R);
    %Extract the output from the eigenvectors
    n = V(:,1)/norm(V(:,1));
end

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
function [y,MSE] = delta2static(x, delta_order, dimS, scale)
[nFr, dim] = size(x);
nDelta = dim/dimS;

if length(delta_order)==1 && nDelta > 1
    delta_order = ones(nDelta, 1) * delta_order;
end

% Initialize the static features
y = x(:,1:dimS);

D = genDeltaTransform(nFr, delta_order(1));
A = D*D;

scaleDim = size(scale, 1);
if scaleDim == 1
    ScaleD = eye(nFr)*scale(2);
    ScaleA = eye(nFr)*scale(3);
elseif scaleDim == nFr
    ScaleD = diag(scale(:,2));
    ScaleA = diag(scale(:,3));
else
    fprintf('Scale is neigher scalar or have the same dim as number of frame!\n', );
end

R = eye(nFr) + D'*ScaleD*D + A'*ScaleA*A;
p = y + D' * ScaleD * x(:,dimS+1:dimS*2) + A' * ScaleA * x(:,dimS*2+1:dimS*3);
y2 = inv(R) * p;

y = y2;

MSE = [];
if 1
    x2 = comp_dynamic_feature(x(:,1:dimS), 2, 2);
    x3 = comp_dynamic_feature(y, 2, 2);
    MSE(1,1) = sum(sum((x2(:,1:dimS)-x(:,1:dimS)).^2));
    MSE(1,2) = sum(sum((x3(:,1:dimS)-x(:,1:dimS)).^2));
    MSE(2,1) = sum(sum((x2(:,dimS+1:dimS*2)-x(:,dimS+1:dimS*2)).^2));
    MSE(2,2) = sum(sum((x3(:,dimS+1:dimS*2)-x(:,dimS+1:dimS*2)).^2));
    MSE(3,1) = sum(sum((x2(:,dimS*2+1:dimS*3)-x(:,dimS*2+1:dimS*3)).^2));
    MSE(3,2) = sum(sum((x3(:,dimS*2+1:dimS*3)-x(:,dimS*2+1:dimS*3)).^2));
end

author: singaxiong

repo: SignalGraph


text:
% Author: Housam Binous

% Forced Duffing Oscillator

% National Institute of Applied Sciences and Technology, Tunis, TUNISIA

% Email: binoushousam@yahoo.com

function xdot=duffing(t,x)

global gamma omega epsilon GAM OMEG

xdot(1)=-gamma*x(1)+omega^2*x(2)-epsilon*x(2)^3+GAM*cos(OMEG*t);
xdot(2)=x(1);

xdot=xdot';

end

author: Sable

repo: mcbench-benchmarks


text:
function [ xmin, ynewlo, icount, numres, ifault ] = nelmin ( fn, n, start, ...
  reqmin, step, konvge, kcount )

%*****************************************************************************80
%
%% NELMIN minimizes a function using the Nelder-Mead algorithm.
%
%  Discussion:
%
%    This routine seeks the minimum value of a user-specified function.
%
%    Simplex function minimisation procedure due to Nelder+Mead(1965),
%    as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
%    subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
%    25, 97) and Hill(1978, 27, 380-2)
%
%    The function to be minimized must have the form:
%
%      function value = fn ( x )
%
%    where X is a vector, and VALUE is the (scalar) function value.
%    The name of this function must be passed as the argument FN.
%
%    This routine does not include a termination test using the
%    fitting of a quadratic surface.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by R ONeill.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    John Nelder, Roger Mead,
%    A simplex method for function minimization,
%    Computer Journal,
%    Volume 7, 1965, pages 308-313.
%
%    R ONeill,
%    Algorithm AS 47:
%    Function Minimization Using a Simplex Procedure,
%    Applied Statistics,
%    Volume 20, Number 3, 1971, pages 338-345.
%
%  Parameters:
%
%    Input, real value = FN ( X ), the name of the MATLAB function which
%    evaluates the function to be minimized, preceded by an "@" sign.
%
%    Input, integer N, the number of variables.
%
%    Input, real START(N).  On input, a starting point
%    for the iteration.  On output, this data may have been overwritten.
%
%    Input, real REQMIN, the terminating limit for the variance
%    of function values.
%
%    Input, real STEP(N), determines the size and shape of the
%    initial simplex.  The relative magnitudes of its elements should reflect
%    the units of the variables.
%
%    Input, integer KONVGE, the convergence check is carried out
%    every KONVGE iterations.
%
%    Input, integer KCOUNT, the maximum number of function
%    evaluations.
%
%    Output, real XMIN(N), the coordinates of the point which
%    is estimated to minimize the function.
%
%    Output, real YNEWLO, the minimum value of the function.
%
%    Output, integer ICOUNT, the number of function evaluations.
%
%    Output, integer NUMRES, the number of restarts.
%
%    Output, integer IFAULT, error indicator.
%    0, no errors detected.
%    1, REQMIN, N, or KONVGE has an illegal value.
%    2, iteration terminated because KCOUNT was exceeded without convergence.
%
  xmin = [];
  ynewlo = [];
  icount = 0;
  numres = 0;

  ccoeff = 0.5;
  ecoeff = 2.0;
  eps = 0.001;
  rcoeff = 1.0;
%
%  Check the input parameters.
%
  if ( reqmin <= 0.0 )
    ifault = 1;
    return
  end

  if ( n < 1 )
    ifault = 1;
    return
  end

  if ( konvge < 1 )
    ifault = 1;
    return
  end

  jcount = konvge;
  dn = n;
  nn = n + 1;
  dnn = nn;
  del = 1.0;
  rq = reqmin * dn;
%
%  Initial or restarted loop.
%
  while ( 1 )

    for i = 1 : n
      p(i,nn) = start(i);
    end

    y(nn) = fn ( start );
    icount = icount + 1;

    for j = 1 : n
      x = start(j);
      start(j) = start(j) + step(j) * del;
      for i = 1 : n
        p(i,j) = start(i);
      end
      y(j) = fn ( start );
      icount = icount + 1;
      start(j) = x;
    end
%
%  The simplex construction is complete.
%
%  Find highest and lowest Y values.  YNEWLO = Y(IHI) indicates
%  the vertex of the simplex to be replaced.
%
    ylo = y(1);
    ilo = 1;

    for i = 2 : nn
      if ( y(i) < ylo )
        ylo = y(i);
        ilo = i;
      end
    end
%
%  Inner loop.
%
    while ( 1 )

      if ( kcount <= icount )
        break
      end

      ynewlo = y(1);
      ihi = 1;

      for i = 2 : nn
        if ( ynewlo < y(i) )
          ynewlo = y(i);
          ihi = i;
        end
      end
%
%  Calculate PBAR, the centroid of the simplex vertices
%  excepting the vertex with Y value YNEWLO.
%
      for i = 1 : n
        z = 0.0;
        for j = 1 : nn
          z = z + p(i,j);
        end
        z = z - p(i,ihi);
        pbar(i) = z / dn;
      end
%
%  Reflection through the centroid.
%
      for i = 1 : n
        pstar(i) = pbar(i) + rcoeff * ( pbar(i) - p(i,ihi) );
      end

      ystar = fn ( pstar );
      icount = icount + 1;
%
%  Successful reflection, so extension.
%
      if ( ystar < ylo )

        for i = 1 : n
          p2star(i) = pbar(i) + ecoeff * ( pstar(i) - pbar(i) );
        end

        y2star = fn ( p2star );
        icount = icount + 1;
%
%  Check extension.
%
        if ( ystar < y2star )

          for i = 1 : n
            p(i,ihi) = pstar(i);
          end

          y(ihi) = ystar;
%
%  Retain extension or contraction.
%
        else

          for i = 1 : n
            p(i,ihi) = p2star(i);
          end

          y(ihi) = y2star;

        end
%
%  No extension.
%
      else

        l = 0;
        for i = 1 : nn
          if ( ystar < y(i) )
            l = l + 1;
          end
        end

        if ( 1 < l )

          for i = 1 : n
            p(i,ihi) = pstar(i);
          end

          y(ihi) = ystar;
%
%  Contraction on the Y(IHI) side of the centroid.
%
        elseif ( l == 0 )

          for i = 1 : n
            p2star(i) = pbar(i) + ccoeff * ( p(i,ihi) - pbar(i) );
          end

          y2star = fn ( p2star );
          icount = icount + 1;
%
%  Contract the whole simplex.
%
          if ( y(ihi) < y2star )

            for j = 1 : nn
              for i = 1 : n
                p(i,j) = ( p(i,j) + p(i,ilo) ) * 0.5;
                xmin(i) = p(i,j);
              end
              y(j) = fn ( xmin );
              icount = icount + 1;
            end

            ylo = y(1);
            ilo = 1;

            for i = 2 : nn
              if ( y(i) < ylo )
                ylo = y(i);
                ilo = i;
              end
            end

            continue
%
%  Retain contraction.
%
          else

            for i = 1 : n
              p(i,ihi) = p2star(i);
            end
            y(ihi) = y2star;

          end
%
%  Contraction on the reflection side of the centroid.
%
        elseif ( l == 1 )

          for i = 1 : n
            p2star(i) = pbar(i) + ccoeff * ( pstar(i) - pbar(i) );
          end

          y2star = fn ( p2star );
          icount = icount + 1;
%
%  Retain reflection?
%
          if ( y2star <= ystar )

            for i = 1 : n
              p(i,ihi) = p2star(i);
            end
            y(ihi) = y2star;

          else

            for i = 1 : n
              p(i,ihi) = pstar(i);
            end
            y(ihi) = ystar;

          end

        end

      end
%
%  Check if YLO improved.
%
      if ( y(ihi) < ylo )
        ylo = y(ihi);
        ilo = ihi;
      end

      jcount = jcount - 1;

      if ( 0 < jcount )
        continue
      end
%
%  Check to see if minimum reached.
%
      if ( icount <= kcount )

        jcount = konvge;

        z = 0.0;
        for i = 1 : nn
          z = z + y(i);
        end
        x = z / dnn;

        z = 0.0;
        for i = 1 : nn
          z = z + ( y(i) - x )^2;
        end

        if ( z <= rq )
          break
        end

      end

    end
%
%  Factorial tests to check that YNEWLO is a local minimum.
%
    for i = 1 : n
      xmin(i) = p(i,ilo);
    end

    ynewlo = y(ilo);

    if ( kcount < icount )
      ifault = 2;
      break
    end

    ifault = 0;

    for i = 1 : n
      del = step(i) * eps;
      xmin(i) = xmin(i) + del;
      z = fn ( xmin );
      icount = icount + 1;
      if ( z < ynewlo )
        ifault = 2;
        break
      end
      xmin(i) = xmin(i) - del - del;
      z = fn ( xmin );
      icount = icount + 1;
      if ( z < ynewlo )
        ifault = 2;
        break
      end
      xmin(i) = xmin(i) + del;
    end

    if ( ifault == 0 )
      break
    end
%
%  Restart the procedure.
%
    for i = 1 : n
      start(i) = xmin(i);
    end

    del = eps;
    numres = numres + 1;

  end

  return
end


text:
function B = getbase(F)

B = getbase(lmi(F));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function r = roots(varargin)
%ROOTS   Zero contours of a DISKFUN.
%   R = ROOTS(F) returns the zero contours of F as a quasimatrix of
%   array-valued chebfuns. Each column of R is one zero contour. This
%   command only finds contours when there is a change of sign and it may
%   group intersecting contours in a non-optimal way. Contours are computed
%   to, roughly, four digits of precision. In particular, this command
%   cannot reliably compute isolated real roots of F or zero curves lying
%   close to the boundary of the domain.
%
%   R = ROOTS(F, G) returns the isolated points of F and G.
%
%
%   In the special case when F is of length 1 then the zero contours are found
%   to full precision.
%
% See also DISKFUN2V/ROOTS, CHEBFUN2/ROOTS, CHEBFUN2V/ROOTS

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Check for empty:
if ( isempty(varargin{ 1 }) )
    r = [];
    return
end

% convert f to a polar chebfun2
f = varargin{1};
f = cart2pol(f);

%check for a second diskfun
if ( nargin > 1 )
    if isa(varargin{2}, 'diskfun')
        g = varargin{2};
        dom = [-pi pi 0 1];
        % we require f and g to be chebfuns. Cart2pol isn't enough because
        % it keeps the rows represented in the Fourier basis.
        f = chebfun2(@(t,r) feval(f, t, r), dom);
        g = chebfun2(@(t,r) feval(g, t, r, 'polar'), dom);
        r = roots@separableApprox(f,g, varargin{3:end});

        %convert to cartesian
        if ~isempty(r)
            [x,y] = pol2cart(r(:,1),r(:,2));
            r = [x,y];
        end
        return
    end
end

rts = roots@separableApprox(f, varargin{2:end});

%convert from polar coords to complex-valued chebfuns:
x = chebpts(max(length(rts),101));

r = feval(rts,x);
rvals = imag(r).*exp(1i*real(r)); %c = re^(i*x)

r = chebfun(rvals);

% The code below creates a real-valued parametrization of the curves.
%
% Go through each component and make it an array-valued chebfun:
% for k = 1:size(vals, 2)
%     comp = feval(rts(:, k), x);
%     AX = imag(comp) .* cos(real(comp));
%     AY = imag(comp) .* sin(real(comp));
%     r{k} = chebfun([AX, AY ]);
% end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function hb_io_test09 ( input_file )

%*****************************************************************************80
%
%% HB_IO_TEST09 tests HB_FILE_READ;
%
%  Modified:
%
%    04 February 2005
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_IO_TEST09\n' );
  fprintf ( 1, '  HB_FILE_READ reads all the data in an HB file.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reading the file "%s".\n', input_file );

  input_unit = fopen ( input_file );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_IO_TEST09 - Fatal error!\n' );
    fprintf ( 1, '  Error opening the file.\n' );
    error ( 'HB_IO_TEST09 - Fatal error!' );
  end

  [ title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd, mxtype, ...
    nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, rhsfmt, rhstyp, ...
    nrhs, nrhsix, colptr, rowind, values, rhsval, rhsptr, rhsind, guess, ...
    exact ] = hb_file_read ( input_unit );
%
%  Print out the header information.
%
  hb_header_print ( title, key, totcrd, ptrcrd, indcrd, valcrd, ...
    rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, ...
    rhsfmt, rhstyp, nrhs, nrhsix );
%
%  Print the structure information.
%
  hb_structure_print ( ncol, mxtype, nnzero, neltvl, colptr, rowind );
%
%  Print the values.
%
  hb_values_print ( ncol, colptr, mxtype, nnzero, neltvl, values );

  if ( 0 < rhscrd )
%
%  Print a bit of the right hand sides.
%
    if ( rhstyp(1) == 'F' )

      r8mat_print_some ( nrow, nrhs, rhsval, 1, 1, 5, 5, '  Part of RHS' );

    elseif ( rhstyp(1) == 'M' & mxtype(3) == 'A' )

      i4vec_print_some ( nrhs+1, rhsptr, 10, '  Part of RHSPTR' );
      i4vec_print_some ( nrhsix, rhsind, 10, '  Part of RHSIND' );
      r8vec_print_some ( nrhsix, rhsvec, 10, '  Part of RHSVEC' );

    elseif ( rhstyp(1) == 'M' & mxtype(3) == 'E' )

      r8mat_print_some ( nnzero, nrhs, rhsval, 1, 1, 5, 5, '  Part of RHS' );

    end
%
%  Print a bit of the starting guesses.
%
    if ( rhstyp(2) == 'G' )

      r8mat_print_some ( nrow, nrhs, guess, 1, 1, 5, 5, '  Part of GUESS' )

    end
%
%  Print a bit of the exact solutions.
%
    if ( rhstyp(3) == 'X' )

      r8mat_print_some ( nrow, nrhs, exact, 1, 1, 5, 5, '  Part of EXACT' );

    end

  end

  fclose ( input_unit );

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To run:  Edit user options in npscalc.m and execute the code.
%
% This code is intended to accompany the paper:
%
% S. N. Friedman, G. S. K. Fung, J. H. Siewerdsen, and B. M. W. Tsui.  "A
% simple approach to measure computed tomography (CT) modulation transfer
% function (MTF) and noise-power spectrum (NPS) using the American College
% of Radiology (ACR) accreditation phantom," Med. Phys. 40, 051907-1 -
% 051907-9 (2013).
% http://dx.doi.org/10.1118/1.4800795
%
% This code is free to distribute (see below).
%
% This program requires the following 4 files:
%   npscalc.m
%   SliceBrowser.m
%   SliceBrowser.fig
%   license.txt
%
% Purpose:  To calculate the 3D NPS using CT data of the American College
%           of Radiology (ACR) accreditation phantom.
%
% Input:    Two consecutive scans of the phantom are needed.  The program
%           requires a data directory to be selected in which there are
%           only two subdirectories containing only CT slices corresponding
%           to the third module of the phantom.  Be careful of partial
%           volume effects with surrounding modules.
%
%           i.e.,  datadir
%                      |-> scan 1 dir
%                      |           | -> only module 3 slices
%                      |
%                      |-> scan 2 dir
%                                  |-> only module 3 slices
%
% CCopyright 2012 Saul N. Friedman
%   Distributed under the terms of the "New BSD License."  Please see
%   license.txt.
%
% N.B.:  SliceBrowser.m is an altered version of code by Marian Uhercik.
%        Please see license.txt.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;
clc;

fprintf('****************************************************************************\n');
fprintf('CT NPS Calculation\n\n')
fprintf('This code is intended to accompany the paper:\n');
fprintf('S. N. Friedman, G. S. K. Fung, J. H. Siewerdsen, and B. M. W. Tsui.  "A\n')
fprintf('simple approach to measure computed tomography (CT) modulation transfer \n')
fprintf('function (MTF) and noise-power spectrum (NPS) using the American College \n')
fprintf('of Radiology (ACR) accreditation phantom," Med. Phys. 40, 051907-1 - \n')
fprintf('051907-9 (2013).\n')
fprintf('http://dx.doi.org/10.1118/1.4800795\n\n')
fprintf('This code can be freely distributed.  Please see license.txt for details.\n');
fprintf('****************************************************************************\n\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User selected options
%

selectdir = 1;  % 1 = prompted for directly locations via dialogue box, 0 = hardcode locations below
    datadir = fullfile('data');
    resultsdir = fullfile('results',datadir);

saveresults = 1;  % Choose whether to save the calculated 3D NPS to disk as a series of DICOM files.

calczero = 1;  % when set to 0, the NPS(0,0,0) value is estimated based on the surrounding 26 voxel values to try remove any zero-frequency artifacts

rstep = 0.025;  % This will be the effective bin size after forming the radially averaged NPS.

showfig = 3;  % decide how important a figure should be in order to be plotted, 1 = no figures, 4 = all figures

% specify size of ROI (pixels/slices)
deltax = 128;
deltay = deltax;

calcdeltaz = 0; % 1 = calculate size of ROI in z (slices) to be same physical size as x dimension, 0 = use value below
    deltaz = 128;% % specify number of slices to use in ROI

% specify overlap fraction for ROIs in z direction
%  N.B. value is typically [0,0.5] and represents overlap with EACH
%  neighboring ROI (i.e., at 0.5 all pixels overlap with at least one
%  neighbor)

overlapz = 0.5;

% Parameters pertaining to placement of ROI in the axial plane.
%  N.B. overlap of ROIs will depend on a combination of the radius chosen
%  in rcenter and the distance between the centers of ROIs chosen.

rcenter = 25; % radius in mm along which to center ROIS, a value of 25 mm corresponds to half the radius of the phantom
deltar = 20;  % distance along the radial arc in pixels between centers of ROIs in the axial plane
deltaangle = 2*pi/180;  % the range of angles in radians over which to look for a suitable pixel on which to center an ROI

%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start of analysis code
%

homedr = cd;

if selectdir ==1
    fprintf('Select the data directory.\n\n');
    datadir = uigetdir(homedr,'Pick data directory');
    if saveresults == 1
        fprintf('Select the results directory.\n\n');
        resultsdir = uigetdir(homedr,'Pick results directory');
    end
else
    datadir = fullfile(homedr,datadir);
    resultsdir = fullfile(homedr,resultsdir);
end

if saveresults == 1
    mkdir(resultsdir)
end

cd(datadir)

% get list of scan directories

scandir = dir;

cd(homedr)

while or(strcmp(scandir(1).name,'.'),strcmp(scandir(1).name,'..'))
    scandir = scandir(2:end);
end

Nscan = size(scandir,1);

if Nscan < 2
    fprintf('At least 2 scans are needed to perform analysis.\n')
    return
end

if Nscan > 2 == 1
    Nscan = 2;
    fprintf('Only using the first 2 directories of data.\n')
end

nfilescheck = 0;

for i = 1:Nscan
    cd(fullfile(datadir,scandir(i).name))

    % get list of filenames

    ffiles = dir;

    cd(homedr)

    while or(strcmp(ffiles(1).name,'.'),strcmp(ffiles(1).name,'..'))
        ffiles = ffiles(2:end);
    end

    Nz = size(ffiles,1);

     if nfilescheck ~= Nz && nfilescheck ~=0
        fprintf('Error.  Scan directories do not contain the same number of slices.\n')
        return
    end

    nfilescheck = Nz;

    if i == 1

        % Read DICOM header info and determine correct mapping function to get HU
        % values as well as voxel sizes.

        im = single(dicomread(fullfile(datadir,scandir(i).name,ffiles(1).name)));

        imheader = dicominfo(fullfile(datadir,scandir(i).name,ffiles(1).name));

        b = imheader.RescaleIntercept;
        m = imheader.RescaleSlope;

        im = im.* m + b;

        pixelx = imheader.PixelSpacing(1);
        pixely = imheader.PixelSpacing(2);
        pixelz = imheader.SliceThickness;

        if showfig > 2
            figure
            imagesc(im)
            colormap gray
            hold on
            axis image;
            impixelinfo;
        end

        Nx = size(im,2);
        Ny = size(im,1);

        imaxisx = [1:Nx] .* pixelx;
        imaxisy = [1:Ny] .* pixely;

        imavg = zeros(Ny,Nx);
        im3D = zeros(Ny,Nx,Nz);
    end

    % Load data into memory and calculate an average 2D image to later
    % determine phantom location

    for z = 1:Nz
        im = single(dicomread(fullfile(datadir,scandir(i).name,ffiles(z).name)));
        im = im.* m + b;
        im3D(:,:,z,i) = im;
        imavg = imavg + im;
    end

end

imavg = imavg ./ Nz;

if deltaz  > Nz
    deltaz = Nz;
    val = sprintf('Selected ROI size in z is larger than the number of slices.  Using value of %d instead.\n\n',Nz);
    fprintf(val)
end

% Detrend the data and plot if desired

fprintf('Detrending data.\n\n');

im3Dzm(:,:,:) = im3D(:,:,:,1) - im3D(:,:,:,2);

if showfig > 2 && size(im3Dzm,3) > 1
    yaxis = [1:size(im3Dzm,1)] .* pixely;
    xaxis = [1:size(im3Dzm)].* pixelx;
    zaxis = [1:size(im3Dzm,3)] .* pixelz;
    axlabel = ' (mm)';

    setappdata(0,'vaxis',yaxis)
    setappdata(0,'uaxis',xaxis)
    setappdata(0,'waxis',zaxis)
    setappdata(0,'axlabel',axlabel)

    SliceBrowser(im3Dzm(:,:,:))
end

% Determine which pixels correspond to phantom cylinder

imlog = zeros(size(imavg));

level = graythresh(abs(imavg));
imlog = im2bw(abs(imavg),level);

% Calculate the center of the cylinder

fprintf('Calculating the center of the phantom and determining phantom-pixel locations.\n\n');

if exist('bwconncomp')==2
    CC = bwconncomp(imlog,26);
    for i = 1:CC.NumObjects;
        Lobj(i) =  length(CC.PixelIdxList{1,i}) ;
    end
else
    [CC CCn] = bwlabel(imlog,8);
    for i = 1:CCn
        Lobj(i) = length(find(CC == i));
    end
end

L = find(Lobj == max(Lobj));

S = regionprops(CC,'Centroid');

c.y = S(L).Centroid(2);
c.x = S(L).Centroid(1);

if showfig>3
    figure
    imagesc([imaxisx(1) imaxisx(end)], [imaxisy(1) imaxisy(end) ], imlog)
    xlabel('mm')
    ylabel('mm')
    title('Cylinder with center marked')
    colormap gray
    hold on
    axis image;
    impixelinfo;
    if exist('bwconncomp')==2
        [I J] = ind2sub(size(imavg),CC.PixelIdxList{1,L});
    else
        [I,J] = find(CC == L);
    end
    plot(J*pixely,I.*pixelx,'.r')
    plot(c.x*pixelx,c.y*pixely,'*y')
    hold off
end

xaxis = ([1:size(imavg,2)] - c.x) .* pixelx;
yaxis = ([1:size(imavg,1)] - c.y) .* pixely;

% calculate polar coordinates for to place ROIs at a specified radial
% direction

for i = 1 :size(imavg,2)
    for j = 1 : size(imavg,1)
        [th(j,i) r(j,i)] = cart2pol(xaxis(i),yaxis(j));
    end
end

rcalc = round(rcenter/pixelx);
anglestep = deltar * pixelx / rcalc;

% Calculate ROI size in z direction such that it is the same physical
% dimension as the x ROI, if desired.

if calcdeltaz == 1
    deltaz = round(deltax * pixelx ./ pixelz);
end

% Determine the number of ROIs in each direction.

Ndeltaxy = floor(2*pi/anglestep);
Ndeltaz = floor(Nz ./ (deltaz .* (1- overlapz)) - 1 / (1 - overlapz) + 1);

% Determine the offset needed to center the ROIs in z direction.
offsetz = round((Nz - deltaz - round((deltaz - overlapz * deltaz)*(Ndeltaz - 1)))/2);

if showfig > 2
    figure
    imagesc(im)
    colormap gray
    hold on
end

ROIavg = zeros(deltay,deltax,deltaz);

nROI = 0;

for i = 1:Ndeltaxy;

    R1 = find( abs(r - rcalc) < pixelx);
    R2 = find( th < -pi + i*anglestep + deltaangle);
    R3 = find( th > -pi + i*anglestep - deltaangle);
    R4 = intersect(R1,R2);
    R5 = intersect(R3,R4);
    R = min(R5);

    [I J] = ind2sub(size(im),R);

    y1 = J-round(deltay/2)+1;
    y2 = J+round(deltay/2);
    x1 = I-round(deltax/2)+1;
    x2 = I+round(deltax/2);

    for z = 1:Ndeltaz
        nROI=nROI+1;
        ROIzm(:,:,:,nROI) =  im3Dzm(y1:y2,x1:x2,(offsetz + 1 + (deltaz * (z-1))):(offsetz + deltaz + deltaz * (z-1))) ;
    end

    if showfig > 2
        plot(I,J,'r.')
        plot([x1,x2,x2,x1,x1],[y1,y1,y2,y2,y1],'r-')
    end
end

% Calculated power spectrum realizations from each ROI.  The ensemble
% average of these is the calculated estimate of the NPS.

fprintf('Calculating the NPS.\n');
fprintf('The may take several minutes. \n');
fprintf('Please be patient.\n\n');

Wavg = zeros(size(ROIzm,1),size(ROIzm,2),size(ROIzm,3));

for i = 1:nROI

    W = abs( fftn((ROIzm(:,:,:,i))) ) .^2;
    W = fftshift(W);

    Wavg = W + Wavg;
end

Wavg = Wavg ./( Ndeltaxy * Ndeltaz .* deltax .* deltay .* deltaz .* 2 ) .* pixelx .* pixely .* pixelz;

% Calculate the axes values in the frequency domain.

yaxis = ( [0:size(ROIzm,1)-1]/size(ROIzm,1) - 0.5 ) ./ pixely;
xaxis = ( [0:size(ROIzm,2)-1]/size(ROIzm,2) - 0.5 ) ./ pixelx;
zaxis = ( [0:size(ROIzm,3)-1]/size(ROIzm,3) - 0.5 ) ./ pixelz;
axlabel = ' (mm^{-1})';

setappdata(0,'vaxis',yaxis)
setappdata(0,'uaxis',xaxis)
setappdata(0,'waxis',zaxis)
setappdata(0,'axlabel',axlabel)

% The calculated zero-frequency value often contains an artifact.  The value
% can be approximated based on an average of the surrounding values.

if calczero == 0;
    J = find(abs(yaxis)==min(abs(yaxis)));
    I = find(abs(xaxis)==min(abs(xaxis)));
    if size(Wavg,3) > 1
        K = find(abs(zaxis)==min(abs(zaxis)));
        Wavg(J,I,K) = (Wavg(J,I,K+1) + Wavg(J,I,K-1) + Wavg(J,I+1,K) + Wavg(J,I+1,K+1) + Wavg(J,I+1,K-1) + Wavg(J,I-1,K) + Wavg(J,I-1,K+1) + Wavg(J,I-1,K-1) + Wavg(J+1,I,K)  + Wavg(J+1,I,K+1) + Wavg(J+1,I,K-1) + Wavg(J+1,I+1,K) + Wavg(J+1,I+1,K+1) + Wavg(J+1,I+1,K-1) + Wavg(J+1,I-1,K)+ Wavg(J+1,I-1,K+1)+ Wavg(J+1,I-1,K-1)+ Wavg(J-1,I,K) + Wavg(J-1,I,K+1)+ Wavg(J-1,I,K-1)+ Wavg(J-1,I+1,K) + Wavg(J-1,I+1,K+1) + Wavg(J-1,I+1,K-1) + Wavg(J-1,I-1,K) + Wavg(J-1,I-1,K+1) + Wavg(J-1,I-1,K-1)) / 26;
    else
        Wavg(J,I) = (Wavg(J,I+1) + Wavg(J,I-1) + Wavg(J+1,I) + Wavg(J+1,I+1) + Wavg(J+1,I-1) + Wavg(J-1,I) + Wavg(J-1,I+1) + Wavg(J-1,I-1)) ./ 8;
    end
end

% Plot a 1D slightly off-centered NPS profile in each direction.

eps = 0;

if showfig > 2
    figure
    subplot(3,1,1)
    plot(xaxis,Wavg(:,round(deltay/2)+eps,round(deltaz/2)+eps));
    xlabel('x-dir (mm^{-1})')
    ylabel('NPS_x (HU^2 mm^3)')
    subplot(3,1,2)
    plot(yaxis,Wavg(round(deltax/2)+eps,:,round(deltaz/2)+eps));
    xlabel('y-dir (mm^{-1})')
    ylabel('NPS_y (HU^2 mm^3)')
    subplot(3,1,3)
    zvalue = squeeze(Wavg(round(deltay/2)+eps,round(deltax/2)+eps,:));
    plot(zaxis,zvalue,'.-');
    xlabel('z-dir (mm^{-1})')
    ylabel('NPS_z (HU^2 mm^3)')
end

% Plot the 3D NPS.

if showfig > 1 && size(Wavg,3) > 3
    SliceBrowser(Wavg)
end

% Calculate radially averaged NPS profile for xy plane

% Switch to polar coords
imslice = Wavg(:,:,round(deltaz/2)+eps);
r = zeros(size(imslice));
th = zeros(size(imslice));

for i = 1 :size(imslice,2)
    for j = 1 : size(imslice,1)
        %         r(j,i) = sqrt(xaxis(i).^2 + yaxis(j).^2);
        %         th(j,i) = atan2(yaxis(j),xaxis(i));
        [th(j,i) r(j,i)] = cart2pol(xaxis(i),yaxis(j));
    end
end

% Create oversampled radial NPS profile

rmax = max(max(r));

npsr = zeros(ceil(rmax/rstep)+1,1);
nsamp = zeros(length(npsr),1);

rbin = 0:rstep:rmax+rstep;

if showfig > 3
    figure
    imagesc([xaxis(1) xaxis(end)], [yaxis(1) yaxis(end)], imslice)
    colormap gray
    hold on
    plot(0,0,'r*')
    hold off
    axis image;
    impixelinfo;

    figure
    imagesc(imslice)
    colormap gray
    hold on
    axis image;
    impixelinfo;
end

% data with radius falling within a given bin are averaged together for a
% low noise approximation of the ESF at the given radius

for i = 1:length(rbin)
    R1 = find(r >= rbin(i));
    R2 = find(r < rbin(i) + rstep);
    R = intersect(R1,R2);
    if showfig > 3
        [X Y] = ind2sub(size(imslice),R);
        plot(Y,X,'.')
        hold all
    end
    npsr(i) = sum(imslice(R));
    nsamp(i) = length(R);
end

i1 = find(nsamp,1,'first');
i2 = find(nsamp,1,'last');
nsamp = nsamp(i1:i2);
npsr = npsr(i1:i2);
rbin = rbin(i1:i2);

I = find(nsamp > 0);
npsr(I) = npsr(I)./nsamp(I);

% Plot radially averaged NPS
if showfig > 2
    figure
    plot(rbin,npsr,'.-')
    xlabel('spatial frequency (mm^{-1})')
    ylabel('NPS_r (HU^2 mm^3)')
    grid on
    axis([rbin(1) rbin(end) -1.2*min(npsr) 1.1*max(npsr)])
end

fprintf('NPS calculated.\n\n');

% Save 3D NPS data as DICOM slices and the 1D profile along w = 0 as text
% files if indicated.

if saveresults == 1

    NPSheader = imheader;

    NPSheader.RescaleIntercept = 0;
    NPSheader.RescaleSlope = 1;
    NPSheader.PixelSpacing(1) = xaxis(2) - xaxis(1);
    NPSheader.PixelSpacing(2) = yaxis(2) - yaxis(1);
    NPSheader.SliceThickness = zaxis(2) - zaxis(1);

    mkdir(fullfile(resultsdir,'3D'));

    for k = 1:size(Wavg,3)
        filename = sprintf('NPS%04d.dcm',k);
        dicomwrite(uint16(Wavg(:,:,k)),fullfile(resultsdir,'3D',filename),NPSheader)
    end
    mkdir(fullfile(resultsdir,'1D'));

    save(fullfile(resultsdir,'1D','axisvalues.txt'),'rbin','-ASCII')
    save(fullfile(resultsdir,'1D','NPSrvalues.txt'),'npsr','-ASCII')

    fprintf('Results saved.\n\n');

end

return;
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
% Usage:
% spheres_out = sphere_3d(sphere_centers, r, [unit_length]);
%
% sphere_centers is row or column vector or matrix with rows equal to a
% point in some 3d space, must have 3 columns or rows. If both, it is
% assumed that X Y Z dimensions are specified by different rows, and that
% columns are points. Output is in the same format as the input.
%
% r is a vector of radii with length equal to the number of rows of
% sphere_centers. A scalar input may also be used, resulting in the same
% radii for all spheres.
%
% spheres_out is a matrix with rows equal to all points in the
% spheres defined by sphere_centers and r
%
% unit_length specifies the unit size of the space in each dimension (or a
% scalar for isotropic sizes). default is 1.
%
% sphere_3d assumes that the point 0, 0, 0 exists in the space, and will only
% report points that are integer multiples of unit length distance from
% 0, 0, 0.

function spheres_out = sphere_3d(sphere_centers, r, varargin)
    flip = 0;
    if isempty(sphere_centers) || isempty(r) || (max(size(sphere_centers)) ~= 3 & min(size(sphere_centers))~= 3), error('Error: incorrect inputs'), end
    if size(sphere_centers, 1) ~= 3
        sphere_centers = sphere_centers';
        flip = 1;
    end

    for i = 1:length(varargin)
        if isnumeric(varargin{i}),
            if length(varargin{i}) == 1, unit_length = [varargin{i} varargin{i} varargin{i}];
            elseif length(varargin{i}) == 3, unit_length = varargin{i};
            else error('Error: unit_length incorrectly specified')
            end
        else disp(['Warning: do not recognize input ' num2str(i+nargin-length(varargin))])
        end
    end
    if ~exist('unit_length', 'var'), unit_length = [1 1 1]; end

    if length(r) == 1, r(1:size(sphere_centers, 2)) = r; end

    spheres_out = [];

    for i = 1:size(sphere_centers, 2)
        boxmin = space_match(sphere_centers(:,i)-r(i), unit_length, 'force');
        boxmax = space_match(sphere_centers(:,i)+r(i), unit_length, 'force');
        candidates = zeros([3 length(boxmin(1):unit_length(1):boxmax(1))*length(boxmin(2):unit_length(2):boxmax(2))*length(boxmin(3):unit_length(3):boxmax(3))]);
        count = 0;
        for m = boxmin(1):unit_length(1):boxmax(1)
            for n = boxmin(2):unit_length(2):boxmax(2)
                for p = boxmin(3):unit_length(3):boxmax(3)
                    count = count+1;
                    candidates(:,count) = [m n p]';
                end
            end
        end
        for k = 1:size(candidates, 2)
            if sqrt(sum((sphere_centers(:,i)-candidates(:,k)).^2))>r(i)
                candidates(:,k) = NaN;
            end
        end
        [row, col] = find(isnan(candidates));
        col = unique(col);
        if ~isempty(col), candidates(:,col) = []; end
        spheres_out(:,end+1:end+size(candidates, 2)) = candidates;
    end
    spheres_out = space_match(spheres_out, unit_length);
    spheres_out = unique(spheres_out', 'rows')';
    if isempty(spheres_out), spheres_out = space_match(sphere_centers, unit_length, 'force'); end
    if flip, spheres_out = spheres_out'; end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeter'};
weights = [1 0.1];
constraint = {'volume'};
optimizer = 'IPOPT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan