text:
function [alpha_log_all,alpha_log_hc,h_0_threshold,n_soundings] = sounding_new(start_date,end_date,plot_flag,hydro,wet)
% This function computes the refractivity and integrated LOS delays.
% In addition it will estimate the power-law decay cofficient and height at which the relative
% tropopsheric delays are approximately zero. Both coefficients are given based on the full heifght range
% and the lower height range up to hc. Note that alpha is dimensionless and
% h0 is the height in m units. Sounding data needs to be stored in the
% correct format! Read below:
%
% ---- Input ----
% All inputs will be retreived from the aps_parms file. Set the inputs by change
% setparm_aps('Parameter_name',new_value)
%
% Variables include:
% look_angle     Mean look angle of the data in rad. By default this is
%                assumed to be 21 degees.
% lambda         Radar wavelength in m. By default this is assumed to be
%                0.056 cm (C-band).
% h_0_threshold  Consider only soundings that go up till the h0 threshold in km.
%                Being the height at which the differential delay is small.
%                When '0' is specified it is estimated from the data. Note that
% 		 soundings should go up high enough to cover the lower tropopshere (up to 16 km).
%                By default the threshold is set to be 15 km.
% h_thres_hc     Height range [0 h_thres_hc] over which the power-law decay
%                coefficient is estimated computed. By default this is set to be 4 km.
% start_date	 Start date of the sounding period, by default [], full period
%                is considered. Specify as a string in 'yyyymmdd' format.
% end_date       End date of the sounding period, by default [], full period
%                is considered. Specify as a string in 'yyyymmdd' format.
% time_stamp	 Include a time stamp to it. This is a column vector with
%                strings e.g. ['00';'12'] for 00Z and 12Z. Only those files
%                ending with this are considered in the computation.
% sounding_dir	 Optional argument giving directly the full path to the soundings.
%                The files on this path should be the YYYYMMDD_HH.mat files.
% error_promp_flag By default ([] or 'y') this is turn on. Can be usefull to turn of
%                when running in a batch mode. Instead NaN values will be outputed.
%
% Output:
% alpha and h0. The computed sounding delays: refractivity, the (mean) LOS delay [m],
% the (mean) LOS phase delay [rad], and corresponding heights are appended to the existing
% sounding .mat files.
%
%
% NOTE on data format:
% sounding data needs to be stored in a folder called sounding_data,
% within your processing directory. Within this folder, each sounding
% aquisition needs to be stored in as 8 digit .mat files, e.g. YYYYMMDD_HH.mat
% format, with the  pressure (hPa), temperature (degree), relative humidity (%)
% and heights (m) as matlab variables P, T, RH and h.
%
% OPTIONAL:
% - visulize intermediate results (refractivity, delay,
%   mean delay and hc relation with H), by putting plot_flag to 1. By
%   default this is not done.
% - Overwrite all the data to recompute the delays by putting recompute flag to 1
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% September 2010 --- Bekaert David
% Modifications:
% 01/11/2011:  DB    Modify from initial version
% 01/12/2011:  DB    Making compatible to other datasets
% 01/11/2012:  DB    Made program more efficient
% 05/11/2012:  DB    Fix bug in loading of the files
% 10/02/2013:  DB    Integrate with other part of toolbox
% 18/03/2013:  DB    Include a time stamp to it.
% 18/03/2013:  DB    Include a start and end time.
% 18/03/2013:  DB    Account for NaN values by removing them for all data columns
% 19/03/2013:  DB    Set the min height range of the delay curve by comparing it to
%                    the median of the height range
% 19/03/2013:  DB    Include a flag to suppress error messages by lack of sounding data and
%                    directly pass NaN value through it.
% 25/03/2013:  DB    Compute the power law coefficient from loglog relation
% 03/04/2013:  DB    Include h0 estimation option from the net delays
% 24/04/2013:  DB    Incorporate in the bigger aps toolbox. Change towards
%                    loglog powerlaw computation.
% 10/05/2013:  DB    Loading the default parameters from the parms_aps file.
% 02/10/2013:  DB    Change filename to sounding and clean the syntax
% 29/10/2013:  DB    Allow the look angle to be a specified as a file
% 29/11/2013:  DB    Get stamps processing flag before lookangle  loading
% 06/11/2014:  DB    Add hydrostatic and wet delay options

% --------- VARIABELS ---------- %
fontsize = 20;                      % fontsize of the figures
save_fig =0;                        % when 1 save figures
recompute = 0;                      % Use earlier computed data unless it is missing
if nargin<3 || isempty(plot_flag)
    plot_flag =0;                       % Control plots of the refractivity, delay and net delay.
end
if nargin<4
   hydro=1;
end
if nargin<5
    wet=1;
end
netdelay_color = [];
curdir = pwd;

hydro
wet

%% Setting the defaults were needed
look_angle = getparm_aps('look_angle');
lambda = getparm_aps('lambda');
h_0_threshold= getparm_aps('sounding_h0');
h_thres_hc= getparm_aps('sounding_h_alpha_thres');

if nargin<1 || isempty(start_date)
    start_date= getparm_aps('sounding_start_date');
end
if nargin<2 || isempty(end_date)
    end_date= getparm_aps('sounding_end_date');
end
time_stamp= getparm_aps('sounding_time_stamp');
sounding_dir= getparm_aps('sounding_dir');
error_promp_flag= getparm_aps('sounding_error_promp');

clear h_0_threshold_default h_thres_hc_default look_angle_default lambda_default

stamps_processed = getparm_aps('stamps_processed');

% checking if the look angle is a file or not
if ischar(look_angle)
	% look angle is specified as a file [rad]
	if strcmp(stamps_processed,'y')
           look_angle = load(look_angle);
	       look_angle = look_angle.la;
           look_angle = mean(look_angle);
	else
	       look_angle = load(look_angle);
	       % use the mean of the look angles
	       look_angle = mean(look_angle);
	end
end

% convert units
theta = look_angle;                       % mean angle of incidence [rad]
h_0_threshold = h_0_threshold*1000;       % put threshold to [m]
h_thres_hc = h_thres_hc*1000;             % put threshold to [m]
clear look_angle

% get height information
if plot_flag==1
    hgt_matfile = getparm_aps('hgt_matfile');
    hgt = load(hgt_matfile);
    max_height = max(hgt.hgt);
    min_height = min(hgt.hgt);
    clear hgt hgt_matfile
end

% Include a while loop. Depending on error message flag, it will break out of the loop and pass
% NaN values as output.
continue_flag =1;
abord_flag = 0;		% when 1 abnormal termination get NaN values instead.
while continue_flag

%% Getting the file list of the sounding data
if isempty(sounding_dir)~=1
	if exist([sounding_dir filesep],'dir')~=7
        error('myApp:argChk', ['The specified filepath of the sounding data does not exist,...  \nAbort,... \n'])
    end
	cd(sounding_dir)
else
	if exist('sounding_data/','dir')~=7
    		error('myApp:argChk', ['There is no sounding_data directory,...  \nAbort,... \n'])
	end
	cd  sounding_data
end
if exist('sounding.list','file')~=2
    % making a list of all the sounding files
    [dummy dummy2] = system('echo sounding_list > sounding.list');
    clear dummy dummy2
    for k=1:size(time_stamp,1)
        command_str = ['ls [0-9]???????_' time_stamp(k,:) '.mat >> sounding.list'];
        [dummy dummy2] = system(command_str);
        clear dummy dummy2
    end
end
temp = tdfread('sounding.list');
[dummy dummy2] = system('rm sounding.list');
clear dummy dummy2
date_list_temp = temp.sounding_list(:,[1:8]);
sounding_list = temp.sounding_list;

% selecting a date range when requested
clear ix

if isempty(start_date)~=1 && isempty(end_date)~=1
	ix = find(datenum(date_list_temp,'yyyymmdd')>=datenum(start_date,'yyyymmdd') & datenum(date_list_temp,'yyyymmdd')<=datenum(end_date,'yyyymmdd'));
	save_name = ['Delay_' start_date '_' end_date];
elseif isempty(start_date)~=1  && isempty(end_date)==1
        ix = find(datenum(date_list_temp,'yyyymmdd')>=datenum(start_date,'yyyymmdd'));
        save_name = ['Delay_' start_date '_end' ];
elseif isempty(start_date)==1 && isempty(end_date)~=1
        ix = find(datenum(date_list_temp,'yyyymmdd')<=datenum(end_date,'yyyymmdd'));
        save_name = ['Delay_start_' end_date];
else
        save_name = ['Delay'];
end
% checking if the save figures folder exist
if save_fig ==1 & exist('figures','dir')~=7
    mkdir('figures')
end

if  isempty(start_date)~=1 || isempty(end_date)~=1
	if isempty(ix)~=1
		date_list_temp = date_list_temp(ix,:);
		sounding_list = sounding_list(ix,:);
	else
		fprintf('No sounding has been acquired in this period \n' )
		continue_flag = 0;
                abord_flag = 1;
                break
	end
	clear temp ix
end

%% Computing refractivity
fprintf('Computing the refractivity \n')
ix_skip_sounding = [];			% in case soundings are rejected along the way
warning('off','MATLAB:load:variableNotFound')
for i=1:size(sounding_list,1)
    temp = load(sounding_list(i,:),'h_min_souding');
    if isfield(temp,'h_min_souding') && recompute==0
        % storing variables for latter processing
        h_range(i,1) = temp.h_min_souding;          % [m]

        if  plot_flag==1
            data = load(sounding_list(i,:),'h_max_souding','Ndelay');
            Ndelay = data.Ndelay;
        else
            data = load(sounding_list(i,:),'h_max_souding');
        end
        h_range(i,2) = data.h_max_souding;          % [m]
        clear data
        skip_sounding = 0;

    else
        load(sounding_list(i,:));
        % This includes variables:
        % - P (pressure in [hPa])
        % - T (temperature in [degrees])
        % - RH (relative humidity in [%])
        % - h (altitude in [m])

        % coping with NaN values in the RH data
        ix1 = find(isnan(RH)==1);
        ix2 = find(isnan(h)==1);
        ix = unique([ix1; ix2]);
        P(ix)=[];
        RH(ix)=[];
        T(ix)=[];
        h(ix)=[];
        clear ix ix1 ix2

        % coping with errors in the data when a double recording was made
        ix_repeat = find(diff(sort(h))==0);
        h(ix_repeat) = [];
        P(ix_repeat)=[];
        RH(ix_repeat)=[];
        T(ix_repeat)=[];
        clear ix_repeat

        % Checking if there is still data left
        if isempty(h)==1 && size(sounding_list,1)==1 && strcmp(error_promp_flag,'y')
             error('myApp:argChk', ['Datafile contains no numeric data. \n'])
        elseif isempty(h)==1
            fprintf([sounding_list(i,:) ' sounding skipped as it containes no numeric data\n'])
            skip_sounding = 1;
            ix_skip_sounding = [ix_skip_sounding ; i];
        else
            skip_sounding = 0;
        end

        % when possible compute the refractivity
        if skip_sounding==1
            % case of no data
            h_range(i,1) = NaN;          % [m]
            h_range(i,2) = NaN;          % [m]
        else
            % saterated pressure es
            Rv = 461.524;                   % [J/K]
            T0 = 273.16;                    % [k]
            L  = 2.5e6;                     % [J/kg]                Latent heat
            e0 = 6.11;                      % [hPa]
            T  = T+273.15;                  % [degree] to [K]       Temperature
            es = e0.*exp(L./Rv.*(1./T0-1./T));

            % partial pressure of water vapour e
            RH = RH./100;                   % [%] to [-]            Relative humidity
            e  = RH.*es;                    % [hPa]

            % Refractivety - Hydrostatic term
            k1 = 77.6;                      % [K/hPa]
            Ndelay.Nhydro = k1.*P./T;              % [K/hPa*hPa/K] = [-]

            % Refractivity - Wet term
            k2_n = 23.3;                    % [K/hPa]
            k3   = 3.75e5;                  % [K^2/hPa]
            Ndelay.Nwet = k2_n.*e./T+k3.*e./T.^2;  % [ppm]

            % refractivety
            Ndelay.N = Ndelay.Nhydro+Ndelay.Nwet;

            % save the heights of the delay as well
            Ndelay.h = h;

            h_min_souding = min(h);
            h_max_souding = max(h);

            % saving the data
            save(sounding_list(i,:),'-append','Ndelay','h_min_souding','h_max_souding')

            % storing variables for latter processing
            h_range(i,1) = h_min_souding;          % [m]
            h_range(i,2) = h_max_souding;          % [m]
            clear h_max_souding h_min_souding

        end
    end

    if floor(i/10)*10 == i
        fprintf([num2str(i) ' refraction dates completed out of ' num2str(size(sounding_list,1)) ' \n']);
    end

    % Test visualizing refractivity with height
    if  plot_flag==1 && skip_sounding~=1
        if i==1
           figure1 = figure('name','Refractivity with height');
           ylabel('Height [km]','fontsize',12)
           xlabel('Refractivity [ppm]','fontsize',12)
           line_colors = hsv(size(sounding_list,1));
           box on
           set(gca,'fontsize',12)
        end
        figure(figure1)
        hold on
        if hydro==1 && wet==0
            if i==1
                title('Refractivity with height','fontsize',12)
            end
            plot(Ndelay.Nhydro,Ndelay.h./1000,'-','color',line_colors(i,:))
        elseif hydro==0 && wet==1
            if i==1
                title('Refractivity with height (hydro)','fontsize',12)
            end
            plot(Ndelay.Nwet,Ndelay.h./1000,'-','color',line_colors(i,:))
        else
            if i==1
                title('Refractivity with height (wet)','fontsize',12)
            end
            plot(Ndelay.N,Ndelay.h./1000,'-','color',line_colors(i,:))
        end
    end
    clear P h RH T Ndelay N
