text:
% SVM_TEST Calculate labels for an SVM model
%
% Usage
%    [labels, votes, feature_labels] = SVM_TEST(db, model, test_set)
%
% Input
%    db (struct): The database containing the feature vector.
%    model (struct): The affine space model obtained from svm_train.
%    test_set (int): The object indices of the testing instances.
%
% Output
%    labels (int): The assigned labels.
%    votes (int): The number of votes for each testing instance and class
%       pair.
%    feature_labels (int): The labels assigned to the individual features.
%
% See also
%    SVM_TRAIN, CLASSIF_ERR, CLASSIF_RECOG

function [labels,votes,feature_labels,K,sv_coef,dec] = svm_test(db,model,test_set)
	% Determine indices of feature vectors belonging to the test objects.
	ind_features = [db.indices{test_set}];

	if ~model.full_test_kernel && ...
		(model.svm.Parameters(2) == 4 || ...
		model.svm.Parameters(2) == 5 || ...
		model.svm.Parameters(2) == 6 || ...
		model.svm.Parameters(2) == 7)
		% If we're working with a precalculated kernel and only some of the testing
		% vectors are available, only use these (unless full_test_kernel is set).

		[kernel_mask,kernel_ind] = ismember(1:size(db.features,2),db.kernel.kernel_set);
		ind_features = ind_features(kernel_mask(ind_features));
	end

	% Obtain the feature vector labels for ind_features.
	[feature_labels,temp,K,sv_coef,dec] = svm_feature_test(db,model.svm,ind_features,model.full_test_kernel);

	for l = 1:length(test_set)
		% For each object in the test set, count the votes assigned to each of its feature vectors.

		ind = find(ismember(ind_features,db.indices{test_set(l)}));
		votes(:,l) = histc(feature_labels(ind),1:length(db.src.classes));
		[temp,labels(l)] = max(votes(:,l));
	end
end

function [labels,votes,K,sv_coef,dec] = svm_feature_test(db,model,ind_features0,full_test_kernel)
	% As in SVM_CALC_KERNEL, it's useful to break down kernel computations into
	% blocks.
	block_size = 8192;

	% How many classes are in this model?
	class_ct = model.nr_class;

	% Prepare the matrix of support vector coefficients for each pair of classes,
	% since each corresponds to a different SVM.
	sv_coef = zeros(model.totalSV,class_ct*(class_ct-1)/2,class(db.features));

	pairs = [];

	% Fill the sv_coef matrix. For more details on the model.sv_coef matrix, see
	% the LIBSVM FAQ. The difference with our sv_coefs is that our stores the
	% coefficients for each class pair, with zeros for the vectors belonging to
	% neither class.
	r = 1;
	for n1 = 1:class_ct
		for n2 = n1+1:class_ct
			pairs = [pairs [n1; n2]];

			class1_SVs = 1+sum(model.nSV(1:n1-1)):sum(model.nSV(1:n1));
			class2_SVs = 1+sum(model.nSV(1:n2-1)):sum(model.nSV(1:n2));

			sv_coef(class1_SVs,r) = model.sv_coef(class1_SVs,n2-1);
			sv_coef(class2_SVs,r) = model.sv_coef(class2_SVs,n1);
			r = r+1;
		end
	end

	% Prepare matrix of decision values for each testing vector and class pair.
	dec = zeros(length(ind_features0),size(sv_coef,2),class(db.features));

	if full_test_kernel && ...
		(model.Parameters(2) == 4 || ...
		model.Parameters(2) == 5 || ...
		model.Parameters(2) == 6 || ...
		model.Parameters(2) == 7)
		% If we have to recalculate kernel, prepare square norms. This should only
		% be necessary in Gaussian kernel case, though (6 & 7)...
		norm2 = sum(abs(db.features(:,db.kernel.kernel_set(model.SVs))).^2,1);
	end

	n = 1;
	while n <= length(ind_features0)
		% Only deal with testing vectors block by block.
		mask = n:min(n+block_size-1,length(ind_features0));
		ind_features = ind_features0(mask);

		if model.Parameters(2) == 0			% linear kernel
			% Linear kernel was calculated in LIBSVM; recalculate it here.
			K = db.features(:,ind_features).'*model.SVs.';
		elseif model.Parameters(2) == 2		% Gaussian kernel
			% Gaussian kernel was calculated in LIBSVM; recalculate it here.
			norm1 = sum(abs(db.features(:,ind_features)).^2,1);
			model.SVs=full(model.SVs);
			norm2 = sum(abs(model.SVs.').^2,1);
			norm2=full(norm2);
			K = bsxfun(@plus,norm1.',norm2)-...
				2*db.features(:,ind_features).'*model.SVs.';
			K = exp(-model.Parameters(4)*K);
		elseif full_test_kernel && ...
			(model.Parameters(2) == 4 || ...
			model.Parameters(2) == 5 || ...
			model.Parameters(2) == 6 || ...
			model.Parameters(2) == 7)
			% Kernel was precalculated, but for testing, we have to recalculate it. This
			% can be necessary if the set of kernel vectors is restricted.

			if strcmp(db.kernel.kernel_type,'linear')
				% Calculate the linear kernel of the testing vector with respect to the
				% support vectors.
				K = db.features(:,ind_features).'*db.features(:,db.kernel.kernel_set(model.SVs));
			elseif strcmp(db.kernel.kernel_type,'gaussian')
				% Same for Gaussian kernel, except that we don't multiply by gamma and
				% exponentiate here. This is done later, as in the precalculated case.
				norm1 = sum(abs(db.features(:,ind_features)).^2,1);
				K = bsxfun(@plus,norm1.',norm2)-...
					2*db.features(:,ind_features).'*db.features(:,db.kernel.kernel_set(model.SVs));
			else
				error('unknown kernel type');
			end
		elseif model.Parameters(2) == 4	|| ...
			model.Parameters(2) == 6
			% Precalculated kernel in square format, so just extract the rows corre-
			% sponding to the support vectors and the columns corresponding to the
			% testing vectors.

			% Translate to the kernel indexation.
			[kernel_mask,kernel_ind] = ismember(1:size(db.features,2),db.kernel.kernel_set);
			K = db.kernel.K(1+model.SVs,kernel_ind(ind_features)).';
		elseif model.Parameters(2) == 5 || ...
			model.Parameters(2) == 7
			% Precalculated kernel in triangular format. Need to unpack the kernel for
			% desired support vectors and testing vectors.

			% Prepare the kernel matrix.
			K = zeros(length(ind_features),length(model.SVs),class(db.kernel.K));

			% Obtain kernel indexation.
			[kernel_mask,kernel_ind] = ismember(1:size(db.features,2),db.kernel.kernel_set);

			for k = 1:length(ind_features)
				% For each of the testing vectors, extract the indices in the triangular
				% kernel.

				% Translate to the kernel indexation.
				idx = kernel_ind(ind_features(k));

				% Separate support vectors into those with index below testing vector index
				% and those above. Need to cast to integers since floating-point indexing
				% causes problems. Take 64 bits to be safe...
				lower_SVs = int64(model.SVs(model.SVs<=idx));
				upper_SVs = int64(model.SVs(model.SVs>idx));

				% Lower SVs are obtained from the SV rows in the idx column.
				K(k,model.SVs<=idx) = ...
					db.kernel.K((idx-1)*(idx-1+3)/2+1+lower_SVs);
				% Upper SVs are obtained from the idx row in the SV columns.
				K(k,model.SVs>idx) = ...
					db.kernel.K((upper_SVs-1).*(upper_SVs-1+3)/2+1+idx);
			end
		end

		if (model.Parameters(2) == 4 && strcmp(db.kernel.kernel_type,'gaussian')) || ...
			model.Parameters(2) == 6 || model.Parameters(2) == 7
			% If we have a Gaussian kernel, we need to multiply by gamma and exponentia-
			% te.
			K = exp(-model.Parameters(4)*K);
		end

		% Calculate the decision values for the current block. See LIBSVM FAQ for
		% the detailed formula.
		dec(mask,:) = bsxfun(@minus,K*sv_coef,model.rho.');

		n = n+length(mask);
	end

	% Prepare the votes matrix, which will contain the votes assigned to each
	% class for each testing vector.
	votes = zeros(size(dec,1),class_ct);

	for r = 1:size(pairs,2)
		% For each pair, calculate the "winners" and add one to their tally.
		n1 = pairs(1,r);
		n2 = pairs(2,r);

		votes(dec(:,r)>=0,n1) = votes(dec(:,r)>=0,n1)+1;
		votes(dec(:,r)<0,n2) = votes(dec(:,r)<0,n2)+1;
	end

	% Assign the class to that with the largest number of votes.
	[temp,labels] = max(votes,[],2);

	% Translate labels to those used in the model.
	labels = model.Label(labels);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function omrun (modelname,datasetfile)  % Filename des *.mo-files, modellname
%
% Execute compiled Modelica-Model with 'datasetfile', containing
% parameters, etc. while the existing datasetfile 'modelname_init.txt' is copied to 'modelname_init_old.txt.
%
% SYNTAX: omrun(modelname,datasetfile) optional: datasetfile
% Beispiel: omrun('package.model','mydasetfile')
%
% Feedback/problems: Christian Schaad, ingenieurbuero@christian-schaad.de

delete([modelname,'_res.mat']); % Delete old result file

linux=isunix;
if exist('datasetfile')
% Originales Datasetfile sichern...
%switch linux
%    case 1
%eval(['!cp ',modelname,'_init.txt ',modelname,'_init_old.txt']);
%    case 0
%eval(['!copy ',modelname,'_init.txt ',modelname,'_init_old.txt']);
%end
% Neues Datasetfile kopieren...
switch linux
    case 1
eval(['!cp ',datasetfile,' ',modelname,'_init.xml']);
    case 0
        eval(['!copy ',datasetfile,' ',modelname,'_init.xml']);
end
end
% Ausf???hren....

switch linux
    case 1
    %    eval(['system(''./',modelname,''') '])      %%% Use this for OM_1.4.4 (linux)
	 eval(['system(''LD_LIBRARY_PATH= ./',modelname, ''') '])    %%% Use this for OM_1.4.3   (linux)
    case 0
        eval(['system(''',modelname,'.exe'')'])
end

author: Sable

repo: mcbench-benchmarks


text:
%RNA Expression Analysis
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 12/22/2011
classdef RNAExpression
    methods (Static = true)
        function run(sim, fileName)
            %import classes
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.PrintUtil;
            import edu.stanford.covert.util.ComputationUtil;

            %load simulation
            g = sim.gene;
            r = sim.state('Rna');

            %get data
            rnaExp = r.expectedGeneExpression(:, 1);

            exp_mRNAWt = 0.041;
            exp_rRNAWt = [0.017; 0.271; 0.525];
            exp_sRNAWt = 0;
            exp_tRNAWt = 0.146;
            assertEqual(1, exp_mRNAWt + sum(exp_rRNAWt) + exp_sRNAWt + exp_tRNAWt);

            %calculate gene molecular weights approximately
            geneMWs = (r.matureRNAGeneComposition * (r.molecularWeights(r.aminoacylatedIndexs) ./ (r.matureRNAGeneComposition' * g.lengths))) ...
                .* g.lengths;

            %check RNA MWs computed correctly
            monocistronicIdxs = find(sum(r.matureRNAGeneComposition(g.mRNAIndexs, :), 1) == 1);
            monocistronicmRNAGeneIdxs = g.mRNAIndexs(any(r.matureRNAGeneComposition(g.mRNAIndexs, monocistronicIdxs), 2));
            assertElementsAlmostEqual(geneMWs(monocistronicmRNAGeneIdxs), r.molecularWeights(r.aminoacylatedIndexs(monocistronicIdxs)), 'relative', 1e-12);
            assertElementsAlmostEqual(geneMWs(g.rRNAIndexs), r.molecularWeights(r.aminoacylatedIndexs(r.matureRRNAIndexs)), 'relative', 1e-12);
            assertElementsAlmostEqual(geneMWs(g.sRNAIndexs), r.molecularWeights(r.aminoacylatedIndexs(r.matureSRNAIndexs)), 'relative', 1e-12);
            assertElementsAlmostEqual(geneMWs(g.tRNAIndexs), r.molecularWeights(r.aminoacylatedIndexs(r.matureTRNAIndexs)), 'relative', 1e-12);
            assertElementsAlmostEqual(r.molecularWeights(r.aminoacylatedIndexs), r.matureRNAGeneComposition' * geneMWs, 'relative', 1e-12);

            %reconstruct gene expression
            mRNAWt = rnaExp(g.mRNAIndexs)' * geneMWs(g.mRNAIndexs);
            rRNAWt = rnaExp(g.ribosomalRRNAIndexs) .* geneMWs(g.ribosomalRRNAIndexs);
            sRNAWt = rnaExp(g.sRNAIndexs)' * geneMWs(g.sRNAIndexs);
            tRNAWt = rnaExp(g.tRNAIndexs)' * geneMWs(g.tRNAIndexs);

            rnaExp(g.mRNAIndexs) = rnaExp(g.mRNAIndexs) * exp_mRNAWt / mRNAWt;
            rnaExp(g.ribosomalRRNAIndexs) = rnaExp(g.ribosomalRRNAIndexs) .* exp_rRNAWt ./ rRNAWt;
            rnaExp(g.sRNAIndexs) = rnaExp(g.sRNAIndexs) * exp_sRNAWt / sRNAWt;
            rnaExp(g.tRNAIndexs) = rnaExp(g.tRNAIndexs) * exp_tRNAWt / tRNAWt;

            rnaExp = rnaExp / sum(rnaExp);

            exp_RNAWt = [exp_mRNAWt; exp_rRNAWt; exp_sRNAWt; exp_tRNAWt];
            RNAwt = [mRNAWt; rRNAWt; sRNAWt; tRNAWt];
            exp_RNAWt = exp_RNAWt / sum(exp_RNAWt);
            RNAwt = RNAwt / sum(RNAwt);

            %plot RNA weight fractions
            [axesHandle, figHandle] = PlotUtil.newAxesHandle();

            cla(axesHandle);
            hold(axesHandle, 'on');
            h = [
                plot(axesHandle, RNAwt(1), exp_RNAWt(1), 'r.')
                plot(axesHandle, RNAwt(2:4), exp_RNAWt(2:4), 'g.')
                plot(axesHandle, RNAwt(5), exp_RNAWt(5), 'b.')
                plot(axesHandle, RNAwt(6), exp_RNAWt(6), 'c.')
                ];
            legend(h, {'mRNA', 'rRNA', 'sRNA', 'tRNA'}, 'Location', 'NorthWest');
            xlabel(axesHandle, 'Reconstruction - Earliest Stage', 'FontSize', 12)
            ylabel(axesHandle, 'Reconstruction - Later Stage', 'FontSize', 12)
            set(axesHandle, 'xscale', 'log', 'yscale', 'log', 'xlimmode', 'auto',  'ylimmode', 'auto', 'xminortick', 'off', 'yminortick', 'off')
            axis(axesHandle, 'equal');
            xlim(axesHandle, [min([xlim(axesHandle) ylim(axesHandle)]) max([xlim(axesHandle) ylim(axesHandle)])]);
            ylim(axesHandle, xlim(axesHandle));
            line(xlim(axesHandle), ylim(axesHandle), 'Parent', axesHandle, 'LineStyle', ':', 'Color', [0.5 0.5 0.5])

            saveas(figHandle, [fileName '-EarlyVsLaterPhaseReconstruction-WeightFractions.pdf']);
            close(figHandle);

            %plot all RNA
            [axesHandle, figHandle] = PlotUtil.newAxesHandle();

            cla(axesHandle);
            hold(axesHandle, 'on');
            h = [
                plot(axesHandle, rnaExp(g.mRNAIndexs), r.expectedGeneExpression(g.mRNAIndexs, 1), 'r.')
                plot(axesHandle, rnaExp(g.rRNAIndexs), r.expectedGeneExpression(g.rRNAIndexs, 1), 'g.')
                plot(axesHandle, rnaExp(g.sRNAIndexs), r.expectedGeneExpression(g.sRNAIndexs, 1), 'b.')
                plot(axesHandle, rnaExp(g.tRNAIndexs), r.expectedGeneExpression(g.tRNAIndexs, 1), 'c.')
                ];
            legend(h, {'mRNA', 'rRNA', 'sRNA', 'tRNA'}, 'Location', 'NorthWest');
            xlabel(axesHandle, 'Reconstruction - Earliest Stage', 'FontSize', 12)
            ylabel(axesHandle, 'Reconstruction - Later Stage', 'FontSize', 12)
            set(axesHandle, 'xscale', 'log', 'yscale', 'log', 'xlimmode', 'auto',  'ylimmode', 'auto', 'xminortick', 'off', 'yminortick', 'off')
            axis(axesHandle, 'equal');
            xlim(axesHandle, [min([xlim(axesHandle) ylim(axesHandle)]) max([xlim(axesHandle) ylim(axesHandle)])]);
            ylim(axesHandle, xlim(axesHandle));
            line(xlim(axesHandle), ylim(axesHandle), 'Parent', axesHandle, 'LineStyle', ':', 'Color', [0.5 0.5 0.5])

            saveas(figHandle, [fileName '-EarlyVsLaterPhaseReconstruction.pdf']);
            close(figHandle);

            %save to excel sheet
            content = [g.wholeCellModelIDs num2cell(rnaExp) num2cell(r.expectedGeneExpression(:, 1))];
            colLabels = {'Gene', 'Expression-Early Reconstruction', 'Expression-Later Reconstruction'};
            PrintUtil.printToFile(content, colLabels, [fileName '.xls'], 'RNA Expression');
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function ideal_vocal_tract_GUI25
% Modifiable runGUI file
clc;
clear all;
fileName = 'ideal_vocal_tract.mat';    %USER - ENTER FILENAME
fileData=load(fileName);
temp=fileData(1).temp;

f = figure('Visible','on',...
'Units','normalized',...
'Position',[0,0,1,1],...
'MenuBar','none',...
'NumberTitle','off');

Callbacks_ideal_vocal_tract_GUI25(f,temp);    %USER - ENTER PROPER CALLBACK FILE
%panelAndButtonEdit(f, temp);       % Easy access to Edit Mode

% Note comment PanelandBUttonCallbacks(f,temp) if panelAndButtonEdit is to
% be uncommented and used
end

% ideal_vocal_tract_GUI25 design
% 2 Panels
%   #1 - input parameters
%   #2 - graphics displays
% 6 Graphic Panels
%   #1 - excitation impulses in time
%   #2 - excitation line frequency spectrum
%   #3 - VT impulse response
%   #4 - VT frequency response
%   #5 - periodic vowel in time
%   #6 - periodic vowel spectrum
% 1 TitleBox
% 6 Buttons
%   #1 - popupmenu - vowel choices
%   #2 - editable button - T: pitch period
%   #3 - editable button - N: analysis frame length
%   #4 - popupmenu - log/linear spectrum plot
%   #5 - pushbutton - Run ideal vocal tract
%   #6 - pushbutton - Close GUI


text:
function so = streamFilter(si,filter,varargin)
% STREAMFILTER  Filtered stream
while ~builtin('feval',filter,head(si),varargin{:})
    si = tail(si);
end
so = {head(si),delayEval(@streamFilter,{tail(si),filter,varargin{:}})};


text:
function recoverScenario(struct, idx_recover, mode, file_name)

% Authors:
%       Salih Guemues
%       Alexander Wischnewski
% Description:
%   function used to recalculate one simulation step to recover the
%   specific QP definition of the (soft) TMPC
% Inputs/parameters:
%   struct:      Log file struct
%   idx_recover: ID of the optimization problem which is plotted (debug_slow.tmpc_cnt)
%   mode:        Choose 'DD' for data dictionary parameters or 'script' for setupTMPCstruct to
%                   obtain the parameters
%   file_name:   if given, store results in seperate .mat file for dedicated visualization purposes

%% -------------------- HOW TO USE ---------------------------- %%
% Run Simulink model controller_dev ->
% execute convertSimLogs(logsout,'SimTest.mat'); ->
% load data via data = load('SimTest.mat');
% call this script via recoverScenario(data, t_recover, mode) and specify t_recover and mode

%% -------------------- DEFINE PARAMETERS ---------------------------- %%
% maximum number of iterations (high to check whether online solution may not have converged)
max_iter = 1000;

%% ---------------------- IMPORT CONTROLLER AND VEHICLE DATA ------------------------- %%
% get access to dictionary with controller parameters
sysDict = Simulink.data.dictionary.open('il_mvdc_mpc.sldd');
dDataSectObj2 = getSection(sysDict,'Design Data');
if(strcmp(mode, 'DD'))
    % load controller parameters from data dictionary
    sys = getValue(getEntry(dDataSectObj2,'sys'));
elseif(strcmp(mode, 'script'))
    % load controller parameters from initialize script
    sys = setupTMPCStruct(false);
else
    disp('No valid mode has been specified. Use DD or script');
end
help = getValue(getEntry(dDataSectObj2,'P_VDC_VirtualController'));
P_VDC_VirtualController = help.Value;
help = getValue(getEntry(dDataSectObj2,'P_VDC_PositiveAxLimScale'));
P_VDC_PositiveAxLimScale = help.Value;

% get access to dictionary with controller parameters
vehDict = Simulink.data.dictionary.open('il_vehicleparameter.sldd');
dDataSectObj3 = getSection(vehDict,'Design Data');
% get some more vehicle parameters
drag_coefficient = getValue(getEntry(dDataSectObj3, 'drag_coefficient'));
help = getValue(getEntry(dDataSectObj3,'roh_air'));
roh_air = help.Value;
vehiclemass_kg = getValue(getEntry(dDataSectObj3, 'vehiclemass_kg'));

%% ------------ IMPORT TRAJECTORY, VEHICLE, PREDICTION DATA AND RECONSTRUCT SCENARIO --------------- %%
% create TargetTrajectory struct and fill with logged data
TargetTrajectory.LapCnt         = struct.debug_slow.debug_slow_tartraj_LapCnt.Data(idx_recover,:);
TargetTrajectory.TrajCnt        = struct.debug_slow.debug_slow_tartraj_TrajCnt.Data(idx_recover,:);
TargetTrajectory.PointIdx       = uint16(0);
TargetTrajectory.s_loc_m        = struct.debug_slow.debug_slow_tartraj_s_loc_m.Data(idx_recover,:)';
TargetTrajectory.s_glob_m       = struct.debug_slow.debug_slow_tartraj_s_glob_m.Data(idx_recover,:)';
TargetTrajectory.x_m            = struct.debug_slow.debug_slow_tartraj_x_m.Data(idx_recover,:)';
TargetTrajectory.y_m            = struct.debug_slow.debug_slow_tartraj_y_m.Data(idx_recover,:)';
TargetTrajectory.psi_rad        = struct.debug_slow.debug_slow_tartraj_psi_rad.Data(idx_recover,:)';
TargetTrajectory.kappa_radpm    = struct.debug_slow.debug_slow_tartraj_kappa_radpm.Data(idx_recover,:)';
TargetTrajectory.v_mps          = struct.debug_slow.debug_slow_tartraj_v_mps.Data(idx_recover,:)';
TargetTrajectory.ax_mps2        = struct.debug_slow.debug_slow_tartraj_ax_mps2.Data(idx_recover,:)';
TargetTrajectory.banking_rad    = struct.debug_slow.debug_slow_tartraj_banking_rad.Data(idx_recover,:)';
TargetTrajectory.ax_lim_mps2    = struct.debug_slow.debug_slow_tartraj_ax_lim_mps2.Data(idx_recover,:)';
TargetTrajectory.ay_lim_mps2    = struct.debug_slow.debug_slow_tartraj_ay_lim_mps2.Data(idx_recover,:)';

% create VehicleDynamicState struct and fill with logged data
VehicleDynamicState.Pos.x_m             = struct.debug_slow.debug_slow_x_real_m.Data(idx_recover,:);
VehicleDynamicState.Pos.y_m             = struct.debug_slow.debug_slow_y_real_m.Data(idx_recover,:);
VehicleDynamicState.Pos.psi_rad         = struct.debug_slow.debug_slow_psi_real_rad.Data(idx_recover,:);
VehicleDynamicState.beta_rad            = struct.debug_slow.debug_slow_beta_real_rad.Data(idx_recover,:);
VehicleDynamicState.v_mps               = struct.debug_slow.debug_slow_vx_real_mps.Data(idx_recover,:);

ax_dist_mps2 = struct.debug_slow.debug_slow_ax_dist_mps2.Data(idx_recover, :)';
ay_dist_mps2 = struct.debug_slow.debug_slow_ay_dist_mps2.Data(idx_recover, :)';
% recalculate uncertainty matrices
UncertaintyTube = calcTubeShapeMatrices(sys.A_d, sys.B_d, P_VDC_VirtualController, sys.M_1, [ax_dist_mps2'; ay_dist_mps2'], sys.N_hor+1);
[PathPos, ~] = localTrajectoryMatching(VehicleDynamicState, TargetTrajectory);
% extract terminal set speed
v_terminal_mps = struct.debug_slow.debug_slow_be_l_abs.Data(idx_recover, 1);
dot_d_numerical = struct.debug_slow.debug_slow_dot_d_numerical_mps.Data(idx_recover, 1);

% get predictions vx_pred, d_psi_pred, s_dot_pred
x_pred_log       = struct.debug_slow.debug_slow_x_pred_m.Data(idx_recover, :)';
y_pred_log       = struct.debug_slow.debug_slow_y_pred_m.Data(idx_recover, :)';
vx_pred_log      = struct.debug_slow.debug_slow_vx_pred_mps.Data(idx_recover, :)';
d_pred_log       = struct.debug_slow.debug_slow_d_pred_m.Data(idx_recover, :)';
dot_d_pred_log   = struct.debug_slow.debug_slow_dot_d_pred_mps.Data(idx_recover, :)';
ax_pred_log      = struct.debug_slow.debug_slow_ax_pred_mps2.Data(idx_recover, :)';
ax_tire_pred_log = struct.debug_slow.debug_slow_ax_tire_pred_mps2.Data(idx_recover, :)';
ay_pred_log      = struct.debug_slow.debug_slow_ay_pred_mps2.Data(idx_recover, :)';

% get linearization and target trajectory
x_traj              = struct.debug_slow.debug_slow_x_traj_m.Data(idx_recover, :)';
y_traj              = struct.debug_slow.debug_slow_y_traj_m.Data(idx_recover, :)';
psi_traj            = struct.debug_slow.debug_slow_psi_traj_rad.Data(idx_recover, :)';
v_traj              = struct.debug_slow.debug_slow_v_traj_mps.Data(idx_recover, :)';
kappa_traj          = struct.debug_slow.debug_slow_kappa_traj_radpm.Data(idx_recover, :)';
ax_diff_traj        = struct.debug_slow.debug_slow_ax_diff_traj_mps2m.Data(idx_recover, :)';
ax_traj             = struct.debug_slow.debug_slow_ax_traj_mps2.Data(idx_recover, :)';
ay_traj             = struct.debug_slow.debug_slow_ay_traj_mps2.Data(idx_recover, :)';
ax_lim_mps2         = struct.debug_slow.debug_slow_ax_lim_mps2.Data(idx_recover, :)';
ay_lim_mps2         = struct.debug_slow.debug_slow_ay_lim_mps2.Data(idx_recover, :)';
d_lim_ub_m          = struct.debug_slow.debug_slow_d_lim_ub_m.Data(idx_recover, :)';
d_lim_lb_m          = struct.debug_slow.debug_slow_d_lim_lb_m.Data(idx_recover, :)';
d_Target_m          = struct.debug_slow.debug_slow_d_Target_m.Data(idx_recover, :)';
dot_d_Target_mps    = struct.debug_slow.debug_slow_dot_d_Target_mps.Data(idx_recover, :)';
ax_traj_old         = struct.debug_slow.debug_slow_ax_traj_mps2.Data(idx_recover-1, :)';
ay_traj_old         = struct.debug_slow.debug_slow_ay_traj_mps2.Data(idx_recover-1, :)';
vx_lin_mps          = struct.debug_slow.debug_slow_vx_lin_mps.Data(idx_recover, :)';
kappa_lin_radpm     = struct.debug_slow.debug_slow_kappa_lin_radpm.Data(idx_recover, :)';
u_opt_total_sim     = struct.debug_slow.debug_slow_u_opt_total.Data(idx_recover, :)';
u_opt_old_log       = struct.debug_slow.debug_slow_u_opt_total.Data(idx_recover-1, :)';

%% --------------------- SETUP INITIAL QP  ------------------------------ %%
% Define problem data by reconstructing it using the sparsity pattern
% init empty matrices
P_help = zeros(sys.osqp_n, sys.osqp_n);
A_help = zeros(sys.osqp_m, sys.osqp_n);
% fill values according to sparsity pattern
P_help(sys.P_i_lin) = sys.P_x_par;
% NOTE THAT P IS STORED AS UPPER TRIANGULAR, however calculations require
% FULL P MATRIX --> add transpose and substract diagonal elements
P_help = P_help + P_help' - diag(diag(P_help));
A_help(sys.A_i_lin) = sys.A_x_par;
q = sys.osqp_qpar;
l = sys.l_par;
u = sys.u_par;
P = sparse(P_help);
A = sparse(A_help);
% Create an OSQP object
prob = osqp;
% change settings
settings = prob.default_settings();
settings.alpha = 1.0;
settings.verbose = true;
settings.scaling = 0;  % number of scaling iterations
settings.max_iter = max_iter;
settings.warm_start = true;
settings.eps_abs = 1e-5;
settings.eps_rel = 1e-5;
% Setup workspace, change alpha parameter and disable prints
prob.setup(P, q, A, l, u, settings);

%% ----------- CALL prepareOptimizationProblem TO RECOVER QP ----------------- %%
[f, lb, ub, A_x, be_u_abs, be_l_abs, ~, error_state] =...
    prepareOptimizationProblem(VehicleDynamicState, dot_d_numerical, PathPos, ...
    v_traj, ax_diff_traj, ax_traj, ay_traj, ax_lim_mps2, ay_lim_mps2, ...
    d_lim_ub_m, d_lim_lb_m, d_Target_m, dot_d_Target_mps, ...
    vx_lin_mps, kappa_lin_radpm, UncertaintyTube, v_terminal_mps, ...
    ax_traj_old, ay_traj_old, u_opt_old_log, ...
    sys, P_VDC_VirtualController, drag_coefficient, roh_air, vehiclemass_kg, ...
    P_VDC_PositiveAxLimScale, true, 0.9);

%% --------------------- UPDATE AND SOLVE QP  ------------------------------ %%
% update QP vectors
prob.update('q', f, 'l', lb, 'u', ub);
prob.update('Ax', A_x);
% solve problem and get solution. two runs necessary to replicate warm start behavior of online
% solution while running the car.
res = prob.solve();
u_opt_total = res.x;
% meaning of status as in https://osqp.org/docs/interfaces/status_values.html
status = res.info.status_val;

% debug output
disp(['The problem has been solved during reoptimization with status: ' num2str(status)]);
disp(['It took ' num2str(res.info.iter) ' iterations and ' num2str(res.info.solve_time) ' seconds']);

%% ----------- CALL transformMPCResult TO RECOVER OPTIMIZED TRAJ. AND PREDICTIONS ----------------- %%
[u_opt_total, x_pred, y_pred, x_pred_left, y_pred_left, x_pred_right, y_pred_right, ...
    vx_pred, d_pred, dot_d_pred, ~, ax_pred, ax_tire_pred, ay_pred, ~, ~]  ...
    = transformMPCResult(status, u_opt_total, error_state,...
    x_traj, y_traj, psi_traj, v_traj, ax_diff_traj, d_Target_m, ax_lim_mps2, ay_lim_mps2, ...
    vx_lin_mps, kappa_lin_radpm, UncertaintyTube, ...
    sys, drag_coefficient, roh_air, vehiclemass_kg, P_VDC_PositiveAxLimScale);

%% plotting of optimization problem result

% recalculate tubes
[lb_v, ub_v] = calcBoundsPlot(zeros(sys.N_hor, 1), UncertaintyTube, false, P_VDC_VirtualController, sys.N_hor, repmat([1; 0; 0], 1, sys.N_hor));
[lb_d, ub_d] = calcBoundsPlot(zeros(sys.N_hor, 1), UncertaintyTube, false, P_VDC_VirtualController, sys.N_hor, repmat([0; 1; 0], 1, sys.N_hor));
[lb_dot_d, ub_dot_d] = calcBoundsPlot(zeros(sys.N_hor, 1), UncertaintyTube, false, P_VDC_VirtualController, sys.N_hor, repmat([0; 0; 1], 1, sys.N_hor));
[~, ub_tire1] = calcBoundsPlot(zeros(sys.N_hor+1, 1), UncertaintyTube, true, ...
    P_VDC_VirtualController, sys.N_hor+1, [1./(P_VDC_PositiveAxLimScale*ax_lim_mps2), 1./ay_lim_mps2]');
[~, ub_tire2] = calcBoundsPlot(zeros(sys.N_hor+1, 1), UncertaintyTube, true, ...
    P_VDC_VirtualController, sys.N_hor+1, [1./ax_lim_mps2, 1./ay_lim_mps2]');
[~, ub_tire3] = calcBoundsPlot(zeros(sys.N_hor+1, 1), UncertaintyTube, true, ...
    P_VDC_VirtualController, sys.N_hor+1, [-1./ax_lim_mps2, 1./ay_lim_mps2]');
[~, ub_tire4] = calcBoundsPlot(zeros(sys.N_hor+1, 1), UncertaintyTube, true, ...
    P_VDC_VirtualController, sys.N_hor+1, [-1./(P_VDC_PositiveAxLimScale*ax_lim_mps2), 1./ay_lim_mps2]');

% calculate limits for ax depending on positive accelerations
idx_pos = ax_traj > 0;
ax_lim_mps2_traj = ax_lim_mps2;
ax_lim_mps2_traj(idx_pos) = P_VDC_PositiveAxLimScale*ax_lim_mps2(idx_pos)';
idx_pos = ax_tire_pred > 0;
ax_lim_mps2_pred = ax_lim_mps2;
ax_lim_mps2_pred(idx_pos) = P_VDC_PositiveAxLimScale*ax_lim_mps2(idx_pos)';
idx_pos = ax_tire_pred_log > 0;
ax_lim_mps2_pred_log = ax_lim_mps2;
ax_lim_mps2_pred_log(idx_pos) = P_VDC_PositiveAxLimScale*ax_lim_mps2(idx_pos)' ;

figure;
grid on; hold on;
plot(x_traj, y_traj);
plot(x_pred, y_pred);
plot(x_pred_log, y_pred_log, '*');
plot(x_pred_left, y_pred_left, 'k--');
plot(x_pred_right, y_pred_right, 'k--', 'HandleVisibility','off');
grid on; axis equal;
legend('Target', 'Optimized', 'Logged', 'Tube');
xlabel('x East in m');
ylabel('y North in m');

% compare target and reoptimized trajectory
figure;
ax11 = subplot(2, 2, 1);
grid on; hold on;
plot(v_traj);
plot(vx_pred);
plot(vx_pred_log, '*');
plot(vx_pred(1:sys.N_hor) + ub_v', 'k--');
plot(vx_pred(1:sys.N_hor) + lb_v', 'k--', 'HandleVisibility','off');
plot([sys.N_hor+1, sys.N_hor+3], [be_u_abs(1), be_u_abs(1)], 'k');
plot([sys.N_hor+1, sys.N_hor+3], [be_l_abs(1), be_l_abs(1)], 'k', 'HandleVisibility','off');
plot(vx_lin_mps);
xlabel('Discretization points');
ylabel('Velocity in mps');
legend('Target', 'Optimized', 'Logged', 'Tube', 'Terminal Set', 'Linearization');

ax12 = subplot(2, 2, 2);
grid on; hold on;
plot(ax_traj./ax_lim_mps2_traj + kappa_traj.*v_traj.^2./ay_lim_mps2);
plot(ax_tire_pred./ax_lim_mps2_pred + ay_pred./ay_lim_mps2);
plot(ax_tire_pred_log./ax_lim_mps2_pred_log + ay_pred_log./ay_lim_mps2, '*');
plot(ax_tire_pred./ax_lim_mps2_pred + ay_pred./ay_lim_mps2 + ub_tire1', 'k--');
plot(ax_tire_pred./ax_lim_mps2_pred + ay_pred./ay_lim_mps2 - ub_tire2', 'k--', 'HandleVisibility','off');
plot([0, sys.N_hor], [1, 1], 'k');
plot([0, sys.N_hor], [-1, -1], 'k');
grid on;
xlabel('Discretization points');
ylabel('Acceleration constraint normalized');
legend('Target', 'Optimized', 'Logged', 'Tube', 'Limit');

ax13 = subplot(2, 2, 3);
grid on; hold on;
plot(d_Target_m)
plot(d_pred);
plot(d_pred_log, '*');
plot(d_pred(1:sys.N_hor) + ub_d', 'k--');
plot(d_pred(1:sys.N_hor) + lb_d', 'k--', 'HandleVisibility','off');
plot(d_lim_ub_m, 'k');
plot(d_lim_lb_m, 'k');
grid on;
xlabel('Discretization points');
ylabel('Lateral deviation in m');
legend('Target', 'Optimized', 'Logged', 'Tube', 'Limit');

ax14 = subplot(2, 2, 4);
grid on; hold on;
plot(-ax_traj./ax_lim_mps2 + kappa_traj.*v_traj.^2./ay_lim_mps2);
plot(-ax_tire_pred./ax_lim_mps2_pred + ay_pred./ay_lim_mps2);
plot(-ax_tire_pred_log./ax_lim_mps2_pred_log + ay_pred_log./ay_lim_mps2, '*');
plot(-ax_tire_pred./ax_lim_mps2_pred + ay_pred./ay_lim_mps2 - ub_tire3', 'k--');
plot(-ax_tire_pred./ax_lim_mps2_pred + ay_pred./ay_lim_mps2 + ub_tire4', 'k--', 'HandleVisibility','off');
plot([0, sys.N_hor], [1, 1], 'k');
plot([0, sys.N_hor], [-1, -1], 'k');
grid on;
xlabel('Discretization points');
ylabel('Acceleration constraint normalized');
legend('Target', 'Optimized', 'Logged', 'Tube', 'Limit');

linkaxes([ax11, ax12, ax13, ax14], 'x');

figure;
ax21 = subplot(2, 3, 1);
grid on; hold on;
plot(ax_traj);
plot(ax_pred);
plot(ax_pred_log, '*');
plot(-u_opt_total(1:2:2*(sys.N_hor+1)));
plot(ax_lim_mps2, 'k');
plot(-ax_lim_mps2, 'k');
grid on;
xlabel('Discretization points');
ylabel('Acceleration x in mps2');
legend('Target', 'Optimized', 'Logged', 'delta ax', 'Limit');

ax22 = subplot(2, 3, 2);
grid on; hold on;
plot(kappa_traj.*v_traj.^2);
plot(ay_pred);
plot(ay_pred_log, '*');
plot(u_opt_total(2:2:2*(sys.N_hor+1)));
plot(ay_lim_mps2, 'k');
plot(-ay_lim_mps2, 'k');
grid on;
xlabel('Discretization points');
ylabel('Acceleration y in mps2');
legend('Target', 'Optimized', 'Logged', 'delta ay', 'Limit');

ax23 = subplot(2, 3, 3);
grid on; hold on;
plot(dot_d_Target_mps)
plot(dot_d_pred);
plot(dot_d_pred_log, '*');
plot(dot_d_pred_log(1:sys.N_hor) + ub_dot_d', 'k--');
plot(dot_d_pred_log(1:sys.N_hor) + lb_dot_d', 'k--', 'HandleVisibility','off');
grid on;
xlabel('Discretization points');
ylabel('Lateral deviation derivative in mps');
legend('Target', 'Optimized', 'Logged', 'Tube');

ax26 = subplot(2, 3, 4);
grid on; hold on;
plot(u_opt_total(2*(sys.N_hor+1)+3:sys.n_slacks:end), 'b');
plot(u_opt_total_sim(2*(sys.N_hor+1)+3:sys.n_slacks:end), 'b*');
plot(u_opt_total(2*(sys.N_hor+1)+4:sys.n_slacks:end), 'c');
plot(u_opt_total_sim(2*(sys.N_hor+1)+4:sys.n_slacks:end), 'c*');
plot([0, sys.N_hor], [0, 0], 'k');
plot([0, sys.N_hor], [1.5, 1.5]/2, 'k');
grid on;
xlabel('Discretization points');
ylabel('Slacks lateral');
legend('Optimized upper', 'Logged upper', 'Optimized lower', 'Optimized lower', 'Limit');

ax27 = subplot(2, 3, 5);
grid on; hold on;
plot(u_opt_total(2*(sys.N_hor+1)+1:sys.n_slacks:end));
plot(u_opt_total_sim(2*(sys.N_hor+1)+1:sys.n_slacks:end), '*');
plot([0, sys.N_hor], [0, 0], 'k');
plot([0, sys.N_hor], [sys.slack_lim_rel, sys.slack_lim_rel], 'k');
grid on;
xlabel('Discretization points');
ylabel('Slacks tire 1 lower');
legend('Optimized upper', 'Logged', 'Limit');

ax28 = subplot(2, 3, 6);
grid on; hold on;
plot(u_opt_total(2*(sys.N_hor+1)+2:sys.n_slacks:end));
plot(u_opt_total_sim(2*(sys.N_hor+1)+2:sys.n_slacks:end), '*');
plot([0, sys.N_hor], [0, 0], 'k');
plot([0, sys.N_hor], [sys.slack_lim_rel, sys.slack_lim_rel], 'k');
grid on;
xlabel('Discretization points');
ylabel('Slacks tire 2 lower');
legend('Optimized', 'Logged', 'Limit');

linkaxes([ax21, ax22, ax23, ax26, ax27, ax28], 'x');

%% store data in external file for further visualization
if nargin > 3
    save(file_name);
end

end

author: TUMFTM


text:
% op_CSICombineCoils.m
%
% Combines the coils of MRSI data. First, coil's phases are adjusted to
% match the phase of the first coil. Then, sensetivity weights for coils are
% calculated by the formula w(i) = S(i)/sqrt(sum(S^2)) where S is intensity of the
% first point of the fids in the center of k space and i is the coil
% number. Coils are then summed by F = sum(W(i)*f(i)) where W(i) is coil's
% sensetivity weight at coil i, f(i) is the fids at coil i and F is the
% final summed signal.
%
% USAGE:
% in = op_CSICombineCoils(in)
%
% INPUT:
% in            = input MRSI object for coils to be combined
% samplePoint(optional)   = samplePoint from fid to be used for phase
%               correction. Default to first point.
% phaseMap(optional)   = 3D matrix of phases for all coils and points. Dimensions are
%                       [num coils, num x, num y]
% weightMap(optional)   = 3D matrix of weights for all coils and points.
%                         Dimensions are the same as above.
%
% OUTPUT:
% in            = MRSI object with combined coils in fids and specs.

function [MRSIStruct, phaseMap, weightMap] = op_CSICombineCoils(MRSIStruct, ...
                                        samplePoint, phaseMap, weightMap)
    arguments
        MRSIStruct (1,1) struct
        samplePoint (1,1) double = 1
        phaseMap double = []
        weightMap double = []
    end
    checkArguments(MRSIStruct);

    [MRSIStruct, pastPermute, pastSize] = reshapeDimensions(MRSIStruct, {'t', 'coils', 'y', 'x'});

    data = getData(MRSIStruct);
    extraDimension = 4;
    if isempty(phaseMap)
        %phase map arranged by coils, y coordinate, x coordinate and extra
        %dimension

        %calculate the angle of each each coordinate for all the coils
        phaseMap = squeeze(angle(data(samplePoint, :, :, :, :)));
        phaseMap = mean(phaseMap, extraDimension);
    end
    %apply phase map to data
    data = applyMap(MRSIStruct, exp(-1i*phaseMap), data);

    if isempty(weightMap)
        %get weights from all positions and coils
        weights = squeeze(data(samplePoint, :, :, :, :));
        %Get the root sum squared value along coil dimension
        weightMap = normalize(weights, 1, 'norm', 2);
        weightMap = mean(weightMap, extraDimension);
    end
    %adding weights to each coil
    data = applyMap(MRSIStruct, weightMap, data);
    %add coils together
    MRSIStruct = combineCoilData(MRSIStruct, data, pastPermute, pastSize);

end

%check arguments to make sure function runs as expected
function checkArguments(in)
    %some pre condition checks and setting default values
    if in.flags.addedrcvrs == 1
        error('coils already combined!')
    end
    checkSpatialFT(in);
%    checkSpectralFT(in);
end

% apply phase or weight map to data. Maps are in dimensiosn (coils, y, x).
% Applies to each fid point and extra dimension.
function data = applyMap(MRSIStruct, map, data)
    % repmat maps to same legnth as time dimension
    tSize = getSizeFromDimensions(MRSIStruct, {'t'});
    mapWithTime = repmat(map, [ones(1, ndims(map)) tSize]);

    dimensions = 1:ndims(mapWithTime);
    timeFirstOrder = circshift(dimensions, 1);
    timeFirstMap = permute(mapWithTime, timeFirstOrder);
    data = data.*timeFirstMap;
end

% sets the weighted data from phase and weight maps to MRSI Struct and sums
% along coil dimension. Finally updates properties of the MRSI struct.
function MRSIStruct = combineCoilData(MRSIStruct, data, pastPermute, pastSize)
    coilDimension = getDimension(MRSIStruct, 'coils');
    % set weighted data back to structure
    MRSIStruct = setData(MRSIStruct, data);
    %permute back to original dimensions
    MRSIStruct = reshapeBack(MRSIStruct, pastPermute, pastSize);

    % sum along coil dimensions
    data = getData(MRSIStruct);
    data = squeeze(sum(data, coilDimension));
    MRSIStruct = setData(MRSIStruct, data);

    % update structure
    MRSIStruct = removeDimension(MRSIStruct, 'coils');
    MRSIStruct = setFlags(MRSIStruct, 'addedrcvrs', true);
end

author: CIC-methods

repo: FID-A

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


text:
function K = intmat(N, varargin)
%INTMAT   Spectral integration matrix.
%   K = INTMAT(N) returns the N x N integration matrix associated with the
%   Chebyshev spectral collocation method at second-kind Chebyshev points.
%
%   K = INTMAT(N, P) returns the N x N integration matrix of order P.
%
%   K = INTMAT(N, P, DOM) scales the integration matrix K to the domain
%   DOM. DOM should be a 1x2 vector.
%
%   K = INTMAT([M N]) returns the M x N first-order rectangular integration
%   matrix which maps from an N-point Chebyshev grid of the second kind to an
%   M-point Chebyshev grid of the same kind.
%
%   K = INTMAT([M N], P) returns an M x N rectangular integration matrix of
%   order P which maps from an N-point to an M-point Chebyshev grid, both of
%   second kind.
%
%   K = INTMAT([M N], P, DOM) returns the same K but scaled to the domain DOM.
%
% See also CUMSUMMAT, DIFFMAT, DIFFROW, INTROW.

% This code has been written for Aurentz and Trefethen,
% "Block operators and spectral discretizations", and is not
% yet fully compliant with standard Chebfun coding practices.

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

% Parse inputs
[m, n, p, dom] = parseInputs(N, varargin{:});

% Build Lagrange basis
K = chebfun(eye(n),dom);

% Integrate p times
K = cumsum(K,p);

% Evaluate at grid of size m
K = feval(K,chebpts(m,dom));

end

function [m, n, p, dom] = parseInputs(N, varargin)
% Parse the inputs to INTMAT.

p = 1;
dom = [-1 1];

if ( isscalar(N) )
    n = N;
    m = n;
else
    m = N(1);
    n = N(2);
end

for j = 1:numel(varargin)
    v = varargin{j};
    if ( isnumeric(v) )
        if ( isscalar(v) )
            p = v;
        else
            dom = v;
        end
    else
        error('CHEBFUN:intmat:unknown', ...
            'Unknown input of type %s.', class(v));
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Con = createControls(Robot)

% CREATECONTROLS  Create controls structure Con.
%   Con = CREATECONTROLS(ROBOT) generates a control structure Con. This
%   structure depends on the robot's motion model, which contains the
%   fields:
%       .u      nominal value of control
%       .uStd   Standard deviation of control noise
%       .U      Covariances matrix
%
%   The resulting structure needs to be updated during execution time with
%   data from one of these origins:
%       1. read from odometry sensors
%       2. generated by a control software
%       3. generated by the simulator.
%
%   See also CREATEROBOTS.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

switch Robot.motion

    case {'constVel'}

        Con.u    = [Robot.dv;deg2rad(Robot.dwDegrees)];
        Con.uStd = [Robot.dvStd;deg2rad(Robot.dwStd)];
        Con.U    = diag(Con.uStd.^2);

    case {'odometry'}

        Con.u    = [Robot.dx;deg2rad(Robot.daDegrees)];
        Con.uStd = [Robot.dxStd;deg2rad(Robot.daStd)];
        Con.U    = diag(Con.uStd.^2);

    otherwise
        error('Unknown motion model %s from robot %d.',Robot.motion,Robot.id);
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function w = load_file(request)
   %LOAD_FILE
   % uses the DATASOURCE class

   w = waveform; w = w();
   if isstruct(request)
      [ds, chanInfo, startTimes, endTimes, ~] = unpackDataRequest(request);
      w = load_objects_from_file(ds, 'waveform', chanInfo, startTimes, endTimes);
      w = w(:);
   else
      filen = request;

      if ~exist(filen,'file')
         disp(['unable to find file: ', filen]);
         return
      end

      %find out which variables exist within the file...
      fcontents = whos('-file',filen);
      variablesToLoad = ismember({fcontents.class},{'waveform', 'cell'});
      if ~any(variablesToLoad),
         %disp(['no items of the this class exist within the file:',filen]);
         return
      end
      stuff =load(filen,fcontents(variablesToLoad).name);

      fieldn = fieldnames(stuff);

      %Look through all variables that were loaded, looking for the desired
      %class.  If a cell is encountered, then look for occurrences of that class
      %within the cell (recursively).

      myObjectMask = false(size(fieldn));
      holder = cell(size(fieldn));
      for fieldidx = 1:numel(fieldn)
         thisfield =  fieldn{fieldidx};
         %grab all myObjects
         if isa(stuff.(thisfield),'waveform')
            holder(fieldidx) = {reshape(stuff.(thisfield),1,numel(stuff.(thisfield)))};
            myObjectMask(fieldidx) = true;
            stuff.(thisfield) = {};
         elseif isa(stuff.(thisfield),'cell')
            holder(fieldidx) = {myObjectsFromCells('waveform',stuff.(thisfield))};
            myObjectMask(fieldidx) = true;
         else
            stuff.(thisfield) = {};
         end
      end
      w = [holder{myObjectMask}];

   end
end

function mObj = myObjectsFromCells(searchClass, mycell)
   myObjectMask = false(size(mycell));
   for i=1:numel(mycell);
      if isa(mycell{i},searchClass),
         myObjectMask(i) = true;
         mycell(i) = {reshape(mycell{i},1,numel(mycell{i}))}; %make all myObjects 1xN
      elseif isa(mycell{i},'cell') %it's a cell, let's recurse
         mycell(i) = {myObjectsFromCells(mycell{i})};  %pull myObjects from the cell and bring them to this level.
      end
   end
   mObj= [mycell{myObjectMask}];
end
%}

function mObj = getFromCells(searchClass, mycell)
   % returns an 1xN array of objects
   searchFn = @(x) isa (x, searchClass);
   makeRows = @(x) reshape(x,1,numel(x));

   target = cellfun(searchFn, mycell);
   objs = cellfun(makeRows, mycell(target), 'uniformoutput', false);

   holdsCell = cellfun(@iscell, mycell);
   if any(holdsCell)
      mObj = getFromCells(searchClass, mycell{holdsCell}); %recurse
   else
      mObj = {};
   end
   mObj = [mObj{:} objs{:}];
end

author: geoscience-community-codes

repo: GISMO


text:
function dt6 = dtiTalNormalize(dt6file,destinationDir)
% dt6 = dtiTalNormalize(dt6file,destinationDir)
%
% Talirach Normalization:
% Takes in a dt6 filename and scales according to the talScale factors
% stored within.
%
% Outputs new file of form originalFileName_tal.mat in destination
% directory
%
% dtiTalNormalize uses mrAnatResliceSpm
%

load(dt6file);

% Get a bounding box in talairach space that is the same size as the input
% image.
origin  = anat.xformToAcPc\[0 0 0 1]';
origin  = origin(1:3)';
bb = [-anat.mmPerVox .* (origin-1) ; anat.mmPerVox .*(size(anat.img)-origin)]';

ts = anat.talScale;
ts.talScaleDir = 'tal2acpc';
% warp T1 to tal space
ts.outMat = inv(anat.xformToAcPc);
[anat.img,newT1Xform] = mrAnatResliceSpm(anat.img, ts, bb', [1 1 1], [7 7 7 0 0 0]);

ts.outMat = inv(anat.xformToAcPc*xformToAnat);
[b0,newB0Xform]  = mrAnatResliceSpm(b0, ts, bb', mmPerVox, [7 7 7 0 0 0]);
[dt6, newDt6Xform, def] = mrAnatResliceSpm(dt6, ts, bb', mmPerVox, [1 1 1 0 0 0]);
dt6(isnan(dt6)) = 0;

%PPD Correction
dt6 = dtiXformTensorsPPD(dt6,def,1);

anat.xformToAcPc = newT1Xform;
xformToAcPc = anat.xformToAcPc*xformToAnat;

b0(isnan(b0)) = 0;
b0 = int16(b0);
[subDir subFile x xx] = fileparts(dt6file);
save([destinationDir,subFile,'_tal'],'b0','xformToAcPc','xformToAnat','anat','notes','mmPerVox','dt6');
if(nargout==0)
    clear dt6;
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function h = polarerrorpatch(theta,r,r_err,color,varargin)

[isclosed,leftovers] = myparse_nocheck(varargin,'isclosed',true);

if isclosed,
  r = [r,r(1)];
  r_err = [r_err,r_err(1)];
  theta = [theta,theta(1)];
end

% bounds
plus_x = max(0,(r+r_err)).*cos(theta);
plus_y = max(0,(r+r_err)).*sin(theta);
minus_x = max(0,(r-r_err)).*cos(theta);
minus_y = max(0,(r-r_err)).*sin(theta);

% edge of patch
err_x = [plus_x,fliplr(minus_x)];
err_y = [plus_y,fliplr(minus_y)];

h = patch(err_x,err_y,color,leftovers{:});

author: kristinbranson

repo: JAABA


text:
%   CLASS GUI_About
% =========================================================================
%
% DESCRIPTION
%   class to manage the about window of goGPSz
%
% EXAMPLE
%   ui = GUI_About.getInstance();
%
% FOR A LIST OF CONSTANTs and METHODS use doc Core_UI

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, Giulio Tagliaferro, Eugenio Realini ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef GUI_About < handle

    properties (Constant, Access = 'protected')
        BG_COLOR = Core_UI.DARK_GREY_BG;
    end

    %% PROPERTIES SINGLETON POINTERS
    % ==================================================================================================================================================
    properties % Utility Pointers to Singletons
        log
        state
    end

    %% PROPERTIES GUI
    % ==================================================================================================================================================
    properties
        w_main      % Handle of the main window
        win         % Handle to this window
    end

    %% PROPERTIES STATUS
    % ==================================================================================================================================================
    properties (GetAccess = private, SetAccess = private)
    end

    %% METHOD CREATOR
    % ==================================================================================================================================================
    methods (Static)
        function this = GUI_About(w_main)
            % GUI_MAIN object creator
            this.init();
            this.openGUI();
            if nargin == 1
                this.w_main = w_main;
            end
        end
    end
    %% METHODS INIT
    % ==================================================================================================================================================
    methods
        function init(this)
            this.log = Core.getLogger();
            this.state = Core.getState();
        end

        function openGUI(this)
            % Main Window ----------------------------------------------------------------------------------------------

            win = figure( 'Name', 'About', ...
                'Visible', 'on', ...
                'DockControls', 'off', ...
                'MenuBar', 'none', ...
                'ToolBar', 'none', ...
                'NumberTitle', 'off', ...
                'Position', [0 0 680 540], ...
                'Resize', 'off');

            this.win = win;

            if isunix && not(ismac())
                win.Position(1) = round((win.Parent.ScreenSize(3) - win.Position(3)) / 2);
                win.Position(2) = round((win.Parent.ScreenSize(4) - win.Position(4)) / 2);
            else
                win.OuterPosition(1) = round((win.Parent.ScreenSize(3) - win.OuterPosition(3)) / 2);
                win.OuterPosition(2) = round((win.Parent.ScreenSize(4) - win.OuterPosition(4)) / 2);
            end

            try
                main_vb = uix.VBox('Parent', win, ...
                    'Padding', 5, ...
                    'BackgroundColor', Core_UI.DARKER_GREY_BG);
            catch
                log = Core.getLogger;
                log.setOutMode(1,[],0); % to plot a Warning I need to disable GUI and enable
                log.addError('Please install GUI Layout Toolbox (https://it.mathworks.com/matlabcentral/fileexchange/47982-gui-layout-toolbox)');
                open('GUI Layout Toolbox 2.3.4.mltbx');
                log.newLine();
                log.addWarning('After installation re-run goGPS');
                close(win);
                return;
            end
            top_bh = uix.HBox('Parent', main_vb);

            logo_GUI_About.BG_COLOR = Core_UI.DARK_GREY_BG;
            left_tbv = uix.VBox('Parent', top_bh, ...
                'BackgroundColor', logo_GUI_About.BG_COLOR, ...
                'Padding', 5);

            % Logo/title box -------------------------------------------------------------------------------------------

            logo_g = uix.Grid('Parent', left_tbv, ...
                'Padding', 5, ...
                'BackgroundColor', logo_GUI_About.BG_COLOR);

            logo_ax = axes( 'Parent', logo_g);
            logo_g.Widths = 64;
            logo_g.Heights = 64;
            [logo, transparency] = Core_UI.getLogo();
            logo(repmat(sum(logo,3) == 0,1,1,3)) = 0;
            logo = logo - 20;
            image(logo_ax, ones(size(logo)), 'AlphaData', transparency);
            logo_ax.XTickLabel = [];
            logo_ax.YTickLabel = [];
            axis off;

            Core_UI.insertEmpty(left_tbv, logo_GUI_About.BG_COLOR);
            left_tbv.Heights = [82 -1];

            % Title Panel -----------------------------------------------------------------------------------------------
            right_tvb = uix.VBox('Parent', top_bh, ...
                'Padding', 5, ...
                'BackgroundColor', logo_GUI_About.BG_COLOR);

            top_bh.Widths = [106 -1];

            title = uix.HBox('Parent', right_tvb, ...
                'BackgroundColor', logo_GUI_About.BG_COLOR);

            txt = this.insertBoldText(title, 'goGPS', 10, Core_UI.LBLUE, 'left');
            txt.BackgroundColor = logo_GUI_About.BG_COLOR;
            title_l = uix.VBox('Parent', title, 'BackgroundColor', GUI_About.BG_COLOR);
            title.Widths = [54 -1];
            Core_UI.insertEmpty(title_l, logo_GUI_About.BG_COLOR)
            txt = this.insertBoldText(title_l, ['- software V' Core.GO_GPS_VERSION], 8, [], 'left');
            txt.BackgroundColor = logo_GUI_About.BG_COLOR;
            title_l.Heights = [2, -1];

            % Disclaimer Panel -----------------------------------------------------------------------------------------------
            Core_UI.insertEmpty(right_tvb, logo_GUI_About.BG_COLOR)
            txt = this.insertText(right_tvb, {['This release is loosely based on the original goGPS software, most of the code has '...
                'been rewritten for better performances, automation, flexibility, but with the primary goal of quasi static processing. ' ...
                'Legacy version can still be found at:'], ...
                'https://github.com/goGPS-Project/goGPS_MATLAB/tree/legacy'}, 7, [], 'left');
            txt.BackgroundColor = logo_GUI_About.BG_COLOR;

            right_tvb.Heights = [20 3 -1];

            credits_bh = uix.HBox('Parent', main_vb, ...
                'Padding', 5, ...
                'BackgroundColor', GUI_About.BG_COLOR);
            credits_legacy = uix.VBox('Parent', credits_bh, ...
                'BackgroundColor', GUI_About.BG_COLOR);
            credits_new = uix.VBox('Parent', credits_bh, ...
                'BackgroundColor', GUI_About.BG_COLOR);

            % Old Credits -----------------------------------------------------------------------------------------------

            this.insertBoldText(credits_legacy, 'Legacy version', 9);
            bar_cont = uix.HBox('Parent', credits_legacy, ...
                'BackgroundColor', GUI_About.BG_COLOR);

            Core_UI.insertEmpty(bar_cont, GUI_About.BG_COLOR)
            Core_UI.insertHBarDark(bar_cont);
            Core_UI.insertEmpty(bar_cont, GUI_About.BG_COLOR)
            bar_cont.Widths = [-1 220 -1];

            this.insertBoldText(credits_legacy, 'Developers of the original core', 8);
            this.insertText(credits_legacy, {'Mirko Reguzzoni', 'Eugenio Realini'}, 7);
            Core_UI.insertEmpty(credits_legacy, GUI_About.BG_COLOR)

            this.insertBoldText(credits_legacy, 'Founding fathers', 8);
            this.insertText(credits_legacy, {'Alessio Dominioni', ...
                'Roberto Terruzzi', ...
                'Marco Colla', ...
                'Marcello Palmulli', ...
                'Michele Tettamanti'}, 7);

            Core_UI.insertEmpty(credits_legacy, GUI_About.BG_COLOR)
            this.insertBoldText(credits_legacy, 'Contributors (in temporal order)', 8);
            this.insertText(credits_legacy, {'Sara Lucca', ...
                'Lisa Pertusini', ...
                'Daniele Sampietro', ...
                'Luana Valentini', ...
                'Alemu Befkadu Nigussie', ...
                'Ivan Reguzzoni', ...
                'Stefano Caldera', ...
                'Damiano Triglione', ...
                'Giuliano Sirioni', ...
                'Antonio Herrera Olmo', ...
                'Andrea Gatti', ...
                'Hendy F. Suhandri', ...
                'Andrea Nardo', ...
                'Daisuke Yoshida', ...
                'Marco Negretti', ...
                'Yuya Iwaki', ...
                'Louis Osei-Poku', ...
                'Paiyuan Zhou'}, 7);

            credits_legacy.Heights = [22 5 18 14*2 5 18 14*5 5 18 14*17];

            % New Credits -----------------------------------------------------------------------------------------------

            this.insertBoldText(credits_new, 'Object oriented new version', 9);
            bar_cont = uix.HBox('Parent', credits_new, ...
                'BackgroundColor', GUI_About.BG_COLOR);

            Core_UI.insertEmpty(bar_cont, GUI_About.BG_COLOR)
            Core_UI.insertHBarDark(bar_cont);
            Core_UI.insertEmpty(bar_cont, GUI_About.BG_COLOR)
            bar_cont.Widths = [-1 220 -1];

            this.insertBoldText(credits_new, 'Developers of the new goGPS core', 8);
            this.insertText(credits_new, {'Andrea Gatti', 'Giulio Tagliaferro'}, 7);
            Core_UI.insertEmpty(credits_new, GUI_About.BG_COLOR)

            this.insertBoldText(credits_new, 'Other Contributors', 8);
            this.insertText(credits_new, {'Eugenio Realini', ...
                'Stefano Caldera', ...
                'Lisa Pertusini', ...
                'Andreas Krietemeyer', ...
                'Alice Bonfiglio', ...
                'Stefano Barindelli',...
                'Alessandra Mascitelli', ...
                }, 7);
            Core_UI.insertEmpty(credits_new, GUI_About.BG_COLOR)

            % IMPORTANT NOTE:
            % If you are a contributor and your name is not in this list feel free to add your name

            this.insertBoldText(credits_new, 'Acknoledgments', 8);
            this.insertText(credits_new, {'This release is funded by', 'GReD srl (www.g-red.eu)'}, 7);
            txt = this.insertText(credits_new, {'I dedicate my work on this software to my father, Tiziano Gatti, ', 'who taught me, among countless things, the importance of putting passion in everything you choose to do. I will never forget you,'}, 7, [], 'right'); txt.FontAngle = 'italic';
            txt = this.insertText(credits_new, {'Andrea Gatti'}, 7, [20 183 242]/255, 'right'); txt.FontAngle = 'italic';

            credits_new.Heights = [22 5 18 14*2 5 18 14*7 5 18 -1 14*3 28];

            % Manage dimension -------------------------------------------------------------------------------------------

            main_vb.Heights = [84 -1];

            %session_height = sum(left_bv.Children(2).Children(1).Heights);

            this.win.Visible = 'on';
            %uiwait(this.win);
        end
    end
    %% METHODS INSERT
    % ==================================================================================================================================================
    methods (Static)
        function txt = insertBoldText(parent, title, font_size, color, alignment)
            if nargin < 4 || isempty(color)
                color = Core_UI.WHITE;
            end
            if nargin < 5 || isempty(alignment)
                alignment = 'center';
            end
            txt = uicontrol('Parent', parent, ...
                'Style', 'Text', ...
                'String', title, ...
                'ForegroundColor', color, ...
                'HorizontalAlignment', alignment, ...
                'FontSize', Core_UI.getFontSize(font_size), ...
                'FontWeight', 'bold', ...
                'BackgroundColor', GUI_About.BG_COLOR);
        end

        function txt = insertText(parent, title, font_size, color, alignment)
            if nargin < 4 || isempty(color)
                color = Core_UI.WHITE;
            end
            if nargin < 5 || isempty(alignment)
                alignment = 'center';
            end
            txt = uicontrol('Parent', parent, ...
                'Style', 'Text', ...
                'String', title, ...
                'ForegroundColor', color, ...
                'HorizontalAlignment', alignment, ...
                'FontSize', Core_UI.getFontSize(font_size), ...
                'BackgroundColor', GUI_About.BG_COLOR);
        end

    end
    %% METHODS getters
    % ==================================================================================================================================================
    methods
    end

    %% METHODS EVENTS
    % ==================================================================================================================================================
    methods (Access = public)
    end
end

author: goGPS-Project


text:
function [A, C, b, f] = HALS_2d(Y, A, C, b, f, params)
%% Hierarchical alternating least square method for solving NMF problem
% Y = A*C + b*f

%input:
%   Y:      d1 X d2 X T, raw data. It will be reshaped to (d1*d2) X T in this
%   function
%   A:      (d1*d2) X K, initial value of spatial components
%   C:      K X T, initial value of temporal components
%   b:      (d1*d2) X nb, initial value of background spatial component
%   f:      nb X T, initial value of background temporal component
%   params: parameters used in this function.
%       bSiz:   blur size. A box kernel (bSiz X bSiz) will be convolved
%       with each neuron's initial spatial component, then all nonzero
%       pixels will be picked as pixels to be updated, and the rest will be
%       forced to be 0.
%       maxIter: maximum iteration of iterating HALS.

% Author: Pengcheng Zhou, Columbia University, based on a python
% implementation from Johannes Friedrich, Columbia University, 2015.

%% parameters
if isfield(params, 'maxIter'), maxIter = params.maxIter; else maxIter=5; end
if isfield(params, 'search_method'); method=params.search_method; else method='ellipse'; end
if and(isfield(params, 'bSiz'), strcmpi(method, 'dilate'))
    params.se = strel('disk', params.bSiz);
end
% search locations
IND = determine_search_location(A, method, params);

%% update spatial and temporal components neuron by neurons

for miter=1:maxIter
    %% update neurons
    Yac = Y - b*f;
    ind_del = find(std(A,0,1)==0);
    A(:, ind_del) = [];
    C(ind_del, :) = [];
    IND(:, ind_del) = [];
    %   temporal
    C = HALS_temporal(Yac, A, C, 5);

    ind_del = find(std(C,0,2)==0);
    A(:, ind_del) = [];
    C(ind_del, :) = [];
    IND(:, ind_del) = [];
    %   spatial
    A = HALS_spatial(Yac, A, C, IND, 5);

    %% update background
    Ybg = Y-A*C;
    % temporal
    f = HALS_temporal(Ybg, b, f, 5);
    % spatial
    b = HALS_spatial(Ybg, b, f, [], 5);
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [R2,R2_states] = tudaR2(X,Y,T,tuda,Gamma)
% Training explained variance per time point per each state
%
% R2_states is per state, R2 is for the entire model
% X and Y must be in matrix format (2D)
% All trials must have equal size, ie all elements of T must be equal
%
% Diego Vidaurre

if any(T(1)~=T)
    error(['tudaR2 requires that all trials have the same length, ' ...
        'i.e. that all elements in T are equal']);
end

[X,Y,T] = preproc4hmm(X,Y,T,tuda.train);

intercept = all(X(:,1)==1);

N = length(T); ttrial = T(1); p = size(X,2);
K = length(tuda.state); q = size(Y,2);

if intercept
    mY = repmat(mean(Y),size(Y,1),1); % one single mean value for the entire data
    mY = reshape(mY,[ttrial N q]);
else
    mY = zeros(ttrial,N,q);
end

mat1 = ones(ttrial,q);
Y = reshape(Y,[ttrial N q]);
Yhat = zeros(ttrial,N,q);

R2_states = zeros(ttrial,q,K);
beta = tudabeta(tuda);
for k = 1:K
    Yhatk = X * beta(:,:,k); %tuda.state(k).W.Mu_W(1:p,p+1:end);
    Yhatk = reshape(Yhatk,[ttrial N q]);
    if k > size(Gamma,2) % NESS model baseline state
        break%
        noGamma = prod(1-Gamma,2);
        Yhat = Yhat + Yhatk .* reshape(repmat(noGamma,1,q),[ttrial,N,q]);
    else
        Yhat = Yhat + Yhatk .* reshape(repmat(Gamma(:,k),1,q),[ttrial,N,q]);
    end
    if intercept
        mYk = mean(reshape(Y,[ttrial*N q]) .* repmat(Gamma(:,k),1,q));
        mYk = reshape(repmat(mYk,size(Y,1),1),[ttrial N q]);
    else
        mYk = zeros(ttrial,N,q);
    end
    ek = permute(sum((Yhatk - Y).^2,2),[1 3 2]); % ttrial x q
    e0k = permute(sum((mYk - Y).^2,2),[1 3 2]); % ttrial x q
    R2_states(:,:,k) = mat1 - ek ./ e0k ; % we do not compute the
end

e = permute(sum((Yhat - Y).^2,2),[1 3 2]); % ttrial x q
e0 = permute(sum((mY - Y).^2,2),[1 3 2]); % ttrial x q
R2 = mat1 - e ./ e0 ;

end

author: OHBA-analysis

repo: HMM-MAR


text:
function y = tapas_condhalluc_obs_sim(r, infStates, p)
% Simulates responses according to the condhalluc_obs model
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Inverse decision temparature beta is the only parameter
be = p;

% Prediction trajectory
mu1hat = infStates(:,1,1);

% Get true-positive rate corresponding to stimuli
tp = r.u(:,2);

% Calculate belief x using Bayes' theorem
x = tp.*mu1hat./(tp.*mu1hat + (1-mu1hat).^2);

% Belief is mu1hat in trials where there is no tone
x(find(tp==0)) = mu1hat(find(tp==0));

% Apply the logistic sigmoid to the inferred beliefs
prob = tapas_sgm(be.*(2.*x-1),1);

% Initialize random number generator
if isnan(r.c_sim.seed)
    rng('shuffle');
else
    rng(r.c_sim.seed);
end

% Simulate
y = binornd(1, prob);

return;

author: translationalneuromodeling

repo: tapas


text:
function opts=addCelListToStruct(opts,args)
%%ADDCELLLISTTOSTRUCT Given a cell array full of string-value pairs, such
%          as one might obtain using a varargin input to a function, take
%          the structure opts and set opts.stringName=value for all string
%          names and value pairs in args.
%
%INPUTS: opts A structure to which the string name-value pairs in args are
%             to be added. opts.stringName must already exist for all
%             possible string names. Typically, these will be set with
%             default values.
%        args A cell array where ars{2*i+1} is a string and args{2*i+2} is
%             the value associated with the string.
%
%OUTPUTS: opts The structure with the elements corresponding to string
%              names in args changed to the values that are also provided.
%
%This function is most useful if one stores default values for a function
%in a structure opts and then lets the user change individual entries with
%a varargin input.
%
%EXAMPLE:
%Consider the following structure, which holds default values where the
%second one is to be changed to 2:
% opts=[];
% opts.one=1;
% opts.two=20;
% args=cell(2,1);
% args{1}='two';
% args{2}=2;
% opts=addCelListToStruct(opts,args)
%The new opts variable has opts.one=1 and opts.two=2.
%
%February 2016 Hatim F. Alqadah, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if nargin < 2
        error('Invalid number of arguments');
    end

    if(isempty(args))
       return;
    end

    assert(isstruct(opts),'opts should be a struct');
    assert(iscell(args),'args should be a cell array');
    assert(isvector(args),'args should be a vector');
    n = length(args);
    if(mod(n,2)~=0)
        error('args should be given as field-value pairs');
    end

    %Get the option names as specifed in opts
    optnames=fieldnames(opts);
    for arg=1:2:n
        name = args{arg};
        if(~ischar(name))
             warning('MATLAB:invalid_fieldname','encountered a non-string option name');
        end
        %Look for a match in the field names
        idx = find(strcmpi(optnames,name),1);
        if(~isempty(idx))
            opts.(optnames{idx})=args{arg+1};
        else
            warning('MATLAB:invalid_option',['unrecognized option ',name,' skipping...']);
        end
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [pred, accuracy] = cosmo_crossvalidate(ds, classifier, partitions, opt)
% performs cross-validation using a classifier
%
% [pred, accuracy] = cosmo_crossvalidate(dataset, classifier, partitions, opt)
%
% Inputs
%   ds                  struct with fields .samples (PxQ for P samples and
%                       Q features) and .sa.targets (Px1 labels of samples)
%   classifier          function handle to classifier, e.g.
%                       @classify_naive_baysian
%   partitions          For example the output from nfold_partition
%   opt                 optional struct with options for classifier
%     .normalization    optional, one of 'zscore','demean','scale_unit'
%                       to normalize the data prior to classification using
%                       zscoring, demeaning or scaling to [-1,1] along the
%                       first dimension of ds. Normalization
%                       parameters are estimated using the training data
%                       and applied to the testing data.
%     .pca_explained_count   optional, transform the data with PCA prior to
%                            classification, and retain this number of
%                            components
%     .pca_explained_ratio   optional, transform the data with PCA prior to
%                            classification, and retain the components that
%                            explain this percentage of the variance
%                            (value between 0-1)
%    .check_partitions  optional (default: true). If set to false then
%                       partitions are not checked for being set properly.
%    .average_train_X   average the samples in the train set using
%                       cosmo_average_samples. For X, use any parameter
%                       supported by cosmo_average_samples, i.e. either
%                       'count' or 'ratio', and optionally, 'resamplings'
%                       or 'repeats'.
%
% Output
%   pred                Qx1 array with predicted class labels.
%                       elements with no predictions have the value NaN.
%   accuracy            scalar classification accuracy
%   test_chunks         Qx1 array with chunks of input dataset, if each
%                       prediction was based using a single classification
%                       step. Predictions with no or more than one
%                       classification step are set to NaN
%
% Examples:
%     % generate dataset with 3 targets and 4 chunks, first target is 3
%     ds=cosmo_synthetic_dataset('ntargets',3,'nchunks',4,'target1',3);
%     % use take-1-chunk for testing crossvalidation
%     partitions=cosmo_nfold_partitioner(ds);
%     classifier=@cosmo_classify_naive_bayes;
%     % run crossvalidation
%     [pred,accuracy]=cosmo_crossvalidate(ds, classifier, ...
%                                                       partitions);
%     % show targets, chunks, and predictions labels for each of the
%     % four folds
%     cosmo_disp({ds.sa.targets,ds.sa.chunks,pred},'threshold',inf)
%     %|| { [ 3    [ 1    [   3   NaN   NaN   NaN
%     %||     4      1        4   NaN   NaN   NaN
%     %||     5      1        5   NaN   NaN   NaN
%     %||     3      2      NaN     3   NaN   NaN
%     %||     4      2      NaN     5   NaN   NaN
%     %||     5      2      NaN     5   NaN   NaN
%     %||     3      3      NaN   NaN     3   NaN
%     %||     4      3      NaN   NaN     4   NaN
%     %||     5      3      NaN   NaN     5   NaN
%     %||     3      4      NaN   NaN   NaN     3
%     %||     4      4      NaN   NaN   NaN     4
%     %||     5 ]    4 ]    NaN   NaN   NaN     5 ] }
%     cosmo_disp(accuracy)
%     %||  0.917
%     %
%     % use take-2-chunks out for testing crossvalidation, LDA classifier
%     partitions=cosmo_nchoosek_partitioner(ds,2);
%     classifier=@cosmo_classify_lda;
%     % run crossvalidation
%     [pred,accuracy]=cosmo_crossvalidate(ds, classifier, ...
%                                                           partitions);
%     % show targets, chunks, and predictions labels for each of the
%     % four folds
%     cosmo_disp({ds.sa.targets,ds.sa.chunks,pred},'threshold',inf)
%     %|| { [ 3    [ 1    [   5     5     3   NaN   NaN   NaN
%     %||     4      1        4     4     4   NaN   NaN   NaN
%     %||     5      1        5     5     4   NaN   NaN   NaN
%     %||     3      2        3   NaN   NaN     3     3   NaN
%     %||     4      2        4   NaN   NaN     4     4   NaN
%     %||     5      2        5   NaN   NaN     4     5   NaN
%     %||     3      3      NaN     5   NaN     3   NaN     3
%     %||     4      3      NaN     4   NaN     4   NaN     4
%     %||     5      3      NaN     5   NaN     5   NaN     5
%     %||     3      4      NaN   NaN     3   NaN     3     3
%     %||     4      4      NaN   NaN     4   NaN     4     5
%     %||     5 ]    4 ]    NaN   NaN     5   NaN     3     3 ] }
%     cosmo_disp(accuracy)
%     %||   0.778
%     %
%     % as the example above, but (1) use z-scoring on each training set
%     % and apply the estimated mean and std to the test set, and (2)
%     % use odd-even partitioner
%     opt=struct();
%     opt.normalization='zscore';
%     partitions=cosmo_oddeven_partitioner(ds);
%     % run crossvalidation
%     [pred,accuracy]=cosmo_crossvalidate(ds, classifier, ...
%                                                        partitions, opt);
%     % show targets, predicted labels, and accuracy
%     cosmo_disp({ds.sa.targets,ds.sa.chunks,pred},'threshold',inf)
%     %|| { [ 3    [ 1    [ NaN     5
%     %||     4      1      NaN     4
%     %||     5      1      NaN     5
%     %||     3      2        3   NaN
%     %||     4      2        4   NaN
%     %||     5      2        5   NaN
%     %||     3      3      NaN     5
%     %||     4      3      NaN     4
%     %||     5      3      NaN     5
%     %||     3      4        3   NaN
%     %||     4      4        4   NaN
%     %||     5 ]    4 ]      3   NaN ] }
%     cosmo_disp(accuracy)
%     %||   0.75
%
% Notes:
%   - to apply this to a dataset struct as a measure (for searchlights),
%     consider using cosmo_crossvalidation_measure
%   - to average samples in the training set prior to training, use the
%     options provided by cosmo_average_samples prefixed by
%     'average_train_'. For example, to take averages of 5 samples, and use
%     each sample in the input approximately 4 times, use:
%           opt.average_train_count=5;
%           opt.average_train_resamplings=4;
%
% See also: cosmo_crossvalidation_measure, cosmo_average_samples
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<4,opt=struct(); end
    if ~isfield(opt, 'normalization'),
        opt.normalization=[];
    end
    if ~isfield(opt, 'check_partitions'),
        opt.check_partitions=true;
    end

    if ~isempty(opt.normalization);
        normalization=opt.normalization;
        opt.autoscale=false; % disable for {matlab,lib}svm classifiers
    else
        normalization=[];
    end

    if all(isfield(opt, {'pca_explained_count','pca_explained_ratio'}))
            error(['pca_explained_count and pca_explained_ratio are ' ...
                'mutually exclusive'])
    elseif isfield(opt, 'pca_explained_count');
        arg_pca='pca_explained_count';
        arg_pca_value=opt.pca_explained_count;
    elseif isfield(opt, 'pca_explained_ratio');
        arg_pca='pca_explained_ratio';
        arg_pca_value=opt.pca_explained_ratio;
    else
        arg_pca=[];
    end

    if opt.check_partitions
        cosmo_check_partitions(partitions, ds);
    end

    % see if samples in training set are to be averaged
    [do_average_train, average_train_opt]=get_average_train_opt(opt);

    train_indices = partitions.train_indices;
    test_indices = partitions.test_indices;

    npartitions=numel(train_indices);

    nsamples=size(ds.samples,1);

    % space for output (one column per partition)
    % the k-th column contains predictions for the k-th fold
    % (with values of NaNs if there was no prediction)
    pred=NaN(nsamples,npartitions);

    targets=ds.sa.targets;

    % process each fold
    for fold=1:npartitions
        train_idxs=train_indices{fold};
        test_idxs=test_indices{fold};
        % for each partition get the training and test data, store in
        % train_data and test_data
        train_data = ds.samples(train_idxs,:);
        train_targets = targets(train_idxs);

        if do_average_train
            % make a minimal dataset, then use cosmo_average_samples to
            % compute averages
            n_train=numel(train_idxs);
            ds_train=struct();
            ds_train.samples=train_data;
            ds_train.sa.chunks=ones(n_train,1);
            ds_train.sa.targets=train_targets;

            ds_train_avg=cosmo_average_samples(ds_train,average_train_opt);
            train_data=ds_train_avg.samples;
            train_targets=ds_train_avg.sa.targets;
        end

        test_data = ds.samples(test_idxs,:);

        % apply pca
        if ~isempty(arg_pca)
            [train_data,pca_params]=cosmo_map_pca(train_data,...
                    arg_pca,arg_pca_value);
            test_data=cosmo_map_pca(test_data,'pca_params',pca_params);
        end

        % apply normalization
        if ~isempty(normalization)
            [train_data,params]=cosmo_normalize(train_data,normalization);
            test_data=cosmo_normalize(test_data,params);
        end

        % then get predictions for the training samples using
        % the classifier, and store these in the k-th column of all_pred.
        p = classifier(train_data, train_targets, test_data, opt);

        pred(test_idxs,fold) = p;
    end

    % compute accuracies
    has_prediction_mask=~isnan(pred);
    correct_mask=bsxfun(@eq,targets,pred) & has_prediction_mask;

    accuracy=sum(correct_mask)/sum(has_prediction_mask);

function [do_average_train, average_train_opt]=get_average_train_opt(opt)
    persistent cached_opt;
    persistent cached_do_average_train;
    persistent cached_average_train_opt;

    if ~isequal(cached_opt,opt)
        cached_average_train_opt=struct();
        cached_do_average_train=false;
        prefix='average_train_';
        keys=fieldnames(opt);
        for k=1:numel(keys)
            key=keys{k};
            sp=cosmo_strsplit(key,prefix);
            if isempty(sp{1})
                % starts with average_train, so take the rest of the key
                % and flag cached_do_average_train
                cached_average_train_opt.(sp{2})=opt.(key);
                cached_do_average_train=true;
            end
        end

        cached_opt=opt;
    end

    do_average_train=cached_do_average_train;
    average_train_opt=cached_average_train_opt;

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%% SPARSE_BATCH_LOCAL uses the BATCH command to run the SPARSE code locally.
%
%  Discussion:
%
%    The SPARSE code is a function, so first we must write a script
%    called SPARSE_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2012
%
%  Author:
%
%    John Burkardt
%
  workers = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_BATCH_LOCAL\n' );
  fprintf ( 1, '  Run SPARSE_SCRIPT locally.\n' );
%
%  The BATCH command sends the script for execution.
%
  my_job = batch ( 'sparse_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'sparse_fun' }, ...
    'matlabpool', workers );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  Here is where we can use the LOAD command
%  to make the script's workspace available.
%  The SPARSE function returns a matrix called A_sparse.
%
  load ( my_job );
%
%  Get the size of the matrix, just to verify it is here.
%
  size ( A_sparse );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_BATCH_LOCAL\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
function i4vec_ascends_test ( )

%*****************************************************************************80
%
%% I4VEC_ASCENDS_TEST tests I4VEC_ASCENDS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  test_num = 6;
%
%  Each ROW of this definition is a COLUMN of the matrix.
%
  x_test = [ ...
    1, 3, 2, 4; ...
    2, 2, 2, 2; ...
    1, 2, 2, 4; ...
    1, 2, 3, 4; ...
    4, 4, 3, 1; ...
    9, 7, 3, 0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_ASCENDS_TEST\n' );
  fprintf ( 1, '  I4VEC_ASCENDS determines if an integer vector ascends.\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    x(1:n) = x_test(1:n,test);

    i4vec_print ( n, x, '  Test vector:' );

    fprintf ( 1, '  I4VEC_ASCENDS =  %d\n', i4vec_ascends ( n, x ) );

  end

  return
end


text:
% [y,npadl,npadu] = padgrab(x,padv,l1,u1,l2,u2,...)
function [y,npadl,npadu] = padgrab(x,padv,varargin)

% parse arguments
usagestring = 'Usage: padgrab(x,padv,l1,u1,l2,u2,...)';
nd1 = ndims(x);
if mod(numel(varargin),2) ~= 0 || isempty(varargin),
  error(usagestring);
end
isones = cellfun(@(x) x == 1, varargin(1:2:end-1)) & cellfun(@(x) x == 1, varargin(2:2:end));
sz = size(x);
isones(1:nd1) = isones(1:nd1) & sz == 1;
nd = find(~isones,1,'last');
if isempty(nd),
  nd = 1;
end
varargin = varargin(1:2*nd);
isrowvec = nd == 1 && nd1 == 2 && size(x,1) == 1;
if isrowvec,
  varargin = [{1,1},varargin];
  nd = 2;
end
if mod(nd,1) ~= 0 || nd > nd1,
  error(usagestring);
end
l = cell2mat(varargin(1:2:end));
u = cell2mat(varargin(2:2:end));
%if any(u < l),
%  error(usagestring);
%end
%sz = size(x);
if nd < nd1,
  x = reshape(x,[sz(1:nd-1),prod(sz(nd:end)),1]);
  sz = size(x);
end

l1 = max(min(l,sz),1);
u1 = max(min(u,sz),1);
a = cell(1,nd);
aisempty = false(1,nd);
for i = 1:nd,
  if sz(i) <= 0 || l(i) > u(i) || u(i) < 1 || l(i) > sz(i),
    a{i} = [];
    aisempty(i) = true;
  else
    a{i} = l1(i):u1(i);
  end
end
y = x(a{:});
npadl = l1-l+aisempty;
npadu = u-u1;
npadl = npadl + (npadu < 0).*npadu;

idxlow = u < 1;
if any(idxlow),
  npadl(idxlow) = u(idxlow) - l(idxlow) + 1;
  npadu(idxlow) = 0;
end
idxhigh = l > sz;
if any(idxhigh),
  npadu(idxhigh) = u(idxhigh) - l(idxhigh) + 1;
  npadl(idxhigh) = 0;
end

if any(npadl > 0),
  y = padarray(y,npadl,padv,'pre');
end
if any(npadu > 0),
  y = padarray(y,npadu,padv,'post');
end

% if dotranspose,
%   y = y';
% end


text:
function boundary_segment_num = p05_boundary_segment_num ( )

%*****************************************************************************80
%
%% P05_BOUNDARY_SEGMENT_NUM counts the boundary segments in problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Output, integer BOUNDARY_SEGMENT_NUM, the number of boundary segments.
%
  boundary_segment_num = 1;

  return
end


text:
function pass = test_functionForm

%% Building blocks
dom = [-2 2];
I = operatorBlock.eye(dom);
D = operatorBlock.diff(dom);
Z = operatorBlock.zeros(dom);
x = chebfun('x', dom);
u = sin(x.^2);
U = operatorBlock.mult(u);

%% Operator instantiations
eyeop = toFunction(I);
err(1) = norm( eyeop(u) - u );
diffop = toFunction(D);
TwoX = diffop(x.*x);   % should be the chebfun 2*x
TwoXAgain = D*(x.*x);
err(2) = norm( 2*x - TwoX );
multop = toFunction(U);
err(3) = norm( multop(x+1) - u.*(x+1) );   % should be zero
% shorthand notation, U*chebfun
err(4) = norm( multop(x+1) - U*(x+1) );

%%

pass = err < 1e-14;

end


text:
function ob = gg_Gnufft(om,M,Msp,R)
%function ob = gg_Gnufft(om, M, Msp, R)
%
% Object constructor for Leslie Greengard's 2-dimensional NUFFT using separable
% gaussian interpolator.
%
% in
%	om [Min,2]	"digital" frequencies in radians
%	M [1]		image dimensions (N1,N2,...,Nd)
%	Msp [1]		# of neighbors used
%	R [1]		oversampling ratio
% out
%	ob              the NUFFT object
%
% A type 1 NUFFT is evaluated using x = G' * b
% A type 2 NUFFT is evaluated using b = G * x
%
% Copyright 2008-6-6	Will Grissom and Jeff Fessler	The University of Michigan

% use greengard-supplied tau
tau = 1/M/M*pi/R/(R-0.5)*Msp; % width of gaussian

% initialize nufft
st = gg_Gnufft_init(om,M,Msp,R,tau);

% build an object with it.
ob = Fatrix([size(om,1) M*M], st, 'forw', @gg_Gnufft_forw, 'back', @gg_Gnufft_back, 'caller', mfilename);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function scan_data = dk8rd(ufilename,plotmode,verbose)
%DK8RD - load Dektak 8 surface profilometer data file
%
% SCAN_DATA = DK8RD(UFILENAME,PLOTMODE,VERBOSE)
%
%DK8RD reads the data from file UFILENAME.txt, and and returns a 2-column
% array SCAN_dATA:
%  (:,1) is the horizontal position data in microns
%  (:,2) is the vertical profile data in microns
%
% Options:
% If VERBOSE == 0, no messages will be displayed.
% If VERBOSE == 1, the data summary will be displayed. (default)
% If VERBOSE == 2, the file header and data summary will be displayed.
% If PLOTMODE == 0, the data will not be displayed.
% if PLOTMODE == 1, the data will be displayed in a new figure. (default)
%
% NOTE: The function can only read ASCII text data files, not binary data
%  files. If you supply a filename without an extension, the extension
%  ".txt" is appended.
%
% Examples:
%
%  If the data is stored in a text file called "mydata.txt", then:
%
%      dk8rd mydata;
%
%  will load the file, plot the data, and print a summary. You can also use:
%
%      dk8rd;
%
%  and you will be prompted for the file name. The command:
%
%      profile_data=dk8rd('mydata',0,1);
%
%  will open the file mydata.txt, plot the data,
%   and save the data to the variable profile_data in the workspace
%   without printing anything in the command window.
%
%
% M.A. Hopcroft
%  mhopeng at gmail dot com
%
% The license governing use of this code is in the accompanying file
%  "license.txt".
%
versionstr='dk8rd 1.3';

% MH Jul2013
% v1.3  add uifileget
% MH Sep2012
% v1.2  updates for modern matlab
% MH MAY2007
% v1.1  re-order input (verbose last)
%       filename w/ or w/o .txt
% MH OCT2006
% v1.0 out of beta
%
%#ok<*NBRAK>

% handle input arguments - set defaults
% default to verbose mode 1 (summary only)
if nargin < 3
    verbose = 1;
end

% default to plot mode 1
if nargin < 2
    plotmode = 1;
end

% get interactive if not being called by another program
if (nargout == 0)
    if (verbose >= 1 && nargin == 0)
        clc;
        fprintf(1, '\n\n  %s\n  %s\n\n', versionstr, 'Load Dektak 8 data files');
    end
end

%%%%
% open user's file
if nargin >= 1 % if the user already specified a file, use it
    userfilename=ufilename;

    % check to see if file exists
    cfid=fopen(userfilename);
    if cfid== -1
        % try with ".txt"
        userfilename=strcat(userfilename,'.txt');
        cfid=fopen(userfilename);
        if cfid== -1
            disp(pwd);
            fprintf(1,'ERROR: data file ["%s"] not found in current directory',userfilename);
            scan_data=[];
            return
        end
        fclose(cfid);
    end

else % if user did not specify the file name, ask for one
    [ufilename, ufilepath]=uigetfile('.txt','Select the Dektak 8 data file');
    if isequal(ufilename,0)
        scan_data=[];
        return
    end
    userfilename=fullfile(ufilepath,ufilename);
    %userfilename=input(' Enter the name of the Dektak 8 .txt file: ','s');
end

% open the file
userfile=fopen(userfilename);
[fpath, fname, fext]=fileparts(userfilename);

%%%%
% process the file header
% the first 44 lines are header information
%  read them in and get some data: scan length, number of points, etc.
% display the header for plotmode 2 only
% note the the Dektak 8 file format is very similar to the Dektak 3
if verbose == 2
    fprintf(1,'\n\n %s\n %s %s\n\n %s\n\n', versionstr, 'Filename:', fname, 'File header:');
end

% read header lines, if they are useful, save the info
while 1
    % get the first word on each line and check to see if it is important
    % the line identifier is followed by a ":", so break the line into
    %  two parts based on the :
    %  if the data is useful, scan it from the second part of the line

    % get a line from the file
    linein=fgetl(userfile);
    % handle spurious carriage returns
    if isempty(linein), linein=fgetl(userfile); end
    % figure out where the ":" is so we can divide the line
    id_pos=find(linein==':');
    % some lines may be malformed or empty
    if ~isempty(id_pos)
        id_pos=id_pos(1); % some lines have extra :'s
        % line_id is the identifier
        line_id=linein(1:id_pos);
        % line_data is everything following the identifier
        line_data=linein(id_pos+1:end);
    else
        % line_id is the identifier
        line_id=linein;
    end

    % check to see if the line identifier is interesting, if so, scan
    % the data from line_data. exit when SCALED DATA is reached

    if strcmp(line_id,'Sclen:')==1
        scan_length=sscanf(line_data,'%g',1);
    elseif strcmp(line_id,'Units:')==1
        headerUnits=sscanf(line_data,'%s',1); % check for unusual "units" value
        if ~strcmpi(headerUnits,'SCI')
            fprintf(1,'dk8rd: WARNING: Header Units value is "%s".\n  Data is assumed to have units of Angstroms. You should verify this.\n',headerUnits);
        end
    elseif strcmp(line_id,'NumPts:')==1
        num_points=sscanf(line_data,'%g',1);
    elseif strcmp(line_id,'Hsf:')==1
        stepsize=sscanf(line_data,'%g',1);
        % note re force- the force output appears to be messed up
        %  in the data files I have seen, so the input is restricted
        %  to the first four characters of the line
    elseif strcmp(line_id,'Force:')==1
        stylus_force=sscanf(line_data(1:4),'%d',1);
    elseif strcmp(line_id,'Date:')==1
        scan_date=sscanf(line_data,'%s',1);
        scan_time=sscanf(line_data,'%*s %s',1);
    elseif strcmp(line_id,'Tspeed:')==1
        scan_speed=sscanf(line_data,'%s',1);
    elseif strcmp(line_id,'SCALED DATA:')==1
        break % exit loop
    end

    % print the line to the command window if verbose mode 2 is set
    if verbose == 2
        % avoid printing empty lines
        if ~isempty(linein)
            fprintf(1,'%s\n',linein);
        end
    end
end

% % the position of the data points is now known
% stepsize=scan_length/num_points;
% % if horiz_res ~= stepsize, error('error reading file: horiz_res ~= stepsize'); end
% scan_data(:,1)=[stepsize:stepsize:scan_length]';

%%%%
% get data from file
%  data is a simple list of vertical position in Angstroms
%   (formatted in 5 columns of text)

% create the x axis
scan_data(:,1)=[stepsize:stepsize:scan_length]';

% read the data points
[scan_data(:,2),count]=fscanf(userfile,'%f',inf);
if count~=num_points
    fprintf(1,'dk8rd: WARNING: The number of data points (%g) does not equal the value in the file header (%g).\n  The file may be corrupted.',...
        count,num_points);
end
% scale data to um from Angstroms
scan_data(:,2)=scan_data(:,2)./10000;

% close the file
fclose(userfile);

% plot
if plotmode == 1

    figure;
    plot(scan_data(:,1),scan_data(:,2),'.-');
    title(['dk8rd: ' fname],'Interpreter','none','FontSize',16,'FontWeight','bold');
    xlabel('Scan Length (\mum)','FontSize',14);
    ylabel('Profile (\mum)','FontSize',14);
    grid on;

end

% print summary data for the user
if verbose >= 1

    fprintf(1,'\n');
    fprintf(1,' %s\n', 'Data Summary:');
    fprintf(1,'\n');
    fprintf(1,' %s %s\n', 'File name:', [fname fext]);
    fprintf(1,' %s %s %s\n','Scan Date:',scan_date,scan_time);

    fprintf(1,'\n');

    fprintf(1,' %s %s\n', 'Scan Speed:', scan_speed);
    fprintf(1,' %s %d\n', 'Number of data points read:', count);
    fprintf(1,' %s %d %s\n', 'Scan Length:', scan_length, 'um');
    fprintf(1,' %s %g %s\n', 'Scan Step:', stepsize, 'um');
    fprintf(1,' %s %g %s\n', 'Nominal Stylus Force:', stylus_force, 'mg');
    fprintf(1,'   %s\n', '[Note: stylus force is incorrectly reported by beta versions of Dektak 8 software]');

    fprintf(1,'\n');

end

author: Sable