text:
  function back = my_cbct_back(proj, cg, ig, varargin)
% function back = cbct_back(proj, cg, ig, varargin)
%|
%| cone-beam backprojector for feldkamp.m
%| in
%|	proj	[ns nt na]	cone-beam projection views
%| option
%|	'use_mex' 1|2|3		default: 1 mex with loop in mex
%|					2 mex with loop in matlab
%|					3 mex with "st" data order (slower)
%|					0 no mex; use matlab
%|	'ia_skip'		default: 1
%|	'offset_source'		default: 0
%|	'nthread'		default: jf('ncore')
%| out
%|	back	[nx ny nz]	back projection result
%|
%| Copyright 2004-8-28 Jeff Fessler, University of Michigan

if nargin == 1 && streq(proj, 'test'), cbct_back_test, return, end
if nargin < 3, help(mfilename), error(mfilename), end

arg.ia_skip = 1;
arg.use_mex = 1;
arg.offset_source = 0;
arg.nthread = jf('ncore');
arg.back_call = @fdk_mex_call;
arg.scale_dang = true;
arg = vararg_pair(arg, varargin);

if arg.use_mex
	if cg.pitch ~= 0 || any(cg.source_zs ~= 0)
		fail('helix not yet supported')
	end
	back = cbct_back_mex(proj, cg.ns, cg.nt, cg.na, ...
		cg.ds, cg.dt, cg.offset_s, cg.offset_t, arg.offset_source, ...
		cg.dsd, cg.dso, cg.dfs, cg.orbit, cg.orbit_start, ...
		ig.mask_or, ig.nz, ig.dx, ig.dy, ig.dz, ...
		[ig.offset_x ig.offset_y ig.offset_z], ...
		arg.ia_skip, arg.scale_dang, ...
		arg.use_mex, arg.nthread, arg.back_call);

else

	back = cbct_back_mat(proj, cg, cg.ns, cg.nt, cg.na, ...
		cg.ds, cg.dt, cg.offset_s, cg.offset_t, arg.offset_source, ...
		cg.dsd, cg.dso, cg.dfs, cg.orbit, cg.orbit_start, ...
		cg.source_zs, ...
		ig.mask_or, ig.nz, ig.dx, ig.dy, ig.dz, ...
		[ig.offset_x ig.offset_y ig.offset_z], ...
		arg.ia_skip, arg.scale_dang);
end

end % cbct_back()

%
% cbct_back_mex()
% mex backprojector
%
function img = cbct_back_mex(proj, ns, nt, na, ...
	ds, dt, offset_s, offset_t, offset_source, ...
	dsd, dso, dfs, orbit, orbit_start, ...
	mask, nz, dx, dy, dz, offset_xyz, ia_skip, scale_dang, ...
	use_mex, nthread, back_call)

[nx ny] = size(mask);
ia_list = 1:ia_skip:na;
betas = deg2rad(orbit_start + orbit * [0:na-1] / na); % [na] source angles
proj = single(proj);
nthread = int32(nthread);

switch use_mex
case 1 % loop in mex
	proj = permute(proj, [2 1 3]); % sta -> tsa
	img = back_call('fdk,ts,back', ...
			int32([nx ny nz]), [dx dy dz], ...
			offset_xyz, uint8(mask), ...
			dso, dsd, dfs, [ds dt], [offset_s offset_t], ...
			proj(:,:,ia_list), betas(ia_list), nthread);
%	img = double6(img);
	img = permute(img, [2 3 1]); % zxy -> xyz

case 2 % loop in matlab (for testing)
	proj = permute(proj, [2 1 3]); % sta -> tsa
	img = 0;
	for ia=ia_list
		ticker(mfilename, ia, na)
		% note: 2006-5-30: replaced -dy with dy
		tmp = back_call('fdk,ts,back', ...
			int32([nx ny nz]), [dx dy dz], ...
			offset_xyz, uint8(mask), ...
			dso, dsd, dfs, [ds dt], [offset_s offset_t], ...
			proj(:,:,ia), betas(ia), nthread);
%		tmp = double6(tmp);
		img = img + tmp;
	end
	img = permute(img, [2 3 1]); % zxy -> xyz

case 3 % fdk,st (for testing only - slower!)
	warn 'fdk,st does not handle projection view edges completely'
	img = back_call('fdk,st,back', ...
			int32([nx ny nz]), [dx dy dz], ...
			offset_xyz, uint8(mask), ...
			dso, dsd, dfs, [ds dt], [offset_s offset_t], ...
			proj(:,:,ia_list), betas(ia_list), nthread);
%	img = double6(img);

otherwise
	fail 'bug'
end

if scale_dang % final "\der angle" scale:
	img = (0.5 * deg2rad(abs(orbit)) / (na/ia_skip)) * img;
end

end % cbct_back_mex()

%
% cbct_back_mat()
% matlab backprojector (slower)
%
function img = cbct_back_mat(proj, cg, ns, nt, na, ...
	ds, dt, offset_s, offset_t, offset_source, ...
	dsd, dso, dfs, orbit, orbit_start, ...
	source_zs, ...
	mask, nz, dx, dy, dz, offset_xyz, ia_skip, scale_dang)

if any(source_zs ~= 0)
	warn('helix not yet tested')
end

[nx ny] = size(mask);
betas = deg2rad(orbit_start + orbit * [0:na-1] / na); % [na] source angles

% precompute as much as possible
wx = (nx-1)/2 + offset_xyz(1);
wy = (ny-1)/2 + offset_xyz(2);
wz = (nz-1)/2 + offset_xyz(3);
[xc yc] = ndgrid(([0:nx-1] - wx) * dx, ([0:ny-1] - wy) * dy);
zc = ([0:nz-1] - wz) * dz;

if 0 % limit back-projection to FOV?  removed 2008-10-9
	rr = sqrt(xc.^2 + yc.^2); % [nx,ny]
	smax = ((ns-1)/2-abs(offset_s)) * ds; % maximum detector s coordinate

	if isinf(dfs)
		gamma_max = atan(smax/dsd);
	elseif dfs == 0
		gamma_max = smax / dsd;
	end

	rmax = dso * sin(gamma_max);
	mask = mask & (rr < rmax);
end
clear wx wy wz rr smax rmax

xc = xc(mask); % [np] pixels within mask
yc = yc(mask);

ws = (ns+1)/2 + offset_s; % trick: +1 because matlab starts from 1
wt = (nt+1)/2 + offset_t;

% loop over slices
img = zeros([size(mask) nz]);
sdim = [ns+3 nt+3]; % trick: extra zeros saves indexing in loop
proj1 = zeros(sdim);
ticker reset

% Greg Handy stuff:
num_turns = orbit/360;
halfNumAngles = round((na/num_turns)/2);
myPitch = cg.pitch * cg.nt * dso / dsd * cg.dt;
h = myPitch / (2*pi);
dist = .5 * myPitch;
zindex = 1;

for iz=1:nz

	% Greg Handy stuff:
	% Just trying to limit the number of angles used
	lower_limit = zc(iz) - dist;

	% enter into the acceptable range for the z-slices
	while zindex < size(cg.source_zs,1)+1 && cg.source_zs(zindex) < lower_limit
		zindex = zindex + 1;
	end
	lowerA = zindex - halfNumAngles;
	upperA = zindex + halfNumAngles;

	% a check to prevent error; images without enough views will be poor
	if lowerA <= 0
		lowerA = 1;
	end
	if upperA > na
		upperA = na;
	end

	% loop over each projection angle
	img2 = 0;
	for ia=lowerA:ia_skip:upperA
		ticker(mfilename, [iz ia], [nz na])
		beta = betas(ia);

		x_beta = +xc * cos(beta) + yc * sin(beta);
		y_betas = dso - (-xc * sin(beta) + yc * cos(beta));

		% detector indices
		if isinf(dsd) || isinf(dso) % par
			mag = ones(size(y_betas));
		else
			mag = dsd ./ y_betas;
		end

		if isinf(dfs) ... % flat
			|| isinf(dsd) || isinf(dso) % par
			sprime = mag .* x_beta;
		elseif (dfs == 0) % arc
			r_loop = x_beta - offset_source;
			sprime = dsd * atan2(r_loop, y_betas);
		end

		tprime = mag * (zc(iz) - source_zs(ia));
		bs = sprime / ds + ws;
		bt = tprime / dt + wt;

		% bi-linear interpolation:
		is = floor(bs); % left bin
		it = floor(bt);

		wr = bs - is;	% left weight
		wl = 1 - wr;	% right weight
		wu = bt - it;	% upper weight
		wd = 1 - wu;	% lower weight

		ibad = (is < 0) | (is > ns) | (it < 0) | (it > nt);
		is(ibad) = ns+1; % trick! point at harmless zeros
		it(ibad) = nt+1;

		proj1(1+[1:ns],1+[1:nt]) = proj(:,:,ia); % trick: left side
		p1 =	wl .* proj1(sub2ind(sdim, is+1,it+1)) + ...
			wr .* proj1(sub2ind(sdim, is+2,it+1));
		p2 =	wl .* proj1(sub2ind(sdim, is+1,it+2)) + ...
			wr .* proj1(sub2ind(sdim, is+2,it+2));

		p0 = wu .* p1 + wd .* p2; % vertical interpolation

		if isinf(dfs) ... % flat
			|| isinf(dsd) || isinf(dso) % par
			p0 = p0 .* mag.^2; % back-projection weighting for flat
		elseif dfs == 0 % arc
			p0 = p0 .* (dsd.^2) ./ (r_loop.^2 + y_betas.^2);
		end

		img2 = img2 + p0;
	end % ia

	img(:,:,iz) = embed(img2, mask);
end % iz

if scale_dang % final "\der angle" scale:
	img = (0.5 * deg2rad(abs(orbit)) / (na*ia_skip)) * img;
end

end % cbct_back_mat()

%
% fdk_mex_call()
%
function out = fdk_mex_call(varargin)
if exist('fdk_mex') == 3
	printm 'using fdk_mex'
	out = fdk_mex(varargin{:});
elseif exist('jf_mex') == 3
	out = jf_mex(varargin{:});
else
	fail('bug: neither fdk_mex nor jf_mex found')
end
end % fdk_mex_call()

%
% cbct_back_test()
% compare various versions to check consistency
%
function cbct_back_test

down = 8; % fast test
% todo: test parallel, flat, arc
%cg = ct_geom('ge1', 'nt', 800, 'na', 50*down, 'dsd', inf, 'down', down);
cg = ct_geom('ge1', 'nt', 800, 'na', 50*down, 'down', down);
%cg.dt = -cg.dt; % todo: handle this case
ig = image_geom('nx', 512, 'ny', 496, 'nz', 480, 'fov', 500, ...
	'mask', 'all-but-edge', 'down', down);

ell = [ ... % somewhat realistic phantom object
	[30 10 10	150 150 280	0 0 1000]; % 30cm diam
	[80 10 10	50 50 30	0 0 300]; % bone-like inserts
	[-10 -40 75	40 40 40	0 0 600];
	[-10 80 -20	30 30 30	0 0 600];
];

if 1
	proj = ellipsoid_proj(cg, ell);
	proj = fdk_filter(proj, 'ramp', cg.dsd, cg.dfs, cg.ds);
else
	proj = cg.zeros;
	proj(:,:,9) = 1;
end
% im clf, im(proj, 'true projections'), cbar, return

if 0 % zero outer edges of projection for comparing ,st to ,ts
	proj([1 cg.ns], :, :) = 0;
	proj(:, [1 cg.nt], :) = 0;
	f.compare_st = true;
else
	f.compare_st = false;
end

	ia_skip = 3; % stress test and makes it faster too
	args = {proj, cg, ig, 'ia_skip', ia_skip, 'scale_dang', 0};
if 1
	back0 = cbct_back(args{:}, 'use_mex', 0);
%	im(back0), return
end

if exist('jf_mex') == 3
	for ii = 1:3
		use_mex = ii;
		printm('testing jf_mex with use_mex=%d', use_mex)
		back1{ii} = cbct_back(args{:}, ...
			'use_mex', use_mex, 'back_call', @jf_mex);
	end
	printm('mpd mex1 vs mat: %g%%', max_percent_diff(back1{1}, back0))
	printm('mpd mex1 vs mex2: %g%%', max_percent_diff(back1{1}, back1{2}))
	if f.compare_st % only if edges are zeroed
	printm('mpd mex1 vs mex3: %g%%', max_percent_diff(back1{1}, back1{3}))
	end
end

if exist('fdk_mex') == 3
	printm('found fdk_mex:\n %s', which('fdk_mex'))
	back2 = cbct_back(args{:}, 'use_mex', 1, 'back_call', @fdk_mex);
	printm('mpd jf_mex vs fdk_mex: %g%%', max_percent_diff(back1{1}, back2))

	if 1
		im plc 1 3
		iz = 1:ig.nz;
		im(1, back1{1}(:,:,iz))
		im(2, back2(:,:,iz))
		err = back2(:,:,iz) - back1{1}(:,:,iz);
		im(3, err)
	end
end

end % cbct_back_test()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [cg, psg] = crandn(rgau,m)
% Generate correlated Gaussian sequences by Fourier synthesis.
%
% Input parameters:
% rgau = correlation function - length n/2
% m    = number of realisations
%
% Output:
% cg   = m x n matrix containing m sequences of n correlated variates from
%        a zero mean, unit variance normal distribution
% psg  = input power spectrum (Fourier transform of correlation function)
%
% Note: Since this uses the fast Fourier transform, it will be fastest if
% n is a power of 2.
%
% S Bocquet 12 August 2008  Tested with MATLAB Version 7.6.0.324 (R2008a)
%
if ~isvector(rgau)
  error('crandn:vector','Correlation function must be a vector')
end
if ((rgau(1) ~= 1) || (max(abs(rgau)) > 1))
  error('crandn:cfnorm','Correlation function not properly normalised');
end
if abs(rgau(end)) > 0.01
  warning('crandn:cfdecay','Correlation function does not go to zero')
end
n = 2*length(rgau);
rgau2 = [rgau(1:end),0,rgau(end:-1:2)]; % two sided correlation function
psg = abs(fft(rgau2)); % Fourier transform of correlation function
nrm = mean(psg)/n; % normalisation
h = sqrt(psg/nrm); % normalised filter function
cg = zeros(n,m);
for i=1:m
  if mod(i,2)==1
    g = complex(randn(1,n),randn(1,n));
    cg2 = ifft(g.*h); % correlated Gaussian process
    cg(:,i) = real(cg2);
  else
    cg(:,i) = imag(cg2);
  end
end
end

author: Sable

repo: mcbench-benchmarks


text:
%% housekeeping
clear
clc
close all
%% "rise" the model
m=rise('fs2000_rise');
%% pushing the particular vector (different from the initial estimation point)
% Note we have to declare the regime of the parameter

calibration=struct();
	calibration.alp=0.330;
	calibration.bet=0.990;
	calibration.gam=0.003;
	calibration.mst=1.011;
	calibration.rho=0.700;
	calibration.psi=0.787;
	calibration.del=0.020;
	calibration.sig_a=0.014;
	calibration.sig_m=0.005;
m=set(m,'parameters',calibration);

%% solve the model and
m=solve(m);

%% print the solution

m.print_solution

disp('This should be identical to the dynare solution (see folder dynare_version)')

%% lets do some estimation: we need the data
% get the data from the dynare_version folder

run dynare_version\fsdat_simul

% rise needs data to be passed as a ts object with dates and
% names of the variables. So we proceed to constructing the database

% we use the same start date as in the Schorfheide paper even though the
% data are not the same
startdate='1950q1';
database=ts(startdate,[gp_obs,gy_obs],{'gp_obs','gy_obs'});

%% pass the data to the rise object

% 192 observations are used in dynare
end_date= obs2date(startdate,192);
% the loglinear option of dynare implies:
% 1- we have to take the log of our data.
% 2- we have to exponentiate the corresponding variables in the rise model
% file. Dynare does this by taking the log of the steady state during
% estimation and this does not ring too transparent to me.

% In rise, we can just take the log of the whole database. This is
% what we do in passing the data. But we also need to tell ts
% that the logged variables should have the same names as the original
% ones. In order to do that, we add set the flag to true when taking the
% log.
vnames=database.varnames;
database=log(database);
database.varnames=vnames;
m=set(m,'data',database,'estim_end_date',end_date);

%% estimate the model
profile off
profile on
m=estimate(m);%,'optimizer',@csminwellwrap,'debug',true
profile off
profile viewer

%% do posterior simulation
[objective,lb,ub,x0,SIG]=pull_objective(m);

SIG=utils.cov.nearest(SIG);

draws_mcmc = 1000; % number of parameter draws through MCMC.
ndraws_burnin = floor(0.1*draws_mcmc);
mcmc_options=struct('burnin',ndraws_burnin,'N',draws_mcmc,'thin',1,...
    'nchain',2);
Results=mh_sampler(objective,lb,ub,mcmc_options,x0,SIG);

%% update the description of the parameters
m=set(m,'tex_name',...
    {
    'alp','$\alpha$'
    'bet','$\beta$'
    'gam','$\gamma$'
    'rho','$\rho$'
    'psi','$\psi$'
	'del','$\delta$'
    'sig_a','$\sigma_a$'
    'sig_m','$\sigma_m$'
    });
%% plot priors, posteriors, priors and posteriors
plot_priors(m)
plot_posteriors(m,Results)
plot_priors_and_posteriors(m,Results)
%% plot priors, posteriors, priors and posteriors for a subset of parameters
close all
myparams={'alp','gam','psi'};
plot_priors(m,myparams)
plot_posteriors(m,Results,myparams)
plot_priors_and_posteriors(m,Results,myparams)
%% check curvature at the mode
mode_curvature(m)

%% check curvature at the mode for a subset of parameters
mode_curvature(m,myparams)

%% check curvature at the mode for a subset of parameters
profile off
profile on
mode_curvature(m,myparams)
profile off
profile viewer

author: jmaih


text:
function [system] = LoadStandingSystem(model, load_path)
    % Load the hybrid system model

    if nargin < 2
        load_path = [];
    end

    doubleStance = sys.domains.DoubleStance(model, load_path);
%     rightStance = sys.domains.RightStance(model, load_path);
%     leftStance = sys.domains.LeftStance(model, load_path);
%
%
%
%     leftImpact = RigidImpact('LeftImpact', rightStance, 'leftFootHeight'); % To leftStance
%     leftImpact.addImpactConstraint(struct2array(leftStance.HolonomicConstraints), load_path);
%
%     rightImpact = RigidImpact('RightImpact', leftStance, 'rightFootHeight'); % To rightStance
%     rightImpact.addImpactConstraint(struct2array(rightStance.HolonomicConstraints), load_path);
%
    io_control  = IOFeedback('IO');

    system = HybridSystem('atlas');
    system = addVertex(system, 'DoubleStance', 'Domain', doubleStance, ...
        'Control', io_control);
%     system = addVertex(system, 'RightStance', 'Domain', rightStance, ...
%         'Control', io_control);
%     system = addVertex(system, 'LeftStance', 'Domain', leftStance, ...
%         'Control', io_control);
%
%     srcs = {'RightStance'
%         'LeftStance'};
%
%     tars = {'LeftStance'
%         'RightStance'};
%
%     system = addEdge(system, srcs, tars);
%     system = setEdgeProperties(system, srcs, tars, ...
%         'Guard', {leftImpact, rightImpact});

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function annolist = compute_det_bbox(expidx,annolist,nms_thresh,det_thresh,bbox_offset)

fprintf('compute_det_bbox()\n');

p = exp_params(expidx);
predDir = p.multicutDir;

nmissing = 0;
bVis = false;
imgidxs_missing = [];
pidx = 13; % neck
headSize = p.refHeight/8;

n = 0;

for imgidx = 1:length(annolist)

    fprintf('.');
    keypointsAll(imgidx).imgname = annolist(imgidx).image.name;
    fname = [predDir '/imgidx_' padZeros(num2str(imgidx),4) '_cidx_1_14'];
    img = imread(annolist(imgidx).image.name);
    try
        load(fname,'unLab','unPos','unProb');
    catch
        annolist(imgidx).annorect(1).bbox = [];
        annolist(imgidx).annorect(2:end) = [];
        imgidxs_missing = [imgidxs_missing; imgidx];
        nmissing = nmissing + 1;
        continue;
    end

    bbox = [unPos(:,1) - 1.5*headSize unPos(:,2) - headSize,...
                 unPos(:,1) + 1.5*headSize unPos(:,2) + 7*headSize];

    bbox(:,1) = bbox(:,1) - bbox_offset;
    bbox(:,2) = bbox(:,2) - bbox_offset;
    bbox(:,3) = bbox(:,3) + bbox_offset;
    bbox(:,4) = bbox(:,4) + bbox_offset;

    bbox = [bbox unProb(:,pidx)];
    bbox(bbox(:,5) < det_thresh,:) = [];

    if (isempty(bbox))
        annolist(imgidx).annorect(1).bbox = [];
        annolist(imgidx).annorect(2:end) = [];
        continue;
    end

    I = nms(bbox,nms_thresh);
    bbox = bbox(I,:);
    bbox(:,1) = max(bbox(:,1),1);
    bbox(:,2) = max(bbox(:,2),1);
    bbox(3) = min(bbox(3),size(img,2));
    bbox(4) = min(bbox(4),size(img,1));

    n = n + size(bbox,1);
    for ridx = 1:size(bbox)
        annolist(imgidx).annorect(ridx).bbox = bbox(ridx,1:4);
    end

    if (size(bbox,1) < length(annolist(imgidx).annorect))
        annolist(imgidx).annorect(size(bbox,1)+1:end) = [];
    end
    if (bVis)
        colors = {'r','g','b','c','m','y'};
        markers = {'+','o','s'}; %,'x','.','-','s','d'
        img = imread(annolist(imgidx).image.name);
        figure(100); clf; imagesc(img); axis equal; hold on;
        for ridx = 1:size(bbox,1)
            bb = annolist(imgidx).annorect(ridx).bbox;
            rectangle('Pos',[bb(:,1) bb(:,2) bb(:,3)-bb(:,1) bb(:,4)-bb(:,2)],'edgeColor',colors{mod(ridx,6)+1},'LineWidth',5);
            axis off;
        end
    end

    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(annolist));
    end