end
clear Rv T0 L e0 es e k1 k2_n k3 i figure1
% remove those sounding acqusitions that have no data coverage
if isempty(ix_skip_sounding)~=1
	sounding_list(ix_skip_sounding,:)=[];
	h_range(ix_skip_sounding,:)=[];
end

%% Computing the height range
fprintf('Computing delay \n')

% defining the height range (equal for the whole dataset)
h_max = min(h_range(:,2));
h_max_mean = mean(h_range(:,2));

if h_0_threshold<=0
  % estimate h0 from the data
  h0_estimate_flag = 1;
  h_0_threshold = 0;			% set the height threshold low such all the delays can be used to compute a net delay
else
  h0_estimate_flag = 0;
end

% check if the delays goes high enough
if h_max<h_0_threshold
    % the set threshold is lower than the lowest value.
    % setting the h_max to h0, but checking the number of soundings that need to be rejected for this.
    ix = find(h_range(:,2)<h_0_threshold);
    n = [1:size(sounding_list,1)]';
    n(ix)=[];
    h_range(ix,:)=[];
    if isempty(ix)==1
        if strcmp(error_promp_flag,'y')
            error('myApp:argChk', ['None of the soundings reach higher than ',num2str(h_max/1000),' km. Try setting h_0_threshold lower. \n'])
        else
            fprintf(['None of the soundings reach higher than ',num2str(h_max/1000),' km. Recommend try to set h_0_threshold lower \n'])
            continue_flag = 0;
            abord_flag = 1;
            break
        end
    end
    clear ix
    if size(n,1)<3 || isempty(n)==1
        if strcmp(error_promp_flag,'y')
            error('myApp:argChk', 'Put h_0_threshold to a lower value. \n')
        else
            fprintf('Try to put h_0_threshold to a lower value \n')
            continue_flag = 0;
            abord_flag = 1;
            break
        end
    end
    if size(n,1)<5
        fprintf(['Only ',num2str(size(n,1)),' soundings are used to compute mean delay \n']);
    end
    sounding_list = sounding_list(n,:);
    h_max = h_0_threshold;
    clear n
elseif h_max>h_0_threshold && h0_estimate_flag==0
    % h0 is lower than the h_max.
    h_max = h_0_threshold;
elseif h_max_mean-h_max >1000 && h0_estimate_flag==1
    % Make h_max larger possible as the h0 value needs to be estimated
    % this will be at the cost of some soundings but they will be included
    % latter on
    h_max = h_max_mean;
    ix = find(h_range(:,2)<h_max);
    h_range(ix,:)=[];
    sounding_list(ix,:) = [];
end

% including a check for the minimum range by comparing it to the median value
h_min_median = median(h_range(:,1));
h_min = max(h_range(:,1));

% removing soundings that do not start low enough
if h_min-h_min_median > 1000
    h_min = h_min_median;

    ix = find(h_range(:,1)>h_min);
    n = [1:size(sounding_list,1)]';
    n(ix)=[];
    h_range(ix,:)=[];
    if isempty(ix)==1
         error('myApp:argChk', ['This should not occur.'])
    end
    clear ix
    if size(n,1)<3 || isempty(n)==1
        error('myApp:argChk', 'Too few soundings are left for mean delay computation. Check for incomplete data files or differences in sounding height ranges \n.')
    end
    if size(n,1)<5
        fprintf(['Only ',num2str(size(n,1)),' soundings are used to compute mean delay \n']);
    end
    sounding_list = sounding_list(n,:);
    clear n
end

fprintf([num2str(size(sounding_list,1)),' soundings are used to compute mean delay \n']);
h_delay = [h_min:1:h_max]';                                              % set to 1 [m] interval
clear h_range h_min h_max

