text:
%otyp==0->cont, otyp==1->disc
function [A N]=sys_ned_dcm(llh, vel_n, Cbn, acc, otyp, dt)
% position   (1-3)
% velocity   (4-6)
% attitude   (7-9)

%system disturbance coefs
N=zeros(9,6);
N(7:9,4:6)=-Cbn; %attitude
N(4:6,1:3)=Cbn; %velocity

%system matrix
A=zeros(9);
[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(llh);
tL=sL/cL;
Rn_h=Rn+llh(3);
Re_h=Re+llh(3);
acc_n=Cbn*acc;

A(1,3)=-vel_n(1)/(Rn_h)^2;
A(1,4)=1/(Rn_h);

A(2,1)=vel_n(2)*tL/(Re_h)/cL;
A(2,3)=-vel_n(2)/(Re_h)/(Re_h)/cL;
A(2,5)=1/(Re_h)/cL;

A(3,6)=-1;

A(4,1)=-(2*WIE_E*cL*vel_n(2)+vel_n(2)*vel_n(2)/(Re_h)/cL/cL);
A(4,3)=(-vel_n(1)*vel_n(3)/(Rn_h)/(Rn_h)+vel_n(2)*vel_n(2)*tL/(Re_h)/(Re_h));
A(4,4)=vel_n(3)/(Rn_h);
A(4,5)=-(2*WIE_E*sL+2*vel_n(2)*tL/(Re_h));
A(4,6)=vel_n(1)/(Rn_h);
A(4,8)=-acc_n(3);
A(4,9)=acc_n(2);

A(5,1)=(2*WIE_E*cL*vel_n(1)-2*WIE_E*sL*vel_n(3)+vel_n(2)*vel_n(1)/(Re_h)/cL/cL);
A(5,3)=-(vel_n(2)*vel_n(3)+vel_n(2)*vel_n(1)*tL)/(Re_h)/(Re_h);
A(5,4)=(2*WIE_E*sL+vel_n(2)*tL/(Re_h));
A(5,5)=(vel_n(3)+vel_n(1)*tL)/(Re_h);
A(5,6)=(2*WIE_E*cL+vel_n(2)/(Re_h));
A(5,7)=acc_n(3);
A(5,9)=-acc_n(1);

A(6,1)=2*WIE_E*sL*vel_n(2);
A(6,3)=(vel_n(1)*vel_n(1)/(Rn_h)/(Rn_h)+vel_n(2)*vel_n(2)/(Re_h)/(Re_h));
A(6,4)=-2*vel_n(1)/(Rn_h);
A(6,5)=-2*(WIE_E*cL+vel_n(2)/(Re_h));
A(6,7)=-acc_n(2);
A(6,8)=acc_n(1);

A(7,1)=-WIE_E*sL;
A(7,3)=-vel_n(2)/(Re_h)/(Re_h);
A(7,5)=1/(Re_h);
A(7,8)=-(WIE_E*sL+vel_n(2)*tL/(Re_h));
A(7,9)=vel_n(1)/(Rn_h);

A(8,3)=vel_n(1)/(Rn_h)/(Rn_h);
A(8,4)=-1/(Rn_h);
A(8,7)=(WIE_E*sL+vel_n(2)*tL/(Re_h));
A(8,9)=(WIE_E*cL+vel_n(2)/(Re_h));

A(9,1)=-(WIE_E*cL+vel_n(2)/(Re_h)/cL/cL);
A(9,3)=vel_n(2)*tL/(Re_h)/(Re_h);
A(9,5)=-tL/(Re_h);
A(9,7)=-vel_n(1)/(Rn_h);
A(9,8)=-(WIE_E*cL+vel_n(2)/(Re_h));

%discretize A;
if (otyp==0)
    A=A;
elseif(otyp==1)
    A=expm(A*dt);
end

repo: ekfmonoslam


text:
function a = sparse(a)
%SPARSE       Convert Taylor to sparse
%

% written  05/21/09     S.M. Rump
%

  a.t = sparse(a.t);

repo: OpenMAS


text:
clear Models, clear MM

% User-Specified Parameters
% ===============================================================

% Conditions and stimuli
% ---------------------------------------------------------------
GA.conditions = [1 2];
GA.freqConditions = [.5 .5];
GA.scanLength = 480;
GA.ISI = 2;
GA.TR = 2;
GA.doepochs = 1;        % build epochs instead of events

% Hox optimization parameters
% A hox gene is a master gene.  These numbers control the stimulus
% parameters, which can be optimized within the GA rather than
% pre-specified.
% Hox elements in this GA go in the stimlist as the first several elements.
% hox sequence currently codes, in number order:
% ISI TR cuelen cuerest stimlen stimrest resplen resprest (in s)
% TR does not work right now
% use zeros to use default rather than variable hox parameters
% ---------------------------------------------------------------
GA.numhox = 8;
GA.hoxrange = [0 0; 1 1; 1 10; 1 10; 1 10; 1 10; 1 10; 1 10];
    % rows are each hox gene's allowable range, cols index hox elements

% Genetic algorithm parameters
% ---------------------------------------------------------------
nmodels = 1;
GA.cbalColinPowerWeights = [0 1 1 1];	% 1 = cbal, 2 = eff, 3 = hrf shape, 4 = freq
GA.numGenerations = 100000;
GA.sizeGenerations = 20;
GA.maxTime = 20;						% max time to run in s, or Inf
GA.alph = 2.1;
GA.plotFlag = 0;

% Filtering, counterbalancing, and design tolerance
% ---------------------------------------------------------------
GA.lowerLimit = [];
GA.HPlength = [120];
GA.LPsmooth = [];
GA.maxOrder = 1;
GA.NumStimthresh = [];
GA.maxCbalDevthresh = []; 	% counterbalancing deviation
GA.maxFreqDevthresh = [];

% Contrast setup
% ---------------------------------------------------------------
GA.contrasts = [-1 0 0 1 0 0; 0 -1 0 0 1 0; 1 -1 0 1 -1 0];      % three elements per condition (3 epochs)!
GA.contrastweights = [1 1 1];	                      % or predictor weights, if no contrasts

% Autocorrelation and special options
% ---------------------------------------------------------------
AutocorrelationFileName = 'myscannerxc';
GA.restlength = [];
GA.restevery = [];
GA.trans2switch = 0;
GA.trans2block = 0;
GA.dofirst = 0;
GA.nonlinthreshold = [2];

eval(['load ' AutocorrelationFileName]);
GA.xc = myscannerxc;

% =============================================
% * vary by parameter - set this in this script
% * if running multiple models, nmodels > 1
% ---------------------------------------------

varyparam = 0;
fieldname = 'alph';		% or 'freqConditions(end), etc.
incrementby = 0.3;
incrementevery = 5;		% every n models
contrastsofinterest = 1;% contrasts or predictors, if no cons specified

% =============================================

if varyparam
	eval(['paramvalue = GA.' fieldname ';'])
	disp(' ');disp('* *********************************************************************')
	disp('*  You have selected to vary a parameter across models in ga_gui_script')
	disp(['*  GA.' fieldname ' starting at ' num2str(paramvalue)])
	disp('* *********************************************************************')
end

for nm = 1:nmodels

   eval(['diary model' num2str(nm) '.log'])
   disp(['Starting Model ' num2str(nm)])
    % freqConditions

   M = optimizeGA_epochs(GA);

   Models(:,nm) = M.stimlist;
   MM{nm} = M;
   save GAworkspace

   if varyparam
   	if isfield(M,'consebeta'),
		Fit(1,nm) = mean(M.consebeta(contrastsofinterest));
	else
		Fit(1,nm) = mean(M.sebeta(contrastsofinterest));
	end
   	eval(['Fit(2,nm) = GA.' fieldname ';'])
   	FC(nm,:) = GA.freqConditions;
   	ALPH(nm,:) = GA.alph;
	str = ['if mod(nm,' num2str(incrementevery) ') == 0,GA.' fieldname ' = GA.' fieldname '+' num2str(incrementby) ';,end']
   	eval(str);
	eval(['paramvalue = GA.' fieldname ';'])
	disp('* *********************************************************************')
	disp('*  Varying parameter option on.')
	disp(['*  GA.' fieldname ' is now ' num2str(paramvalue)])
	disp('* *********************************************************************')
   end

   eval(['save model' num2str(nm) '_' M.date ' M'])
   save GAworkspace
   diary off
   fprintf('\n')
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function nmt_sourceoriplot(cfg)
% NMT_SOURCEORIPLOT
% plots functional source orientation data on slices or on
% a surface, optionally as an overlay on anatomical MRI data, where
% statistical data can be used to determine the opacity of the mask. Input
% data comes from FT_SOURCEANALYSIS.
%
% Use as
%   ft_sourceoriplot(cfg, data)
%
% No cfg options are implemented yet.
%
% Requires spm_ov_quivernmt.m to be installed as SPM plugin in
% spm_orthviews subdirectory

global st

% check that spm_ov_quivernmt has been installed in the SPM plug-in directory
if(~exist('spm_orthviews/spm_ov_quivernmt.m','file'))
    error(['Please move or link ' which('spm_ov_quivernmt') ' to: ' fileparts(which('spm_ov_reorient.m'))])
end

if(size(st.nmt.ori,3) > 1)
    ori = st.nmt.ori(:,:,st.nmt.cfg.time_idx(1));
else
    ori = st.nmt.ori;
end
ori = ori ./ sqrt(sum(ori.^2,2)); % normalize

fname{1} = fullfile(tempdir,'oriX.nii');
fname{2} = fullfile(tempdir,'oriY.nii');
fname{3} = fullfile(tempdir,'oriZ.nii');
fname{4} = fullfile(tempdir,'orimask.nii');

for ii=1:3
    oriwrapped{ii} = st.vols{1}.blobs{1}.vol;
    oriwrapped{ii}(st.nmt.cfg.inside_idx) = ori(:,ii);
end
oriwrapped{4} = double(isfinite(oriwrapped{1}));

for ii=1:4
    vol.fname = fname{ii};
    vol.dim = size(st.vols{1}.blobs{1}.vol);
    vol.mat = st.vols{1}.blobs{1}.mat;
    vol.dt = [16 0];
    spm_write_vol(vol,oriwrapped{ii});
end

%%
for ii=1:3
    orivol{ii}=spm_vol(fname{ii});
end
orimaskvol=spm_vol(fname{4});
spm_orthviews('quivernmt','delete',1);
spm_orthviews('quivernmt','init',1,orivol,orimaskvol)
spm_orthviews('redraw')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function net  = merge_batch_norm(net)
% MERGE_BATCH_NORM function merges Batch Normalization layers into preceding Conv layers.
%
%   NET = merge_batch_norm(NET)
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

  % find BatchNorm layers
  names = dagFindLayersOfType(net, 'dagnn.BatchNorm') ;

  % if no BatchNorm return
  if isempty(names), return;  end

  fprintf('>> Merging batch norm layers with the preceding conv layers\n');
  for name = names
    name = char(name) ;
    layer = net.layers(net.getLayerIndex(name)) ;

    % merge into previous conv layer
    playerName = dagFindLayersWithOutput(net, layer.inputs{1}) ;
    playerName = playerName{1} ;
    playerIndex = net.getLayerIndex(playerName) ;
    player = net.layers(playerIndex) ;
    if ~isa(player.block, 'dagnn.Conv')
      error('!!! Batch normalization cannot be merged as it is not preceded by a conv layer.') ;
    end

    % if the convolution layer does not have a bias,
    % recreate it to have one
    if ~player.block.hasBias
      block = player.block ;
      block.hasBias = true ;
      net.renameLayer(playerName, 'tmp') ;
      net.addLayer(playerName, ...
                   block, ...
                   player.inputs, ...
                   player.outputs, ...
                   {player.params{1}, sprintf('%s_b',playerName)}) ;
      net.removeLayer('tmp') ;
      playerIndex = net.getLayerIndex(playerName) ;
      player = net.layers(playerIndex) ;
      biases = net.getParamIndex(player.params{2}) ;
      net.params(biases).value = zeros(block.size(4), 1, 'single') ;
    end

    filters = net.getParamIndex(player.params{1}) ;
    biases = net.getParamIndex(player.params{2}) ;
    multipliers = net.getParamIndex(layer.params{1}) ;
    offsets = net.getParamIndex(layer.params{2}) ;
    moments = net.getParamIndex(layer.params{3}) ;

    [filtersValue, biasesValue] = mergeBatchNorm(...
      net.params(filters).value, ...
      net.params(biases).value, ...
      net.params(multipliers).value, ...
      net.params(offsets).value, ...
      net.params(moments).value) ;

    net.params(filters).value = filtersValue ;
    net.params(biases).value = biasesValue ;
  end

  % Remove Batch Norm layers, they are already merged
  names = dagFindLayersOfType(net, 'dagnn.BatchNorm') ;
  for i = 1:numel(names)
    layer = net.layers(net.getLayerIndex(names{i})) ;
    net.removeLayer(names{i}) ;
    net.renameVar(layer.outputs{1}, layer.inputs{1}, 'quiet', true) ;
  end

  % -------------------------------------------------------------------------
  function layers = dagFindLayersWithOutput(net, outVarName)
  % -------------------------------------------------------------------------
  layers = {} ;
  for l = 1:numel(net.layers)
    if any(strcmp(net.layers(l).outputs, outVarName))
      layers{1,end+1} = net.layers(l).name ;
    end
  end

  % -------------------------------------------------------------------------
  function layers = dagFindLayersOfType(net, type)
  % -------------------------------------------------------------------------
  layers = [] ;
  for l = 1:numel(net.layers)
    if isa(net.layers(l).block, type)
      layers{1,end+1} = net.layers(l).name ;
    end
  end

  % -------------------------------------------------------------------------
  function [filters, biases] = mergeBatchNorm(filters, biases, multipliers, offsets, moments)
  % -------------------------------------------------------------------------
  % wk / sqrt(sigmak^2 + eps)
  % bk - wk muk / sqrt(sigmak^2 + eps)
  a = multipliers(:) ./ moments(:,2) ;
  b = offsets(:) - moments(:,1) .* a ;
  biases(:) = biases(:) + b(:) ;
  sz = size(filters) ;
  numFilters = sz(4) ;
  filters = reshape(bsxfun(@times, reshape(filters, [], numFilters), a'), sz) ;

author: filipradenovic

repo: cnnimageretrieval


text:
function seq = prtUtilTopographicalSort(adj)
% [seq] = toposort(adj)  A Topological ordering of nodes in a directed graph
% INPUT:  adj  -  adjacency matrix
% OUTPUT: seq  -  a topological ordered sequence of nodes
%                 or an empty matrix if graph contains cycles
%
% Taken from Bayesian Network Toolbox by Minka

N = size(adj);
indeg = sum(adj,1);
outdeg = sum(adj,2);
seq = [];
for i = 1:N,
    idx = find(indeg==0);    % Find nodes with indegree 0
    if isempty(idx),   % If can't find than graph contains a cycle
        seq = [];
        break;
    end;
    [dummy, idx2] = max(outdeg(idx)); % Remove the node with the max number of connections
    indx = idx(idx2);
    seq = [seq, indx];
    indeg(indx) = -1;
    idx = find(adj(indx,:));
    indeg(idx) = indeg(idx)-1;
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function p00_story ( prob )

%*****************************************************************************80
%
%% P00_STORY prints the "story" for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  if ( prob == 1 )
    p01_story ( );
  elseif ( prob == 2 )
    p02_story ( );
  elseif ( prob == 3 )
    p03_story ( );
  elseif ( prob == 4 )
    p04_story ( );
  elseif ( prob == 5 )
    p05_story ( );
  elseif ( prob == 6 )
    p06_story ( );
  elseif ( prob == 7 )
    p07_story ( );
  elseif ( prob == 8 )
    p08_story ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_STORY - Fatal error!\n' );
    fprintf ( 1, '  Unexpected input value of PROB.\n' );
    error ( 'P00_STORY - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function EB=EBCM(I_prop1)
[m,n]=size(I_prop1);
h_sob=fspecial('sobel');
I_prop1_norm=double(I_prop1./255);
I_filt_norm=imfilter(I_prop1_norm,h_sob,'replicate');
gij_xij=(I_prop1_norm).*(I_filt_norm);
h_avg=fspecial('average'); %%%default size is 3 by 3
I_filt_avg=imfilter(I_filt_norm,h_avg,'replicate');
gij_xij_avg=imfilter(gij_xij,h_avg,'replicate');
eij=(gij_xij_avg./(I_filt_avg+0.0001));
cij=abs(I_prop1_norm-eij)./abs(I_prop1_norm+eij+0.0001);
cij_u=sum(sum(uint8(round(cij*255))));
EB=cij_u/(m*n);
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = recover(lmi_variables)
%RECOVER Create NCVAR object using variable indicies

if isempty(lmi_variables)
    varargout{1} = [];
else
    if isa(lmi_variables,'ncvar') | isa(lmi_variables,'lmi') | isa(lmi_variables,'constraint')
        varargout{1} = flush(recover(depends(lmi_variables)));
    else
        n = length(lmi_variables);
        i=1:n;
        if nargout <= 1

            lmi_variables = lmi_variables(:)';
            basis = sparse(i,i+1,ones(n,1),n,n+1);
            if any(diff(lmi_variables)<0)
                [i,j]=sort(lmi_variables);
                basis = [basis(:,1) basis(:,j+1)];
                lmi_variables = lmi_variables(j);
            end

            [un_Z_vars2] = uniquestripped(lmi_variables);
            if length(un_Z_vars2) < length(lmi_variables)
                [un_Z_vars,hh,jj] = unique(lmi_variables);
                if length(lmi_variables) ~= length(un_Z_vars)
                    basis = basis*sparse([1 1+jj(:)'],[1 1+(1:length(jj))],ones(1,1+length(jj)))';
                    lmi_variables = un_Z_vars;
                end
            end

            varargout{1} = ncvar(n,1,[],lmi_variables(:)',basis,0);

        else
            x = ncvar(n,1,[],lmi_variables(:)',sparse(i,i+1,ones(n,1),n,n+1),0);
            for i = 1:length(lmi_variables)
                varargout{i} = flush(x(i));
            end
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Johann Diep (johann.diep@esa.int) - November 2021
%
% This script builds the structure for the Katwijk GNSS dataset.

warning off;

%% Conversion Constants

G =  9.80665;       % Gravity constant, m/s^2
G2MSS = G;          % g to m/s^2
MSS2G = (1/G);      % m/s^2 to g
D2R = (pi/180);     % degrees to radians
R2D = (180/pi);     % radians to degrees
KT2MS = 0.514444;   % knot to m/s
MS2KMH = 3.6;       % m/s to km/h

%% GNSS ERROR PROFILE

gnss_data_planetary = readtable('./external-data/katwijk/gnss/gps-latlong.txt');
load('./katwijk-planetary-data/data/imu_planetary.mat');

% GNSS data structure:
%         t: Mx1 time vector (seconds).
%       lat: Mx1 latitude (radians).
%       lon: Mx1 longitude (radians).
%         h: Mx1 altitude (m).
%       vel: Mx3 NED velocities (m/s).
%       std: 1x3 position standard deviations, [lat lon h] (rad, rad, m).
%      stdm: 1x3 position standard deviations, [lat lon h] (m, m, m).
%      stdv: 1x3 velocity standard deviations, [Vn Ve Vd] (m/s).
%      larm: 3x1 lever arm from IMU to GNSS antenna (x-fwd, y-right, z-down) (m).
%      freq: 1x1 sampling frequency (Hz).
%   zupt_th: 1x1 ZUPT threshold (m/s).
%  zupt_win: 1x1 ZUPT time window (seconds).
%       eps: 1x1 time interval to compare IMU time vector to GNSS time vector (seconds).

gnss_planetary.stdm = 0.01 * [1 1 1]; % degradation

% GNSS positions
gnss_planetary.lat = table2array(gnss_data_planetary(:,3)).*D2R;
gnss_planetary.lon = table2array(gnss_data_planetary(:,4)).*D2R;
gnss_planetary.h = table2array(gnss_data_planetary(:,5));

gnss_planetary = gnss_m2r(gnss_planetary.lat(1),gnss_planetary.h(1),gnss_planetary); % convert error from meters to radians

% parsing the timestamps of the measurements
Date = cell2mat(table2array(gnss_data_planetary(:,1)));
Year = str2num(Date(:,1:4));
Month = str2num(Date(:,6:7));
Day = str2num(Date(:,9:10));
Hour = str2num(Date(:,12:13))-1;
Minutes = str2num(Date(:,15:16));
Seconds = str2num(Date(:,18:19));
MilliSeconds = str2num(Date(:,21:23));
DateTime = datetime(Year, Month, Day, Hour, Minutes, Seconds, MilliSeconds);
Timestamp = posixtime(DateTime);
gnss_planetary.t = Timestamp - Timestamp(1);

gnss_planetary.freq = 1/mean(diff(gnss_planetary.t)); % estimating the frequency
[gnss_planetary.vel,~] = pllh2vned(gnss_planetary); % estimating the velocity

gnss_planetary.larm = [0,0,0]'; % GNSS lever arm from IMU to GNSS antenna
gnss_planetary.eps = mean(diff(imu_planetary.t))/3; % rule of thumb for choosing eps
gnss_planetary.stdv = gnss_planetary.stdm/100; % educated guess

% Parameters for ZUPT detection algorithm  copied from above
gnss_planetary.zupt_th = 0; % ZUPT threshold (m/s).
gnss_planetary.zupt_win = 4; % ZUPT time window (seconds).

% Increase frequency by interpolation
t_q = 0:0.05:gnss_planetary.t(end);
lon_q = interp1(gnss_planetary.t,gnss_planetary.lon,t_q);
lat_q = interp1(gnss_planetary.t,gnss_planetary.lat,t_q);
h_q = interp1(gnss_planetary.t,gnss_planetary.h,t_q);
v_q = interp1(gnss_planetary.t,gnss_planetary.vel,t_q);
freq_q = 1/mean(diff(t_q));

gnss_planetary.t = t_q';
gnss_planetary.lon = lon_q';
gnss_planetary.lat = lat_q';
gnss_planetary.h = h_q';
gnss_planetary.vel = v_q;
gnss_planetary.freq = freq_q;

gnss_planetary_r = gnss_planetary; % reference

% remove datapoints
gnss_planetary.t(2000:4000) = [];
gnss_planetary.lon(2000:4000) = [];
gnss_planetary.lat(2000:4000) = [];
gnss_planetary.h(2000:4000) = [];
gnss_planetary.vel(2000:4000,:) = [];

gnss_planetary_sparse_r = gnss_planetary; % reference
[gnss_planetary,~] = gnss_gen(gnss_planetary_sparse_r, gnss_planetary); % degredation

%% Saving

save('./katwijk-planetary-data/data/gnss_planetary_r.mat','gnss_planetary_r');
save('./katwijk-planetary-data/data/gnss_planetary_sparse_r.mat','gnss_planetary_sparse_r');
save('./katwijk-planetary-data/data/gnss_sparse_planetary.mat','gnss_planetary');

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function value = r8mat_maxcol_minrow ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_MAXCOL_MINROW gets the maximum column minimum row of an R8MAT.
%
%  Discussion:
%
%    value = max ( 1 <= I <= N ) ( min ( 1 <= J <= M ) A(I,J) )
%
%    For a given matrix, MAT_MAXCOL_MINROW <= MAT_MINROW_MAXCOL
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the matrix.
%
%    Output, real VALUE, the maximum column
%    minimum row entry of A.
%
  value = max ( min ( a(1:m,1:n)' ) );

  return
end


text:
function [exchModel,unusedMets] = setExchangeBounds(model,mets,lb,ub,closeOthers,mediaOnly)
% setExchangeBounds
%   Define the exchange flux bounds for a given set of metabolites.
%
% Input:
%   model         a model structure
%   mets          a cell array of metabolite names (case insensitive) or
%                 metabolite IDs, or a vector of metabolite indices
%                 (opt, default all exchanged metabolites)
%   lb            lower bound of exchange flux. Can be either a vector of
%                 bounds corresponding to each of the provided metabolites,
%                 or a single value that will be applied to all.
%                 (opt, default to model.annotation.defaultLB if it exists,
%                 otherwise -1000)
%   ub            upper bound of exchange flux. Can be either a vector of
%                 bounds corresponding to each of the provided metabolites,
%                 or a single value that will be applied to all.
%                 (opt, default to model.annotation.defaultUB if it exists,
%                 otherwise 1000)
%   closeOthers   close exchange reactions for all other exchanged
%                 metabolites not present in the provided list. This will
%                 prevent IMPORT of the metabolites, but their EXPORT will
%                 not be modified.
%                 (opt, default true)
%   mediaOnly     only consider exchange reactions involving exchange to or
%                 from the extracellular (media) compartment. Reactions
%                 such as "sink" reactions that exchange metabolites
%                 directly with an intracellular compartment will therefore
%                 be ignored even though "getExchangeRxns" identifies such
%                 such reactions as exchange reactions.
%                 Note: The function will attempt to identify the
%                 extracellular compartment by the "compNames" field, and
%                 also requires the "metComps" field to be present,
%                 otherwise the mediaOnly flag will be ignored.
%                 (opt, default false)
%
% Output:
%   exchModel     a model structure with updated exchange flux bounds for
%                 the provided set of metabolites
%   unusedMets    metabolites provided by the user that were not used
%                 because they are not involved in any exchange reactions
%                 in the model
%
% NOTE: Exchange reactions involving more than one metabolite will be
% ignored.
%
% Usage: exchModel = setExchangeBounds(model,mets,lb,ub,closeOthers,mediaOnly);

% handle input arguments
if nargin < 2
    mets = [];
elseif ~islogical(mets) || ~isnumeric(mets)
    mets=convertCharArray(mets);
end

if nargin < 3 || isempty(lb)
    if isfield(model,'annotation') && isfield(model.annotation,'defaultLB')
        lb = model.annotation.defaultLB;
    else
        lb = -1000;
    end
end

if nargin < 4 || isempty(ub)
    if isfield(model,'annotation') && isfield(model.annotation,'defaultUB')
        ub = model.annotation.defaultUB;
    else
        ub = 1000;
    end
end

if nargin < 5 || isempty(closeOthers)
    closeOthers = true;
end

if nargin < 6
    mediaOnly = false;
elseif mediaOnly
    if ~all(isfield(model,{'compNames','metComps'}))
        error('mediaOnly option requires the "compNames" and "metComps" model fields.');
    end
end

% for models with an "unconstrained" field, generate a version of the model
% where all the unconstrained metabolites have stoichiometric coeffs set to
% zero, but remain in the model to avoid index changes
model_temp = model;
if isfield(model,'unconstrained')
    model_temp.S(model_temp.unconstrained == 1,:) = 0;
end

% find exchange rxns, ignoring those involving more than one metabolite
[~,exchRxnInd] = getExchangeRxns(model);
exchRxnInd(sum(model_temp.S(:,exchRxnInd) ~= 0,1) > 1) = [];

% find all exchanged metabolites
exchMetInd = find(any(model_temp.S(:,exchRxnInd) ~= 0,2));

if ( mediaOnly )
    % ignore exchanged metabolites in non-extracellular compartments and
    % any exchange reactions involving these metabolites
    [~,extCompInd] = ismember('extracellular',lower(model.compNames));
    if extCompInd > 0
        ignoreMet = (model.metComps(exchMetInd) ~= extCompInd);
    else
        error('Could not find any compartments named "extracellular".');
    end
    ignoreRxn = any(model.S(exchMetInd(ignoreMet),exchRxnInd) ~= 0,1);
    exchMetInd(ignoreMet) = [];
    exchRxnInd(ignoreRxn) = [];
end

% Check that all exchange reactions are formulated in the same direction.
% If not, this means that negative flux indicates import for some exchange
% reactions, but indicates export for others. Therefore, the LB and UB
% would need to be specified differently depending on the exchange reaction
% direction, which is error-prone.
if all(all(model.S(exchMetInd,exchRxnInd) <= 0))
    importDir = 'backward';
elseif all(all(model.S(exchMetInd,exchRxnInd) >= 0))
    importDir = 'forward';
else
    fprintf('WARNING: Some exchange reactions differ in direction, and therefore have opposite meanings of LB and UB.\n');
    if closeOthers
        fprintf('         Therefore, the "closeOthers" option will be set to FALSE.\n');
    end
    closeOthers = false;
end

% prepare exchanged metabolites and bounds
if ~isempty(mets)

    % prepare lb and ub variables
    if numel(lb) == 1
        lb = lb*ones(size(mets));
    elseif numel(lb) ~= numel(mets)
        error('lb must be a single value or a vector of equal length as mets');
    end
    if numel(ub) == 1
        ub = ub*ones(size(mets));
    elseif numel(ub) ~= numel(mets)
        error('ub must be a single value or a vector of equal length as mets');
    end

    % map provided mets to exchanged metabolites
    if isnumeric(mets)
        % mets are provided as indices
        [isExch,metInd] = ismember(mets,exchMetInd);
    elseif sum(ismember(lower(mets),lower(model.metNames))) > sum(ismember(mets,model.mets))
        % assume that mets are provided as names
        [isExch,metInd] = ismember(lower(mets),lower(model.metNames(exchMetInd)));
    else
        % assume that mets are provided as met IDs
        [isExch,metInd] = ismember(mets,model.mets(exchMetInd));
    end

    % get exchanged met indices and determine unused (non-exchanged) mets
    unusedMets = mets(~isExch);
    exchMetInd = exchMetInd(metInd(isExch));

    % update bound vectors
    lb = lb(isExch);
    ub = ub(isExch);

else
    % if no mets were provided, use all exchanged mets
    if numel(lb) > 1 || numel(ub) > 1
        error('Only one upper and one lower bound may be provide if metabolites are not specified.');
    else
        lb = lb*ones(size(exchMetInd));
        ub = ub*ones(size(exchMetInd));
    end
    unusedMets = {};
end

% check that at least one exchanged metabolite matches
if isempty(exchMetInd)
    fprintf('None of the provided metabolites were found in any exchange reactions.\n');
    exchModel = model;
    return
end

% determine which metabolite is exchanged in each exchange reaction
[metInd,rxnInd] = find(model_temp.S(exchMetInd,exchRxnInd) ~= 0);

% check for any metabolites that are exchanged in more than one reaction
[tbl,i,~] = unique(metInd,'first');
repeatedInds = find(not(ismember(1:numel(tbl),i)));
multiMetInd = exchMetInd(metInd(repeatedInds));
if ~isempty(multiMetInd)
    fprintf('WARNING: The following metabolites are involved in more than one exchange reaction:\n');
    fprintf('\t%s\n',model.metNames{multiMetInd(1:min(numel(multiMetInd),10))});
    if numel(multiMetInd) > 10
        fprintf('\t... and %u more.\n',numel(multiMetInd)-10);
    end
end

% set exchange reaction bounds
model = setParam(model,'lb',exchRxnInd(rxnInd),lb(metInd));
model = setParam(model,'ub',exchRxnInd(rxnInd),ub(metInd));

if closeOthers
    % constrain import of all other exchange reactions to zero
    constrainInd = setdiff(exchRxnInd,exchRxnInd(rxnInd));
    if strcmp(importDir,'backward')
        model = setParam(model,'lb',constrainInd,0);
    else
        model = setParam(model,'ub',constrainInd,0);
    end
end

% assign output model
exchModel = model;

author: SysBioChalmers

repo: RAVEN


text:
function mean = log_normal_mean ( a, b )

%*****************************************************************************80
%
%% LOG_NORMAL_MEAN returns the mean of the Lognormal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = exp ( a + 0.5 * b * b );

  return
end


text:
function y = VBA_sample(form,suffStat,N)
% legacy code
s = warning ('on');
warning ('*** The function `VBA_sample` is now deprecated. Please see `VBA_random` for an alternative.')
warning (s);

if nargin < 3
    N = 1;
end

switch form
    case 'gaussian'
        if isscalar(suffstat.mu)
            N = {1, N};
        else
            N = {N};
        end

        y = VBA_random ('Gaussian', suffStat.mu, suffStat.Sigma, N{:});

    case 'gamma'
        y = VBA_random ('Gamma', suffStat.a, suffStat.b, 1, N);

    case 'dirichlet'
        y = VBA_random ('Dirichlet', suffStat.d, N);

     case 'bernoulli'
        y = VBA_random ('Bernoulli', suffStat.p, 1, N);

    case 'binomial'
        y = VBA_random ('Binomial', suffStat.n, suffstat.p, 1, N);

    case 'multinomial'
        y = VBA_random ('Multinomial', suffStat.n, suffstat.p, N);

end

author: MBB-team

repo: VBA-toolbox


text:
function [outputMovie, movieDims, nPixels, nFrames] = loadMovieList(movieList, varargin)
	% Load movies, automatically detects type (avi, tif, or hdf5) and concatenates if multiple movies in a list.
		% NOTE:
			% The function assumes input is 2D time series movies with [x y frames] as dimensions
			% If movies are different sizes, use largest dimensions and align all movies to top-left corner.
	% Biafra Ahanonu
	% started: 2013.11.01

	[outputMovie, movieDims, nPixels, nFrames] = ciapkg.io.loadMovieList(movieList, 'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
classdef SSC_OMP < handle
    % Implements sparse subspace clustering algorithm using OMP algorithm
    properties
        Quiet = false
        RepresentationMethod
        % Options to be passed on to the solver which is computing the representations
        RepSolverOptions = struct
    end

    properties(SetAccess=private)
        % Dimension of signal space
        N
        % Number of signals in the space
        S
        % A matrix of size NxS where each column is one signal vector
        Data
        % The sparsity level or the largest dimension of the sparse subspaces
        K
        % The expected number of subspaces
        NumSubspaces
        % Labels obtained through clustering
        Labels
        % Representation matrix (each column is a representation vector for one data vector)
        Representation
        % Adjacency matrix
        Adjacency
        % The spectral clusterer used in the algorithm
        Clusterer
        % Threshold for residual norm
        ResidualNormThreshold
        % Number of iterations for each vector
        Iterations
    end

    methods
        function self = SSC_OMP(X, K, NumSubspaces, ...
            ResidualNormThreshold, RepresentationMethod)
            % Constructor
            self.Data = X;
            self.K = K;
            if nargin < 3
                NumSubspaces = -1;
            end
            self.NumSubspaces = NumSubspaces;
            if nargin < 4
                ResidualNormThreshold = 1e-3;
            end
            self.ResidualNormThreshold = ResidualNormThreshold;
            if nargin < 5
                RepresentationMethod = spx.cluster.ssc.OMP_REPR_METHOD.FLIPPED_OMP_MATLAB;
            end
            self.RepresentationMethod = RepresentationMethod;
            [n, s] = size(X);
            self.N = n;
            self.S = s;
            self.Labels = ones(s, 1);
            % Each data vector is represented using other data vectors in same space
            self.Representation = zeros(s, s);
        end

        function result = solve(self)
            % prepare sparse representations
            tstart = tic;
            self.recover_coefficients();
            representation_time = toc(tstart);
            self.build_adjacency();

            % conduct spectral clustering
            tstart = tic;
            result = spx.cluster.spectral.simple.normalized_symmetric_sparse(self.Adjacency, self.NumSubspaces);
            cluster_labels = result.labels;
            result.clustering_time = toc(tstart);
            % We are disabling our version of spectral clustering for now.
            % clusterer = spx.cluster.spectral.Clustering(self.Adjacency);
            % keep reference for debugging purposes.
            % self.Clusterer = clusterer;
            % clusterer.NumClusters = self.NumSubspaces;
            % cluster_labels = clusterer.cluster_random_walk();
            self.Labels = cluster_labels;
            result.Labels = self.Labels;
            result.Z = self.Representation;
            result.W = self.Adjacency;
            result.representation_time = representation_time;
        end

    end

    methods(Access=private)

        function recover_coefficients(self)
            % Computes sparse representations of the data vectors
            data_matrix = self.Data;
            % Number of data vectors
            ns = self.S;
            % sparsity level
            nk = self.K;
            quiet = self.Quiet;
            rnorm_thr = self.ResidualNormThreshold;

            if self.RepresentationMethod.isClassicOMP_C()
                % Classic OMP C version
                C = spx.fast.omp_spr(data_matrix, nk, rnorm_thr);
                self.Representation = C;
            elseif self.RepresentationMethod.isBatchOMP_C()
                % Batch OMP C version
                C = spx.fast.batch_omp_spr(data_matrix, nk, rnorm_thr);
                self.Representation = C;
            elseif self.RepresentationMethod.isFlippedOMP_MATLAB()
                % flipped OMP MATLAB version
                [representations, iterations] = spx.cluster.ssc.flipped_omp(...
                    data_matrix, nk, rnorm_thr, quiet);
                self.Representation = representations;
                self.Iterations = iterations;
            elseif self.RepresentationMethod.isBatchFlippedOMP_MATLAB()
                C = spx.cluster.ssc.batch_flipped_omp(data_matrix, nk, rnorm_thr, quiet);
                self.Representation = C;
            elseif self.RepresentationMethod.isBatchFlippedOMP_C()
                C = spx.fast.batch_flipped_omp_spr(data_matrix, nk, rnorm_thr);
                self.Representation = C;
            elseif self.RepresentationMethod.isGOMP_C()
                nl  = 2;
                options.verbose = 0;
                C = spx.fast.gomp_spr(data_matrix, nk, nl, rnorm_thr, options);
                self.Representation = C;
            elseif self.RepresentationMethod.isMC_OMP()
                options = self.RepSolverOptions;
                options.quiet = quiet;
                C = spx.cluster.ssc.mc_omp(...
                    data_matrix, nk, rnorm_thr, options);
                self.Representation = C;
            else
                error('Invalid representation method.');
            end
            if ~quiet
                fprintf('\n');
            end
        end

        function detect_outliers(self)
            % Identifies outliers in the representations
        end

        function build_adjacency(self)
            C = abs(self.Representation);
            % Normalize the matrix by column wise maximums
            C = spx.norm.normalize_linf(C);
            % disp(C(:, 1));
            % Make it symmetric
            C = C + C';
            % Keep it
            self.Adjacency = C;
            % disp(self.Adjacency(:, 1));
        end

    end

end

author: indigits

repo: sparse-plex


text:
function M = grassmannfactory(n, p, k)
% Returns a manifold struct to optimize over the space of vector subspaces.
%
% function M = grassmannfactory(n, p)
% function M = grassmannfactory(n, p, k)
%
% Grassmann manifold: each point on this manifold is a collection of k
% vector subspaces of dimension p embedded in R^n.
%
% The metric is obtained by making the Grassmannian a Riemannian quotient
% manifold of the Stiefel manifold, i.e., the manifold of orthonormal
% matrices, itself endowed with a metric by making it a Riemannian
% submanifold of the Euclidean space, endowed with the usual inner product.
% In short: it is the usual metric used in most cases.
%
% This structure deals with matrices X of size n x p x k (or n x p if
% k = 1, which is the default) such that each n x p matrix is orthonormal,
% i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for
% i = 1 : k if k > 1. Each n x p matrix is a numerical representation of
% the vector subspace its columns span.
%
% By default, k = 1.
%
% See also: stiefelfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%   March 22, 2013 (NB) : Implemented geodesic distance.
%   April 17, 2013 (NB) : Retraction changed to the polar decomposition, so
%                         that the vector transport is now correct, in the
%                         sense that it is compatible with the retraction,
%                         i.e., transporting a tangent vector G from U to V
%                         where V = Retr(U, H) will give Z, and
%                         transporting GQ from UQ to VQ will give ZQ: there
%                         is no dependence on the representation, which is
%                         as it should be. Notice that the polar
%                         factorization requires an SVD whereas the qfactor
%                         retraction requires a QR decomposition, which is
%                         cheaper. Hence, if the retraction happens to be a
%                         bottleneck in your application and you are not
%                         using vector transports, you may want to replace
%                         the retraction with a qfactor.
%   July  4, 2013 (NB)  : Added support for the logarithmic map 'log'.
%   July  5, 2013 (NB)  : Added support for ehess2rhess.
%   June 24, 2014 (NB)  : Small bug fix in the retraction, and added final
%                         re-orthonormalization at the end of the
%                         exponential map. This follows discussions on the
%                         forum where it appeared there is a significant
%                         loss in orthonormality without that extra step.
%                         Also changed the randvec function so that it now
%                         returns a globally normalized vector, not a
%                         vector where each component is normalized (this
%                         only matters if k>1).

    assert(n >= p, ...
           ['The dimension n of the ambient space must be larger ' ...
	        'than the dimension p of the subspaces.']);

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

    if k == 1
        M.name = @() sprintf('Grassmann manifold Gr(%d, %d)', n, p);
    elseif k > 1
        M.name = @() sprintf('Multi Grassmann manifold Gr(%d, %d)^%d', ...
                             n, p, k);
    else
        error('k must be an integer no less than 1.');
    end

    M.dim = @() k*p*(n-p);

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

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

    M.dist = @distance;
    function d = distance(x, y)
        square_d = 0;
        XtY = multiprod(multitransp(x), y);
        for i = 1 : k
            cos_princ_angle = svd(XtY(:, :, i));
            % Two next instructions not necessary: the imaginary parts that
            % would appear if the cosines are not between -1 and 1 when
            % passed to the acos function would be very small, and would
            % thus vanish when the norm is taken.
            % cos_princ_angle = min(cos_princ_angle,  1);
            % cos_princ_angle = max(cos_princ_angle, -1);
            square_d = square_d + norm(acos(cos_princ_angle))^2;
        end
        d = sqrt(square_d);
    end

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

    % Orthogonal projection of an ambient vector U to the horizontal space
    % at X.
    M.proj = @projection;
    function Up = projection(X, U)

        XtU = multiprod(multitransp(X), U);
        Up = U - multiprod(X, XtU);

    end

    M.tangent = M.proj;

	M.egrad2rgrad = M.proj;

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, H)
        PXehess = projection(X, ehess);
        XtG = multiprod(multitransp(X), egrad);
        HXtG = multiprod(H, XtG);
        rhess = PXehess - HXtG;
    end

    M.retr = @retraction;
    function Y = retraction(X, U, t)
        if nargin < 3
            t = 1.0;
        end
        Y = X + t*U;
        for i = 1 : k
            % We do not need to worry about flipping signs of columns here,
            % since only the column space is important, not the actual
            % columns. Compare this with the Stiefel manifold.
            % [Q, unused] = qr(Y(:, :, i), 0); %#ok
            % Y(:, :, i) = Q;

            % Compute the polar factorization of Y = X+tU
            [u, s, v] = svd(Y(:, :, i), 'econ'); %#ok
            Y(:, :, i) = u*v';
        end
    end

    M.exp = @exponential;
    function Y = exponential(X, U, t)
        if nargin == 3
            tU = t*U;
        else
            tU = U;
        end
        Y = zeros(size(X));
        for i = 1 : k
            [u s v] = svd(tU(:, :, i), 0);
            cos_s = diag(cos(diag(s)));
            sin_s = diag(sin(diag(s)));
            Y(:, :, i) = X(:, :, i)*v*cos_s*v' + u*sin_s*v';
            % From numerical experiments, it seems necessary to
            % re-orthonormalize. This is overall quite expensive.
            [q, unused] = qr(Y(:, :, i), 0); %#ok
            Y(:, :, i) = q;
        end
    end

    % Test code for the logarithm:
    % Gr = grassmannfactory(5, 2, 3);
    % x = Gr.rand()
    % y = Gr.rand()
    % u = Gr.log(x, y)
    % Gr.dist(x, y) % These two numbers should
    % Gr.norm(x, u) % be the same.
    % z = Gr.exp(x, u) % z needs not be the same matrix as y, but it should
    % v = Gr.log(x, z) % be the same point as y on Grassmann: dist almost 0.
    M.log = @logarithm;
    function U = logarithm(X, Y)
        U = zeros(n, p, k);
        for i = 1 : k
            x = X(:, :, i);
            y = Y(:, :, i);
            ytx = y.'*x;
            At = y.'-ytx*x.';
            Bt = ytx\At;
            [u, s, v] = svd(Bt.', 'econ');

            u = u(:, 1:p);
            s = diag(s);
            s = s(1:p);
            v = v(:, 1:p);

            U(:, :, i) = u*diag(atan(s))*v.';
        end
    end

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

    M.rand = @random;
    function X = random()
        X = zeros(n, p, k);
        for i = 1 : k
            [Q, unused] = qr(randn(n, p), 0); %#ok<NASGU>
            X(:, :, i) = Q;
        end
    end

    M.randvec = @randomvec;
    function U = randomvec(X)
        U = projection(X, randn(n, p, k));
        U = U / norm(U(:));
    end

    M.lincomb = @lincomb;

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

    % This transport is compatible with the polar retraction.
    M.transp = @(x1, x2, d) projection(x2, d);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [n, p, k]);
    M.vecmatareisometries = @() true;

end

% Linear combination of tangent vectors
function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>

    if nargin == 3
        d = a1*d1;
    elseif nargin == 5
        d = a1*d1 + a2*d2;
    else
        error('Bad use of grassmann.lincomb.');
    end

end

author: andrewssobral

repo: lrslibrary


text:
classdef DNNSGAII < ALGORITHM
% <multi> <real/integer> <multimodal>
% Decision space based niching NSGA-II

%------------------------------- Reference --------------------------------
% J. Liang, C. Yue, and B. Qu. Multimodal multi-objective optimization: A
% preliminary study, Proceedings of the IEEE Congress on Evolutionary
% Computation, 2016, 2454-2461.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection_Mod(round(Problem.N/2),round(Problem.N/2),Population.decs,FrontNo,-CrowdDis);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                [Population,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef LibrationPointEnum < matlab.mixin.SetGet
    %LibrationPointEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        L1('L1')
        L2('L2');
    end

    properties
        name char = '';
    end

    methods
        function obj = LibrationPointEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function [listBoxStr, enums] = getListBoxStr()
            m = enumeration('LibrationPointEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
            enums = m;
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('LibrationPointEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('LibrationPointEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


function [samples] = to_qpsk(bitvec)
    assert(iscolumn(bitvec) || isrow(bitvec), "Bit vector must be a row or column vector");
    assert(mod(length(bitvec), 2) == 0, "Length of bit vector must be an even number")

    % Keep the power of each point at 1.0
    amplitude = 1/sqrt(2);

    samples = zeros(length(bitvec) / 2, 1);

    for idx=1:2:length(bitvec)
        bits = bitvec(idx:idx+1);

        if (bits == [0, 0])
            sample = amplitude + amplitude * 1j;
        elseif (bits == [0 1])
            sample = amplitude - amplitude * 1j;
        elseif (bits == [1 0])
            sample = -amplitude + amplitude * 1j;
        else
            sample = -amplitude - amplitude * 1j;
        end

        samples((idx + 1) / 2) = sample;
    end
end


text:
% kaiserbeta() - Estimate Kaiser window beta
%
% Usage:
%   >> beta = pop_kaiserbeta(dev);
%
% Inputs:
%   dev       - scalar maximum passband deviation/ripple
%
% Output:
%   beta      - scalar Kaiser window beta
%
% References:
%   [1] Proakis, J. G., & Manolakis, D. G. (1996). Digital Signal
%       Processing: Principles, Algorithms, and Applications (3rd ed.).
%       Englewood Cliffs, NJ: Prentice-Hall
%
% Author: Andreas Widmann, University of Leipzig, 2005

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2005-2014 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
% $Id$

function [ beta ] = kaiserbeta(dev)

devdb = -20 * log10(dev);

if devdb > 50
    beta = 0.1102 * (devdb - 8.7);
elseif devdb >= 21
    beta = 0.5842 * (devdb - 21)^0.4 + 0.07886 * (devdb - 21);
else
    beta = 0;
end

end


text:
function result = GR_PF_DMtest(y1,bench,qn)

% INPUT: squared forecast errors from model (y1) and benchmark (bench); qn is bandwidth in NW
% OUTPUT: pval, the p-value of the test
% This is valid with NaN values

a=isfinite(y1); c=isfinite(bench);
ac=a+c; ac=find(ac==2);
y1=y1(ac); bench=bench(ac);
kk1=size(y1,2);
teststatv=[]; pval=[];
for i1=1:kk1
        P = length(y1);
%         y=(y1(:,i1)-true).^2;
        y=y1(:,i1)-bench;       % 2018_06_18: loss difference
        variance = nw(y,qn)/P;
        teststat = mean(y)/sqrt(variance);
        teststatv = [teststatv; teststat];
        pval = [pval; 1-cdf('chi2',teststat^2,1)];
end
result.teststat=teststatv;
result.pval=pval;

function result = nw(y,qn)
%input: y is a T*k vector and qn is the truncation lag
%output: the newey west HAC covariance estimator
%Formulas are from Hayashi
[T,k]=size(y); ybar=ones(T,1)*((sum(y))/T);
dy=y-ybar;
G0=dy'*dy/T;
for j=1:qn-1
   gamma=(dy(j+1:T,:)'*dy(1:T-j,:))./(T-1);
   G0=G0+(gamma+gamma').*(1-abs(j/qn));
end
result=G0;

author: european-central-bank

repo: BEAR-toolbox


text:
%{
Run this file to compile the mex files if you need to
(or if the binaries work for your computer, then you can skip this)
This file also contains two simple test problems

Stephen Becker, Feb 14, 2012  srbecker@alumni.caltech.edu
Updated         May  3, 2012  adding the f2c version for Windows userss
Updated         Jun 24, 2012  adding a lot of trouble-shooting support mainly for 64-bit linux
MAJOR Update    Feb 21, 2015  stephen.becker@colorado.edu
    completely redid mex file and lbfgs library
    (converted LBFGSB fortran code to C). Many changes. Most of the Matlab
    interface remains the same, but mex interface has changed.
    Should not affect end-user. Compilation is SO MUCH EASIER than it
    was with Fortran
    Most of the C code is the same as the Fortran, but a few fcn signatures
    different, and chanted "task" from str to int, and modified print functions
    No longer able to print to file

See also lbfgsb.m and lbfgsb_wrapper.c
%}

%% -- Compile --

% run mex -setup if you haven't already...

SRC_DIR = fullfile('..','src');
% Note: make sure to run be in the directory containing this file
% when you run this, otherwise above path is wrong
INCLUDES = SRC_DIR;
SRC = {'lbfgsb.c','linesearch.c','subalgorithms.c','print.c',...
    'linpack.c','miniCBLAS.c','timer.c'};
for i = 1:length(SRC)
    SRC{i} = fullfile(SRC_DIR,SRC{i});
end

if ispc
    % do not specify -lm flag
    mex('lbfgsb_wrapper.c','-largeArrayDims','-UDEBUG',...
        ['-I',SRC_DIR], SRC{:} );
else
    mex('lbfgsb_wrapper.c','-largeArrayDims','-lm','-UDEBUG',...
        ['-I',SRC_DIR], SRC{:} );
end

%% test the new function
disp('=== lbfgsb "driver1" test problem (Rosenbrock, 25 dimensions) === ');
% Here's the test problem included with lbfgsb called 'driver1'
% (It's a version of the Rosenbrock test function)

n   = 25;

l   = ones(n,1); u = l;
odd = 1:2:n;
even= 2:2:n;
l(odd) = 1.0;
u(odd) = 1.0e2;
l(even)= -1.0e2;
u(even)=  1.0e2;

opts    = struct( 'x0', 3*ones(n,1) );
opts.printEvery     = 2; % controls how often we print output from .m file wrapper
opts.m  = 5;
%opts.maxIts = 10;
opts.errFcn = @(x) norm(x-1); % for now just an arbitrary fcn
% opts.verbose = -1; % default is -1, i.e., no output from mex

[x,f,info] = lbfgsb( @driver1, l, u, opts );

% The true objective value is 0.
if abs(f) < 1e-8
    disp('Success!');
    semilogy( abs(info.err(:,1)-f),'o-' );
    xlabel('iteration');
    ylabel('error in objective function');
else
    disp('Something didn''t work right :-(  ');
end

% the structure info.err contains the objective function (1st column)
%   and norm(gradient,Inf) (2nd column)

%% another test function, the 2D Rosenbrock function
disp('=== Rosenbrock test function, 2D === ');
n = 2;

fxy = @(x,y) 100*( y-x.^2).^2  +  (1-x ).^2 ;
f   = @(x)   fxy( x(1,:), x(2,:) );
gxy = @(x,y) [100*(4*x.^3-4*x.*y)+2*x-2; 100*(2*y-2*x.^2)];
g   = @(x)   gxy( x(1,:), x(2,:) );

% There are no constraints
l   = -inf(n,1);
u   = inf(n,1);

opts    = struct( 'x0', [-1.9;2] );
opts.printEvery     = 1;
opts.m  = 5;

% Here's an example of using an error function. For Rosenbrock,
%   we know the true solution, so we can measure the error at every
%   iteration:
trueSoln = [1;1];
% "errFcn" will be printed to the screen
opts.errFcn     = @(x) norm(x-trueSoln)/max(norm(trueSoln),1);
% "outputFcn" will save values in the "info" output
opts.outputFcn  = opts.errFcn;

% Ask for very high accuracy
opts.pgtol      = 1e-10;
opts.factr      = 1e3;

% The {f,g} is another way to call it
[x,f,info] = lbfgsb( {f,g} , l, u, opts );

if abs(f) < 1e-8
    disp('Success!');
% since we included opts.outputFcn, the info.err now has 3 columns.
%   The first 2 columns are the same as before; the 3rd column
%   is the output of our outputFcn
semilogy( info.err(:,3)-f,'o-' ); xlabel('iteration'); ylabel('relative error in iterate function');
else
    disp('Something didn''t work right :-(  ');
end

repo: Deep-ADMM-Net


function [YYact,XXact] = YXB_(YY,lags,constant_timetrend)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'YXB_' organizes the data in the form
% of Y = XB+E
% NO dummy observations
% constant and trends are at the end

% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018
% Revised, 9/11/2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<3
    constant = 1;
    timetrend = 0;
else
    constant  = constant_timetrend(1);
    timetrend = constant_timetrend(2);
end

nlags_   = lags;                 % number of lags   */
T0       = lags;                 % size of pre-sample */

nv      = size(YY,2);            %* number of variables */
nobs    = size(YY,1)-T0;         %* number of observations */

% Actual observations

YYact = YY(T0+1:T0+nobs,:);
XXact = zeros(nobs,nv*nlags_);
i = 1;

while (i <= nlags_)
    XXact(:,(i-1)*nv+1:i*nv) = YY(T0-(i-1):T0+nobs-i,:);
    i = i+1;
end

if constant
    % last column of XXact = constant
    XXact = [XXact ones(nobs,1)];
end

if timetrend
    % last column of XXact = constant
    XXact = [XXact (1:nobs)'];
end


text:
function [EccA] = vect_solveKepler(meanALL, eccALL)
%solveKepler Summary of this function goes here
%   Detailed explanation goes here.
    EccA = zeros(size(eccALL));
    tol = 1E-12;

    if(any(eccALL < 1))
        ecc = eccALL(eccALL < 1);
        mean = meanALL(eccALL < 1);

        mean = AngleZero2Pi(mean);

        E0 = zeros(size(mean));
        E1 = ones(size(mean)) * Inf;
        if(any(mean > pi))
            E0(mean > pi) = mean(mean > pi) - ecc(mean > pi);
        end
        if(any(mean <= pi))
            E0(mean <= pi) = mean(mean <= pi) + ecc(mean <= pi);
        end

        first = true;
        while(any(abs(E1 - E0) > tol))
            if(~first)
                E0 = E1;
            end
            E1 = E0 + (mean - E0 + ecc.*sin(E0))./(1-ecc.*cos(E0));
            if(first)
                first = false;
            end
        end
        EccA(eccALL < 1) = E1;
    end

    if(any(eccALL >= 1))
        ecc = eccALL(eccALL >= 1);
        mean = meanALL(eccALL >= 1);

        H0 = zeros(size(mean));
        H1 = ones(size(mean)) * Inf;

        if(any(ecc < 1.6))
            bool = (-pi < mean & mean < 0) | mean > pi;
            H0(bool) = mean(bool) - ecc(bool);
            H0(not(bool)) = mean(not(bool)) + ecc(not(bool));
        end
        if(any(ecc >= 1.6))
            bool = ecc<3.6 & abs(mean) > pi;
            H0(bool) = mean(bool) - sign(mean(bool)).*ecc(bool);
            H0(not(bool)) = mean(not(bool))./(ecc(not(bool))-1);
        end

        first = true;
        while(any(abs(H1 - H0) > tol))
            if(~first)
                H0 = H1;
            end
            H1 = H0 + (mean - ecc.*sinh(H0) + H0)./(ecc.*cosh(H0) - 1);

            if(any(isnan(H1))) %this code necessary to prevent Inf/Inf = NaN cases where H0 is very large and the sinh/cosh functions give Inf as an answer
                bool1 = isnan(H1) & -ecc.*sinh(H0) == Inf;
                bool2 = isnan(H1) & -ecc.*sinh(H0) == -Inf;

                H1(bool1) = H0(bool1) + ones(size(H0(bool1)));
                H1(bool2) = H0(bool2) - ones(size(H0(bool2)));
            end

            if(first)
                first = false;
            end
        end
        EccA(eccALL >= 1) = H1;
    end
end


text:
function spline_test17 ( )

%*****************************************************************************80
%
%% TEST17 tests SPLINE_CUBIC_SET, SPLINE_CUBIC_VAL.
%
%  Discussion:
%
%    For boundary condition 0, the spline should come very close within
%    the interpolation interval.
%
%    For conditions 1 and 2, the spline should be essentially exactly equal
%    to the data, inside and outside the interpolation interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST17\n' );
  fprintf ( 1, '  SPLINE_CUBIC_SET sets up a cubic spline;\n' );
  fprintf ( 1, '  SPLINE_CUBIC_VAL evaluates it.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cubic data, unevenly spaced knots.\n' );

  for i = 1 : n
    t(i) = ( ( i - 1 ) / ( n - 1 ) )^2;
  end

  for i = 1 : n
    y(i) =  fcube ( t(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data to be interpolated:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of data values = %d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       T             Y\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '%14f  %14f\n', t(i), y(i) );
  end
%
%  Try boundary condition types 0, 1 and 2.
%
  for k = 0 : 2

    if ( k == 0 )

      ibcbeg = 0;
      ybcbeg = 0.0E+00;

      ibcend = 0;
      ybcend = 0.0E+00;

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Boundary condition 0 at both ends:\n' );
      fprintf ( 1, '  Spline is quadratic in boundary intervals.\n' );

    elseif ( k == 1 )

      ibcbeg = 1;
      ybcbeg = fpcube ( t(1) );

      ibcend = 1;
      ybcend = fpcube ( t(n) );

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Boundary condition 1 at both ends:\n' );
      fprintf ( 1, '  Y''(left) =  %f\n', ybcbeg );
      fprintf ( 1, '  Y''(right) = %f\n', ybcend );

    elseif ( k == 2 )

      ibcbeg = 2;
      ybcbeg = fppcube ( t(1) );

      ibcend = 2;
      ybcend = fppcube ( t(n) );

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Boundary condition 2 at both ends:\n' );
      fprintf ( 1, '  YP"(left) =  %f\n', ybcbeg );
      fprintf ( 1, '  YP"(right) = %f\n', ybcend );

    end

    ypp = spline_cubic_set ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SPLINE"(T), F"(T):\n' );
    fprintf ( 1, '\n' );
    for i = 1: n
      fprintf ( 1, '%14f  %14f\n', ypp(i), fppcube(t(i)) );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  T, SPLINE(T), F(T)\n' );
    fprintf ( 1, '\n' );

    for i = 0 : n

      if ( i == 0 )
        jhi = 1;
      elseif ( i < n )
        jhi = 2;
      else
        jhi = 2;
      end

      for j = 1 : jhi

        if ( i == 0 )
          tval = t(1) - 1.0E+00;
        elseif ( i < n )
          tval = ( ( jhi - j + 1 ) * t(i)     ...
                 + (       j - 1 ) * t(i+1) ) ...
                 / ( jhi );
        else
          if ( j == 1 )
            tval = t(n);
          else
            tval = t(n) + 1.0E+00;
          end
        end

        [ yval, ypval, yppval ] = spline_cubic_val ( n, t, y, ypp, tval );

        fprintf ( 1, '%14f  %14f  %14f\n', tval, yval, fcube ( tval ) );

      end
    end

  end

  return
end


text:
%TRANSL Create or unpack an SE3 translational transform
%
% Create a translational transformation matrix::
%
% T = TRANSL(X, Y, Z) is an SE(3) homogeneous transform (4x4) representing
% a pure translation of X, Y and Z.
%
% T = TRANSL(P) is an SE(3) homogeneous transform (4x4) representing a
% translation of P=[X,Y,Z]. If P (Mx3) it represents a sequence and T
% (4x4xM) is a sequence of homogeneous transforms such that T(:,:,i)
% corresponds to the i'th row of P.
%
% Unpack the translational part of a transformation matrix::
%
% P = TRANSL(T) is the translational part of a homogeneous transform T as a
% 3-element column vector.  If T (4x4xM) is a homogeneous transform
% sequence the rows of P (Mx3) are the translational component of the
% corresponding transform in the sequence.
%
% [X,Y,Z] = TRANSL(T) is the translational part of a homogeneous transform
% T as three components.  If T (4x4xM) is a homogeneous transform sequence
% then X,Y,Z (1xM) are the translational components of the corresponding
% transform in the sequence.
%
% Notes::
% - Somewhat unusually this function performs a function and its inverse.  An
%   historical anomaly.
%
% See also CTRAJ.

% Copyright (C) 1993-2015, 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

function [t1,t2,t3] = transl(x, y, z)
    if nargin == 1
        if ishomog(x)
            if ndims(x) == 3
                % transl(T)  -> P, trajectory case
                if nargout == 1
                    t1 = squeeze(x(1:3,4,:))';
                elseif nargout == 3
                    t1 = squeeze(x(1,4,:))';
                    t2 = squeeze(x(2,4,:))';
                    t3 = squeeze(x(3,4,:))';
                end
            else
                % transl(T)  -> P
                if nargout == 1 || nargout == 0
                    t1 = x(1:3,4);
                elseif nargout == 3
                    t1 = x(1,4);
                    t2 = x(2,4);
                    t3 = x(3,4);
                end

            end
        elseif length(x) == 3
            % transl(P) -> T
            t = x(:);
            t1 =    [eye(3)          t(:);
                0   0   0   1];
        else
            % transl(P) -> T, trajectory case
            n = numrows(x);
            t1 = repmat(eye(4,4), [1 1 n]);
            t1(1:3,4,:) = x';
        end
    elseif nargin == 3
        % transl(x,y,z) -> T
        t = [x; y; z];
        t1 =    rt2tr( eye(3), t);
    end

author: Allopart

repo: rbpf-gmapping


text:
function [ r, s1, s2, s3 ] = r8_random ( s1, s2, s3 )

%*****************************************************************************80
%
%% R8_RANDOM returns a pseudorandom number between 0 and 1.
%
%  Discussion:
%
%    This function returns a pseudo-random number rectangularly distributed
%    between 0 and 1.   The cycle length is 6.95E+12.  (See page 123
%    of Applied Statistics (1984) volume 33), not as claimed in the
%    original article.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 July 2008
%
%  Author:
%
%    Original FORTRAN77 version by Brian Wichman, David Hill.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Brian Wichman, David Hill,
%    Algorithm AS 183: An Efficient and Portable Pseudo-Random
%    Number Generator,
%    Applied Statistics,
%    Volume 31, Number 2, 1982, pages 188-190.
%
%  Parameters:
%
%    Input, integer S1, S2, S3, three values used as the
%    seed for the sequence.  These values should be positive
%    integers between 1 and 30,000.
%
%    Output, real R, the next value in the sequence.
%
%    Output, integer S1, S2, S3, updated seed values.
%
  s1 = mod ( 171 * s1, 30269 );
  s2 = mod ( 172 * s2, 30307 );
  s3 = mod ( 170 * s3, 30323 );

  r = mod ( s1 / 30269.0 ...
          + s2 / 30307.0 ...
          + s3 / 30323.0, 1.0 );

  return
end


text:
function [solution, LPProblem] = solveCobraLPCPLEX(LPProblem, printLevel, basisReuse, conflictResolve, contFunctName, minNorm, interface)
% Calls CPLEX to solve an LP problem
% By default, use the matlab interface to cplex written by TOMLAB, in
% preference to the one written by ILOG.
%
% USAGE:
%
%    [solution, LPProblem] = solveCobraLPCPLEX(LPProblem, printLevel, basisReuse, conflictResolve, contFunctName, minNorm, interface)
%
% INPUT:
%
%    LPProblem:          Structure containing the following fields describing the LP problem to be solved
%
%                          * .A - LHS matrix
%                          * .b - RHS vector
%                          * .c - Objective coeff vector
%                          * .lb - Lower bound vector
%                          * .ub - Upper bound vector
%                          * .osense - Objective sense (-1 max, +1 min)
%                          * .rxns - (optional) cell array of reaction abbreviations (necessary for
%                            making a readable confilict resolution file).
%                          * .csense - (optional) Constraint senses, a string containting the constraint sense for
%                            each row in `A` ('E', equality, 'G' greater than, 'L' less than).
%                          * .LPBasis - (optional) Basis from previous solution of similar LP problem.
%                            See `basisReuse`
% OPTIONAL INPUTS:
%    printLevel:         Printing level in the CPLEX m-file and CPLEX C-interface.
%
%                          * 0 - Silent
%                          * 1 - Warnings and Errors
%                          * 2 - Summary information (Default)
%                          * 3 - More detailed information
%                          * > 10 - Pause statements, and maximal printing (debug mode)
%    basisReuse:         0 - Use this for one of solution of an LP (Default);
%                        1 - Returns a basis for reuse in the next LP i.e. outputs `model.LPBasis`
%    conflictResolve:    0 (Default);
%                        1 If LP problem is proven to be infeasible by CPLEX,
%                        it will print out a 'conflict resolution file',
%                        which indicates the irreducible infeasible set of
%                        equaltiy & inequality constraints that together,
%                        combine to make the problem infeasible. This is
%                        useful for debugging an LP problem if you want to
%                        try to resolve a constraint conflict
%    contFunctName:      structure or function with parameters (only for `tomlab_cplex` or `ILOGcomplex`)
%
%                        - when using the `tomlab_cplex` interface
%
%                              1. contFunctName = [] Use all default CLPEX control parameters, (Default);
%                              2. contFunctName = someString e.g. 'someFunctionName'
%                                 uses the user specified control parameters defined
%                                 in `someFunctionName.m` (see template function CPLEXParamSet for details).
%                              3. contFunctName = `cpxControl` structure (output from a file like `CPLEXParamSet.m`)
%
%                        - when using the `ILOGcomplex` interface (parameter structure for Cplex). The full set of parameters can be obtained by calling `Cplex().Param`. For example:
%
%                              - `[solverParams.simplex.display, solverParams.tune.display, solverParams.barrier.display, solverParams.sifting.display, solverParams.conflict.display] = deal(0);`
%                              - `[solverParams.simplex.tolerances.optimality, solverParams.simplex.tolerances.feasibility] = deal(1e-9, 1e-8);`
%
%    minNorm:            {(0), 1 , `n x 1` vector} If not zero then, minimise the Euclidean length
%                        of the solution to the LP problem. Gives the same objective,
%                        but minimises the square of flux. `minNorm` ~1e-6 should be
%                        high enough for regularisation yet keep the same objective
%    interface:          {'ILOGcomplex', 'ILOGsimple', 'tomlab_cplex'}
%                        Default is the `tomlab_cplex` interface
%
% OUTPUT:
%    solution:           Structure containing the following fields describing a LP solution:
%
%                          * .full:               Full LP solution vector
%                          * .obj:                Objective value
%                          * .rcost:              Lagrangian multipliers to the simple inequalties (Reduced costs)
%                          * .dual:               Lagrangian multipliers to the equalities
%                          * .nInfeas:            Number of infeasible constraints
%                          * .sumInfeas:          Sum of constraint violation
%                          * .stat:               COBRA Standardized solver status code:
%
%                            * 1 - Optimal solution
%                            * 2 - Unbounded solution
%                            * 0 - Infeasible
%                            * -1 - No solution reported (timelimit, numerical problem etc)
%                          * .origStat:           CPLEX status code. Use `cplexStatus(solution.origStat)` for more information from the CPLEX solver
%                          * .solver              solver used by `cplex`
%                          * .time                time taken to solve the optimization problemtime taken to solve the optimization problem
%
% OPTIONAL OUTPUT:
%    LPProblem:          with field:
%
%                          * .LPBasis:            When input `basisReuse = 1`, we return a basis for reuse in the next LP
%
% CPLEX consists of 4 different LP solvers which can be used to solve sysbio optimization problems
% you can control which of the solvers, e.g. simplex vs interior point solver using the
% CPLEX control parameter cpxControl.LPMETHOD. At the moment, the solver is
% automatically chosen for you
%
% .. Note:
%       ILOG CPLEX parameters
%       https://www.ibm.com/support/knowledgecenter/SSSA5P_12.6.3/ilog.odms.studio.help/pdf/paramcplex.pdf
%
%       TOMLAB Cplex parameters
%       http://tomwiki.com/CPLEX_Parameter_Table
%
% .. Author: - Ronan Fleming

if ~exist('printLevel','var')
    printLevel=0;
end
if ~exist('basisReuse','var')
    basisReuse=0;
end
if ~exist('conflictResolve','var')
    conflictResolve=0;
end
if ~exist('interface','var')
    interface='tomlab_cplex';
end
if strcmp(interface,'tomlab_cplex')
    if ~exist('contFunctName','var')
        cpxControl=[];
    else
        if isstruct(contFunctName)
            cpxControl=contFunctName;
        else
            if ~isempty(contFunctName)
                %calls a user specified function to create a CPLEX control structure
                %specific to the users problem. A TEMPLATE for one such function is
                %CPLEXParamSet
                cpxControl=eval(contFunctName);
            else
                cpxControl=[];
            end
        end
    end
end
if ~exist('minNorm','var')
    minNorm=0;
end

if basisReuse
    if isfield(LPProblem, 'LPBasis')
        basis = LPProblem.LPBasis;
        % use advanced starting information when optimization is initiated.
        cpxControl.advance = 1;
        cpxControl.ADVIND = 1;
    else
        basis=[];
    end
else
    basis=[];
    % do not use advanced starting information when optimization is initiated.
    cpxControl.advance = 0;
    cpxControl.ADVIND = 0;
end

if ~isfield(LPProblem,'A')
    if ~isfield(LPProblem,'S')
        error('Equality constraint matrix must either be a field denoted A or S.')
    end
    LPProblem.A=LPProblem.S;
end

if ~isfield(LPProblem,'csense')
    nMet=size(LPProblem.A);
    if printLevel>0
        fprintf('%s\n','Assuming equality constraints, i.e. S*v=b');
    end
    %assuming equality constraints
    LPProblem.csense(1:nMet,1)='E';
end

if ~isfield(LPProblem,'osense')
    %assuming maximisation
    LPProblem.osense=-1;
    if printLevel>0
        fprintf('%s\n','Assuming maximisation of objective');
    end
end

if size(LPProblem.A,2)~=length(LPProblem.c)
    error('dimensions of A & c are inconsistent');
end

if size(LPProblem.A,2)~=length(LPProblem.lb) || size(LPProblem.A,2)~=length(LPProblem.ub)
    error('dimensions of A & bounds are inconsistent');
end

%get data
[c,x_L,x_U,b,csense,osense] = deal(LPProblem.c,LPProblem.lb,LPProblem.ub,LPProblem.b,LPProblem.csense,LPProblem.osense);
%modify objective to correspond to osense
c=full(c*osense);

%cplex expects it dense
b=full(b);
%Conflict groups descriptor (cpxBuildConflict can be used to generate the input). Set this if
%conflict refinement is desired in the case that infeasibility is detected
%by CPLEX.
if conflictResolve
    % set to deterministic mode to get reproducible conflict resolve file
    if (isfield(cpxControl,'PARALLEL') && cpxControl.PARALLEL ~=1) || (isfield(cpxControl,'PARALLELMODE') && cpxControl.PARALLELMODE ~=1)
        fprintf('PARALLEL / PARALLELMODE Parameter was changed to 1 to ensure a reproducible log file\n');
        cpxControl.PARALLEL = 1;
        cpxControl.PARALLELMODE = 1;
    end
    [m_lin,n]=size(LPProblem.A);
    m_quad=0;
    m_sos=0;
    m_log=0;
    %determines how elaborate the output is
    mode='full';%'minimal';
    fprintf('%s\n%s\n','Building Structure for Conflict Resolution...','...this slows CPLEX down so should not be used for repeated LP');
    confgrps = cpxBuildConflict(n,m_lin,m_quad,m_sos,m_log,mode);
    prefix=pwd;
    suffix='LP_CPLEX_conflict_file.txt';
    conflictFile=[prefix filesep suffix];
else
    confgrps=[]; conflictFile=[];
end

%Name of file to write the CPLEX log information to. If empty, no log is
%written.
logfile=[];

%Name of a file to save the CPLEX problem object (Used for submitting
%possible bugs in CPLEX to ILOG)
savefile=[]; savemode=[];
% savefile='C:\CPLEX_possible_bug.txt';

% vector defining which callbacks to use in CPLEX. If the ith entry of the logical vector
% callback is set, the corresponding callback is defined. The callback calls the m-file specified
% in Table 7 below. The user may edit this file, or make a new copy, which is put in a directory
% that is searched before the cplex directory in the Matlab path.
callback=[]; %I'm not really sure what this option means as yet

%this is not a tomlab problem so this is not needed
Prob=[];

% variables not used in LP problems
IntVars=[]; PI=[]; SC=[]; SI=[]; sos1=[]; sos2=[];

%quadratic constraint matrix, size n x n
if sum(minNorm)~=0
    if length(minNorm)==1
        % same weighting of min norm for all variables
        F=speye(length(c))*minNorm;
    else
        if length(minNorm)~=length(c)
            error('Either minNorm is a scalar, or is an n x 1 vector')
        else
            % individual weighting of min norm for all variables
            F=spdiags(minNorm,0,length(c),length(c));
        end
    end
else
    F=[];
end
%Structure array defining quadratic constraints
qc=[];

%Structure telling whether and how you want CPLEX to perform a sensitivity analysis (SA).
%This may be useful in future but probably will have more meaning with an
%additional term in the objective
saRequest =[];

%Vector with MIP starting solution, if known
xIP=[];

%Logical constraints, i.e. an additional set of single-sided linear constraints that are controlled
%by a binary variable (switch) in the problem
logcon=[];

%Report of incompatibility R2016b - ILOGcomplex interface
verMATLAB = version('-release');
if str2num(verMATLAB(1:end-1)) >= 2016 && strcmp(interface, 'ILOGcomplex')
    error(['MATLAB ',verMATLAB, ' and the ILOGcomplex interface are not compatible. Select ILOGsimple or tomlab_cplex as a CPLEX interface.'])
end

%call cplex
tic;
switch interface
    case 'ILOGcomplex'
        %complex ibm ilog cplex interface
        if ~isempty(csense)
            %set up constant vectors for CPLEX
            b_L(csense == 'E',1) = b(csense == 'E');
            b_U(csense == 'E',1) = b(csense == 'E');
            b_L(csense == 'G',1) = b(csense == 'G');
            b_U(csense == 'G',1) = Inf;
            b_L(csense == 'L',1) = -Inf;
            b_U(csense == 'L',1) = b(csense == 'L');
        else
            b_L = b;
            b_U = b;
        end

        % Initialize the CPLEX object
        try
            ILOGcplex = Cplex('fba');
        catch ME
            error('CPLEX not installed or licence server not up')
        end

        ILOGcplex.Model.sense = 'minimize';

        % Now populate the problem with the data
        ILOGcplex.Model.obj   = c;
        ILOGcplex.Model.lb    = x_L;
        ILOGcplex.Model.ub    = x_U;
        ILOGcplex.Model.A     = LPProblem.A;
        ILOGcplex.Model.lhs   = b_L;
        ILOGcplex.Model.rhs   = b_U;

        if ~isempty(F)
            %quadratic constraint matrix, size n x n
            ILOGcplex.Model.Q=F;
        end

        if ~exist('contFunctName','var')
            cpxControl= struct();
        else
            %Read ILOG cplex parameters
            ILOGcplex = setCplexParam(ILOGcplex, cpxControl, printLevel);
        end

        if printLevel==0
            ILOGcplex.DisplayFunc=[];
        else
            %print level
            ILOGcplex.Param.barrier.display.Cur = printLevel;
            ILOGcplex.Param.simplex.display.Cur = printLevel;
            ILOGcplex.Param.sifting.display.Cur = printLevel;
        end

        % Optimize the problem
        ILOGcplex.solve();
        %http://www-01.ibm.com/support/knowledgecenter/SSSA5P_12.2.0/ilog.odms.cplex.help/Content/Optimization/Documentation/CPLEX/_pubskel/CPLEX1210.html
        if ILOGcplex.Solution.status == 1
            solution.obj        = osense*ILOGcplex.Solution.objval;
            solution.full       = ILOGcplex.Solution.x;
            solution.rcost      = ILOGcplex.Solution.reducedcost;
            solution.dual       = ILOGcplex.Solution.dual;
            solution.nInfeas    = NaN;
            solution.sumInfeas  = NaN;
            %solution.stat       = ILOGcplex.Solution.
            solution.origStat   = ILOGcplex.Solution.status;
            solution.solver     = ILOGcplex.Solution.method;
            solution.time       = ILOGcplex.Solution.time;
            solution.kappa      = ILOGcplex.Solution.quality.kappa.value;
        else
            warning(['IBM CPLEX STATUS = ' int2str(ILOGcplex.Solution.status) ', see: http://www-01.ibm.com/support/knowledgecenter/SSSA5P_12.2.0/ilog.odms.cplex.help/Content/Optimization/Documentation/CPLEX/_pubskel/CPLEX1210.html'])
            solution.origStat   = ILOGcplex.Solution.status;
            solution.full       = NaN;
            solution.obj        = NaN;
            solution.rcost      = NaN;
            solution.dual       = NaN;
            solution.nInfeas    = NaN;
            solution.sumInfeas  = NaN;
            solution.solver     = NaN;
            solution.time       = NaN;
        end
    case 'ILOGsimple'
        try
            ILOGcplex = Cplex('fba');
        catch ME
            error('CPLEX not installed or licence server not up')
        end
        %simple ibm ilog cplex interface
        options = cplexoptimset;
        if printLevel == 0
            options = cplexoptimset(options,'Display','off');
        else
            options = cplexoptimset(options,'Display','on');
        end

        if ~isempty(csense)
            if sum(minNorm)~=0
                Aineq = [LPProblem.A(csense == 'L',:); - LPProblem.A(csense == 'G',:)];
                bineq = [b(csense == 'L',:); - b(csense == 'G',:)];
                %             min      0.5*x'*H*x+f*x or f*x
                %             st.      Aineq*x     <= bineq
                %             Aeq*x    = beq
                %             lb <= x <= ub
                [x,fval,exitflag,output,lambda] = cplexqp(F,c,Aineq,bineq,LPProblem.A(csense == 'E',:),b(csense == 'E',1),x_L,x_U,[],options);
            else
                Aineq = [LPProblem.A(csense == 'L',:); - LPProblem.A(csense == 'G',:)];
                bineq = [b(csense == 'L',:); - b(csense == 'G',:)];
                %        min      c*x
                %        st.      Aineq*x <= bineq
                %                 Aeq*x    = beq
                %                 lb <= x <= ub
                [x,fval,exitflag,output,lambda] = cplexlp(c,Aineq,bineq,LPProblem.A(csense == 'E',:),b(csense == 'E',1),x_L,x_U,[],options);
            end
            %primal
            solution.obj=osense*fval;
            solution.full=x;
            %this is the dual to the equality constraints but it's not the chemical potential
            solution.dual=lambda.eqlin;
        else
            Aineq=[];
            bineq=[];
            if sum(minNorm)~=0
                [x,fval,exitflag,output,lambda] = cplexqp(F,c,Aineq,bineq,LPProblem.A,b,x_L,x_U,[],options);
            else
                [x,fval,exitflag,output,lambda] = cplexlp(c,Aineq,bineq,LPProblem.A,b,x_L,x_U,[],options);
            end
            solution.obj=osense*fval;
            solution.full=x;
            %this is the dual to the equality constraints but it's not the chemical potential
            solution.dual=sparse(size(LPProblem.A,1),1);
            solution.dual(csense == 'E')=lambda.eqlin;
            %this is the dual to the inequality constraints but it's not the chemical potential
            solution.dual(csense == 'L')=lambda.ineqlin(1:nnz(csense == 'L'),1);
            solution.dual(csense == 'G')=lambda.ineqlin(nnz(csense == 'L')+1:end,1);
        end
        %this is the dual to the simple ineequality constraints : reduced costs
        solution.rcost=lambda.lower-lambda.upper;
        solution.nInfeas = [];
        solution.sumInfeas = [];
        solution.origStat = output.cplexstatus;
    case 'tomlab_cplex'
        %tomlab cplex interface
        if ~isempty(csense)
            %set up constant vectors for CPLEX
            b_L(csense == 'E',1) = b(csense == 'E');
            b_U(csense == 'E',1) = b(csense == 'E');
            b_L(csense == 'G',1) = b(csense == 'G');
            b_U(csense == 'G',1) = Inf;
            b_L(csense == 'L',1) = -Inf;
            b_U(csense == 'L',1) = b(csense == 'L');
        else
            b_L = b;
            b_U = b;
        end

        %tomlab cplex interface
        %   minimize   0.5 * x'*F*x + c'x     subject to:
        %      x             x_L <=    x   <= x_U
        %                    b_L <=   Ax   <= b_U
        [x, slack, v, rc, f_k, ninf, sinf, Inform, basis] = cplex(c, LPProblem.A, x_L, x_U, b_L, b_U, ...
            cpxControl, callback, printLevel, Prob, IntVars, PI, SC, SI, ...
            sos1, sos2, F, logfile, savefile, savemode, qc, ...
            confgrps, conflictFile, saRequest, basis, xIP, logcon);

        solution.full=x;
        %this is the dual to the equality constraints but it's not the chemical potential
        solution.dual=v*osense;%negative sign Jan 25th
        %this is the dual to the simple ineequality constraints : reduced costs
        solution.rcost=rc*osense;%negative sign Jan 25th
        if Inform~=1
            solution.obj = NaN;
        else
            if minNorm==0
                solution.obj=f_k*osense;
            else
                solution.obj=c'*x*osense;
            end
            %     solution.obj
            %     norm(x)
        end
        solution.nInfeas = ninf;
        solution.sumInfeas = sinf;
        solution.origStat = Inform;
    otherwise
        error([interface ' is not a recognised solveCobraLPCPLEX interface'])
end
solution.time=toc;
Inform = solution.origStat;

if Inform~=1 && conflictResolve ==1
    switch interface
        case {'ILOGcomplex','ILOGsimple'}
            if isfield(LPProblem,'mets') && isfield(LPProblem,'rxns')
                %this code reads the conflict resolution file and replaces the
                %arbitrary names with the abbreviations of metabolites and reactions
                [nMet,nRxn]=size(LPProblem.A);
                totAbbr=nMet+nRxn;
                conStrFind=cell(nMet+nRxn,1);
                conStrReplace=cell(nMet+nRxn,1);
                %only equality constraint rows
                for m=1:nMet
                    conStrFind{m,1}=['c' int2str(m) ':'];
                    conStrReplace{m,1}=[LPProblem.mets{m} ':  '];
                end
                %reactions
                for n=1:nRxn
                    conStrFind{nMet+n,1}=['x' int2str(n) ' '];
                    conStrReplace{nMet+n,1}=[LPProblem.rxns{n} ' '];
                end
                fid1 = fopen(suffix);
                fid2 = fopen(['COBRA_' suffix], 'w');
                while ~feof(fid1)
                    tline{1}=fgetl(fid1);
                    %replaces all occurrences of the string str2 within string str1
                    %with the string str3.
                    %str= strrep(str1, str2, str3)
                    for t=1:totAbbr
                        tline= strrep(tline, conStrFind{t}, conStrReplace{t});
                    end
                    fprintf(fid2,'%s\n', tline{1});
                end
                fclose(fid1);
                fclose(fid2);
                %delete other file without replacements
                %         delete(suffix)
                fprintf('%s\n',['Conflict resolution file written to: ' prefix '\COBRA_' suffix]);
                fprintf('%s\n%s\n','The Conflict resolution file gives an irreducible infeasible subset ','of constraints which are making this LP Problem infeasible');
            else
                warning('Need reaction and metabolite abbreviations in order to make a readable conflict resolution file');
            end
    end
else
    if printLevel>0 && Inform~=1
        fprintf('%s\n','No conflict resolution file. Consider to set conflictResolve = 1 next time.');
    end
end

if strcmp(interface, 'tomlab_cplex')
    % Try to give back COBRA Standardized solver status:
    %           1   Optimal solution
    %           2   Unbounded solution
    %           0   Infeasible
    %           -1  No solution reported (timelimit, numerical problem etc)
    if Inform==1
        solution.stat = 1;
        if printLevel>0
            %use tomlab code to print out exit meassage
            [ExitText,ExitFlag] = cplexStatus(Inform);
            solution.ExitText=ExitText;
            solution.ExitFlag=ExitFlag;
            fprintf('\n%s%g\n',[ExitText ', Objective '],  c'*solution.full*osense);
        end
    else
        if Inform==2
            solution.stat = 2;
            %use tomlab code to print out exit meassage
            [ExitText,ExitFlag] = cplexStatus(Inform);
            solution.ExitText=ExitText;
            solution.ExitFlag=ExitFlag;
            fprintf('\n%s%g\n',[ExitText ', Objective '],  c'*solution.full*osense);
        else
            if Inform==3
                solution.stat = 0;
            else
                %this is a conservative view
                solution.stat = -1;
                %use tomlab code to print out exit meassage
                [ExitText,ExitFlag] = cplexStatus(Inform);
                solution.ExitText=ExitText;
                solution.ExitFlag=ExitFlag;
                fprintf('\n%s%g\n',[ExitText ', Objective '],  c'*solution.full*osense);
            end
        end
    end
end

%return basis
if basisReuse
    LPProblem.LPBasis=basis;
end

if sum(minNorm)~=0
    fprintf('%s\n','This objective corresponds to a flux with minimum Euclidean norm.');
    fprintf('%s%d%s\n','The weighting for minimising the norm was ',minNorm,'.');
    fprintf('%s\n','Check that the objective is the same without minimising the norm.');
end

author: opencobra

repo: cobratoolbox


text:
%kinfo 'Info'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros info.pane file
%
% Parameters:
%
% Example: kinfo( {})
%
% Khoros helpfile follows below:

function varargout = kinfo(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kinfo(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={};
maxval={};
minval={};
istoggle=[];
was_set=istoggle * 0;
paramtype={};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'info" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%Sensor Refs=[use (0,1), val, obs noise std, ini std]
%u is assumed to be a random constant (with infinite initial variance)
function [u, ximu, Pimu, Pu, Pux]=selfcalib1(imu_data, Aimu, Qimu, Cimu, Rimu, PimuInv, M, sensor_ref)

DBG=0;  %Debug

nu=size(M,2); %# of kinematic variables
ns=size(M,1); %# of sensors
nimu=size(Aimu,1);    %total # of imu states

if (isempty(sensor_ref)) %not a compulsory argument
    sensor_ref=zeros(size(M,2),4);
end

%System model (augment with kv)
A=diagmat_v000(Aimu, eye(nu),[]);
Q=diagmat_v000(Qimu, zeros(nu), []);
C=[M Cimu];
R=Rimu;
Pinv=diagmat_v000(PimuInv, zeros(nu), []);

%External body acc/rot rate observations
next=0;
Cext=[];
Rext=[];
yext=[];
for in=1:nu
    if (sensor_ref(in,1)==1)     %there is an external obs
        next=next+1;
        Cext=[Cext;zeros(1,size(A,1))];
        Cext(end,in)=1;
        Rext=diagmat_v000(sensor_ref(in,3)^2,Rext,[]);
        yext=[yext;sensor_ref(in,2)];
        if (sensor_ref(in,4)~=0)    %external obs have random bias
            A=diagmat_v000(1, A, []);
            Q=diagmat_v000(0, Q, []);
            Pinv=diagmat_v000(1/sensor_ref(in,4)^2, Pinv, []);
            C=[C zeros(ns,1)];
            Cext=[Cext zeros(next,1)];
            Cext(end,end)=1;
        end
    end
end
nst=size(A,1);

%%%%% Start the estimation routine
%%Sensor outputs
%%For the first sensor data use the innovation form (Pinv is probably singular)
y=imu_data(:,1);
Rinv=inv(R);
P=inv(Pinv+C'*Rinv*C);
xest=P*(C'*Rinv*y);
%%External obs
if (next>0)
    K=P*Cext'*inv(Cext*P*Cext'+Rext);
    xest=xest+K*(yext-Cext*xest);
    P=(eye(nst)-K*Cext)*P*(eye(nst)-K*Cext)'+K*Rext*K';
end

if (DBG)
    ndat=size(imu_data,2);
    debmn=zeros(nu,ndat);
    debop=zeros(size(A,1),ndat);

    WLS=lsmat_v001(M,R,0);
    mean_dat=WLS*y;
    debmn(:,1)=mean_dat;
    debop(:,1)=xest;
end

%For the rest use standard KF
for in=2:size(imu_data,2)
    %%%prediction
    xest=A*xest;
    P=A*P*A'+Q;

    %%%filter
    %%Sensor data
    y=imu_data(:,in);
    K=P*C'*inv(C*P*C'+R);
    %K(1:nimu,:)=0;
    xest=xest+K*(y-C*xest);
    P=(eye(nst)-K*C)*P*(eye(nst)-K*C)'+K*R*K';
    %%External obs
    if (next>0)
        K=P*Cext'*inv(Cext*P*Cext'+Rext);
        xest=xest+K*(yext-Cext*xest);
        P=(eye(nst)-K*Cext)*P*(eye(nst)-K*Cext)'+K*Rext*K';
    end

    if (DBG)
        mean_dat=((in-1)*mean_dat+(WLS*y))/in;
        debmn(:,in)=mean_dat;
        debop(:,in)=xest;
    end
end

%Outputs (combine deterministic and stochastic results)
ximu=xest(nu+1:nu+nimu,1);
Pimu=P(nu+1:nu+nimu,nu+1:nu+nimu);
u=xest(1:nu,:);
Pu=P(1:nu,1:nu);
Pux=P(1:nu,nu+1:nu+nimu);

if (DBG)
    in=1;
    figure;
    plot(imu_data(in,:))
    hold on;
    plot(debmn(in,:),'r')
    plot(debop(in,:),'g')
end

repo: ekfmonoslam


text:
function a = r8ri_to_r8ge ( nz, ija, sa, n )

%*****************************************************************************80
%
%% R8RI_TO_R8GE converts an R8RI matrix to R8GE form.
%
%  Discussion:
%
%    An R8RI matrix is in row indexed sparse storage form.
%
%    A R8GE matrix is in general storage.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
%    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
%    Third Edition,
%    Cambridge University Press, 2007,
%    ISBN13: 978-0-521-88068-8,
%    LC: QA297.N866.
%
%  Parameters:
%
%    Input, integer NZ, the size required for the RI
%    or "row indexed" sparse storage.
%
%    Input, integer IJA(NZ), the index vector.
%
%    Input, real SA(NZ), the value vector.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the matrix stored in GE
%    or "general" format.
%
  a = zeros ( n, n );

  for k = 1 : n
    i = k;
    j = k;
    a(i,j) = sa(k);
  end

  for i = 1 : n
    for k = ija(i) : ija(i+1) - 1
      j = ija(k);
      a(i,j) = sa(k);
    end
  end

  return
end

repo: jburkardt-m