end
fprintf(' done\n');
fprintf('# bbox: %d\n',n);
fprintf('nmissing: %d\n',nmissing);

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% Whole cell simulation class.
% - Runs simulations
% - Downsamples simulations
% - Stores and loads simulation to/from local disk or a database
% - Simulations are stored both on local disk and in the database with the
%   meta data, parameters, and time courses list below. Several dependent
%   time courses are also stored in the database.
% - Opens simulations in dashboard
%
% Options
% - lengthSec
% - stepSizeSec
% - verbosity
% - seed
% - geneticKnockouts
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 9/7/2010
classdef Simulation < handle
    properties (Constant)
        optionNames = {
            'lengthSec'
            'stepSizeSec'
            'verbosity'
            'seed'
            'macromoleculeStateInitialization'
            'geneticKnockouts'
            'stimulus'
            'media'
            };
    end

    properties (SetAccess = private)
        compartment             %compartments
        gene                    %indices of genes

        stateMetadata           %all state IDs, classes, names, etc.
        states                  %cell array of states in alphabetical order
        state_time              %reference to time state
        state_metabolite        %reference to metabolite state
        state_stimulus          %reference to stimulus state

        processMetadata         %all process IDs, classes, names, etc.
        processes               %cell array of processes in alphabetical order
        processInitOrderIndexs  %indices of processes in initialization order
        processEvalOrderIndexs  %indices of processes in state evolution order
        processesInInitOrder    %processes in initialization order
        processesInEvalOrder    %processes in state evolution order
        processIndex_translation        %index of translation process
        processIndex_tRNAAminoacylation %index of tRNA aminoacylation process
    end

    properties (Access = private)
        randStream              %random number stream
    end

    %options. Use applyOptions(name1, value1, ...) to change one or more
    %options before a run.
    properties (SetAccess = private)
        lengthSec   = 50000; %length of simulation (s); ~1.5x 9 h doubling time
        stepSizeSec = 1;     %time scale of simulation (s)
        verbosity   = 1;     %0 = no output, 5 = maximum output
        seed        = [];    %set to any number for reproducible random streams
        macromoleculeStateInitialization = 'multinomial';  %Toggles how inital state is calculated; see initializeState

        geneticKnockouts = {}; %Whole Cell Model IDs of genes to be knocked out
        stimulus         = []; %stimulus set values
        media            = []; %media set values
    end

    methods
        function this = Simulation(kbStates, kbProcesses)
            %construct components
            this.constructRandStream();

            %construct constants
            this.compartment = edu.stanford.covert.cell.sim.constant.Compartment();
            this.gene = edu.stanford.covert.cell.sim.constant.Gene();

            %construct states
            this.states = {};
            if exist('kbStates', 'var')
                sMetadata.wholeCellModelIDs = {kbStates.wholeCellModelID}';
                sMetadata.names             = {kbStates.name}';
                sMetadata.classes           = {kbStates.class}';

                this.stateMetadata = sMetadata;
                this.stateMetadata.indexs = ...
                    this.buildNameIndex(sMetadata.wholeCellModelIDs);
                this.constructStates(sMetadata.wholeCellModelIDs);
            end

            %construct processes
            this.processes = {};
            if exist('kbProcesses', 'var')
                mMetadata.wholeCellModelIDs = {kbProcesses.wholeCellModelID}';
                mMetadata.names             = {kbProcesses.name}';
                mMetadata.classes           = {kbProcesses.class}';
                mMetadata.initOrder         = [kbProcesses.initializationOrder]';
                mMetadata.evalOrder         = [kbProcesses.evaluationOrder]';

                this.processMetadata = mMetadata;
                this.processMetadata.indexs = ...
                    this.buildNameIndex(mMetadata.wholeCellModelIDs);

                this.constructProcesses(mMetadata.wholeCellModelIDs);
            end

            %link processes-states
            for i = 1:numel(this.states)
                this.states{i}.storeObjectReferences(this);
            end

            %link processes-states
            for i = 1:numel(this.processes)
                this.processes{i}.storeObjectReferences(this);
            end
        end
    end

    methods
        function state = state(this, wholeCellModelID)
            i = this.stateIndex(wholeCellModelID);
            if i
                state = this.states{i};
            else
                state = [];
            end
        end

        function process = process(this, wholeCellModelID)
            i = this.processIndex(wholeCellModelID);
            if i
                process = this.processes{i};
            else
                process = [];
            end
        end

        function i = stateIndex(this, wholeCellModelID)
            [~, i] = ismember(wholeCellModelID, this.stateWholeCellModelIDs);
            if i == 0
                [~, i] = ismember(['State_' wholeCellModelID], this.stateWholeCellModelIDs);
            end
        end

        function i = processIndex(this, wholeCellModelID)
            [~,i] = ismember(wholeCellModelID, this.processWholeCellModelIDs);
            if i == 0
                [~, i] = ismember(['Process_' wholeCellModelID], this.processWholeCellModelIDs);
            end
        end

        function result = stateWholeCellModelIDs(this)
            result = cellfun(@(s) s.wholeCellModelID, this.states, 'UniformOutput', false);
        end

        function result = processWholeCellModelIDs(this)
            result = cellfun(@(m) m.wholeCellModelID, this.processes, 'UniformOutput', false);
        end

        function this = applyAllParameters(this, varargin)
            if nargin >= 2 && isstruct(varargin{1})
                values = varargin{1};
            else
                values = cell2struct(varargin(2:2:end), varargin(1:2:end-1), 2);
            end

            this.applyOptions(values);
            this.applyParameters(values);
            this.applyFittedConstants(values);
            this.applyFixedConstants(values);
        end

        function this = applyOptions(this, varargin)
            if nargin >= 2 && isstruct(varargin{1})
                options = varargin{1};
            else
                options = cell2struct(varargin(2:2:end), varargin(1:2:end-1), 2);
            end
            this.setFromStruct(options, 'options', this.optionNames);
            if mod(this.lengthSec, this.stepSizeSec)
                throw(MException('Simulation:error',...
                    sprintf('Simulation length (%d) must be a multiple of step size (%d)',...
                    this.lengthSec, this.stepSizeSec)));
            end
        end

        function this = applyParameters(this, parameters)
            this.setFromStruct(parameters, 'parameters', {});
        end

        %apply perturbations
        %- genetic: adjust constants so that effectively the genes'
        %  products won't be expressed. This is done by setting the gene
        %  products' decay rates to be very high
        function this = applyPerturbations(this)
            %modify constants
            this.applyPerturbationsToConstants();

            %modify initial conditions
            this.applyPerturbationsToState()
        end

        function applyPerturbationsToConstants(this)
            %% genetic knockouts
            %references
            g = this.gene;
            r = this.state('Rna');
            m = this.state('ProteinMonomer');
            c = this.state('ProteinComplex');

            %indices
            [~, geneIdxs] = ismember(this.geneticKnockouts, this.gene.wholeCellModelIDs);
            if ~all(geneIdxs)
                throw(MException('Simulation:error', 'undefined genes'));
            end
            matureRNAIdxs = find(any(r.matureRNAGeneComposition(geneIdxs(~ismember(geneIdxs, g.mRNAIndexs)), :), 1));
            monIdxs = find(ismember(this.gene.mRNAIndexs, geneIdxs));
            cpxIdxs = find(any(any(c.proteinComplexComposition(geneIdxs, :, :), 3), 1));

            %set constants to that effectively the knocked out genes aren't
            %expressed. More specifically we allow genes to be expressed,
            %but then decay their gene products quickly.
            decayRate = 1e6;

            r.decayRates(r.processedIndexs(matureRNAIdxs), :)     = decayRate;
            r.decayRates(r.matureIndexs(matureRNAIdxs), :)        = decayRate;
            r.decayRates(r.boundIndexs(matureRNAIdxs), :)         = decayRate;
            r.decayRates(r.misfoldedIndexs(matureRNAIdxs), :)     = decayRate;
            r.decayRates(r.damagedIndexs(matureRNAIdxs), :)       = decayRate;
            r.decayRates(r.aminoacylatedIndexs(matureRNAIdxs), :) = decayRate;

            m.decayRates(m.nascentIndexs(monIdxs), :)             = decayRate;
            m.decayRates(m.signalSequenceIndexs(monIdxs), :)      = decayRate;
            m.decayRates(m.processedIIndexs(monIdxs), :)          = decayRate;
            m.decayRates(m.processedIIIndexs(monIdxs), :)         = decayRate;
            m.decayRates(m.foldedIndexs(monIdxs), :)              = decayRate;
            m.decayRates(m.matureIndexs(monIdxs), :)              = decayRate;
            m.decayRates(m.inactivatedIndexs(monIdxs), :)         = decayRate;
            m.decayRates(m.boundIndexs(monIdxs), :)               = decayRate;
            m.decayRates(m.misfoldedIndexs(monIdxs), :)           = decayRate;
            m.decayRates(m.damagedIndexs(monIdxs), :)             = decayRate;

            c.decayRates(c.nascentIndexs(cpxIdxs), :)             = decayRate;
            c.decayRates(c.matureIndexs(cpxIdxs), :)              = decayRate;
            c.decayRates(c.inactivatedIndexs(cpxIdxs), :)         = decayRate;
            c.decayRates(c.boundIndexs(cpxIdxs), :)               = decayRate;
            c.decayRates(c.misfoldedIndexs(cpxIdxs), :)           = decayRate;
            c.decayRates(c.damagedIndexs(cpxIdxs), :)             = decayRate;
        end

        function applyPerturbationsToState(this)
            import edu.stanford.covert.cell.sim.constant.Condition;

            %% genetic knockouts
            %references
            g = this.gene;
            r = this.state('Rna');
            m = this.state('ProteinMonomer');
            c = this.state('ProteinComplex');

            %indices
            [~, geneIdxs] = ismember(this.geneticKnockouts, this.gene.wholeCellModelIDs);
            if ~all(geneIdxs)
                throw(MException('Simulation:error', 'undefined genes'));
            end
            matureRNAIdxs = find(any(r.matureRNAGeneComposition(geneIdxs(~ismember(geneIdxs, g.mRNAIndexs)), :), 1));
            monIdxs = find(ismember(this.gene.mRNAIndexs, geneIdxs));
            cpxIdxs = find(any(any(c.proteinComplexComposition(geneIdxs, :, :), 3), 1));

            %modify initial conditions:
            %- set RNAs, protein counts to zero
            %- synchronize molecule counts with other states
            %  - Chromosome: bound proteins
            %  - RNA polmyerase, transcript: No active transcription
            %  - Ribosome, polypeptide: No active translation
            %  - FtsZ ring: No formed ring
            initRNACounts = r.counts;
            r.counts(r.processedIndexs(matureRNAIdxs), :)     = 0;
            r.counts(r.matureIndexs(matureRNAIdxs), :)        = 0;
            r.counts(r.boundIndexs(matureRNAIdxs), :)         = 0;
            r.counts(r.misfoldedIndexs(matureRNAIdxs), :)     = 0;
            r.counts(r.damagedIndexs(matureRNAIdxs), :)       = 0;
            r.counts(r.aminoacylatedIndexs(matureRNAIdxs), :) = 0;
            if any(any(r.updateExternalState(-(initRNACounts - r.counts), true)))
                throw(MException('Simulation:error', 'Unable to properly knockout gene'));
            end

            initMonCounts = m.counts;
            m.counts(m.nascentIndexs(monIdxs), :)        = 0;
            m.counts(m.signalSequenceIndexs(monIdxs), :) = 0;
            m.counts(m.processedIIndexs(monIdxs), :)     = 0;
            m.counts(m.processedIIIndexs(monIdxs), :)    = 0;
            m.counts(m.foldedIndexs(monIdxs), :)         = 0;
            m.counts(m.matureIndexs(monIdxs), :)         = 0;
            m.counts(m.inactivatedIndexs(monIdxs), :)    = 0;
            m.counts(m.boundIndexs(monIdxs), :)          = 0;
            m.counts(m.misfoldedIndexs(monIdxs), :)      = 0;
            m.counts(m.damagedIndexs(monIdxs), :)        = 0;
            if any(any(m.updateExternalState(-(initMonCounts - m.counts), true)))
                throw(MException('Simulation:error', 'Unable to properly knockout gene'));
            end

            initCpxCounts = c.counts;
            c.counts(c.nascentIndexs(cpxIdxs), :)     = 0;
            c.counts(c.matureIndexs(cpxIdxs), :)      = 0;
            c.counts(c.inactivatedIndexs(cpxIdxs), :) = 0;
            c.counts(c.boundIndexs(cpxIdxs), :)       = 0;
            c.counts(c.misfoldedIndexs(cpxIdxs), :)   = 0;
            c.counts(c.damagedIndexs(cpxIdxs), :)     = 0;
            if any(any(c.updateExternalState(-(initCpxCounts - c.counts), true)))
                throw(MException('Simulation:error', 'Unable to properly knockout gene'));
            end

            %% stimulus
            stim = this.state('Stimulus');
            newVals = this.stimulus;
            vals = stim.setValues;
            tfs = true(size(vals, 1), 1);
            for i = 1:size(newVals, 1)
                tfs = tfs & ~any(...
                    ismember(vals(:, Condition.objectCompartmentIndexs), newVals(:, Condition.objectCompartmentIndexs)) & ...
                    ((vals(:, Condition.initialTimeIndexs) >= newVals(:, Condition.initialTimeIndexs) & ...
                    vals(:, Condition.initialTimeIndexs) <= newVals(:, Condition.finalTimeIndexs)) | ...
                    (vals(:, Condition.finalTimeIndexs) >= newVals(:, Condition.initialTimeIndexs) & ...
                    vals(:, Condition.finalTimeIndexs) <= newVals(:, Condition.finalTimeIndexs)) | ...
                    (vals(:, Condition.finalTimeIndexs) <= newVals(:, Condition.initialTimeIndexs) & ...
                    vals(:, Condition.finalTimeIndexs) >= newVals(:, Condition.finalTimeIndexs))));
            end
            stim.setValues = [vals(tfs, :); newVals];

            %% media
            m = this.state('Metabolite');
            newVals = this.media;
            vals = m.setCounts;
            tfs = true(size(vals, 1), 1);
            for i = 1:size(newVals, 1)
                tfs = tfs & ~any(...
                    ismember(vals(:, Condition.objectCompartmentIndexs), newVals(:, Condition.objectCompartmentIndexs)) & ...
                    ((vals(:, Condition.initialTimeIndexs) >= newVals(:, Condition.initialTimeIndexs) & ...
                    vals(:, Condition.initialTimeIndexs) <= newVals(:, Condition.finalTimeIndexs)) | ...
                    (vals(:, Condition.finalTimeIndexs) >= newVals(:, Condition.initialTimeIndexs) & ...
                    vals(:, Condition.finalTimeIndexs) <= newVals(:, Condition.finalTimeIndexs)) | ...
                    (vals(:, Condition.finalTimeIndexs) <= newVals(:, Condition.initialTimeIndexs) & ...
                    vals(:, Condition.finalTimeIndexs) >= newVals(:, Condition.finalTimeIndexs))));
            end
            m.setCounts = [vals(tfs, :); newVals];
        end

        function this = applyFittedConstants(this, value)
            this.setFromStruct(value, 'fittedConstants', {});
        end

        function this = applyFixedConstants(this, value)
            this.setFromStruct(value, 'fixedConstants', {});
        end

        function this = setMetabolicReactionKinetics(this, value)
            %get handle to metabolism submodel
            met = this.process('Metabolism');
            if isempty(met)
                if isempty(value)
                    return;
                else
                    throw(MException('Simulation:error', 'Cannot set kinetic rates'));
                end
            end

            %set kinetics
            rxnIds = fields(value);
            [rxnTfs, rxnIdxs] = ismember(rxnIds, met.reactionWholeCellModelIDs);
            if ~all(rxnTfs)
                throw(MException('Simulation:error', 'Cannot set kinetics of undefined reactions:\n- %s', ...
                    strjoin(sprintf('\n- '), rxnIdxs{~rxnTfs})));
            end

            for i = 1:numel(rxnIds)
                if isfield(value.(rxnIds{i}), 'rev')
                    met.enzymeBounds(rxnIdxs(i), 1) = value.(rxnIds{i}).rev;
                end
                if isfield(value.(rxnIds{i}), 'for')
                    met.enzymeBounds(rxnIdxs(i), 2) = value.(rxnIds{i}).for;
                end
            end

            %copy kinetics over to FBA setup
            met.fbaEnzymeBounds(met.fbaReactionIndexs_metabolicConversion, :) = ...
                met.enzymeBounds(met.reactionIndexs_fba, :);
        end

        function this = applyRandStreamStates(this, value)
            if size(value.simulation(:, 1), 2) > 1
                warning('WholeCell:warning', 'Applying first rand stream state');
            end

            this.randStream.state = value.simulation(:, 1);
            for i = 1:numel(this.states)
                o = this.states{i};
                o.randStream.state = value.states.(o.wholeCellModelID(7:end))(:, 1);
            end
            for i = 1:numel(this.processes)
                o = this.processes{i};
                o.randStream.state = value.processes.(o.wholeCellModelID(9:end))(:, 1);
            end
        end

        function this = applyTimeCourses(this, value)
            this.setFromStruct(value, 'timeCourses', {});
        end

        function constructRandStream(this)
            %seed private stream
            this.randStream = edu.stanford.covert.util.RandStream('mcg16807');

            %seed state private streams
            for i = 1:length(this.states)
                o = this.states{i};
                o.constructRandStream();
            end

            %seed process private streams
            for i = 1:length(this.processes)
                o = this.processes{i};
                o.constructRandStream();
            end
        end

        function this = seedRandStream(this)
            %pick a seed
            if isempty(this.seed)
                this.seed = round(mod(now, 1) * 1e7);
            end

            %seed default stream
            if verLessThan('matlab', '7.12')
                globalStream = RandStream.getDefaultStream(); %#ok<GETRS>
            else
                globalStream = RandStream.getGlobalStream();
            end
            globalStream.reset(1);

            %seed private stream
            this.randStream.reset(this.seed)

            %seed state private streams
            for i = 1:length(this.states)
                o = this.states{i};
                o.seed = this.seed;
                o.seedRandStream();
            end

            %seed process private streams
            for i = 1:length(this.processes)
                o = this.processes{i};
                o.seed = this.seed;
                o.seedRandStream();
            end
        end

        function value = getAllParameters(this)
            import edu.stanford.covert.util.StructUtil;

            value = struct;
            value = StructUtil.catstruct(value, this.getOptions());
            value = StructUtil.catstruct(value, this.getParameters());
            value = StructUtil.catstruct(value, this.getFittedConstants());
            value = StructUtil.catstruct(value, this.getFixedConstants());
        end

        function value = getOptions(this)
            value = this.getAsStruct('options', this.optionNames);
        end

        function value = getParameters(this)
            value = this.getAsStruct('parameters', {});
        end

        function value = getFittedConstants(this)
            value = this.getAsStruct('fittedConstants', {});
        end

        function value = getFixedConstants(this)
            value = this.getAsStruct('fixedConstants', {});
        end

        function value = getMetabolicReactionKinetics(this)
            met = this.process('Metabolism');
            if isempty(met)
                throw(MException('Simulation:error', 'No metabolism process. Cannot get kinetics.'))
            end

            value = struct();
            for i = 1:numel(met.reactionWholeCellModelIDs)
                value.(met.reactionWholeCellModelIDs{i}) = struct(...
                    'rev', met.enzymeBounds(i, 1), ...
                    'for', met.enzymeBounds(i, 2) ...
                    );
            end
        end

        function value = getTimeCourses(this)
            value = this.getAsStruct('timeCourses', {}, true, false);
        end

        function value = getRandStreamStates(this)
            value = struct('simulation', [], 'states', struct(), 'processes', struct());
            value.simulation = this.randStream.state;
            for i = 1:numel(this.states)
                o = this.states{i};
                value.states.(o.wholeCellModelID(7:end)) = o.randStream.state;
            end
            for i = 1:numel(this.processes)
                o = this.processes{i};
                value.processes.(o.wholeCellModelID(9:end)) = o.randStream.state;
            end
        end

        %number of simulation steps
        function value = getNumSteps(this)
            value = this.lengthSec / this.stepSizeSec;
        end

        % Allocate arrays to store values of physical properties over course of
        % simulation. Each of these arrays has one column.
        function this = allocateMemoryForState(this, numTimePoints)
            %states
            for i = 1:numel(this.states)
                this.states{i}.allocateMemory(numTimePoints);
            end

            %processes
            for i = 1:numel(this.processes)
                this.processes{i}.allocateMemoryForState(numTimePoints);
            end
        end
    end

    %loading, allocation, initialization
    methods
        this = initializeConstants(this, knowledgeBase)
        this = initializeState(this)
        this = run(this, varargin)
        [this, requirements, allocations, usages] = evolveState(this, outputDirectory)
        [daughters, state_daughts] = divideState(this)
    end

    %protected helper methods
    methods (Access = protected)
        function this = constructStates(this, stateWholeCellModelIDs)
            if ~exist('stateWholeCellModelIDs', 'var')
                stateWholeCellModelIDs = this.stateMetadata.wholeCellModelIDs;
            end

            this.states = cell(numel(stateWholeCellModelIDs), 1);
            [~, map1] = ismember(stateWholeCellModelIDs, this.stateMetadata.wholeCellModelIDs);
            [~, map2] = ismember(cellfun(@(x) ['State_' x], stateWholeCellModelIDs, 'UniformOutput', false), this.stateMetadata.wholeCellModelIDs);
            map = map1 + map2;

            %construct states
            for i = 1:numel(stateWholeCellModelIDs)
                constructor = str2func(...
                    ['edu.stanford.covert.cell.sim.state.' this.stateMetadata.classes{map(i)}]);
                this.states{i} = constructor(...
                    this.stateMetadata.wholeCellModelIDs{map(i)}, this.stateMetadata.names{map(i)});
            end

            %construct references to states
            this.state_time       = this.state('Time');
            this.state_metabolite = this.state('Metabolite');
            this.state_stimulus   = this.state('Stimulus');
        end

        function this = constructProcesses(this, processWholeCellModelIDs)
            if ~exist('processWholeCellModelIDs', 'var')
                processWholeCellModelIDs = this.processMetadata.wholeCellModelIDs;
            end

            [~, map1] = ismember(processWholeCellModelIDs, this.processMetadata.wholeCellModelIDs);
            [~, map2] = ismember(cellfun(@(x) ['Process_' x], processWholeCellModelIDs, 'UniformOutput', false), this.processMetadata.wholeCellModelIDs);
            map = map1 + map2;

            if ~isempty(map)
                this.processes{length(map),1} = [];
                for i = 1:length(map)
                    constructor = str2func(...
                        ['edu.stanford.covert.cell.sim.process.' this.processMetadata.classes{map(i)}]);
                    this.processes{i} = constructor(...
                        this.processMetadata.wholeCellModelIDs{map(i)}, this.processMetadata.names{map(i)});
                end
            else
                this.processes(:) = [];
            end

            [this.processesInInitOrder, this.processInitOrderIndexs] = this.createProcessOrdering(this.processMetadata.initOrder);
            [this.processesInEvalOrder, this.processEvalOrderIndexs] = this.createProcessOrdering(this.processMetadata.evalOrder);

            this.processIndex_translation = this.processIndex('Translation');
            this.processIndex_tRNAAminoacylation = this.processIndex('tRNAAminoacylation');
        end

        function [processes, indexs] = createProcessOrdering(this, completeOrdering)
            ids = this.processWholeCellModelIDs;
            [~,order] = sort(completeOrdering);
            indexs = zeros(length(this.processes),1,'int8');
            j = 1;
            for i = 1:length(order)
                if completeOrdering(order(i))
                    [~,k] = ismember(this.processMetadata.wholeCellModelIDs(order(i)), ids);
                    if k
                        indexs(j) = k;
                        j = j + 1;
                    end
                end
            end
            indexs = indexs(1:j-1);

            processes = this.processes(indexs);
        end

        function result = buildNameIndex(~, names)
            result = struct;
            for i = 1:length(names)
                result.(names{i}) = i;
            end
        end

        function value = getAsStruct(this, field, fields, statesFlag, processesFlag)
            if nargin < 4
                statesFlag = true;
            end
            if nargin < 5
                processesFlag = true;
            end

            value = struct();
            if statesFlag
                value.states = struct;
            end
            if processesFlag
                value.processes = struct;
            end

            propertyNames = fieldnames(this);
            for i = 1:length(fields)
                if ismember(fields{i}, propertyNames)
                    value.(fields{i}) = this.(fields{i});
                end
            end

            if statesFlag
                for i = 1:length(this.states)
                    id = strrep(this.states{i}.wholeCellModelID, 'State_', '');
                    value.states.(id) = struct;
                    mValue = this.states{i}.(['get' upper(field(1)) field(2:end)]);
                    if ~isempty(fieldnames(mValue))
                        value.states.(id) = mValue;
                    end
                end
            end

            if processesFlag
                for i = 1:length(this.processes)
                    id = strrep(this.processes{i}.wholeCellModelID, 'Process_', '');
                    value.processes.(id) = struct;
                    mValue = this.processes{i}.(['get' upper(field(1)) field(2:end)]);
                    if ~isempty(fieldnames(mValue))
                        value.processes.(id) = mValue;
                    end
                end
            end
        end

        function this = setFromStruct(this, value, field, fields)
            propertyNames = fieldnames(this);
            for i = 1:length(fields)
                if isfield(value, fields{i}) && ismember(fields{i}, propertyNames)
                    try %#ok<TRYNC>
                        this.(fields{i}) = value.(fields{i});
                    end
                end
            end

            if isfield(value, 'states')
                for i = 1:length(this.states)
                    state = this.states{i};
                    id = strrep(state.wholeCellModelID, 'State_', '');
                    if isfield(value.states, id)
                        state.(['set' upper(field(1)) field(2:end)])(value.states.(id));
                    end
                end
            end

            if isfield(value, 'processes')
                for i = 1:length(this.processes)
                    process = this.processes{i};
                    id =  strrep(process.wholeCellModelID, 'Process_', '');
                    if isfield(value.processes, id)
                        process.(['set' upper(field(1)) field(2:end)])(value.processes.(id));
                    end
                end
            end
        end
    end

    %helper methods for testing
    methods
        function value = getForTest(this, propName)
            value = this.(propName);
        end

        function this = setForTest(this, propName, value)
            this.(propName) = value;
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function varargout = asec(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('callback');
        operator.bounds = @bounds;
        operator.derivative = @(x)(1./(x.*(x.^2-1).^0.5));

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

function [L,U] = bounds(xL,xU)
if xU <= -1 || xL >= 1
    L = asec(xL);
    U = asec(xU);
elseif xL < 0 & xU > 0
    L = 0;
    U = pi;
elseif xU < 0 || xL > 0
    L = real(asec(xL));
    U = real(asec(xU));
else
    L = 0;
    U = pi;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function X = uminus(X)
% UMINUS (overloaded)

X.basis = -X.basis;


text:
function [res] = tt_dot(tt1,tt2)
%Scalar product of two TT-tensors in the TT1.0 format
%   [RES]=TT_DOT(TT1,TT2) Scalar product of two TT-tensors in the TT1.0
%   format. Please avoid its usage, use dot function from the
%   object-oriented version
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
d=size(tt1,1);
g0=(tt1{1})'*(tt2{1});
gd=(tt1{d})'*(tt2{d});
r1=size(gd,1); r2=size(gd,2);
gd=reshape(gd,1,r1*r2);
for i=2:d-1
    core1=tt1{i};
    core2=tt2{i};
    core1=ten_conv(core1,2,g0);
    ncur=size(core1,1); r21=size(core1,2); r31=size(core1,3);
    r22=size(core2,2); r32=size(core2,3);
    core1=reshape(core1,[ncur*r21,r31]);
    core2=reshape(core2,[ncur*r22,r32]);
    g0=core1'*core2;
end
r1=size(g0,1); r2=size(g0,2);
g0=reshape(g0,1,r1*r2);
res=g0*gd';
return

end

author: oseledets

repo: TT-Toolbox


text:
% try first with Recon3D
load Recon3DModel_301.mat
start = 'xol7aone[r]';
stop = 'cholate[c]';
s=' ';
osenseStr = 'max';

% add sink reaction for start
model = addReaction(model,['start_', start],'reactionFormula',[start,s,'-->']);
model = addReaction(model,['stop_', stop],'reactionFormula',['-->',s,stop]);
printRxnFormula(model,model.rxns(end-1:end))

%% maximise flux through this dummy stop reaction
model = changeObjective(model,['stop_', stop]);
%force flux through this dummy start reaction
model = changeRxnBounds(model,['start_', start],-1000,'b');
model = changeRxnBounds(model,['stop_', stop],1000,'b');
model = changeRxnBounds(model,'sink_cholate[c]',0,'b');
%%
model0=model;

%%
[vSparse, sparseRxnBool1, essentialRxnBool]  = sparseFBA(model, osenseStr);
model.rxns(sparseRxnBool1)
nnz(sparseRxnBool1)
nnz(vSparse)

%%
rxnPenalty = ones(length(model.rxns),1);
param.printLevel = 1;
param.regularizeOuter = 0;
param.theta=0.1;
[solution,sparseRxnBool2] = findSparsePathway(model,rxnPenalty,param);
model.rxns(sparseRxnBool2)
nnz(sparseRxnBool2)
nnz(solution.v)

%check that both approaches give a sparse solution
assert(nnz(sparseRxnBool1)==nnz(sparseRxnBool2))

author: opencobra

repo: cobratoolbox


text:
function ctranspose(x)
%CTRANSPOSE is not defined for tensors.
%
%   See also TENSOR/PERMUTE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

error('Transpose on tensor is not defined');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function y = cat( dim, varargin )

%Disciplined convex/geometric programming information for CAT:
%   CAT imposes no convexity restrictions on its arguments.

if ~isnumeric( dim ) || any( size( dim ) ~= 1 ) || dim <= 0 || dim ~= floor( dim ),
    error( 'First argument must be a dimension.' );
end

%
% Quick exit
%

if nargin == 2,
    y = varargin{1};
    return
end

%
% Determine the final size and check consistency
%

sz    = [];
nz    = 1;
nzer  = 0;
nargs = 0;
isr   = true;
for k = 1 : nargin - 1,
    x = varargin{k};
    sx = size( x );
    if any( sx ),
        x = cvx( x );
        bx = x.basis_;
        nz = max( nz, size( bx, 1 ) );
        nzer = nzer + nnz( bx );
        sx( end + 1 : dim ) = 1;
        if isempty( sz ),
            sz = sx;
        elseif length( sx ) ~= length( sz ) || nnz( sx - sz ) > 1,
            error( 'All dimensions but the one being concatenated (%d) must be equal.', dim );
        else
            sz( dim ) = sz( dim ) + sx( dim );
        end
        if ~isreal( bx ),
            isr = false;
        end
        if all( sx ),
            nargs = nargs + 1;
            varargin{nargs} = x;
        end
    end
end

%
% Simple cases
%

if nargs == 0,

    if isempty( sz ),
        sz( dim ) = nargin;
    end
    y = cvx( sz, [] );
    return

elseif nargs == 1,

    y = varargin{1};
    return

end

%
% Harder cases
%

msiz = sz( dim );
lsiz = prod( sz( 1 : dim - 1 ) );
rsiz = prod( sz( dim + 1 : end ) );
psz  = lsiz * msiz * rsiz;
issp = cvx_use_sparse( [ nz, psz ], nzer, isr );
for k = 1 : nargs,
    x = varargin{k}.basis_;
    if issp ~= issparse(x),
        if issp,
            x = sparse(x);
        else
            x = full(x);
        end
    end
    x( end + 1 : nz, end ) = 0;
    if rsiz > 1,
        x = reshape( x, numel(x) / rsiz, rsiz );
    end
    varargin{k} = x;
end
if rsiz > 1,
    yb = builtin( 'cat', 1, varargin{1:nargs} );
else
    yb = builtin( 'cat', 2, varargin{1:nargs} );
end
yb = reshape( yb, nz, psz );

%
% Create object
%

y = cvx( sz, yb );

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

repo: REST


text:
% SAP2MATLAB parses Siemens ASCII protocol data and generates a
% corresponding MATLAB data structure.
%
% This function is used for de-serialising the header information
% from a FieldTrip buffer containing fMRI data from a Siemens scanner.

% Copyright (C) 2010, Stefan Klanke,
% 	Modified by Tim van Mourik, 2014
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

warning('Trying to compile MEX file')
oldDir = pwd;
try
    %This is where the header files are located
    options = '-I../include';
    %These two C-files are copiled
    fileNames = '../src/sap2matlab.c ../src/siemensap.c ';
    %the compilation does not depend on external libraries
    libraries = [];
    eval(['mex ' fileNames, libraries, options]);
catch
    rethrow(lasterror)
end

%% Test
%  This file contains an example string that will be parsed by sap2matlab
load('mrprotString.mat', 'apstr');
S = sap2matlab(apstr);

%%
file = '../Projects/TestReadMrprot/mrprot_triotim.txt';
f = fopen(file, 'r');
textData = fread(f);
fclose(f);

S = sap2matlab(sprintf('%s', textData));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function MatingPool = MatingSelection(Population,RefPoint,Range)
% The mating selection of AR-MOEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Calculate the degree of violation of each solution
    CV = sum(max(0,Population.cons),2);

    %% Calculate the fitness of each feasible solution based on IGD-NS
    if sum(CV==0) > 1
        % Calculate the distance between each solution and point
        N = sum(CV==0);
        Distance    = CalDistance(Population(CV==0).objs-repmat(Range(1,:),N,1),RefPoint);
        Convergence = min(Distance,[],2);
        [dis,rank]  = sort(Distance,1);
        % Calculate the fitness of noncontributing solutions
        Noncontributing = true(1,N);
        Noncontributing(rank(1,:)) = false;
        METRIC   = sum(dis(1,:)) + sum(Convergence(Noncontributing));
        fitness  = inf(1,N);
        fitness(Noncontributing) = METRIC - Convergence(Noncontributing);
        % Calculate the fitness of contributing solutions
        for p = find(~Noncontributing)
            temp = rank(1,:) == p;
            noncontributing = false(1,N);
            noncontributing(rank(2,temp)) = true;
            noncontributing = noncontributing & Noncontributing;
            fitness(p) = METRIC - sum(dis(1,temp)) + sum(dis(2,temp)) - sum(Convergence(noncontributing));
        end
    else
        fitness = zeros(1,sum(CV==0));
    end

    %% Combine the fitness of feasible solutions with the fitness of infeasible solutions
    Fitness = -inf(1,length(Population));
    Fitness(CV==0) = fitness;

    %% Binary tournament selection
    MatingPool = TournamentSelection(2,length(Population),CV,-Fitness);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% PLOT   Linear plot.
%    PLOT(X,Y) plots vector Y versus vector X. If X or Y is a matrix,
%    then the vector is plotted versus the rows or columns of the matrix,
%    whichever line up.  If X is a scalar and Y is a vector, disconnected
%    line objects are created and plotted as discrete points vertically at
%    X.
%
%    PLOT(Y) plots the columns of Y versus their index.
%    If Y is complex, PLOT(Y) is equivalent to PLOT(real(Y),imag(Y)).
%    In all other uses of PLOT, the imaginary part is ignored.
%
%    Various line types, plot symbols and colors may be obtained with
%    PLOT(X,Y,S) where S is a character string made from one element
%    from any or all the following 3 columns:
%
%           b     blue          .     point              -     solid
%           g     green         o     circle             :     dotted
%           r     red           x     x-mark             -.    dashdot
%           c     cyan          +     plus               --    dashed
%           m     magenta       *     star             (none)  no line
%           y     yellow        s     square
%           k     black         d     diamond
%           w     white         v     triangle (down)
%                               ^     triangle (up)
%                               <     triangle (left)
%                               >     triangle (right)
%                               p     pentagram
%                               h     hexagram
%
%    For example, PLOT(X,Y,'c+:') plots a cyan dotted line with a plus
%    at each data point; PLOT(X,Y,'bd') plots blue diamond at each data
%    point but does not draw any line.
%
%    PLOT(X1,Y1,S1,X2,Y2,S2,X3,Y3,S3,...) combines the plots defined by
%    the (X,Y,S) triples, where the X's and Y's are vectors or matrices
%    and the S's are strings.
%
%    For example, PLOT(X,Y,'y-',X,Y,'go') plots the data twice, with a
%    solid yellow line interpolating green circles at the data points.
%
%    The PLOT command, if no color is specified, makes automatic use of
%    the colors specified by the axes ColorOrder property.  By default,
%    PLOT cycles through the colors in the ColorOrder property.  For
%    monochrome systems, PLOT cycles over the axes LineStyleOrder property.
%
%    Note that RGB colors in the ColorOrder property may differ from
%    similarly-named colors in the (X,Y,S) triples.  For example, the
%    second axes ColorOrder property is medium green with RGB [0 .5 0],
%    while PLOT(X,Y,'g') plots a green line with RGB [0 1 0].
%
%    If you do not specify a marker type, PLOT uses no marker.
%    If you do not specify a line style, PLOT uses a solid line.
%
%    PLOT(AX,...) plots into the axes with handle AX.
%
%    PLOT returns a column vector of handles to lineseries objects, one
%    handle per plotted line.
%
%    The X,Y pairs, or X,Y,S triples, can be followed by
%    parameter/value pairs to specify additional properties
%    of the lines. For example, PLOT(X,Y,'LineWidth',2,'Color',[.6 0 0])
%    will create a plot with a dark red line width of 2 points.
%
%    Example
%       x = -pi:pi/10:pi;
%       y = tan(sin(x)) - sin(tan(x));
%       plot(x,y,'--rs','LineWidth',2,...
%                       'MarkerEdgeColor','k',...
%                       'MarkerFaceColor','g',...
%                       'MarkerSize',10)
%
%    See also PLOTTOOLS, SEMILOGX, SEMILOGY, LOGLOG, PLOTYY, PLOT3, GRID,
%    TITLE, XLABEL, YLABEL, AXIS, AXES, HOLD, LEGEND, SUBPLOT, SCATTER.
%
%    Reference page in Doc Center
%       doc plot
%
%    Other functions named plot
%
%       alphaShape/plot      empiricalblm/plot    semiconjugateblm/plot
%       blm/plot             fints/plot           tabular/plot
%       conjugateblm/plot    graph/plot           tall/plot
%       customblm/plot       LinearModel/plot     timeseries/plot
%       diffuseblm/plot      polyshape/plot       ts/plot
%       digraph/plot
%


text:
function point_num = sparse_grid_f2s_size ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_F2S_SIZE sizes a sparse grid using Fejer Type 2 Slow rules.
%
%  Discussion:
%
%    The grid is defined as the sum of the product rules whose LEVEL
%    satisfies:
%
%      0 <= LEVEL <= LEVEL_MAX.
%
%    This calculation is much faster than a previous method.  It simply
%    computes the number of new points that are added at each level in the
%    1D rule, and then counts the new points at a given DIM_NUM dimensional
%    level vector as the product of the new points added in each dimension.
%
%    This approach will work for nested families, and may be extensible
%    to other families, and to mixed rules.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the total number of unique
%    points in the grids.
%

%
%  Special case.
%
  if ( level_max < 0 )
    point_num = 0;
    return
  end

  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  Construct the vector that counts the new points in the 1D rule.
%
  new_1d = zeros ( level_max+1, 1 );

  new_1d(0+1) = 1;

  p = 1;
  o = 1;

  for l = 1 : level_max
    p = 2 * l + 1;
    if ( o < p )
      new_1d(l+1) = o + 1;
      o = 2 * o + 1;
    else
      new_1d(l+1) = 0;
    end
  end
%
%  Count the number of points by counting the number of new points
%  associated with each level vector.
%
  level_1d = zeros ( dim_num, 1 );

  point_num = 0;

  for level = 0 : level_max

    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );

      point_num = point_num + prod ( new_1d(level_1d(1:dim_num)+1) );

      if ( ~more )
        break
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [dg_rad, dg_display, dg_range] = default_doa_grid(n, unit, dim)
%DEFAULT_DOA_GRID Creates the search grid for DOA estimation.
%Inputs:
%   n - Number of grid points. For 2D DOAs, n can be 2-element vector,
%       specifying the number of grid points for azimuth and elevation
%       angles, respectively.
%   unit - Can be 'radian', 'degree', or 'sin'.
%   dim - DOA dimension.
%Outputs:
%   dg_rad - Grid of candidate DOAs converted to radians.
%   dg_display - Grid of candidate DOAs in the specified unit.
%   dg_range - Range of the DOA candidates. For the 1D case, this is an 1x2
%              vector. For the 2D case, this is a 2x2 matrix, where the
%              first row represents the range of the azimuth angle, and the
%              second row represents that of the elevation angle.
if isempty(dim)
    if design.dim > 1
        dim = 2;
    else
        dim = 1;
    end
else
    if dim ~= 1 && dim ~= 2
        error('Incorrect DOA dimension.');
    end
end
if dim == 2
    if isscalar(n)
        n = [n n];
    end
end
switch lower(unit)
    case 'radian'
        if dim == 1
            dg_range = [-pi/2 pi/2];
            dg_rad = -pi/2:pi/n:(pi/2-pi/n);
            dg_display = dg_rad;
        else
            dg_range = [0 2*pi;0 pi/2];
            dg_rad = grid2(0, 2*pi, 0, pi/2, n(1), n(2));
            dg_display = dg_rad;
        end
    case 'degree'
        if dim == 1
            dg_range = [-pi/2 pi/2];
            dg_rad = -pi/2:pi/n:(pi/2-pi/n);
            dg_display = rad2deg(dg_rad);
        else
            dg_range = [0 2*pi;0 pi/2];
            dg_rad = grid2(0, 2*pi, 0, pi/2, n(1), n(2));
            dg_display = rad2deg(dg_rad);
        end
    case 'sin'
        if dim == 1
            dg_range = [-1 1];
            dg_display = -1:2/n:(1-2/n);
            dg_rad = asin(dg_display);
        else
            dg_range = [-1 1;0 1];
            dg_display = grid2(-1, 1, 0, 1, n(1), n(2));
            dg_rad = asin(dg_display);
            dg_rad(1,:) = dg_rad(1,:) + pi/2; % azimuth -> [0, pi]
        end
    otherwise
        error('Unknown unit ''%s''.', unit);
end
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function [path P] = path_from_pred(pred,d,varargin)
% PATH_FROM_PRED Convert a predecessor array into a path to a vertex.
%
% path = path_from_pred(pred,d) returns the list of vertices on the path
% from a source vertex to a vertex d.  The predecessor array is a row
% vector such that
%   pred(i) = 0 if vertex i is a source vertex u,
%   pred(i) = 0 if vertex i has no predecessor associated with it, or
%   pred(i) = j where j preceeds vertex i on the shortest path from u to i.
% The vertex d is a number from 1 to n, where n is the length of the
% predecessor array.
%
% The returned path is a 1 by k+1 array where the path from the source to d
% has k edges path lists the order of visting the vertices.  In the case
% that the vertex is the source or unreachable from the source, the path
% just has the single starting vertex.  This reflects one potential problem
% with this function, the source is not uniquely encoded in the predecessor
% array and the source and unreachable vertices are treated equally.
%
% [path P] = path_from_pred(pred,d) also returns the adjacency matrix of the
% directed graph corresponding to the path.  Building this matrix takes
% additional time.
%
% Example:
%    load('graphs/bfs_example.mat');
%    [d dt pred] = bfs(A,1,struct('target', 3));
%    path = path_from_pred(pred,3); % sequence of vertices to vertex 3
%
% See also BFS, DFS, SHORTEST_PATHS

% David Gleich
% Copyright, Stanford University, 2007-2008

%% History
%  2007-04-17: Initial coding
%%

% TODO handle the case of a matrix of predecessors from all_shortest_paths

path = path_from_pred_mex(pred,d);

if nargout == 2
    n = length(pred);
    P = sparse(path(1:end-1), path(2:end), 1, n, n);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function writescores(pyra, model, info, fid)

% Write block scores to cascade score stats file.
%
% pyra     feature pyramid
% model    object model
% info     detection info from gdetect.m
% fid      file descriptor

% indexes into info from getdetections.cc
DET_USE = 1;    % current symbol is used
DET_IND = 2;    % rule index
DET_X   = 3;    % x coord (filter and deformation)
DET_Y   = 4;    % y coord (filter and deformation)
DET_L   = 5;    % level (filter)
DET_DS  = 6;    % # of 2x scalings relative to the start symbol location
DET_PX  = 7;    % x coord of "probe" (deformation)
DET_PY  = 8;    % y coord of "probe" (deformation)
DET_VAL = 9;    % score of current symbol
DET_SZ  = 10;   % <count number of constants above>

for i = 1:size(info,3)
  scores = zeros(model.numblocks, 1);

  for j = 1:model.numsymbols
    % skip unused symbols
    if info(DET_USE, j, i) == 0
      continue;
    end

    if model.symbols(j).type == 'T'
      scores = addfilterfeat(model, scores,           ...
                             info(DET_X, j, i),       ...
                             info(DET_Y, j, i),       ...
                             pyra.padx, pyra.pady,    ...
                             info(DET_DS, j, i),      ...
                             model.symbols(j).filter, ...
                             pyra.feat{info(DET_L, j, i)});
    else
      ruleind = info(DET_IND, j, i);
      if model.rules{j}(ruleind).type == 'D'
        bl = model.rules{j}(ruleind).def.blocklabel;
        dx = info(DET_PX, j, i) - info(DET_X, j, i);
        dy = info(DET_PY, j, i) - info(DET_Y, j, i);
        def = [-(dx^2); -dx; -(dy^2); -dy];
        scores(bl) = scores(bl) + model.rules{j}(ruleind).def.w*def;
      end
      bl = model.rules{j}(ruleind).offset.blocklabel;
      scores(bl) = scores(bl) + model.rules{j}(ruleind).offset.w;
    end
  end
  if abs(sum(scores) - info(DET_VAL, model.start, i)) > 1e-10
    fprintf('%f ~= %f\n', sum(scores), info(DET_VAL, model.start, i));
    error('wrong score');
  end
  dowrite(scores, fid);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stores the filter feature vector in the example ex
function scores = addfilterfeat(model, scores, x, y, padx, pady, ds, fi, feat)
% model object model
% x, y  location of filter in feat (with virtual padding)
% padx  number of cols of padding
% pady  number of rows of padding
% ds    number of 2x scalings (0 => root level, 1 => first part level, ...)
% fi    filter index
% feat  padded feature map

fsz = model.filters(fi).size;
% remove virtual padding
fy = y - pady*(2^ds-1);
fx = x - padx*(2^ds-1);
f = feat(fy:fy+fsz(1)-1, fx:fx+fsz(2)-1, :);

bl = model.filters(fi).blocklabel;
scores(bl) = scores(bl) + model.filters(fi).w(:)' * f(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write scores to fid
function dowrite(scores, fid)

% write total score and block scores
fwrite(fid, [sum(scores); scores], 'double');


text:
function [sourcemodel] = patchsvd(cfg, sourcemodel)

% This function does something like Limpiti et al.
% IEEE trans biomed eng 2006;53(9);1740-54

% Copyright (c) 2006, Jan-Mathijs Schoffelen & Robert Oostenveld, F.C. Donders Centre
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% set the defaults
if ~isfield(cfg, 'patchsvd'),    cfg.patchsvd    = 3;         end
if ~isfield(cfg, 'patchsvdnum'), cfg.patchsvdnum = 5;         end
if ~isfield(cfg, 'feedback'),    cfg.feedback    = 'textbar'; end

if isnumeric(cfg.patchsvd),
  Ndipoles = size(sourcemodel.pos,1);
else
  Ndipoles = 1;
end
Ninside  = length(sourcemodel.inside);
Nchans   = size(sourcemodel.leadfield{sourcemodel.inside(1)}, 1);
lfall    = cell(1,Ndipoles);
coeff    = cell(1,Ndipoles);
nghbr    = cell(1,Ndipoles);