%% computing the integrated refractivity
for i=1:size(sounding_list,1)
    compute_flag=0;
    % checking if it has already been computed. If so use this data at once
    % else compute it
    temp = load(sounding_list(i,:),'h_delay');
    if isfield(temp,'h_delay') && recompute==0
        h_delay_temp=temp.h_delay;
        if h_delay(1)==h_delay_temp(1) && h_delay(end)==h_delay_temp(end)
            if hydro==1 && wet==0
                data = load(sounding_list(i,:),'delay_hydro','phase_delay_hydro');
                if ~isfield(data,'delay_hydro')
                   compute_flag=1;
                else
                    delay = data.delay_hydro;
                    phase_delay = data.phase_delay_hydro;
                end
            elseif hydro==0 && wet==1
                data = load(sounding_list(i,:),'delay_wet','phase_delay_wet');
                if ~isfield(data,'delay_wet')
                    compute_flag=1;
                else
                    delay = data.delay_wet;
                    phase_delay = data.phase_delay_wet;
                end
            else
                data = load(sounding_list(i,:),'delay','phase_delay');
                if ~isfield(data,'delay')
                    compute_flag=1;
                else
                    delay = data.delay;
                    phase_delay = data.phase_delay;
                end
            end

            if compute_flag==0
                % store all delays for later mean LOS delay computation
                delay_matrix(1:length(delay),i) = delay;
                % do the same for the phase delay
                phase_delay_matrix(1:length(delay),i) = phase_delay;
                compute_flag=0;
            end
        else
           % delay needs to be computed
           compute_flag=1;
        end
    else
        compute_flag=1;
    end

    if compute_flag==1
        load(sounding_list(i,:),'Ndelay')

        % interpolation to 1 m interval
        if hydro==1 && wet==0
           N_regular = interp1(Ndelay.h,Ndelay.Nhydro,h_delay,'linear');
        elseif hydro==0 && wet==1
           N_regular = interp1(Ndelay.h,Ndelay.Nwet,h_delay,'linear');
        else
           N_regular = interp1(Ndelay.h,Ndelay.N,h_delay,'linear');
        end

        % integration of refractivity with height (= delay) and projection on the LOS
        delay = zeros([length(N_regular) 1]);                                   % initialize
        delay = flipud(cumsum(flipud(N_regular))./10^6.*cos(theta));            % 10^6 factor is for correction of ppm, theta for LOS projection
        phase_delay = delay*4*pi./lambda;
        clear scale N_regular

        % store all delays for later mean LOS delay computation
        delay_matrix(1:length(delay),i) = delay;
        % do the same for the phase delay
        phase_delay_matrix(1:length(delay),i) = phase_delay;

        % saving the data
        if hydro==1 && wet==0
            delay_hydro = delay;
            phase_delay_hydro = phase_delay;
            save(sounding_list(i,:),'-append','delay_hydro','phase_delay_hydro','h_delay')
        elseif hydro==0 && wet==1
            phase_delay_wet=phase_delay;
            delay_wet = delay;
            save(sounding_list(i,:),'-append','delay_wet','phase_delay_wet','h_delay')
        else
            save(sounding_list(i,:),'-append','delay','phase_delay','h_delay')
        end

    end

    % Test visualizing delay with height
    if plot_flag==1
        if i==1
           figure2 = figure('name','Delay with height');
           xlabel('LOS phase delay [rad]','fontsize',fontsize)
           ylabel('Height [km]','fontsize',fontsize)
           set(gca,'fontsize',fontsize)
        end
        figure(figure2)
        hold on
        plot(phase_delay,h_delay./1000,'-','color',[0.7 0.7 0.7])
        box on

        if i==size(sounding_list,1) && h0_estimate_flag==1
            Ax1 = gca;
            ylimits = get(Ax1,'ylim');
            set(Ax1,'ylim',[0 ylimits(2)])
            % adding a second axis
            box on
            set(gca,'fontsize',fontsize)
            if hydro==1 && wet==0
                title({'LOS hydro delay [m]',' '},'fontsize',fontsize)
                fig_save_name = ['figures' filesep 'Delay_hydro_curves_all.eps'];
            elseif hydro==0 && wet==1
                title({'LOS wet delay [m]',' '},'fontsize',fontsize)
                fig_save_name = ['figures' filesep 'Delay_wet_curves_all.eps'];
            else
                title({'LOS delay [m]',' '},'fontsize',fontsize)
                fig_save_name = ['figures' filesep 'Delay_curves_all.eps'];
            end
            Ax2 = axes('Position',get(Ax1,'Position'),'XAxisLocation','top');
            set(Ax2,'ylim',get(Ax1,'ylim'),'ytick',[],'yticklabel','');
            xlims = get(Ax1,'xlim');
            xtick_new = get(Ax1,'xtick')./xlims(2);
            new_x_label = get(Ax1,'xtick').*lambda./(4*pi);
            xtick_labels = round((new_x_label)*100)./100;
            xtick_labels = num2str(xtick_labels');
            set(Ax2,'xtick',xtick_new,'xticklabel',xtick_labels);
            set(Ax2,'color','none');
            set(Ax2,'fontsize',fontsize)
            box on
            % saving of the figure when requested
            if save_fig==1
                set(figure2,'PaperPositionMode','auto')
                print(figure2,'-depsc','-r150',fig_save_name)
                clear fig_save_name
            end
            clear figure2
        end
    end
    clear Ndelay.N Ndelay.h delay
end
clear i theta

n_soundings = size(delay_matrix,2);

if n_soundings <=3
   continue_flag = 0;
   abord_flag = 1;
   break
end

%% estimating h0 from the delay curves, when required.
% this is based on the net delays and used the std
% as criteria for setting h0
if h0_estimate_flag==1
    n_netdelays_max = 400;
    if (n_soundings*n_soundings-n_soundings)/2<=n_netdelays_max
       	fprintf(['Estimating h0 from all (', num2str((n_soundings*n_soundings-n_soundings)/2),') possible net delay combinations \n'])
    	ix1_random=repmat(1:n_soundings,n_soundings,1);
        ix1_random=single(ix1_random);
    	ix2_random=ix1_random';
    	ix_temp=logical(tril(ones(n_soundings)))&~eye(n_soundings);
    	ix_random=[ix1_random(ix_temp),ix2_random(ix_temp)];
        clear ix1_random ix2_random
    else
        fprintf(['Estimating h0 from (', num2str(n_netdelays_max),') net delay combinations \n'])
        % take 2 random sets of n_netdelays_max
        new_set_search = 1;
        loop_counter = 0;
        while new_set_search
            ix1_random = ceil(rand(n_netdelays_max*5,1)*n_soundings);
            ix2_random = ceil(rand(n_netdelays_max*5,1)*n_soundings);
            ix_random = [ix1_random ix2_random];
            clear ix1_random ix2_random

            % remove those net delays of the same date
            ix_temp = find(ix_random(:,1)-ix_random(:,2)==0);
            ix_random(ix_temp,:)=[];
            clear ix_temp

            % remove repetition of pair combination
            ix_random = unique(ix_random,'rows');

            % check if enough combinations are left otherwize rerun
            if size(ix_random,1)>=n_netdelays_max
               ix_random(n_netdelays_max+1:end,:)=[];
               new_set_search=0;

            end
            loop_counter = loop_counter+1;
            if loop_counter == 50
                fprintf('Having difficulty determining net delay combinations \n')
                keyboard
            end
        end
    end

    % computation of the net delay
    netdelays = delay_matrix(:,ix_random(:,1)) - delay_matrix(:,ix_random(:,2));
    netphasedelays = phase_delay_matrix(:,ix_random(:,1)) - phase_delay_matrix(:,ix_random(:,2));

    if plot_flag==1
        fig_netdelay = figure('name','Net delays');
        if isempty(netdelay_color)
            plot(netdelays,h_delay/1000)
        else
            plot(netdelays,h_delay/1000,'color',netdelay_color)
        end
        max_spacing = max([abs([min(min(netdelays)) max(max(netdelays))])]);
        xlim([-1.1*max_spacing max_spacing*1.1]);
        xlabel('Net LOS delay [m]','fontsize',fontsize)
        ylabel('Height [km]','fontsize',fontsize)
        set(gca,'fontsize',fontsize)

        fig_netphasedelay = figure('name','Net phase delays');
        if isempty(netdelay_color)
            plot(netphasedelays,h_delay/1000)
        else
            plot(netphasedelays,h_delay/1000,'color',netdelay_color)
        end
        max_spacing = max([abs([min(min(netphasedelays)) max(max(netphasedelays))])]);
        xlim([-1.1*max_spacing max_spacing*1.1]);
        xlabel('\Delta\phi_{tropo} [rad]','fontsize',fontsize)
        ylabel('Height [km]','fontsize',fontsize)
        set(gca,'fontsize',fontsize)

    end

    % computation of the netdelay standard deviation
    netdelays_std = std(netdelays,[],2);
    clear netdelays

    % find the height at which the lowest standard deviation becomes larger
    % than the 0.5 cm threshold.
    fprintf(['Estimating h0 from the height at which the std varies more than ' num2str(0.05) ' cm \n'])
    ix = (netdelays_std<=0.0005);
    h_0_threshold = min(h_delay(ix));
    clear ix
    ix = find(h_delay>h_0_threshold);

    % make some sanity checks
    if h_0_threshold<6000
       fprintf(['h0 is estimated to be ' num2str(h_0_threshold/1000) ' km. This is low abord. \n'])
       if strcmp(error_promp_flag,'y')
           error('myApp:argChk', 'Put h_0_threshold to a lower value. \n')
       else
           fprintf('Try setting h_0_threshold manually \n')
           continue_flag = 0;
           abord_flag = 1;
           break
       end
    else
        % plotting the estimate h0 on the netdelay curve
        if plot_flag==1
            figure(fig_netdelay)
            hold on
            plot(get(gca,'xlim'),[h_0_threshold/1000 h_0_threshold/1000],'k--')

            figure(fig_netphasedelay)
            hold on
            plot(get(gca,'xlim'),[h_0_threshold/1000 h_0_threshold/1000],'k--')

            % saving of the figure when requested
            if save_fig==1
                fig_save_name = ['figures' filesep 'NetDelay_curves.eps'];
                set(fig_netdelay,'PaperPositionMode','auto')
                print(fig_netdelay,'-depsc','-r150',fig_save_name)
                clear fig_save_name

                fig_save_name = ['figures' filesep 'NetPhaseDelay_curves.eps'];
                set(fig_netphasedelay,'PaperPositionMode','auto')
                print(fig_netphasedelay,'-depsc','-r150',fig_save_name)
                clear fig_save_name
            end
        end

        fprintf(['h0 is estimated to be ' num2str(h_0_threshold/1000) ' km. \n'])

        % removing all the other variables
        h_delay(ix:end)=[];
        delay_matrix(ix:end,:)=[];
        phase_delay_matrix(ix:end,:)=[];

        % setting the start of the integration to h_0
        delay_matrix = delay_matrix - repmat(delay_matrix(end,:),length(h_delay),1);
        phase_delay_matrix = phase_delay_matrix - repmat(phase_delay_matrix(end,:),length(h_delay),1);

        if plot_flag==1
            figure2 = figure('name','Delay with height');
            plot(phase_delay_matrix,h_delay./1000,'-','color',[0.7 0.7 0.7] )
            xlabel('LOS phase delay [rad]','fontsize',fontsize)
            ylabel('Height [km]','fontsize',fontsize)
            set(gca,'fontsize',fontsize)
            % saving of the figure when requested
            if save_fig==1
                fig_save_name = ['figures' filesep 'Delay_curves.eps'];
                set(figure2,'PaperPositionMode','auto')
                print(figure2,'-depsc','-r150',fig_save_name)
                clear fig_save_name
            end

        end
    end
end

%% Computing mean delay
fprintf('Computing mean delay \n')
delay_mean = mean(delay_matrix,2);                                        % given in [m]
phase_delay_mean = mean(phase_delay_matrix,2);
if plot_flag==1
    % plot mean delay on top of the existing figure
    figure(figure2)
    hold on
    plot(phase_delay_mean,h_delay./1000,'k-','linewidth',2)
    Ax1 = gca;
    ylimits = get(Ax1,'ylim');
    set(Ax1,'ylim',[0 ylimits(2)])
    % adding a second axis
    box on
    set(gca,'fontsize',fontsize)
    title({'LOS delay [m]',' '},'fontsize',fontsize)
    Ax2 = axes('Position',get(Ax1,'Position'),'XAxisLocation','top');
    set(Ax2,'ylim',get(Ax1,'ylim'),'ytick',[],'yticklabel','');
    xlims = get(Ax1,'xlim');
    xtick_new = get(Ax1,'xtick')./xlims(2);
    new_x_label = get(Ax1,'xtick').*lambda./(4*pi);
    xtick_labels = round((new_x_label)*100)./100;
    xtick_labels = num2str(xtick_labels');
    set(Ax2,'xtick',xtick_new,'xticklabel',xtick_labels);
    set(Ax2,'color','none');
    set(Ax2,'fontsize',fontsize)
    % saving of the figure when requested
    if save_fig==1
        fig_save_name = ['figures' filesep  'Delay_curves.eps'];
        set(figure2,'PaperPositionMode','auto')
        print(figure2,'-depsc','-r150',fig_save_name)
        clear fig_save_name
    end
    clear figure2

    % an individual figure for the mean delay
    figure4 = figure('name','Mean delay with height');
    xlabel('LOS phase delay [rad]','fontsize',fontsize)
    ylabel('Height [km]','fontsize',fontsize)
    set(gca,'fontsize',fontsize)
    hold on
    plot(phase_delay_mean,h_delay./1000,'k-','linewidth',2)
    set(gca,'fontsize',fontsize)
end
save([save_name '.mat'],'delay_mean','h_delay','n_soundings','delay_matrix')

%% estimating the power law directly from the log log plot
% computing the loglog powerlaw

h_log = log10(h_0_threshold-h_delay);
phi_delay_log = log10(phase_delay_mean);
delay_log = log10(delay_mean);
A = [h_log ones(size(h_log))];

% based on the full height range
ix_all = find(h_log<=0.001);
A_all = A;
A_all(ix_all,:)=[];
delay_log_temp_all = delay_log;
delay_log_temp_all(ix_all,:)=[];
phase_delay_log_temp_all = phi_delay_log;
phase_delay_log_temp_all(ix_all,:)=[];
coeff_all = inv(A_all'*A_all)*A_all'*delay_log_temp_all;
coeff_all_phase = inv(A_all'*A_all)*A_all'*phase_delay_log_temp_all;
alpha_log_all = coeff_all(1);
K_log_all = exp(coeff_all(2));

% based on the height range till hc_threshold
ix_hc = find(h_log<=0.001 |  h_delay>h_thres_hc);
A_hc = A;
A_hc(ix_hc,:)=[];
delay_log_temp_hc = delay_log;
delay_log_temp_hc(ix_hc,:)=[];
phase_delay_log_temp_hc = phi_delay_log;
phase_delay_log_temp_hc(ix_hc,:)=[];
coeff_hc = inv(A_hc'*A_hc)*A_hc'*delay_log_temp_hc;
coeff_hc_phase = inv(A_hc'*A_hc)*A_hc'*phase_delay_log_temp_hc;
alpha_log_hc = coeff_hc(1);
K_log_hc = exp(coeff_hc(2));

if plot_flag==1
    xlimits = [log10(h_0_threshold-max_height) log10(h_0_threshold+1)];
    % an individual figure for the mean delay
    figure4_1 = figure('name','Log-Log plot of the Mean delay with height');
    ylabel('Log(tropopsheric LOS delay [m])','fontsize',fontsize)
    xlabel('Height [km])','fontsize',fontsize)
    set(gca,'fontsize',fontsize)
    hold on
    plot(h_log,delay_log,'k-','linewidth',2)
    hold on
    plot(A_hc(:,1),A_hc*coeff_hc,'r--','linewidth',2)
    set(gca,'fontsize',fontsize)
    % plotting the powerlaw based on the lower height range
    hold on
    xlim(xlimits)
    ylimits = get(gca,'ylim');
    legend('Mean delay','Power law',2)
    legend boxoff
    set(gca,'xtick',[get(gca,'xtick') log10(h_0_threshold)],'xticklabel',num2str([get(gca,'xtick') log10(h_0_threshold)]'))

    % updating the lables to be heights
    temp_loc = get(gca,'xtick');
    temp = round((h_0_threshold-10.^temp_loc)/1000*100)/100;
    temp_str = num2str(temp');
    set(gca,'xtick',temp_loc,'xticklabel',temp_str)
    box on

    if save_fig==1
        fig_save_name = ['figures' filesep 'loglog_mean_delay.eps'];
        set(figure4_1,'PaperPositionMode','auto')
        print(figure4_1,'-depsc','-r150',fig_save_name)
        clear fig_save_name
    end

    % an individual figure for the mean delay
    figure4_2 = figure('name','Log-Log plot of the Mean phase delay with height');
    ylabel('Log(\phi_{tropo} [rad])','fontsize',fontsize)
    xlabel('Height [km])','fontsize',fontsize)
    set(gca,'fontsize',fontsize)
    hold on
    plot(h_log,phi_delay_log,'k-','linewidth',2)
    hold on
    plot(A_hc(:,1),A_hc*coeff_hc_phase,'r--','linewidth',2)
    set(gca,'fontsize',fontsize)
    % plotting the powerlaw based on the lower height range
    hold on
    xlim(xlimits)
    ylimits = get(gca,'ylim');
    legend('Mean delay','Power law',2)
    legend boxoff
    set(gca,'xtick',[get(gca,'xtick') log10(h_0_threshold)],'xticklabel',num2str([get(gca,'xtick') log10(h_0_threshold)]'))

    % updating the lables to be heights
    temp_loc = get(gca,'xtick');
    temp = round((h_0_threshold-10.^temp_loc)/1000*100)/100;
    temp_str = num2str(temp');
    set(gca,'xtick',temp_loc,'xticklabel',temp_str)
    box on

    if save_fig==1
        fig_save_name = ['figures' filesep 'loglog_mean_phasedelay.eps'];
        set(figure4_2,'PaperPositionMode','auto')
        print(figure4_2,'-depsc','-r150',fig_save_name)
        clear fig_save_name
    end

end

save([save_name '.mat'],'-append','alpha_log_all','alpha_log_hc')
clear ix_all A_all delay_log_temp_all coeff_all ix_hc A_hc delay_log_temp_hc coeff_hc delay_log phi_delay_log h_log

%% evaluating the performance of the estimation by plotting the estimated powerlaw wrt to the mean delay from sounding
if plot_flag==1
    % powerlaw estimated from log log plot full height range
    delay_log_all = K_log_all*(h_0_threshold-h_delay).^alpha_log_all;

    % powerlaw delay from log log plot estimated from the 0 till hc height range
    delay_log_hc = K_log_hc*(h_0_threshold-h_delay).^alpha_log_hc;
    % plot mean delay on top of the existing figure
    hfig_powerlawdelay = figure('name','Fitted power law delay to mean delay curve');
    plot(delay_mean,h_delay./1000,'k-','linewidth',2)
    hold on
    plot(delay_log_hc,h_delay./1000,'r-','linewidth',2)
    hold on
    xlimits = get(gca,'xlim');
    plot([xlimits(1) xlimits(2)],[max_height max_height]./1000,'k--')
    legend('Mean delay','Fitted power law','Max height of the region')
    ylabel('Height [km]','fontsize',fontsize)
    xlabel('LOS delay [m]','fontsize',fontsize)
    set(gca,'fontsize',fontsize)
    box on

    if save_fig==1
        fig_save_name = ['figures' filesep  'MeanDelay_PowerLawDelay.eps'];
        set(hfig_powerlawdelay,'PaperPositionMode','auto')
        print(hfig_powerlawdelay,'-depsc','-r150',fig_save_name)
        clear fig_save_name
    end
    clear hfig_powerlawdelay

end

continue_flag=0;
abord_flag=0;

end

% The while loop was existed because of an error statement.
if abord_flag == 1
	alpha_log_all = NaN;
	alpha_log_hc = NaN;
	n_soundings=NaN;
    h_0_threshold = NaN;
	fprintf('Early termination \n')
else
    h_0_threshold = h_0_threshold/1000;       % put threshold to [km]
	fprintf('DONE \n')
end

cd(curdir)

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
img = which('keuken_2014_enhanced_for_underlay.img');

% sup colliculus
xyz = [5.1767  -32.3069   -4.2592
-5.1767  -32.3069  -4.2592];

myradius = 3.5;

cluster_orthviews();

cl = sphere_roi_tool_2008(img, myradius, xyz, 'useexisting');

r = cluster2region(cl);
obj = region2imagevec(r);

% Mask with brainstem
bstem = fmri_data(which('brainstem_mask_tight_2018.img'));
obj = apply_mask(obj, bstem);
r = region(obj);
orthviews(r)

r(1).shorttitle = 'R_SC';
r(2).shorttitle = 'L_SC';

sc_regions = r;
sc_obj = obj;

%% inf coll

xyz = [4.8186  -37.2304   -9.9917
      -4.8186  -37.2304   -9.9917];

  myradius = 3;

%cluster_orthviews();

cl = sphere_roi_tool_2008(img, myradius, xyz, 'useexisting');

r = cluster2region(cl);
obj = region2imagevec(r);

% Mask with brainstem
bstem = fmri_data(which('brainstem_mask_tight_2018.img'));
obj = apply_mask(obj, bstem);
r = region(obj);
orthviews(r)

r(1).shorttitle = 'R_IC';
r(2).shorttitle = 'L_IC';

ic_regions = r;
ic_obj = obj;

%% Save

cd('/Users/tor/Documents/Code_Repositories/CanlabCore/CanlabCore/canlab_canonical_brains/Thalamus_brainstem_ROIs_surfaces')
savename = 'sup_inf_colliculus_roi_2018_tor';

pag_regions = r;
pag_obj = obj;

save(savename, 'sc_regions', 'sc_obj', 'ic_regions', 'ic_obj');

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [x,fval,exitflag,info] = opti_mkltrnls(fun,grad,x0,ydata,lb,ub,opts)
%OPTI_MKLTRNLS Solve a NLS using MKLTRNLS (Intel MKL Trust Region Solver)
%
%   min sum[ (F(x) - ydata)^2 ]       subject to:   lb <= x <= ub
%    x
%
%   x = opti_mkltrnls(fun,grad,x0,ydata) solves a NLS where fun is the
%   fitting function. grad is an optional gradient of the fitting function
%   and x0 is a starting guess. ydata is the data to fit the function to.
%
%   x = opti_mkltrnls(fun,grad,x0,ydata,lb,ub) solves subject to decision
%   variables bounds lb <= x <= ub. Avoid Infinite bounds.
%
%   x = opti_mkltrnls(fun,...,ub,opts) uses opts to pass optiset options to
%   the solver.
%
%   [x,fval,exitflag,info] = opti_mkltrnls(...) returns the objective value
%   at the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR DTRNLS, DTRNLSBC and DJACOBI

%   Copyright (C) 2011 Jonathan Currie (IPL)

if(nargin < 7), opts = optiset; end
if(nargin < 6), ub = []; end
if(nargin < 5), lb = []; end
if(nargin < 4), error('MKLTRNLS requires at least 4 arguments'); end

%Setup display level
opts.display = dispLevel(opts.display);
opts.optiver = optiver;

%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('MKLTRNLS requires an initial guess, x0!');
end

%Addin mkltrnls settings if specified
if(isfield(opts,'solverOpts') && ~isempty(opts.solverOpts))
    mopts = mkltrnlsset(opts.solverOpts);
else
    mopts = [];
end
%Add OPTI Options
mopts.maxiter = opts.maxiter;
mopts.maxtime = opts.maxtime;
mopts.display = dispLevel(opts.display);
mopts.tolrfun = opts.tolrfun;
mopts.tolafun = opts.tolafun;
mopts.iterfun = opts.iterfun;

t = tic;
%MKLTRNLS requires x0 within bounds
x0 = movex0(lb,ub,x0);
% Run MKLTRNLS
[x, fval, exitflag, iter, feval, geval] = mkltrnls(fun,grad,x0,ydata,lb,ub,mopts);

%Collect Results
info.Iterations = iter;
info.FuncEvals = feval;
info.GradEvals = geval;
info.Time = toc(t);
info.Algorithm = 'MKLTRNLS: Intel MKL Trust Region NLS';

switch(exitflag)
    case 1
        info.Status = 'Optimal';
    case 0
        info.Status = 'Exceeded Iterations';
    case -1
        info.Status = 'Infeasible / Could not Converge';
    case -2
        info.Status = 'Singular / Error';
    case -5
        info.Status = 'User Exit';
    otherwise
        info.Status = 'MKLTRNLS Error';
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function gamma_values_test ( )

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

  n_data = 0;

  while ( 1 )

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

    if ( n_data == 0 )
      break
    end

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

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = imagesc(DEM,varargin)

%IMAGESC Scale data in GRIDobj and display as image object
%
% Syntax
%
%     h = imagesc(DEM,varargin)
%
% Description
%
%     see imagesc
%
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 8. January, 2013

[x,y] = refmat2XY(DEM.refmat,DEM.size);
ht = imagesc(x,y,DEM.Z,varargin{:});

axis xy
axis image

if nargout == 1;
    h = ht;
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function savebinstack(bin, fname)
% Save an binary stack to a file, as an RGB Image.
%
%  If file name contains '??', then image is saved into a serie
%  of files, with increasing index.
%
%   See also:
%   savestack, imread
%

%   ---------
%   author: David Legland, david.legland(at)inrae.fr
%   INRA - Cepia Software Platform
%   created the 10/09/2003.
%   http://www.pfl-cepia.inra.fr/index.php?page=slicer

dim = size(bin);

if contains(fname, '??')
    % save a serie file, one file per slice of images
    disp('save slices - not yet implemented');

else
    % save one file containing all slices of image
    disp('save a stack');
    if length(dim)==3
        img = uint8(zeros([dim(1) dim(2) 3]));
        for i=1:dim(3)
            img(:,:,1) = uint8(bin(:,:,i)*255);
            img(:,:,2) = uint8(bin(:,:,i)*255);
            img(:,:,3) = uint8(bin(:,:,i)*255);
            imwrite(img, fname, 'tif', 'Compression', 'none', ...
                'WriteMode', 'append');
        end
    else
        for i=1:dim(4)
            imwrite(bin(:,:,:,i), fname, 'tif', 'Compression', 'none', ...
                'WriteMode', 'append');
        end
    end
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [csd,Hz] = spm_ccf2csd(ccf,Hz)
% Converts cross covariance function to cross spectral density
% FORMAT [csd,Hz] = spm_ccf2csd(ccf,Hz)
%
% ccf  (N,:,:)          - cross covariance functions
% Hz   (n x 1)          - vector of frequencies (Hz)
%
% csd  (n,:,:)          - cross spectral density (cf, mar.P)
%
% See also:
%  spm_ccf2csd.m, spm_ccf2mar, spm_csd2ccf.m, spm_csd2mar.m, spm_mar2csd.m,
%  spm_csd2coh.m and spm_Q
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_ccf2csd.m 5895 2014-02-26 14:28:23Z karl $

% unpack cells
%--------------------------------------------------------------------------
if iscell(ccf)
    for i = 1:length(ccf)
       csd{i}    = spm_ccf2csd(ccf{i},Hz);
    end
    return
end

% unpack time bins (for time-frequency responses)
%--------------------------------------------------------------------------
if ndims(ccf) == 4
    for i = 1:size(ccf,1)
       csd(i,:,:,:) = spm_ccf2csd(squeeze(ccf(i,:,:,:)),Hz);
    end
    return
end

% Frequencies
%--------------------------------------------------------------------------
ds    = Hz(2) - Hz(1);
gi    = ceil(Hz/ds) + 1;
N     = Hz(end);

% Fourier transform cross-spectral density
%==========================================================================
for i = 1:size(ccf,2)
    if ismatrix(ccf)
        g          = fft(ccf(:,i));
        csd(:,i)   = g(gi)/N/ds;
    else
        for j = 1:size(ccf,3)
            g          = fft(ccf(:,i));
            csd(:,i,j) = g(gi)/N/ds;
        end
    end
end

author: spm


text:
function [X, W, mu_X] = prewhiten(X)
%PREWHITEN Performs prewhitening of a dataset X
%
%   [X, W, mu_X] = prewhiten(X)
%
% Performs prewhitening of the dataset X. Prewhitening concentrates the main
% variance in the data in a relatively small number of dimensions, and
% removes all first-order structure from the data. In other words, after
% the prewhitening, the covariance matrix of the data is the identity
% matrix. The function returns the subtracted data mean in mu_X, and the
% applied linear mapping in W.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    welcome;

    % Compute and apply the ZCA mapping
    mu_X = mean(X, 1);
    X = bsxfun(@minus, X, mu_X);
    mappedX = X / sqrtm(cov(X));
    if nargout > 1
        W = X \ mappedX;
    end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function rgb = yellow

% returns a predefined color as [red green blue] values
%
% skin_surface         = [255 213 119]/255;
% outer_skull_surface  = [140  85  85]/255;
% inner_skull_surface  = [202 100 100]/255;
% cortex = [255 213 119]/255;
% black  = [0   0   0  ]/255;
% white  = [255 255 255]/255;
% red    = [255 0   0  ]/255;
% green  = [0   192 0  ]/255;
% blue   = [0   0   255]/255;
% yellow = [255 255 0  ]/255;
% cortex_light = [199 194 169]/255;
% cortex_dark  = [100 97 85]/255;

rgb = [255 255 0  ]/255;


text:
function desc = calHsvHist(im, seg, numRegion)

    binNum = 50;
    binVal = 0:1/(binNum):1;
    desc = zeros([numRegion binNum*3]);

    cnt = 0;
    ind={};
    for iReg=1:numRegion
        ind{iReg} = seg(:)==iReg;
    end

    for ch=1:3
        for bin=1:binNum
            cnt = cnt + 1;
            I = im(:,:,ch);
            I = ( (I>=binVal(bin)) & (I<binVal(bin+1)) );

            for iReg=1:numRegion

                desc(iReg, cnt) = sum(I(ind{iReg}))/sum(ind{iReg});
            end
        end
    end

%     tmp = sum(desc, 2);
%     desc = (desc ./ repmat(tmp(:), [1 size(desc,2)]))*3;
end

author: kittenish

repo: Image-Shadow-Detection-and-Removal


text:
function b = r8ci_mxv ( n, a, x )

%*****************************************************************************80
%
%% R8CI_MXV multiplies a R8CI matrix times a vector.
%
%  Discussion:
%
%    The R8CI storage format is used for an N by N circulant matrix.
%    An N by N circulant matrix A has the property that the entries on
%    row I appear again on row I+1, shifted one position to the right,
%    with the final entry of row I appearing as the first of row I+1.
%    The format simply records the first row of the matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N), the R8CI matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the product A * x.
%
  for i = 1 : n
    b(i) = sum ( a(n+2-i:n) .* x(1:i-1) ) + sum ( a(1:n+1-i) .* x(i:n) );
  end

  return
end


text:
function collectObjects(fig)
% Helper function for browseFigure.m

% Get handles of already present objects
objHandle = get(fig, 'UserData');
% Find handles of new objects
newObjects  = findobj(fig,'Visible','on','-not','Type','figure');
% Make new objects invisible
set(newObjects, 'Visible', 'off');
% Add new objects handles
objHandle{end+1} = newObjects;
% Pass object handles to figure's UserData
set(fig, 'UserData', objHandle);
axHandle = findobj(newObjects, 'Type', 'axes');
% Make handles invisible to avoid overwriting of the axes
set(axHandle, 'HandleVisibility', 'off');

author: Sable

repo: mcbench-benchmarks


text:
function order = hexagon_unit_set ( rule )

%*****************************************************************************80
%
%% HEXAGON_UNIT_SIZE sizes a quadrature rule inside the unit hexagon in 2D.
%
%  Integration region:
%
%    The definition is given in terms of THETA, the angle in degrees of the
%    vector (X,Y).  The following six conditions apply, respectively,
%    between the bracketing values of THETA of 0, 60, 120, 180, 240,
%    300, and 360.
%
%                              0 <= Y <= - SQRT(3) * X + SQRT(3)
%                              0 <= Y <=                 SQRT(3)/2
%                              0 <= Y <=   SQRT(3) * X + SQRT(3)
%      - SQRT(3) * X - SQRT(3)   <= Y <= 0
%                    - SQRT(3)/2 <= Y <= 0
%        SQRT(3) * X - SQRT(3)   <= Y <= 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Abramowitz and Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964.
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, integer RULE, the rule desired.
%      1, 1 point,  degree 1;
%      2, 4 points, degree 3;
%      3, 7 points, degree 3;
%      4, 7 points, degree 5;
%
%    Output, integer ORDER, the order of the desired rule.
%    If RULE is not legal, then ORDER is returned as -1.
%
  if ( rule == 1 )

    order = 1;
%
%  Stroud rule H2:3-1.
%
  elseif ( rule == 2 )

    order = 4;
%
%  Stroud rule H2:3-2.
%
  elseif ( rule == 3 )

    order = 7;
%
%  Stroud rule H2:5-1.
%
  elseif ( rule == 4 )

    order = 7;

  else

    order = -1;

  end

  return
end


text:
function clipMovie(indir,outdir,startFrame,endFrame,varargin)
% function clipMovie(indir,outdir,startFrame,endFrame,...)

[moviename,trxfilename,removeInvalid] = myparse(varargin,...
  'moviename','movie.ufmf','trxfilename','trx.mat','removeInvalid',false);

if ~exist(outdir,'dir')
  mkdir(outdir,'perframe');
end

if ~exist(fullfile(outdir,moviename),'file')
  copyfile(fullfile(indir,moviename),fullfile(outdir,moviename));
end

if ~strcmp(trxfilename(end-3:end),'.mat')
  trxfilename = [trxfilename '.mat'];
end

Q = load(fullfile(indir,trxfilename));
selflds = {'x','y','theta','a','b','timestamps',...
  'x_mm','y_mm','theta_mm','a_mm','b_mm'};
pfstart = [];
pfend = [];
for ndx = numel(Q.trx):-1:1
  if Q.trx(ndx).firstframe > endFrame || Q.trx(ndx).endframe < startFrame
    if removeInvalid
      Q.trx(ndx) = [];
    else
      Q.trx(ndx).firstframe = startFrame;
      Q.trx(ndx).endframe = startFrame;
      for fnum = 1:numel(selflds)
        Q.trx(ndx).(selflds{fnum}) = Q.trx(ndx).(selflds{fnum})(1);
      end
      Q.trx(ndx).dt = 1;
      Q.trx(ndx).nframes = 1;
      pfstart(ndx) = 1; pfend(ndx) = 1;
    end
  else
    curStart = max(Q.trx(ndx).firstframe,startFrame);
    curEnd = min(Q.trx(ndx).endframe,endFrame);
    for fnum = 1:numel(selflds)
      Q.trx(ndx).(selflds{fnum}) = Q.trx(ndx).(selflds{fnum})( (curStart:curEnd)-Q.trx(ndx).firstframe+1);
    end
    Q.trx(ndx).dt = Q.trx(ndx).dt( (curStart:curEnd-1)-Q.trx(ndx).firstframe+1);
    Q.trx(ndx).endframe = curEnd;
    Q.trx(ndx).nframes = curEnd-curStart+1;
    Q.trx(ndx).firstframe = curStart;
    Q.trx(ndx).off = -curStart+1;
    pfstart(ndx) = curStart-Q.trx(ndx).firstframe+1;
    pfend(ndx) = curEnd-Q.trx(ndx).firstframe+1;
  end
end
if ~isempty(Q.timestamps),
  Q.timestamps = Q.timestamps(startFrame:endFrame);
end
save(fullfile(outdir,trxfilename),'-struct','Q');

if ~exist(fullfile(outdir,'perframe'),'dir'),
  mkdir(fullfile(outdir,'perframe'));
end

if ~removeInvalid
  dd = dir(fullfile(indir,'perframe','*.mat'));
  for pf = 1:numel(dd)
    Q = load(fullfile(indir,'perframe',dd(pf).name));
    for ndx = 1:numel(Q.data)
      if isnan( pfstart(ndx))
        Q.data{ndx} = [];
      else
        Q.data{ndx} = Q.data{ndx}(pfstart(ndx):pfend(end));
      end
    end
    save(fullfile(outdir,'perframe',dd(pf).name),'-struct','Q');

  end
end

author: kristinbranson

repo: JAABA


text:
function sensor_out=la_sen(acc, gyro, gyro_der, sen_la, sen_or, sen_typ)
nsen=size(sen_typ,1);
sensor_out=zeros(nsen,1);
for in=1:nsen;
    %lever arm
    Csm=sen_or(:,:,in);
    la=sen_la(:,in);

    %Acceleration
    acc_sen=Csm'*(acc+cross(gyro,cross(gyro,la))+cross(gyro_der,la));

    %Rotation rate
    gyro_sen=Csm'*gyro;

    if (sen_typ(in)==1)
        sensor_out(in)=acc_sen(1);
    elseif (sen_typ(in)==2)
        sensor_out(in)=gyro_sen(1);
    end
end

repo: ekfmonoslam


text:
function [fv, normals] = MimCreateSurfaceFromSegmentation(segmentation, smoothing_size, small_structures, label, coordinate_system, template_image, limit_to_one_component_per_index, minimum_component_volume_mm3, reporting)
    % MimCreateSurfaceFromSegmentation. Creates a surface mesh from a segmentation volume.
    %
    %     Inputs
    %     ------
    %
    %     filepath, filename - path and filename for STL file
    %
    %     segmentation - a binary 3D PTKImage containing 1s for the voxels to
    %         visualise
    %
    %     smoothing_size - the amount of smoothing to perform. Higher
    %         smoothing makes the image look better but removes small
    %         structures such as airways and vessels
    %
    %     small_structures - set to true for improved visualisation of
    %         narrow structures such as airways and vessels
    %
    %     label - Only voxels of this colour are considered
    %
    %     coordinate_system  a MimCoordinateSystem enumeration
    %         specifying the coordinate system to use
    %
    %     template_image  A PTKImage providing voxel size and image size
    %         parameters, which may be required depending on the coordinate
    %         system
    %
    %     limit_to_one_component_per_index - set this flag to true if each colour
    %         index in the image represents only a single, large object.
    %         This will prevent the appearance of orphaned 'island' components caused by the
    %         image smoothing and surface rendering approximations.
    %
    %     minimum_component_volume_mm3 - if two or more separate components in the
    %         image could share the same colour index, but you want to prevent the
    %         appearance of orphaned 'island' components caused by the
    %         image smoothing and surface rendering approximations, then set this to
    %         an appropriate minimum volume value. Any components smaller than this
    %         value will not be rendered.
    %
    %     reporting - an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if ~isa(segmentation, 'PTKImage')
        error('Requires a PTKImage as input');
    end

    % Isolate this colour component
    sub_seg = segmentation.BlankCopy;
    raw_image = segmentation.GetMappedRawImage;
    sub_seg.ChangeRawImage(raw_image == label);

    sub_seg.CropToFit;

    % Perform a closing operation and filter the image to create a smoother appearance
    if small_structures
        morph_size = 1;
        threshold = 0.2;
    else
        morph_size = 3;
        threshold = 0.5;
    end

    required_padding_mm = 2*threshold + morph_size;
    required_padding = ceil(required_padding_mm/min(segmentation.VoxelSize)) + 1;
    sub_seg.AddBorder(required_padding);

    if morph_size > 0
        sub_seg.BinaryMorph(@imclose, morph_size);
    end

    if smoothing_size > 0
        sub_seg = MimGaussianFilter(sub_seg, smoothing_size);
    end

    if limit_to_one_component_per_index

        % The smoothing may split the segmentation, so we consider only the main component,
        % and remove any other components by setting their value to zero
        surviving_components = sub_seg.RawImage >= threshold;
        surviving_components = MimImageUtilities.GetLargestConnectedComponent(surviving_components);
    else

        cc = bwconncomp(sub_seg.RawImage >= threshold);
        num_pixels = cellfun(@numel, cc.PixelIdxList);
        voxel_threshold = minimum_component_volume_mm3/prod(segmentation.VoxelSize);
        larger_elements = num_pixels > voxel_threshold;
        surviving_components = false(sub_seg.ImageSize);
        for index = 1 : numel(larger_elements)
            if larger_elements(index)
                surviving_components(cc.PixelIdxList{index}) = true;
            end
        end
    end

    sub_seg_raw = sub_seg.RawImage >= threshold;
    sub_seg_raw(~surviving_components) = 0;
    sub_seg.ChangeRawImage(sub_seg_raw);

    [xc, yc, zc] = sub_seg.GetPTKCoordinates;
    [xc, yc, zc] = MimImageCoordinateUtilities.ConvertFromPTKCoordinatesCoordwise(xc, yc, zc, coordinate_system, template_image);

    fv = isosurface(xc, yc, zc, sub_seg.RawImage, threshold);

    if nargout > 1
        normals = isonormals(xc, yc, zc, sub_seg.RawImage, fv.vertices);
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% ----------------------------------------------------------------------------
% function hfssSetColor(fid, Objects, Color)
%
% Description :
% -------------
% Creates the VBScript necessary to set an object's color in HFSS.
%
% Parameters :
% ------------
% fid     - file identifier of the HFSS script file.
% Objects - cell array with name(s) of the object(s) whose color require(s)
%           to be changed.
% Color   - [3x1 vector] represents the [R, G, B] components of the color.
%
% Note :
% ------
%
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ...
% hfssSetColor(fid, 'Substrate', [0, 64, 0]);
% hfssSetColor(fid, {'Sub1', 'Sub2'}, [0, 64, 0]);
%
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% CHANGELOG
%
% ??-???-????: *Initial release.
% 09-Sep-2020: *It can set a color to several objects simultaneously.
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% This file is part of HFSS-MATLAB-API.
%
% HFSS-MATLAB-API 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.
%
% HFSS-MATLAB-API 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
% Foobar; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place, Suite 330, Boston, MA  02111-1307  USA
%
% Copyright 2004, Vijay Ramasami (rvc@ku.edu)
% ----------------------------------------------------------------------------
function hfssSetColor(fid, Objects, Color)

if (~iscell(Objects))
    Objects = cellstr(Objects);
end

fprintf(fid, '\n');
fprintf(fid, 'oEditor.ChangeProperty _\n');
fprintf(fid, '\tArray("NAME:AllTabs", _\n');
fprintf(fid, '\t\tArray("NAME:Geometry3DAttributeTab", _\n');
fprintf(fid, '\t\t\tArray("NAME:PropServers", ');

for n = 1:numel(Objects)-1
    fprintf(fid, '"%s", ', Objects{n});
end
fprintf(fid, '"%s"), _\n', Objects{end});

fprintf(fid, '\t\t\tArray("NAME:ChangedProps",  _\n');
fprintf(fid, '\t\t\t\tArray("NAME:Color", "R:=", %d, "G:=", %d, "B:=", %d) _\n', ...
        Color(1), Color(2), Color(3));
fprintf(fid, '\t\t\t) _\n');
fprintf(fid, '\t\t) _\n');
fprintf(fid, '\t) \n');

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function [P2,D2]=sortem(P,D)
% this function takes in two matrices P and D, presumably the output
% from Matlab's eig function, and then sorts the columns of P to
% match the sorted columns of D (going from largest to smallest)
%
% EXAMPLE:
%
% D =
%    -90     0     0
%      0   -30     0
%      0     0   -60
% P =
%      1     2     3
%      1     2     3
%      1     2     3
%
% [P,D]=sortem(P,D)
% P =
%      2     3     1
%      2     3     1
%      2     3     1
% D =
%    -30     0     0
%      0   -60     0
%      0     0   -90

D2=diag(sort(diag(D),'descend')); % make diagonal matrix out of sorted diagonal values of input D
[c, ind]=sort(diag(D),'descend'); % store the indices of which columns the sorted eigenvalues come from
P2=P(:,ind); % arrange the columns in this order


text:
report_this_filefun(mfilename('fullpath'));

clf
[X,Y,Z] = meshgrid(gx,gy,gz);
v = [];
i = 0.05
ii = 0.05
for i3 = 1:2:length(gy)
    i3 = gy(i3);
    l = loc(2,:) == i3;
    re3 = reshape(as(l),length(gz),length(gx));
    r = reshape(loc(4,l),length(gz),length(gx));
    l = r > tresh;
    re3(l) = zeros(1,length(find(l)))*nan;
    rect = [i ii 0.15 0.15];
    i = i+0.20;
    if i > 0.9 ; i = 0.05 ; ii = ii + 0.20; end
    axes('position',rect)
    pcolor(gx,-gz,re3)
    caxis([min(as) max(as)])
    colormap(jet)
    title([ 'Lat=' num2str(i3) ' km'],'FontSize',8)
    set(gca,'FontSize',8)
    hold on
    l = a.Latitude > i3-dy/2 & a.Latitude < i3+dy/2;
    plot3(a(l,1),-a(l,7),a(l,7)*0,'k.')
    shading interp
    % axis('off')
    hold on
end

author: CelsoReyes


text:
function area = polySgnArea(x,y)
% area of an orientation with sign depending of orientation

area = 0.5*sum((y(2:end)-y(1:end-1)) .* (x(2:end)+x(1:end-1)));

end

author: mtex-toolbox

repo: mtex

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


text:
function [A1, B1, C1, A2, B2, C2] = sss(a, b, c)
%SSS   gives both solutions to the side-side-side problem, in radians.
%
%   SSS(a, b, c) results in NaNs for those indices where the existence
%   condition |pi - a| - |pi - b| <= |pi - c| <= |pi - a| + |pi -b| is not
%   met.
%
%   See also SSSD.

% Rody P.S. Oldenhuis
% Delft University of Technology
% oldenhuis@gmail.com
%
% Created    : 23/Feb/2009
% Last edited: 30/Nov/2012

    % first solution
    A1 = acos2( (cos(a) - cos(b).*cos(c))./(sin(b).*sin(c)), a);
    B1 = acos2( (cos(b) - cos(a).*cos(c))./(sin(a).*sin(c)), b);
    C1 = acos2( (cos(c) - cos(a).*cos(b))./(sin(a).*sin(b)), c);

    % second solution
    A2 = 2*pi - A1;
    B2 = 2*pi - B1;
    C2 = 2*pi - C1;

    % check constraints
    indices = ( ...
        (abs(pi-a) - abs(pi-b)) > abs(pi-c) | ...
        abs(pi-c) > (abs(pi-a) + abs(pi-b)) );
    A1(indices) = NaN;   B1(indices) = NaN;   C1(indices) = NaN;
    A2(indices) = NaN;   B2(indices) = NaN;   C2(indices) = NaN;

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function spheremap(status)
%SPHEREMAP - Map planar points to unit sphere
%    This function takes points drawn on x-y axis and maps them to a sphere
%    in real time.  Click and drag on the planar area to draw.
%
% Michael Hanchak, Dayton, OH, USA  2000-2011

if nargin <1
    status = 'build';
end

switch status
    case 'build'
        clf;set(gcf,'position',[100 100 500 500]);
        axes('units','normalized','position',[.05 .05 .9 .45],'xlim',[-180 180],...
            'ylim',[-90 90],'tag','plane','xtick',[-180 -135 -90 -45 0 45 90 135 180],...
            'ytick',[-90 -45 0 45 90]);hold on;grid
        line(190,0,'marker','.','linestyle','none','markersize',5,...
            'tag','lineplane','color','black');
        axes('units','normalized','position',[.05 .5 .9 .45],'xlim',[-1 1],...
            'ylim',[-1 1],'zlim',[-1 1],'tag','sphere','visible','off');hold on;

        axis vis3d
        [X,Y,Z] = sphere(24);
        h0 = mesh(X,Y,Z);
        set(h0,'facecolor',[1 1 1],'edgecolor',[0 0 0]);
        line(0,0,0,'marker','.','linestyle','none','markersize',5,...
            'tag','linesphere','color','black');
        view(100,30);
        set(gcf,'windowbuttondownfcn',...
            'spheremap(''mouse'');set(gcf,''windowbuttonmotionfcn'',''spheremap(''''mouse'''')'')',...
            'windowbuttonupfcn','set(gcf,''windowbuttonmotionfcn'','' '')');
        clbk = ['h = findobj(''tag'',''lineplane'');set(h,''xdata'',[],''ydata'',[]);',...
            'h = findobj(''tag'',''linesphere'');set(h,''xdata'',[],''ydata'',[],''zdata'',[]);'];
        uicontrol('style','pushbutton','units','normalized','position',[.02 .85 .1 .05],...
            'callback',clbk,'string','Clear');
        set(gcf,'toolbar','figure')

    case 'mouse'
        whichclick = get(gcf,'selectiontype');
        switch whichclick
            case 'normal'
                temp = get(findobj('tag','plane'),'currentpoint');
                point = temp(1,1:2)*pi/180;
                if point(1) <= pi && point(1) >= -pi && point(2) <= pi/2 &&...
                        point(2) >= -pi/2

                    h = findobj('tag','lineplane');
                    x = get(h,'xdata');
                    y = get(h,'ydata');
                    set(h,'xdata',[x, point(1)*180/pi],'ydata',[y, point(2)*180/pi]);

                    h = findobj('tag','linesphere');
                    x = get(h(1),'xdata');
                    y = get(h(1),'ydata');
                    z = get(h(1),'zdata');
                    set(h,'xdata',[x, cos(point(1))*cos(point(2))],...
                        'ydata',[y, sin(point(1))*cos(point(2))],...
                        'zdata',[z, sin(point(2))]);

                    drawnow
                end

            %case 'alt'
            %    set(gcf,'windowbuttondownfcn','','windowbuttonmotionfcn','');
            %    axes(findobj('tag','sphere'));
            %    rot3d;view(100,30);
        end

end


text:
function f = dec2fix(x, nfracbits, nbits)

% DEC2FIX Convert decimal integer to binary string fixed point.
%
% Usage: F = DEC2FIX(X, NFRACBITS, NBITS)
%
% Converts the signed decimal integer given by X (either a scalar, vector,
% or matrix) to the two's complement representation as a string. If X is a
% vector or matrix, F(i, :) is the representation for X(i) (the shape of X
% is not preserved). Note that many fractional numbers that can be
% represented with a finite number of fractional digits cannot be
% represented by a finite number of fractional bits (specifically
% non-powers-of-two like 0.3), so input NFRACBITS is required to specify
% the precision for the number of fractional bits. Even powers-of-two
% fractional decimal numbers (like 0.5) require this input.
%
% Example:
%     >> dec2fix([2.3 2.4 -2.3 -2.4], 3)
%
%     ans =
%
%     010.010
%     010.011
%     101.110
%     101.101
%
% Inputs:
%   -X: decimal integers to convert to two's complement.
%   -NFRACBITS: number of bits to represent fractional part.
%   -NBITS: total number of bits in the representation including the
%   fractional bits (optional, default is the fewest number of bits
%   necessary to represent the integer portion).
%
% Outputs:
%   -F: fixed point representation of X as a string.
%
% See also: FIX2DEC, TWOS2DEC, DEC2TWOS, DEC2BIN, DEC2HEX, DEC2BASE.

error(nargchk(2, 3, nargin));
x = x(:);
maxx = max(abs(x));
nbits_min = max([nextpow2(maxx + (any(x == maxx))) + 1 + nfracbits, ...
    nfracbits]);

% Default number of bits.
if nargin < 3
    nbits = nbits_min;
elseif nbits < nbits_min
      warning('dec2twos:nbitsTooSmall', ['Minimum number of bits to ' ...
        'represent maximum input x is %i, which is greater than ' ...
        'input nbits = %i. Setting nbits = %i.'], ...
        nbits_min, nbits, nbits_min)
    nbits = nbits_min;
end

% Convert to two's complement string.
f = dec2twos(round(x * 2.^nfracbits), nbits);

% Insert binary point.
f(:, end+1) = '.';
f = f(:, [1:(nbits-nfracbits), nbits+1, (nbits-nfracbits+1):nbits]);


function [dyc, dyt]=hamfilter(X,h,d,ck,fig,timeplot, nameplot)

% local projections (direct  forecast) to  perform detrending
% regressions
% y(t+h)= a*y(t)+b*y(t-1)+c*y(t-2)+...+d*y(t-d)+ e(t+h)
% cycle  is  e(t+h)
% trend  is  hat(a)*y(t)+hat(b)*y(t-1)+hat(c)*y(t-2)+...+hat(d)*y(t-d)

% h= horizon of  the  projection
% d= number  of  lags  used
% ck if=1 constant  if =0 no constant in  the projection

% dyc = estimated cycle
% dyt = estimated trend

[enddT,QQ]=size(X);
yc=zeros(enddT,QQ); yt=zeros(enddT,QQ);
R=[];
T = 1:1:enddT;
time = T;
if nargin > 5
    time =timeplot;
end
if nargin > 6
    titleplot = nameplot;
    if length(nameplot) ~= size(X,2)
        error('nameplot size shold be the same as the column of Y')
    end
else
    for v = 1 : size(X,2)
        eval(['titleplot{'   num2str(v) '} =  ''Var' num2str(v) ''';'])
    end
end

for qq=1:QQ
    yh=squeeze(X(d+h:enddT,qq)); %    independent  variable
    if  ck==1
        R=ones(enddT-d-h+1,1);     % constant
    end

    for  jj=1:d
        r=squeeze(X(d+1-jj:enddT-h+1-jj,qq));
        R=[R r];  % dependent  variables
    end

    yc(d+h:enddT,qq) = yh - R*((R'*R)\(R'*yh)); % cycle
    yt(d+h:enddT,qq) = R*((R'*R)\(R'*yh));        % trend

    if  fig==1
        subplot(2,1,1)
        %plot(time(d+h:enddT),yh(1:enddT-d-h,1),'r', 'linewidth',2);hold  on;
        plot(time(d+h:enddT),X(d+h:enddT,qq),'r', 'linewidth',2);hold  on;
        plot(time(d+h:enddT),yt(d+h:enddT,qq),'k--','linewidth',2);hold off; axis  tight;
        legend('data', 'Hamil trend')
        title(titleplot(qq))
        subplot(2,1,2)
        plot(time(d+h:enddT),yc(d+h:enddT,qq),'b', 'linewidth',2);
        legend('Hamil cycle')
        pause
    end
end
dyc=yc;
dyt=yt;
end


text:
function varargout = subsref(f, index)
%SUBSREF   DISKFUN subsref.
%( )
%   F(TH, R, 'polar') returns the values of the DISKFUN F evaluated at the
%   value(s) (TH, R) in polar coordinates. See CHEBFUN/FEVAL for
%   further details.
%
%   F(:, R) returns a chebfun representing the function F along a radial
%   slice fixed by R, and F(TH, :) returns a chebfun representing F along
%   an angular slice.
%
%   F(:, :) returns F.
%
%   F(X, Y) or F(X,Y, 'cart') returns the values of the DISKFUN F evaluated
%   at the value(s) (X,Y) in Cartesian coordinates.
%   The colon operator for Cartesian coordinates is not supported,
%   except for F(:, :), which just returns F.
%
%   F(c) where c is a complex-valued chebfun evaluates F along the contour
%   parametrized by the real and imaginary parts of c, and returns a
%   chebfun representing the values of F along the contour.
%
%   F.PROP returns the property PROP of F as defined by GET(F, 'PROP').
%
% {}
%   F{S1, S2, S3, S4} restricts F to the domain [S1, S2, S3, S4]. See
%   SEPARABLEAPPROX/RESTRICT for further details. Note that
%   F{[S1,S2, S3, S4]} is not supported due to the behaviour of the MATLAB
%   subsref() command.
%
% See also DISKFUN/FEVAL, DISKFUN/GET and DISKFUN/RESTRICT.

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

idx = index(1).subs;

% Check for the cases that must be explicitly handled here, which are
% 1. f(x, y), x, y are numeric, cartesian flag.
% 2. f(t,r) t, r are numeric, polar flag.
% Pass the other cases off to separableApprox/subsref.

if ( strcmp(index(1).type, '()') )
    x = idx{1};
    if ( length(idx) == 3 ) && (strcmp(idx{3}, 'polar'))
        y = idx{2};
        out = feval(f, x, y, 'polar');
        varargout = { out };

    elseif ( length(idx) == 3 ) && (strcmp(idx{3}, 'cart'))
        y = idx{2};
        out = feval(f, x, y, 'cart');
        varargout = { out };

    else  % Pass to separableApprox.
        out = subsref@separableApprox(f, index);
        varargout = { out };
    %else
      % error('CHEBFUN:DISKFUN:subsref:inputs', ...
               % 'Can only evaluate diskfuns at (X,Y) or (TH,R)');
    end

elseif ( strcmp(index(1).type,'()') )
    if ( numel(idx)==4 )
        % This intentionally fails:
        varargout = { restrict(f,[idx{1}, idx{2}, idx{3}, idx{4}]) };

    else
        error('CHEBFUN:DISKFUN:SUBSREF:restrict',...
              'Restriction domain should be given by four corners');
    end

else
    [varargout{1:nargout}] = subsref@separableApprox(f, index);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function p = cleanuppermodel(p)

% We might have created a bilinear model from an original polynomial model.
% We should use the original model when we solve the upper bound problem.
p_bilinear = p;
p = p.originalModel;

p = removeCuts(p);

n_start = length(p.c);

% Quadratic mode, and quadratic aware solver?
bilinear_variables = find(p.variabletype == 1 | p.variabletype == 2);
if ~isempty(bilinear_variables)
    used_in_c = find(p.c);
    quadraticterms = used_in_c(find(ismember(used_in_c,bilinear_variables)));
    if ~isempty(quadraticterms) & p.solver.uppersolver.objective.quadratic.nonconvex
        usedinquadratic = zeros(1,length(p.c));
        for i = 1:length(quadraticterms)
            Qij = p.c(quadraticterms(i));
            power_index = find(p.monomtable(quadraticterms(i),:));
            if length(power_index) == 1
                p.Q(power_index,power_index) = Qij;
            else
                p.Q(power_index(1),power_index(2)) = Qij/2;
                p.Q(power_index(2),power_index(1)) = Qij/2;
            end
            p.c(quadraticterms(i)) = 0;
        end
    end
end

p.lb = p_bilinear.lb(1:length(p.c));
p.ub = p_bilinear.ub(1:length(p.c));
p.bilinears = [];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function res = fieldnames(this, varargin)
% Returns names of the fields in .other
% FORMAT res = fieldnames(this)
%
% An overloaded function...
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: fieldnames.m 2720 2009-02-09 19:50:46Z vladimir $

res = fieldnames(this.other, varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function neigh = grNeighborEdges(edges, node)
%GRNEIGHBOREDGES Find adjacent edges of a given node.
%
%   Deprecated, use 'grAdjacentEdges' instead.
%
%   NEIGHS = grNeighborEdges(EDGES, NODE)
%   EDGES  : the complete edges list
%   NODE   : index of the node
%   NEIGHS : the indices of edges containing the node index
%
%   See also
%     grAdjacentEdges

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2003-08-13
% Copyright 2003 INRA - TPV URPOI - BIA IMASTE

warning('MatGeom:graphs:deprecated', ...
    'function grNeighborEdges is obsolete, use grAdjacentEdges instead');

neigh = find(edges(:,1) == node | edges(:,2) == node);

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function setFigureDefaults()
	% Sets preferred figure defaults.
	% Biafra Ahanonu
	% started: 2013.10.29
	% inputs
		%
	% outputs
		%

	% changelog
		%
	% TODO
		%

	ciapkg.view.setFigureDefaults();
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
close all; clear all; clc;

wave = spx.wavelet.daubechies.wavelet_function(4, 8, 1024, 4);
subplot(221);
t = (1:1024)./1024;
plot(t(300:800),wave(300:800));
ylim([-.25 +.25]);
title(' D4 Wavelet ');

wave = spx.wavelet.daubechies.scaling_function(4, 8, 1024, 4);
subplot(222);
t = (1:1024)./1024;
plot(t(300:800),wave(300:800));
ylim([-.25 +.25]);
title(' D4 Scaling function ');

author: indigits

repo: sparse-plex


text:
%% Data Sturcture for Degree of Freedom (DOF)
%
% For various elements, dofs could be associated to vertices, edges, or
% faces with apporipate orientation. In this doc, we collect data
% structures for commonly used elements.
%
% * <dofedgedoc.html Dof on edges: 2-D>
% * <dof3edgedoc.html Dof on edges: 3-D>
% * <dofP2doc.html Dof of P2 elements>

repo: Aero-matlab


text:
function varargout = ldivide(varargin)
%.\   Pointwise CHEBFUN2 left array divide.
%   F.\G if G is a CHEBFUN2 and F is a double this returns (1/F)*G
%
%   F.\G if G is a double and F is a CHEBFUN2 this returns G\F, but this does
%   not work if F becomes numerically close to zero.
%
%   F.\G we do not allow F and G to both be CHEBFUN2 objects.
%
%   F.\G is the same as the command ldivide(F,G)

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

[varargout{1:nargout}] = ldivide@separableApprox(varargin{:});

end


text:
function r = isinf(a)
%ISINF        Array of 1's for inf components
%
%   r = isinf(a)
%

% written  08/07/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    Matlab flaw for empty sparse corrected
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  r = reshape( isinf(a.x(:)) | any(isinf(a.dx),2) , size(a.x) );

repo: OpenMAS


text:
%% t_installSegmentation
%
% Illustrates how to align the inplane volume from an fMRI session to the
% 3D volume anatomy using sample data set <erniePRF>
%
% Dependencies:
%   Remote Data Toolbox
%
% This tutorial is part of a sequence. Run
%   t_initAnatomyFromFreesurfer
%   t_initVistaSession
%   t_alignInplaneToVolume
% prior to running this tutorial.
%
% Summary
%
% - Specify alignment matrix linking inplane anatomy to volume anatomy
% - Save mrSESSION with alignment matrix
%
% Tested 07/21/2016 - MATLAB r2015a, Mac OS 10.11.6
%
%  See also: t_initAnatomyFromFreesurfer t_initVistaSession t_alignInplaneToVolume
%
% Winawer lab (NYU)

%% Start
% Clean start in case we have a vista session open
mrvCleanWorkspace();

% Remember where we are
curdir = pwd();

%% Organize functional data

% Find ernie PRF session in scratch directory
erniePathTemp      = fullfile(vistaRootPath, 'local', 'scratch', 'erniePRF');

if ~exist(erniePathTemp, 'dir')
    help(mfilename)
    error('Please run  pre-requisite tutorials')
end

% Navigate and create a directory
cd(erniePathTemp)

%% Align inplane to t1 and install Gray/white segmentation

% Open a hidden view
vw = initHiddenInplane();

% Segmentation inputs
%   use command line, not dialog
query = false;

%   keep all gray nodes, including those outside functional FOV
keepAllNodes = true;

%   path to class file
filePaths = fullfile('3DAnatomy', 't1_class.nii.gz');

%   number of layers in gray graph along surface (3 layers for 1 mm voxels)
numGrayLayers = 3;

% Do it
installSegmentation(query, keepAllNodes, filePaths, numGrayLayers);

%% Visualize

% open a UI
vw = mrVista;

% Define an ROI that is the entire functional slab
vw = makeGrayROI(vw);
vw = refreshScreen(vw,0);

% add the local variable, vw, to the UI in case you would like to interact
% with the GUI
updateGlobal(vw)

%% Clean up
close(viewGet(vw, 'figure number'));
mrvCleanWorkspace
cd(curdir)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [S,KL,KH] = use_spm_filter(TR,dims,LChoice,HChoice,HParam,varargin)
% :Usage:
% ::
%
%     function [S,KL,KH] = use_spm_filter(TR,dim of filter,LChoice,HChoice,HP filter in s,[LP Gauss len in s])
%
% :Inputs:
%
%   **K{s}.LChoice:**
%        Low-pass  filtering {'hrf' 'Gaussian' 'none'}
%   **K{s}.LParam:**
%        Gaussian parameter in seconds
%   **K{s}.HChoice:**
%        High-pass filtering {'specify' 'none'}
% ..
%    05/22/01 Tor Wager
% ..

K{1}.RT = TR;
K{1}.LChoice = LChoice;
K{1}.HChoice = HChoice;
K{1}.HParam = HParam;
K{1}.row = ones(dims, 1);

if length(varargin) > 0
        K{1}.LParam = varargin{1};
end

KL = []; KH = [];

spmS = spm_filter('set', K);

S = eye(length(K{1}.row));

if ~strcmp(HChoice,'none')

	KH = full(spmS{1}.KH);
	S = S - KH * pinv(KH);

end

if ~strcmp(LChoice,'none')
	KL = full(spmS{1}.KL);
	S = KL * S;         	% lowpass * highpass; hp = I - res forming mtx of S.KH
end

return

function [vargout] = spm_filter(Action,K,Y)
% filter routine
% FORMAT [K] = spm_filter('set',K)
% FORMAT [Y] = spm_filter('apply',K,Y)
%
% Action    - 'set'   fills in filter structure K
% Action    - 'apply' applies K to Y = K*Y
% K         - filter convolution matrix or:
% K{s}      - cell of structs containing session-specific specifications
%
% K{s}.RT       - repeat time in seconds
% K{s}.row      - row of Y constituting session s
% K{s}.LChoice  - Low-pass  filtering {'hrf' 'Gaussian' 'none'}
% K{s}.LParam   - Gaussian parameter in seconds
% K{s}.HChoice  - High-pass filtering {'specify' 'none'}
% K{s}.HParam   - cut-off period in seconds
%
% K{s}.HP       - low frequencies to be removed
% K{s}.LP       - sparse toepltz low-pass convolution matrix
%
% Y         - data matrix
%
% K         - filter structure
% Y         - filtered data K.K*Y
%___________________________________________________________________________
%
% spm_filter implements band pass filtering in an efficient way by
% using explicitly the projector matrix form of the High pass
% component.  spm_filter also configures the filter structure in
% accord with the specification fields if required
%___________________________________________________________________________
% @(#)spm_filter.m	2.4 Karl Friston 99/08/31

% set or apply
%---------------------------------------------------------------------------
switch Action

	case 'set'
	%-------------------------------------------------------------------
	for s = 1:length(K)

		% matrix order
		%-----------------------------------------------------------
		k     = length(K{s}.row);

		% make low pass filter
		%-----------------------------------------------------------
		switch K{s}.LChoice

			case 'none'
			%---------------------------------------------------
			h       = 1;
			d       = 0;

			case 'hrf'
			%---------------------------------------------------
			h       = spm_hrf(K{s}.RT);
			h       = [h; zeros(size(h))];
			g       = abs(fft(h));
			h       = real(ifft(g));
			h       = fftshift(h)';
			n       = length(h);
			d       = [1:n] - n/2 - 1;

			case 'Gaussian'
			%---------------------------------------------------
			sigma   = K{s}.LParam/K{s}.RT;
			h       = round(4*sigma);
			h       = exp(-[-h:h].^2/(2*sigma^2));
			n       = length(h);
			d       = [1:n] - (n + 1)/2;
			if      n == 1, h = 1; end

			otherwise
			%---------------------------------------------------
			error('Low pass Filter option unknown');
			return

		end

		% create and normalize low pass filter
		%-----------------------------------------------------------
		K{s}.KL = spdiags(ones(k,1)*h,d,k,k);
		K{s}.KL = spdiags(1./sum(K{s}.KL')',0,k,k)*K{s}.KL;

		% make high pass filter
		%-----------------------------------------------------------
		switch K{s}.HChoice

			case 'none'
			%---------------------------------------------------
			K{s}.KH = [];

			case 'specify'
			%---------------------------------------------------
			n       = fix(2*(k*K{s}.RT)/K{s}.HParam + 1);
			X       = spm_dctmtx(k,n);
			K{s}.KH = sparse(X(:,[2:n]));

			otherwise
			%---------------------------------------------------
			error('High pass Filter option unknown');
			return

		end

	end

	% return structure
	%-------------------------------------------------------------------
	vargout = K;

	case 'apply'
	%-------------------------------------------------------------------
	if iscell(K)

		% ensure requisite feild are present
		%-----------------------------------------------------------
		if ~isfield(K{1},'KL')
			K = spm_filter('set',K);
		end

		for s = 1:length(K)

			% select data
			%---------------------------------------------------
			y = Y(K{s}.row,:);

			% apply low pass filter
			%---------------------------------------------------
			y = K{s}.KL*y;

			% apply high pass filter
			%---------------------------------------------------
			if ~isempty(K{s}.KH)
				y = y - K{s}.KH*(K{s}.KH'*y);
			end

			% reset filtered data in Y
			%---------------------------------------------------
			Y(K{s}.row,:) = y;

		end

	% K is simply a convolution matrix
	%-------------------------------------------------------------------
	else
		Y = K*Y;
	end

	% return filtered data
	%-------------------------------------------------------------------
	vargout   = Y;

	otherwise
	%-------------------------------------------------------------------
	warning('Filter option unknown');

end

function C = spm_dctmtx(N,K,n,f)
% Creates basis functions for Discrete Cosine Transform.
% FORMAT C = spm_dctmtx(N,K,n)
%     OR C = spm_dctmtx(N,K)
%     OR D = spm_dctmtx(N,K,n,'diff')
%     OR D = spm_dctmtx(N,K,'diff')
% N - dimension
% K - order
% n - optional points to sample
%____________________________________________________________________________
% spm_dctmtx creates a matrix for the first few basis functions of a one
% dimensional discrete cosine transform.
% With the 'diff' argument, spm_dctmtx produces the derivatives of the
% DCT.
%
% See:    Fundamentals of Digital Image Processing (p 150-154).
%         Anil K. Jain 1989.
%____________________________________________________________________________
% @(#)spm_dctmtx.m	1.3 John Ashburner MRCCU/FIL 96/08/14

d = 0;

if (nargin == 2)
	n = (0:(N-1))';
	if (nargin == 3)
		d = 1;
	end
elseif (nargin == 3)
	if (strcmp(n,'diff'))
		d = 1;
		n = (0:(N-1))';
	else
		n = n(:);
	end
elseif (nargin == 4)
	n = n(:);
	if (strcmp(f,'diff'))
		d = 1;
	else
		error('Incorrect Usage');
	end
else
	error('Incorrect Usage');
end

C = zeros(size(n,1),K);

if (d == 0)
	C(:,1)=ones(size(n,1),1)/sqrt(N);
	for k=2:K
		C(:,k) = sqrt(2/N)*cos(pi*(2*n+1)*(k-1)/(2*N));
	end
else
	for k=2:K
		C(:,k) = -2^(1/2)*(1/N)^(1/2)*sin(1/2*pi*(2*n*k-2*n+k-1)/N)*pi*(k-1)/N;
	end
end

author: canlab

repo: CanlabCore


text:
function life3D(action)
%LIFE3D   MATLAB's version of Conway's Game of Life.
%   "Life" is a cellular automaton invented by John
%   Conway that involves live and dead cells in a
%   rectangular, two-dimensional universe. In
%   MATLAB, the universe is a sparse matrix that
%   is initially all zero.
%
%   Whether cells stay alive, die, or generate new
%   cells depends upon how many of their eight
%   possible neighbors are alive. By using sparse
%   matrices, the calculations required become
%   astonishingly simple. We use periodic (torus)
%   boundary conditions at the edges of the
%   universe. Pressing the "Start" button
%   automatically seeds this universe with several
%   small random communities. Some will succeed
%   and some will fail.
%
%   Expanded to 3D by:
%       Leandro Barajas 06-20-2002 L.G.Barajas@ieee.org
%
%   C. Moler, 7-11-92, 8-7-92.
%   Adapted by Ned Gulley, 6-21-93
%

%   Copyright 1984-2001 The MathWorks, Inc.
%   $Revision: 5.9 $  $Date: 2001/04/15 12:03:02 $

% Possible actions:
% initialize
% start

% Information regarding the play status will be held in
% the axis user data according to the following table:
play= 1;
stop=-1;

if nargin<1,
   action='initialize';
end;

if strcmp(action,'initialize'),
   figNumber=figure( ...
      'Name','Life3D: Conway''s Game of Life in 3-Dimesions', ...
      'NumberTitle','off', ...
      'DoubleBuffer','on', ...
      'Visible','off', ...
      'Color','white', ...
      'BackingStore','off');
   axes( ...
      'Units','normalized', ...
      'Position',[0.05 0.05 0.75 0.90], ...
      'Visible','off', ...
      'DrawMode','fast', ...
      'Color','none',...
      'NextPlot','add');
%      'NextPlot','replace' );

   text(0,0,{'Press the "Start" button to see the Game of Life demo' 'Use the slider to change the number of initial cells.'}, ...
      'HorizontalAlignment','center');
   axis([-1 1 -1 1]);

   %===================================
   % Information for all buttons
   labelColor=[0.8 0.8 0.8];
   yInitPos=0.90;
   xPos=0.85;
   btnLen=0.10;
   btnWid=0.10;
   % Spacing between the button and the next command's label
   spacing=0.05;

   %====================================
   % The CONSOLE frame
   frmBorder=0.02;
   yPos=0.05-frmBorder;
   frmPos=[xPos-frmBorder yPos btnLen+2*frmBorder 0.9+2*frmBorder];
   h=uicontrol( ...
      'Style','frame', ...
      'Units','normalized', ...
      'Position',frmPos, ...
      'BackgroundColor',[0.50 0.50 0.50]);

   %====================================
   % The START button
   btnNumber=1;
   yPos=0.90-(btnNumber-1)*(btnWid+spacing);
   labelStr='Start';
   cmdStr='start';
   callbackStr='life3D(''start'');';

   % Generic button information
   btnPos=[xPos yPos-spacing btnLen btnWid];
   startHndl=uicontrol( ...
      'Style','pushbutton', ...
      'Units','normalized', ...
      'Position',btnPos, ...
      'String',labelStr, ...
      'Interruptible','on', ...
      'Callback',callbackStr);

   %====================================
   % The STOP button
   btnNumber=2;
   yPos=0.90-(btnNumber-1)*(btnWid+spacing);
   labelStr='Stop';
   % Setting userdata to -1 (=stop) will stop the demo.
   callbackStr='set(gca,''Userdata'',-1)';

   % Generic button information
   btnPos=[xPos yPos-spacing btnLen btnWid];
   stopHndl=uicontrol( ...
      'Style','pushbutton', ...
      'Units','normalized', ...
      'Position',btnPos, ...
      'Enable','off', ...
      'String',labelStr, ...
      'Callback',callbackStr);

   %====================================
   % The NumberOfCell Slider
   labelStr='# Cells';

   callbackStr='h=findobj(gcf,''Tag'',''StartCells'');set(h,''Tooltip'',sprintf(''Initial Cells: %4d'',floor(get(h,''Value''))));';
   infoHndl=uicontrol( ...
      'Style','slider', ...
      'Units','normalized', ...
      'Position',[xPos btnWid+btnWid*2+0.05 btnLen/4 0.10*3], ...
      'String',labelStr, ...
      'SliderStep',[.01 .1], ...
      'Max',100, ...
      'Min',1, ...
      'Value',20, ...
      'Tooltip','Initial Cells:  20',...
      'Tag','StartCells', ...
      'Callback',callbackStr);

  %====================================
   % The INFO button
   labelStr='Info';
   callbackStr='life3D(''info'')';
   infoHndl=uicontrol( ...
      'Style','push', ...
      'Units','normalized', ...
      'Position',[xPos 0.20 btnLen 0.10], ...
      'String',labelStr, ...
      'Callback',callbackStr);

   %====================================
   % The CLOSE button
   labelStr='Close';
   callbackStr='close(gcf)';
   closeHndl=uicontrol( ...
      'Style','push', ...
      'Units','normalized', ...
      'Position',[xPos 0.05 btnLen 0.10], ...
      'String',labelStr, ...
      'Callback',callbackStr);

   % Uncover the figure
   hndlList=[startHndl stopHndl infoHndl closeHndl];
   set(figNumber,'Visible','on', ...
      'UserData',hndlList);

   view(3);
   StartCells = 20;
   m = 19;
   colormap(jet(m));
   colorbar;

elseif strcmp(action,'start'),
   objh=findobj(gcf,'Tag','StartCells');
   StartCells = get(objh,'Value');
   m = 19;
   cla;
   axHndl=gca;
   figNumber=gcf;
   hndlList=get(figNumber,'Userdata');
   startHndl=hndlList(1);
   stopHndl=hndlList(2);
   infoHndl=hndlList(3);
   closeHndl=hndlList(4);
   set([startHndl closeHndl infoHndl],'Enable','off');
   set(stopHndl,'Enable','on');

   % ====== Start of Demo
   set(axHndl, ...
      'UserData',play, ...
      'DrawMode','fast', ...
      'Visible','on');
   box off

   X = zeros(m,m,m);

   p = -1:1;
   for count=1:StartCells,
      kx=floor(rand*(m-4))+2;
      ky=floor(rand*(m-4))+2;
      kz=floor(rand*(m-4))+2;
      X(kx+p,ky+p,kz+p)=(rand(3,3,3)>0.5);
   end;

   % The following statements plot the initial configuration.
   % The "find" function returns the indices of the nonzero elements.
   plothandle=[];
   figure(gcf);
   hold on
   colour=jet(m);
   for k=1:m
       [i,j] = find(X(:,:,k));
%       kd = sub2ind(size(X(i,j,k)),i,j)
       kd = k(ones(size(i)));
       if isempty(i)
           i = 1;
           j = 1;
           kd = NaN;
       end
       plothandle(k) = line(i,j,kd);
       set(plothandle(k),'linestyle','none',...
                      'Marker','.',...
                      'MarkerSize',20*2,...
                      'MarkerFaceColor',colour(k,:),...
                      'MarkerEdgeColor',colour(k,:),...
                      'EraseMode','normal');

   end
   drawnow
   axis([0 m+1 0 m+1 0 m+1]);

   hold off

   % Whether cells stay alive, die, or generate new cells depends
   % upon how many of their eight possible neighbors are alive.
   % Here we generate index vectors for four of the eight neighbors.
   % We use periodic (torus) boundary conditions at the edges of the universe.

   n = [m 1:m-1];
   e = [2:m 1];
   s = [2:m 1];
   w = [m 1:m-1];
   u = [m 1:m-1];
   d = [2:m 1];

   rotate3d on
   while get(axHndl,'UserData')==play,
    % How many of eight+5+5 neighbors are alive. (only the ones that share at least one border
    %      N = X(n,:,:) + X(s,:,:) + X(:,e,:) + X(:,w,:) + ...
    %          X(n,e,:) + X(n,w,:) + X(s,e,:) + X(s,w,:);

    % 3D Version

    % Use Euclidean distance
    d1 = 1/1;           % Adjacent cells (Share 4 vertix)
    d2 = 1/sqrt(2);     % Diagonal cells (Share 2 vertix)
    d3 = 1/sqrt(3);     % Double Diagonal cells (Share 1 vertix)

    % How many of 9+9+8 neighbors are alive. (Only the ones that share at least one vertix)
      N = X(n,:,:)*d1 + X(s,:,:)*d1 + X(:,e,:)*d1 + X(:,w,:)*d1 + ...
          X(n,e,:)*d2 + X(n,w,:)*d2 + X(s,e,:)*d2 + X(s,w,:)*d2 + ...
          X(n,:,u)*d2 + X(s,:,u)*d2 + X(:,e,u)*d2 + X(:,w,u)*d2 + ...
          X(n,:,d)*d2 + X(s,:,d)*d2 + X(:,e,d)*d2 + X(:,w,d)*d2 + ...
          X(n,e,u)*d3 + X(s,e,u)*d3 + X(s,w,u)*d3 + X(n,w,u)*d3 + ...
          X(n,e,d)*d3 + X(s,e,d)*d3 + X(s,w,d)*d3 + X(n,w,d)*d3 + ...
          X(:,:,u)*d1 + X(:,:,d)*d1;

      % A live cell with two live neighbors, or any cell with three
      % neigbhors, is alive at the next time step.
      Xold = X;

      CL = 3.0;             % Minimun # cells to create a new one
      CH = 4.0;             % Maximum # cells to create a new one
      KL = 4.0;             % Minimun # cells to kill one

      X = ( (    ((N >= CL) & (N <= CH)) ) &...    % Create cells
         ~(  X & ((N >  KL)           )) ) ;       % Kill cells by overpopulation

      Xsum = sum(X(:));
      if (Xsum>prod(size(X))/8);
         X = ( X & (rand(size(X))>0.1) ); % Expontaneous cell annihilation
      end

      if Xold==X   % if no change the exit
         break;
      end

      % Update plot.
      for k=1:m
          [i,j] = find(X(:,:,k));
          set(plothandle(k),'xdata',i,'ydata',j,'zdata',k(ones(size(i))));
      end

      set(plothandle,'Visible','on');
      box on
      title( sprintf('Total Cells: %d',Xsum))
      xlabel(sprintf('Density: %5.2f%%',mean(X(:))*100))
      ylabel(sprintf('Rate: %5.2f%%',(mean(X(:))-mean(Xold(:)))*100))

      drawnow
      pause(0.2)
   end

   % ====== End of Demo

   set([startHndl closeHndl infoHndl],'Enable','on');
   set(stopHndl,'Enable','off');

elseif strcmp(action,'info');
   helpwin(mfilename);

end;    % if strcmp(action, ...

repo: mcbench-benchmarks


text:
function model = mogUpdateCovariance(model)

% MOGUPDATECOVARIANCE Update the covariances of an MOG model.
% FORMAT
% DESC updates the covariance matrices of a mixtures of
% Gaussians model. The implementation currently uses an
% eigenvalue based update.
% ARG model : the model which is to be updated.
% RETURN model : the model with updated covariances.
%
% SEEALSO : mogCreate, mogUpdateMean, mogEstep
%
% COPYRIGHT : Neil D. Lawrence, 2006

% MLTOOLS

for i = 1:model.m

  centredY = model.Y - repmat(model.mean(i,:), model.N, 1);
  centredY = centredY.*repmat(sqrt(model.posterior(:,i)), 1, model.d);
  switch model.covtype
    case 'ppca'
     C = (centredY'*centredY+0.001*eye(model.d))/sum(model.posterior(:, i)+.001);
     [vec, val] = eig(C);
     val = diag(val);
     [val, ind] = sort(val);
     ind = ind(end:-1:1);
     val = val(end:-1:1);
     vec = vec(:, ind(1:model.q));
     sigma2 = mean(val(model.q+1:end));
     if sigma2<eps
       sigma2 = eps;
     end
     lambda = val(1:model.q) - sigma2;
     %[sigma2, eigVec, lambda] = ppca(C, model.q);
     if length(lambda) ~= model.q
       % Something's wrong here ...
       sigma2 = 1e-6;
       warning('Not enough eigenvalues extracted.')
       lambdaTemp = lambda;
       lambda = zeros(model.q, 1);
       lambda(1:length(lambdaTemp)) = lambdaTemp;
     end

     model.sigma2(i) = sigma2;
     model.W{i} = vec*diag(sqrt(lambda));
     model.U{i} = sqrt(sigma2)*eye(model.d);
     for j = 1:model.q
       model.U{i} = cholupdate(model.U{i}, model.W{i}(:, j));
     end
   case 'spherical'
    model.sigma2(i) = sum(sum(centredY.*centredY))/(model.d*sum(model.posterior(:, i)));
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat