Skip to content

Commit

Permalink
Merge pull request #9 from afvillaverde/dev
Browse files Browse the repository at this point in the history
Dev
  • Loading branch information
afvillaverde authored Jun 20, 2022
2 parents 788af8e + bd8d56b commit 9c67436
Show file tree
Hide file tree
Showing 35 changed files with 1,326 additions and 202 deletions.
52 changes: 30 additions & 22 deletions STRIKE-GOLDD/STRIKE_GOLDD.m
Original file line number Diff line number Diff line change
@@ -1,24 +1,22 @@
%==========================================================================
% StrIkE-GOLDD (Structural Identifiability taken as Extended-Generalized
% Observability with Lie Derivatives and Decomposition)
% Observability with Lie Derivatives and Decomposition):
% a Matlab toolbox for structural identifiability and observability (SIO)
% analysis of nonlinear models.
%--------------------------------------------------------------------------
%
% A Matlab toolbox for structural identifiability and observability
% analysis of nonlinear models
%--------------------------------------------------------------------------
%
% User input can be specified in file 'options.m'
% Documentation available in the 'doc' folder
% User input can be specified in file 'options.m'.
% Documentation available in the 'doc' folder.
%--------------------------------------------------------------------------
%
% Version 3.0, last modified: 19/10/2020
% Alejandro Fernandez Villaverde (afvillaverde@iim.csic.es)
% Version 4.0, last modified: 25/04/2022
% Contact: Alejandro F. Villaverde (afvillaverde@uvigo.gal)
%==========================================================================

function STRIKE_GOLDD(varargin)

fprintf('\n\n -------------------------------- \n');
fprintf(' >>> STRIKE-GOLDD toolbox 3.0 \n');
fprintf(' >>> STRIKE-GOLDD toolbox 4.0 \n');
fprintf(' -------------------------------- \n');

%==========================================================================
Expand All @@ -36,6 +34,7 @@ function STRIKE_GOLDD(varargin)
delete("current_options.m");
nmf = pwd;
case 2
% when STRIKE_GOLDD.m is called by AutoRepar, add extra paths:
[~,paths,opts,submodels,prev_ident_pars] = options;
modelname = varargin{1};
mf = pwd;
Expand All @@ -53,13 +52,28 @@ function STRIKE_GOLDD(varargin)
addpath(genpath(paths.functions));

%==========================================================================
% ORC-DF algorithm (only for affine-in-inputs systems):
if opts.affine==1
ORC_DF(modelname,opts,prev_ident_pars);
return
end % If not, run the FISPO algorithm:
% Read the algorithm chosen in the options file:
switch opts.algorithm
case 1
% Run the FISPO algorithm (remainder of this file)
case 2
prob_obs_test(modelname,opts,prev_ident_pars,nmf);
return
case 3
ORC_DF(modelname,opts,prev_ident_pars,nmf);
return
case 4
Lie_Symmetry
return
case 5
if nargin < 2 % (if nargin == 2, STRIKE-GOLDD is being called by AutoRepar => avoid recursive loop)
AutoRepar
return
end
end

%==========================================================================
% If this line is executed, it means that the FISPO algorithm was chosen.
% Load model:

% Convert model to multi-experiment form (Optional):
Expand All @@ -70,7 +84,7 @@ function STRIKE_GOLDD(varargin)

load(modelname);

fprintf('\n Analyzing the %s model... \n', modelname);
fprintf('\n Analyzing the %s model with the FISPO algorithm...\n',modelname);

tic
%==========================================================================
Expand Down Expand Up @@ -369,10 +383,6 @@ function STRIKE_GOLDD(varargin)
fprintf('\n However, the maximum computation time allowed for calculating each of them has been reached.');
fprintf('\n You can increase it by changing <<opts.maxLietime>> (currently opts.maxLietime = %d)',opts.maxLietime);
unidflag = 0;
if opts.unidentif == 1
skip_elim = 1;
increaseLie = 0;
end
end
end
end
Expand All @@ -394,7 +404,6 @@ function STRIKE_GOLDD(varargin)
end
end
else% If the maxLietime has been reached, but the minimum of Lie derivatives has not been calculated:
if opts.unidentif == 0
if opts.decomp == 0
fprintf('\n More Lie derivatives would be needed to analyse the model.');
fprintf('\n However, the maximum computation time allowed for calculating each of them has been reached.');
Expand All @@ -420,7 +429,6 @@ function STRIKE_GOLDD(varargin)
else
decomp_flag = 1;
end
end
end
end % end of 'if opts.forcedecomp == 0'

Expand Down
Binary file modified STRIKE-GOLDD/doc/STRIKE-GOLDD_manual.pdf
Binary file not shown.
1 change: 0 additions & 1 deletion STRIKE-GOLDD/functions/AutoRepar.m
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@
addpath(resultspath);



%% STRIKE-GOLDD
[modelname,paths,opts,submodels,prev_ident_pars] = options();
if (opts.use_existing_results==0)
Expand Down
15 changes: 10 additions & 5 deletions STRIKE-GOLDD/functions/ORC_DF.m
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
%====ORC-DF: An implementation of the Observability Rank Criterion for====%
%=====================systems with Direct Feedthrough=====================%

function ORC_DF(modelname,opts,prev_ident_pars)
function ORC_DF(modelname,opts,prev_ident_pars,nmf)

tStart=tic;

Expand All @@ -25,7 +25,7 @@ function ORC_DF(modelname,opts,prev_ident_pars)

load(modelname) %#ok<*LOAD>

fprintf('\n >>> Analyzing observability of %s with affine in control system algorithm\n',modelname);
fprintf('\n >>> Analyzing the %s model with the ORC-DF algorithm...\n',modelname);

%======================Dimensions of the problem==========================%

Expand Down Expand Up @@ -829,12 +829,17 @@ function ORC_DF(modelname,opts,prev_ident_pars)
fprintf('\n\n Total execution time: %d \n\n',totaltime);

% Save results:
resultsname = sprintf('id_results_%s',modelname);
fullresultsname = strcat(pwd,filesep,'results',filesep,resultsname,'_',date);

resultsname = sprintf('id_results_%s_%s',modelname,date);
fullresultsname = strcat(nmf,filesep,'results',filesep,resultsname);
warning off 'parallel:cluster:CannotSaveCorrectly'
save(fullresultsname);

% Delete affine model (optional):
if opts.affine_delete_model, delete(affine_model_file);end
if opts.affine_delete_model
if isfile(affine_model_file)
delete(affine_model_file);
end
end

end
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
function [exp_e,exp_p] = expand_lie(pMax,p)
%% Expand, as a Taylor series, the paramter to be removed and exp(epsilon)
%% Expand, as a Taylor series, the parameter to be removed and exp(epsilon)
% Input: pMax --> Maximum degree of Lie series
% p --> Parameter
% Output: exp_e --> Vector with coefficients of Taylor series of exp(epsilon)
Expand Down
Original file line number Diff line number Diff line change
@@ -1,60 +1,60 @@
function [sim_vec] = simp(allVar,vec,eps_var)
%% FUNCTION SIMPLIFY
%
% Function to find some expressions that may have not been simplified
% during the computation of new variables
%
% INPUT: allVar --> All variables
% vec --> vector with final tranformations
% eps_var --> value of exponential(epsilon)
%
% OUTPUT:
% sim_vec --> simplified final transformations vector
%%
syms epsilon
ch=children(vec);
[~,ncol]=size(ch);
if iscell(ch)==0 % For R2020b
ch=num2cell(ch);
end
for i=1:ncol
[~,ncol_2]=size(ch{1,i});
if iscell(ch{1,i})==1 %For R2020b
aux_v=[ch{1,i}{:}];
else
aux_v=[ch{1,i}(:)];
end
for j=1:ncol_2
[n,d]=numden(aux_v(j));
if n==1 %% Only denominator
[c,p]=coeffs(d,allVar);
if (numel(unique(abs(p)))==1 && abs(c)~=1) % Some expression contains exp(ep)
eq=simplify(log(c),'IgnoreAnalyticConstraints',true);
[cc,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc*epsilon),eps_var^(cc));
end
elseif d==1
[c,p]=coeffs(n,allVar);
if (numel(unique(abs(p)))==1 && abs(c)~=1) % Some expression contains exp(ep)
eq=simplify(log(c),'IgnoreAnalyticConstraints',true);
[cc,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc*epsilon),eps_var^(cc));
end
else
[c_n,p_n]=coeffs(n,allVar);
[c_d,p_d]=coeffs(d,allVar);
if numel(unique(abs(p_n)))==1 && abs(c_n)~=1 % Some expression contains exp(ep)
eq=simplify(log(c_n),'IgnoreAnalyticConstraints',true);
[cc_n,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc_n*epsilon),eps_var^(cc_n));
elseif numel(unique(abs(p_d)))==1 && abs(c_d)~=1
eq=simplify(log(c_d),'IgnoreAnalyticConstraints',true);
[cc_d,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc_d*epsilon),eps_var^(cc_d));
end
end
end
end
sim_vec=vec;
end

function [sim_vec] = simp(allVar,vec,eps_var)
%% FUNCTION SIMPLIFY
%
% Function to find some expressions that may have not been simplified
% during the computation of new variables
%
% INPUT: allVar --> All variables
% vec --> vector with final tranformations
% eps_var --> value of exponential(epsilon)
%
% OUTPUT:
% sim_vec --> simplified final transformations vector
%%
syms epsilon
ch=children(vec);
[~,ncol]=size(ch);
if iscell(ch)==0 % For R2020b
ch=num2cell(ch);
end
for i=1:ncol
[~,ncol_2]=size(ch{1,i});
if iscell(ch{1,i})==1 %For R2020b
aux_v=[ch{1,i}{:}];
else
aux_v=[ch{1,i}(:)];
end
for j=1:ncol_2
[n,d]=numden(aux_v(j));
if n==1 %% Only denominator
[c,p]=coeffs(d,allVar);
if (numel(unique(abs(p)))==1 && abs(c)~=1) % Some expression contains exp(ep)
eq=simplify(log(c),'IgnoreAnalyticConstraints',true);
[cc,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc*epsilon),eps_var^(cc));
end
elseif d==1
[c,p]=coeffs(n,allVar);
if (numel(unique(abs(p)))==1 && abs(c)~=1) % Some expression contains exp(ep)
eq=simplify(log(c),'IgnoreAnalyticConstraints',true);
[cc,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc*epsilon),eps_var^(cc));
end
else
[c_n,p_n]=coeffs(n,allVar);
[c_d,p_d]=coeffs(d,allVar);
if numel(unique(abs(p_n)))==1 && abs(c_n)~=1 % Some expression contains exp(ep)
eq=simplify(log(c_n),'IgnoreAnalyticConstraints',true);
[cc_n,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc_n*epsilon),eps_var^(cc_n));
elseif numel(unique(abs(p_d)))==1 && abs(c_d)~=1
eq=simplify(log(c_d),'IgnoreAnalyticConstraints',true);
[cc_d,~]=coeffs(eq,epsilon);
vec(i)=subs(vec(i),exp(cc_d*epsilon),eps_var^(cc_d));
end
end
end
end
sim_vec=vec;
end

File renamed without changes.
File renamed without changes.
41 changes: 41 additions & 0 deletions STRIKE-GOLDD/functions/aux_Prob_Obs_Test/ExpMatrixSeries.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
% -------------------------------------------------------------------------
% computation using a newton operator
%--------------------------------------------------------------------------

function [sol1,sol2]=ExpMatrixSeries(Mat,LocalOrder,Matsize,MyPrime)

if LocalOrder==1
sol1=zeros(1,Matsize*Matsize);
sol1((1:Matsize:Matsize*Matsize)+(0:Matsize-1))=1;
sol2=zeros(1,Matsize*Matsize);
sol2((1:Matsize:Matsize*Matsize)+(0:Matsize-1))=1;
else
NewOrder=2^ceil(log2(LocalOrder)-1);
aux=trunpoly(Mat,LocalOrder);
[aux1,aux2]=ExpMatrixSeries(aux,NewOrder,Matsize,MyPrime);

aux1=trunpoly(aux1,LocalOrder);
aux2=trunpoly(aux2,LocalOrder);

sol2=multmatpolytrun(-aux1,aux2,NewOrder,Matsize);
sol2=mod(sol2,MyPrime);
sol2(end,(1:Matsize:Matsize*Matsize)+(0:Matsize-1))= ...
sol2(end,(1:Matsize:Matsize*Matsize)+(0:Matsize-1))+2;
sol2=multmatpolytrun(aux2,sol2,NewOrder,Matsize);
sol2=mod(sol2,MyPrime);


Matder=dermatpoly(Mat,LocalOrder);
aux1der=dermatpoly(aux1,LocalOrder);
sol1=mod(trunpoly( ...
multmatpolytrun(Matder,aux1,LocalOrder,Matsize)-aux1der ...
,LocalOrder),MyPrime);
sol1=mod(multmatpolytrun(sol2,sol1,LocalOrder,Matsize),MyPrime);
sol1=intmatpoly(sol1,LocalOrder);
sol1=ratmod(sol1,MyPrime);
sol1(end,(1:Matsize:Matsize*Matsize)+(0:Matsize-1))= ...
sol1(end,(1:Matsize:Matsize*Matsize)+(0:Matsize-1))+1;
sol1=mod(multmatpolytrun(aux1,sol1,LocalOrder,Matsize),MyPrime);


end
27 changes: 27 additions & 0 deletions STRIKE-GOLDD/functions/aux_Prob_Obs_Test/InverseMatrixSeries.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
%--------------------------------------------------------------------------
% use of newton operator for computation of InvHomSolInvA
%--------------------------------------------------------------------------

function sol=InverseMatrixSeries(Mat,LocalOrder,MatSize,MyPrime)

if LocalOrder==1
sol=Mat(end,:);
aux=zeros(MatSize,MatSize);
for i=1:MatSize
aux(i,:)=sol((i-1)*MatSize+1:i*MatSize);
end
aux=inv(aux);
for i=1:MatSize
sol((i-1)*MatSize+1:i*MatSize)=aux(i,:);
end
else
NewOrder = 2^ceil(log2(LocalOrder)-1);
aux=trunpoly(Mat,NewOrder);
aux=InverseMatrixSeries(aux,NewOrder,MatSize,MyPrime);
aux=trunpoly(aux,LocalOrder);

sol=mod(multmatpolytrun(-aux,Mat,LocalOrder,MatSize),MyPrime);
sol(end,(1:MatSize:MatSize*MatSize)+(0:MatSize-1))= ...
sol(end,(1:MatSize:MatSize*MatSize)+(0:MatSize-1))+2;
sol=mod(multmatpolytrun(sol,aux,LocalOrder,MatSize),MyPrime);
end
Loading

0 comments on commit 9c67436

Please sign in to comment.