if isnumeric(cfg.patchsvd) && ~isfield(sourcemodel, 'patchindx'),
  fprintf('computing patches in 3D, not taking topology into account\n');
  progress('init', cfg.feedback, 'computing patchsvd');
  for dipindx=1:Ninside
    % renumber the loop-index variable to make it easier to print the progress bar
    i = sourcemodel.inside(dipindx);

    % compute the distance from this dipole to each other dipole
    dist = sqrt(sum((sourcemodel.pos-repmat(sourcemodel.pos(i,:), [Ndipoles 1])).^2, 2));

    % define the region of interest around this dipole
    sel  = find(dist<=cfg.patchsvd);
    sel  = intersect(sel, sourcemodel.inside);
    Nsel = length(sel);

    progress(dipindx/Ninside, 'computing patchsvd %d/%d, Nsel=%d\n', dipindx, Ninside, Nsel);
    % concatenate the leadfield of all dipoles that are inside the ROI into one matrix
    lfr     = cell2mat(sourcemodel.leadfield(sel(:)'));
    % svd of leadfields of dipoles inside the ROI
    [U,S,V] = svd(lfr);

    if cfg.patchsvdnum < 1,
      % Limpiti et al 2006 formula 12
      s          = diag(S).^2;
      s          = cumsum(s)./sum(s);
      n(dipindx) = find(s - cfg.patchsvdnum > 0, 1);
    else
      n(dipindx) = cfg.patchsvdnum;
    end
    lfall{i} = U(:,1:n(dipindx)); %*S(1:n(dipindx),1:n(dipindx)); %klopt dit?
    nghbr{i} = sel;
    coeff{i} = V(1:n(dipindx), :);
  end
  progress('close');
elseif isnumeric(cfg.patchsvd) && isfield(sourcemodel, 'patchindx'),
  fprintf('computing patches in 2D, taking surface topology into account\n');
  progress('init', cfg.feedback, 'computing patchsvd');
  for dipindx=1:Ninside
    % renumber the loop-index variable to make it easier to print the progress bar
    i = sourcemodel.inside(dipindx);

    % define the region of interest around this dipole
    sel  = nearest(sourcemodel.patchsize(i,:), cfg.patchsvd);
    sel  = sourcemodel.patchindx(i,1:sel);
    Nsel = length(sel);

    progress(dipindx/Ninside, 'computing patchsvd %d/%d, Nsel=%d\n', dipindx, Ninside, Nsel);
    % concatenate the leadfield of all dipoles that are inside the ROI into one matrix
    lfr     = cell2mat(sourcemodel.leadfield(sel(:)'));
    % svd of leadfields of dipoles inside the ROI
    [U,S,V] = svd(lfr);

    if cfg.patchsvdnum < 1,
      % Limpiti et al 2006 formula 12
      s          = diag(S).^2;
      s          = cumsum(s)./sum(s);
      n(dipindx) = find(s - cfg.patchsvdnum > 0, 1);
    else
      n(dipindx) = min(cfg.patchsvdnum, size(lfr,2));
    end
    lfall{i} = U(:,1:n(dipindx)); %*S(1:n(dipindx),1:n(dipindx)); %klopt dit?
    nghbr{i} = sel;
    coeff{i} = V(1:n(dipindx), :);
    sv{i}    = s;
  end
  progress('close');
elseif strcmp(cfg.patchsvd, 'all'),
  lfr     = cell2mat(sourcemodel.leadfield(sourcemodel.inside(:)'));
  [U,S,V] = svd(lfr);

  if cfg.patchsvdnum < 1,
    % Limpiti et al 2006 formula 12
    s = diag(S).^2;
    s = cumsum(s)./sum(s);
    n = find(s - cfg.patchsvdnum > 0, 1);
  else
    n = cfg.patchsvdnum;
  end
  lfall{1} = U(:,1:n); %*S(1:n(dipindx),1:n(dipindx)); %klopt dit?
  nghbr{1} = sourcemodel.inside;
  coeff{1} = V(1:n, :);

  %---change output
  sourcemodel.pos    = mean(sourcemodel.pos(sourcemodel.inside,:),1);
  sourcemodel.inside = 1;
  sourcemodel.dim    = [1 1 1];
  sourcemodel.xsourcemodel  = sourcemodel.pos(1);
  sourcemodel.ysourcemodel  = sourcemodel.pos(2);
  sourcemodel.zsourcemodel  = sourcemodel.pos(3);
else
  %do nothing
end

if ~all(n==n(1)),
  nmax = max(n);
  for dipindx = 1:Ninside
    i = sourcemodel.inside(dipindx);
    if n(dipindx)<nmax,
      lfall{i} = [lfall{i} zeros(Nchans, nmax-n(dipindx))];
    end
  end
end

% update the leadfields
sourcemodel.leadfield = lfall;
%sourcemodel.expcoeff  = coeff;
sourcemodel.neighbours= nghbr;

author: spm


text:
function [ x, y, w ] = rule04 ( )

%*****************************************************************************80
%
%% RULE04 returns the rule of degree 4.
%
%  Discussion:
%
%    Order 4 (6 pts)
%    1/6 data for 4-th order quadrature with 2 nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Output, real X(*), Y(*), the coordinates of the nodes.
%
%    Output, real W(*), the weights.
%
  x = [ ...
       0.00000000000000000000000000000000, ...
       0.00000000000000000000000000000000 ];
  y = [ ...
      -0.39011034927098671317419951524681, ...
       0.83747122925185483484415621413487 ];
  w = [ ...
       0.14699335257362380894667320570121, ...
       0.72352316251791601189863277931619E-01 ];

  return
end


text:
function conf = voc_config_inriaperson()
% Configuration used for training the INRIA person model
%
% To use this execute:
%  >> global VOC_CONFIG_OVERRIDE;
%  >> VOC_CONFIG_OVERRIDE = @voc_config_inriaperson;

% SEE: INRIA/README for more details

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf.training.C = 0.006;

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function s = allsubs(x)
%ALLSUBS Generate all possible subscripts for a sparse tensor X.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

%% Generate all possible indicies

% Preallocate (discover any memory issues here!)
s = zeros(prod(x.size),ndims(x));

% Generate appropriately sized ones vectors.
o = cell(ndims(x),1);
for n = 1:ndims(x)
    o{n} = ones(size(x,n),1);
end

% Generate each column of the subscripts in turn
for n = 1:ndims(x)
    i = o;
    i{n} = (1:size(x,n))';
    s(:,n) = khatrirao(i);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function d = det3(A)

% det3 - 3x3 determinant
%
%   d = det3(A);
%
%   A is a [3 3 n] matrix, d is a vector of size d where
%   d(i)=det(A(:,:,i)).
%
%   It computes explicitely the det by expanding along the 1st columns 2x2
%   determinant. Faster than recursive loops.
%
%   Note: works also for 2x2 det.
%   Note: for higher order det, use loop
%
%   Copyright (c) 2008 Gabriel Peyre

if size(A,1)==2 && size(A,2)==2
    %% 2x2 Det %%
    d = A(1,1,:).*A(2,2,:) - A(1,2,:).*A(2,1,:);
elseif size(A,1)==3 && size(A,2)==3
    %% 3x3 Det %%
    d = A(1,1,:).*( A(2,2,:).*A(3,3,:) - A(2,3,:).*A(3,2,:) ) - ...
        A(1,2,:).*( A(2,1,:).*A(3,3,:) - A(2,3,:).*A(3,1,:) ) + ...
        A(1,3,:).*( A(2,1,:).*A(3,2,:) - A(2,2,:).*A(3,1,:) );
else
    n = size(A,3);
    d = zeros(n,1);
    for i=1:n
        d(i) = det(A(:,:,i));
    end
end
d = d(:);

author: gpeyre

repo: matlab-toolboxes


text:
function [out] = evap_23(p1,p2,S,Smax,Ep,dt)
% evap_23 combines evap_5 (evaporation) and evap_6 (transpiration)

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

% Flux function
% ------------------
% Description:  Transpiration from vegetation at the potential rate if
%               storage is above field capacity and scaled by relative
%               storage if not (similar to evap_6), addition of
%               Evaporation from bare soil scaled by relative storage
%               (similar to evap_5)
% Constraints:  Ea <= Ep
%               Ea <= S/dt
% @(Inputs):    p1   - fraction vegetated area [-] (0...1)
%               p2   - field capacity coefficient[-]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               Ep   - potential evapotranspiration rate [mm/d]
%               dt   - time step size [d]

out = min([p1.*Ep+(1-p1).*S./Smax.*Ep, p1*Ep*S./(p2*Smax)+(1-p1).*S./Smax.*Ep,S/dt]);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


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

    methods
        function obj = Rz(varargin)
            obj@ETS2('Rz', varargin{:});
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
%SWIRLFLOW Axisymmetric swirl flow physics mode.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
classdef StudentT2C
%%STUDENTT2C A collection of functions for the bivariate Student-t copulae.
%           Implemented functions include: PDF, tau.
%
%REFERENCES:
%[1] H. Joe and D. Kurowicka, Dependence Modeling: Vine Copula Handbook.
%    World Scientific, 2011.
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    methods(Static)
        function jprob = PDF(u,R,nu)
            %PDF Generates a joint probability for column vector u or
            %    a row vector of joint probabilities for matrix u
            %    with each column as the input vector.
            %
            %INPUTS:
            % u: A 2-by-n matrix where each column represents a
            %    2-dimensional random vector and n is the number of
            %    random vectors. Entries must be in [0,1].
            % R: A 2-by-2 correlation matrix, or a scalar which is the
            %    correlation.
            % nu: Scalar degrees of freedom.
            %
            %OUTPUTS: jprob: A 1-by-n vector of joint probabilities.
            %
            %EXAMPLE 1: Make a contour plot of the density.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            %     for iidx = 1:length(Y)
            %         Z(idx,iidx) = StudentT2C.PDF([X(idx,iidx);Y(idx,iidx)],0.5,3);
            %     end
            % end
            % contourf(X,Y,Z)
            %
            %EXAMPLE 2: Generate an exact joint density plot.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            % for iidx = 1:length(Y)
            % Z(idx,iidx) = StudentT2C.PDF([X(idx,iidx);Y(idx,iidx)],0.95,3);
            % X(idx,iidx) = GumbelD.invCDF(X(idx,iidx),1,2);
            % Y(idx,iidx) = GaussianD.invCDF(Y(idx,iidx));
            % Z(idx,iidx) = Z(idx,iidx).*GumbelD.PDF(X(idx,iidx),1,2).*GaussianD.PDF(Y(idx,iidx));
            % end
            % end
            % surf(X,Y,Z)
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('u','var') || isempty(u)
                jprob = [];
                return
            end
            if ~exist('R','var') || isempty(R)
                rho = 0;
            elseif ~isscalar(R)
                rho = R(2,1);
            else
                rho = R;
            end

            a = StudentTD.invCDF(u(1,:),nu);
            b = StudentTD.invCDF(u(2,:),nu);
            c = (a.^2+b.^2-2*rho*a.*b)/(nu*(1-rho^2));

            jprob = (gamma(0.5*nu+1)/(gamma(nu/2)*nu*pi*sqrt(1-rho^2)))...
                    ./(StudentTD.PDF(a,0,1,nu).*StudentTD.PDF(b,0,1,nu).*(1+c).^(0.5*nu+1));

            %NaN results from 0*inf in the calculation.
            %These should be 0s.
            nanIdxs = isnan(jprob);
            jprob(nanIdxs) = 0;

        end

        function t = tau(R)
            %%TAU Compute Kendall's tau for the copula.
            %
            %INPUTS:
            % R: A 2-by-2 correlation matrix, or a scalar which is the
            %    correlation.
            %
            %OUTPUTS:
            % t: The scalar value of Kendall's tau.
            %
            %Note: This does not depend on the degrees of freedom nu.
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('R','var') || isempty(R)
                rho = 0;
            elseif ~isscalar(R)
                rho = R(2,1);
            else
                rho = R;
            end

            t = 2*arcsin(rho)/pi;
        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 msm_to_mm_test03 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST03 tests MSM_TO_MM_ARRAY_COMPLEX_SKEWSYMMETRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST03\n' );
  fprintf ( 1, '  Convert an MSM to MM array complex skew-symmetric format.\n' );

  output_filename = 'msm_to_mm_test03.mm';
  a = c8mat_indicator ( 4, 4 );
  for j = 1 : 4
    a(j,j) = 0.0;
    for i = j + 1 : 4
      a(i,j) = - a(j,i);
    end
  end
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'array', 'complex', 'skew-symmetric' );

  return
end


text:
function [xUpdate,MUpdate,DUpdate,innov,Pzz,W]=reducedStateUpdate(xPred,PPred,MPred,DPred,z,R,H)
%%REDSTATEFILTERUPDATE Perform the measurement update step in the reduced
%               state estimator. This is a filter that takes measurements
%               in Cartesian coordinates and assumes that the dynamic model
%               includes an unknown parameter that is somehow bounded, but
%               whose contribution is modeled with a mean and covariance
%               matrix. The filter separates contributions due to
%               measurement errors and dynamic model mismatch errors.
%
%INPUTS: xPred The xDimX1 predicted state estimate.
%        PPred The xDimXxDim predicted total state covariance estimate.
%              This is provided by the state prediction step.
%        MPred The xDimXxDim matrix contributing to the total predicted
%              state covariance matrix based solely on measurement errors.
%        DPred The xDimXzDim matrix of bias coefficients that are supposed
%              to relate target state errors to dynamic model parameter
%              uncertainty.
%            z The zDimX1 measurement vector.
%            R The zDimXzDim measurement covariance matrix.
%            H The optional zDimXxDim measurement matrix. The measurement
%              is modeled as z=H*x+noise. If this parameter is omitted or
%              an empty matrix is passed, then H will be taken as a
%              zDimXxDim identity matrix followed by columns of zeros
%              (Assuming that zDim<=xDim. Otherwise, H must be provided).
%
%OUTPUTS: xUpdate The xDimX1 updated state estimate.
%         MUpdate The xDimXxDim updated state covariance matrix
%                 contribution due to measurement errors.
%         DUpdate The xDimXzDim matrix of updated bias coefficients
%                 contributing to the state covariance matrix.
%      innov, Pzz The zDimX1 innovation and a zDimXzDim matrix Pzz that is
%                 akin to an innovation covariance matrix are returned in
%                 case one wishes to analyze the consistency of the
%                 estimator or use those values in gating or likelihood
%                 evaluation.
%               W The gain used in the update. This can be useful when
%                 gating and using the function calcMissedGateCov.
%
%The filter is taken from [1]. Other applications are described in [1].
%Note that the formulation of the dynamic model requires that the
%dimensionality of the measurements does not vary.
%
%In [1] and [2], no clear method of initializing this type of tracking
%filter is provided. A simple way to initialize the filter would be to use
%two Cartesian converted measurements to obtain a state estimate and
%covariance as one would do with a normal Kalman filter (one could, for
%example, use the KalmanFIRSmoother function) and then set MPrev to the
%covariance value obtained while setting DPrev to zero.
%
%REFERENCES:
%[1] P. Mookerjee and F. Reifler, "Reduced state estimator for systems with
%    parametric inputs," IEEE Transactions on Aerospace and Electronic
%    Systems, vol. 40, no. 2, pp. 446-461, Apr. 2004.
%[2] P. Mookerjee and F. Reifler, "Reduced state estimators for consistent
%    tracking of maneuvering targets," IEEE Transactions on Aerospace and
%    Electronic Systems, vol. 41, no. 2, pp. 608-619, Apr. 2005.
%
%July 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(xPred,1);

if(nargin<7||isempty(H))
	zDim=size(z,1);
    H=[eye(zDim,zDim),zeros(zDim,xDim-zDim)];
end

%Equation 20 in [1].
Pzz=H*PPred*H'+R;

%Ensure symmetry
Pzz=(Pzz+Pzz')/2;

W=PPred*H'/Pzz;
L=eye(xDim,xDim)-W*H;

%Equation 23 in [1].
MUpdate=L*MPred*L'+W*R*W';

%Ensure symmetry
MUpdate=(MUpdate+MUpdate')/2;

%Equation 24 in [1].
DUpdate=L*DPred;

innov=z-H*xPred;
%Equation 26 in [1].
xUpdate=xPred+W*innov;

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.


text:
%% Example
% This script is used to initialize the hand guiding functionality on KUKA
% robot remotely

% Mohammad SAFEEA, 9th of June 2017

% TO use this example:
% 1- Start the server on KUKA
% 2- Start the client on Matlab
% 3- Press the white button and hand guide the robot
% 4- To save the coordinates of the robot, long click on the green button
% 5- Repeat three times
% 6- The robot gives you an interval of five seconds, clear the area around
% the robot directly, the robot will reproduce the motion tought.

close all;clear;clc;
warning('off')
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_Touch_pneumatisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
  return;
end
pause(1);

iiwa.startHandGuiding()

p1 = iiwa.getJointsPos();

iiwa.startHandGuiding()

p2 = iiwa.getJointsPos();

iiwa.startHandGuiding()

p3 = iiwa.getJointsPos();

iiwa.startHandGuiding()

fprintf('CLear the area round the robot, the robot is going to move');

pause(5)

relVel=0.1;
iiwa.movePTPJointSpace(p1, relVel)
iiwa.movePTPJointSpace(p2, relVel)
iiwa.movePTPJointSpace(p3, relVel)

%% turn off the server
iiwa.net_turnOffServer();

repo: KST-Kuka-Sunrise-Toolbox


text:
%%**********************************************************************
%% NTscaling: Compute NT scaling matrix
%%
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'.
%% compute W satisfying W*Z*W = X.
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%     important to keep W symmertic.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************************

 function [W,G,sv,gamx,gamz,dd,ee,ff] = ...
          NTscaling(blk,X,Z,Zchol,invZchol);

    numblk = size(blk,1);
    W = cell(numblk,1); G = cell(numblk,1); sv = cell(numblk,1);
    gamx = cell(numblk,1); gamz = cell(numblk,1);
    dd = cell(numblk,1); ee = cell(numblk,1); ff = cell(numblk,1);
%%
    for p = 1:size(blk,1)
        pblk = blk(p,:);
        numblk = length(pblk{2});
        n = sum(pblk{2});
        if strcmp(pblk{1},'l')
           dd{p} = X{p}./Z{p}; %% do not add perturbation, it badly affects cre-a
        elseif strcmp(pblk{1},'q');
           gamx{p} = sqrt(qops(pblk,X{p},X{p},2));
           gamz{p} = sqrt(qops(pblk,Z{p},Z{p},2));
           w2 = gamz{p}./gamx{p};  w = sqrt(w2);
           dd{p} = qops(pblk,1./w2,ones(n,1),4);
           tt = qops(pblk,1./w,Z{p},3) - qops(pblk,w,X{p},4);
           gamtt = sqrt(qops(pblk,tt,tt,2));
           ff{p} = qops(pblk,1./gamtt,tt,3);
           ee{p} = qops(pblk,sqrt(2)./w,ff{p},4);
        elseif strcmp(pblk{1},'s')
           tmp = Prod2(pblk,Zchol{p},X{p},0);
           tmp = Prod2(pblk,tmp,Zchol{p}',1);
           [sv2,V] = blkeig(pblk,tmp);
           sv2 = max(1e-20,sv2);
           sv{p} = sqrt(sv2);
           tmp  = Prod2(pblk,invZchol{p},V);
           G{p} = Prod2(pblk,spdiags(sqrt(sv{p}),0,n,n),tmp');
           W{p} = Prod2(pblk,G{p}',G{p},1);
        end
    end
%%**********************************************************************


text:
function [obs,nobs]=scan_obs_ppk(rtk,obsr)

global gls
nobs0=size(obsr,1); nobs=0; opt=rtk.opt; SNRMASK=32; %#ok
obs=repmat(gls.obs_tmp,nobs0,1);

for i=1:nobs0
    if obsr(i).S(1)~=0&&obsr(i).S(1)/4.0<SNRMASK
        continue;
    end

    obs(nobs+1)=obsr(i);
    nobs=nobs+1;
end

if nobs==0,obs=NaN;return;end
if nobs<nobs0
    obs(nobs+1:end,:)=[];
end

return

repo: GINav