Skip to content

Commit 6a78b73

Browse files
toladnertoladner
committed
CORA v2026.0.1
- contDynamics/conform: added outlier detection - contDynamics/reach: bug fix - contSet/plot: bug fix - global/classes/confPredictor: added, makes set-based predictions with coverage guarantees
1 parent 2415486 commit 6a78b73

File tree

56 files changed

+5513
-134
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

56 files changed

+5513
-134
lines changed

‎.gitignore‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -50,4 +50,7 @@ manual/manual.iml
5050
**/.DS_Store
5151

5252
# exclude folders in workspace
53-
workspace/*
53+
workspace/*
54+
55+
# CORA release
56+
global/functions/release/commitHistory.txt

‎contDynamics/@contDynamics/conform.m‎

Lines changed: 26 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,8 @@
2323
% .sys - updated system
2424
% .p - estimated parameters
2525
% .fval - final optimization cost
26+
% .idzActive - indizes of active / boundary data points (only
27+
% determined if options.cs.determineActive = true, otherwise NaN)
2628
%
2729
% References:
2830
% [1] M. Althoff and J. M. Dolan. Reachability computation of low-order
@@ -80,12 +82,31 @@
8082
% conformance synthesis for dedicated system dynamics
8183
[params,options] = validateOptions(sys,params,options);
8284
[sys_upd,params,options] = aux_updateSys(sys,params,options);
83-
85+
8486
% Identify conformant parameters
8587
if type == "white" || type == "black"
86-
[params,fval,p_opt,union_y_a] = priv_conform_white(sys_upd,params,options);
87-
results.sys = sys;
88+
if options.cs.numOutlier >= 1 && ...
89+
(params.testSuite(1).n_k > 1 || params.testSuite(1).n_s > 1)
90+
throw(CORAerror('CORA:specialError', ['Outlier detection ' ...
91+
'not implemented for n_k > 1 or n_s > 1. ' ...
92+
' Please use options.cs.numOutlier = 0.']))
93+
end
94+
if options.cs.numOutlier >= 1 && any(strcmp(options.cs.outMethod, {'search','searchG'}))
95+
% outlier detection via iterative search
96+
[params,fval,p_opt,union_y_a,idzActive] = priv_conform_iterOD(sys_upd, params, options);
97+
else
98+
if options.cs.numOutlier >= 1 && strcmp(options.cs.outMethod, 'RMSE')
99+
% remove data points with biggest RMSE
100+
deviations = computeOutputDev(params.testSuite,sys_upd);
101+
dev_vec = squeeze(rms(deviations,1));
102+
[~,idzOutlier] = sort(dev_vec,'descend');
103+
options.cs.idzOutlier = idzOutlier(1:options.cs.numOutlier);
104+
end
105+
[params,fval,p_opt,union_y_a,idzActive] = priv_conform_white(sys_upd,params,options);
106+
end
107+
results.sys = sys_upd;
88108
results.unifiedOutputs = union_y_a;
109+
results.idzActive = idzActive;
89110

90111
elseif contains(type,"gray") % "graySim","graySeq","grayLS"
91112
[params,fval,p_opt,sys_opt] = priv_conform_gray(sys_upd,params,options,type);
@@ -119,7 +140,7 @@
119140

120141
% adapt testSuite to the system dynamics (augment nominal input, update
121142
% initial state)
122-
params.testSuite = aux_preprocessTestSuite(sys,params.testSuite,options.cs.recMethod);
143+
params.testSuite = aux_preprocessTestSuite(sys,params.testSuite);
123144

124145
end
125146

@@ -141,7 +162,7 @@
141162

142163
end
143164

144-
function testSuite = aux_preprocessTestSuite(sys,testSuite,recMethod)
165+
function testSuite = aux_preprocessTestSuite(sys,testSuite)
145166
% preprocessTestSuite - preprocess testSuite, i.e.,
146167
% - augment nominal inputs with zeros to match input dimension of sys
147168
% - unify test cases for linear systems,

‎contDynamics/@contDynamics/private/priv_conform_iterOD.m‎

Lines changed: 180 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
function [params_out,fval,p_opt,union_y_a,idzActive] = priv_conform_iterOD(sys, params, options)
2+
% priv_conform_iterOD - conformance synthesis with iterative outlier
3+
% detection based on [1]
4+
%
5+
% Syntax:
6+
% [params_out, results] = priv_conform_iterOD(pred, params, options)
7+
%
8+
% Inputs:
9+
% pred - confPredictor object
10+
% params - Initial test parameters, including the test suite
11+
% options - options for conformance synthesis
12+
%
13+
% Outputs:
14+
% params_out - parameters solving the conformance problem
15+
% fval - conformance cost
16+
% p_opt - estimated parameters containing the scaling factors alpha
17+
% and the change of the center vectors Delta c with
18+
% p_opt = [alpha_x' alpha_u' dc_x' dc_u']'
19+
% union_y_a - n_m x 1 cell array with dim_y x n_k x n_s output arrays
20+
% with union_y_a{m} = testSuite(m).y - y_nom
21+
% m_active - indizes of active data points at the final solution
22+
%
23+
% References:
24+
% [1] Campi et al. 2009, "Interval predictor models: Identification and
25+
% reliability".
26+
% [2] Laura Luetzow, Michael Eichelbeck, Mykel Kochenderfer, and
27+
% Matthias Althoff. "Zono-Conformal Prediction: Zonotope-Based
28+
% Uncertainty Quantification for Regression and Classification
29+
% Tasks," arXiv, 2025.
30+
%
31+
% Other m-files required: none
32+
% Subfunctions: none
33+
% MAT-files required: none
34+
%
35+
% See also: confPredictor
36+
37+
% Authors: Michael Eichelbeck, Laura Luetzow
38+
% Written: 02-October-2025
39+
% Last update: ---
40+
% Last revision: ---
41+
42+
% ------------------------------ BEGIN CODE -------------------------------
43+
44+
if strcmp(options.cs.outMethod,'searchG')
45+
greedy = true;
46+
else
47+
greedy = false; % Default to non-greedy approach if not specified
48+
end
49+
n_out = options.cs.numOutlier;
50+
51+
num_nodes = zeros(n_out+1, 1); % Stores the number of nodes at each level of the search tree
52+
53+
% Data structures to store intermediate results during constraint removal
54+
idxConstr = cell(n_out+1, 1); % Stores reduced test suites at each level of the search tree
55+
fval_all = zeros(n_out+1, 1); % Tracks the cost function values at each level
56+
fval_best = zeros(n_out+1, 1); % Tracks the cost function values at each level
57+
params_best = cell(n_out+1,1);
58+
T_all = zeros(n_out+1, 1);
59+
resOut = cell(n_out, 1); % Stores results for each node in the search tree
60+
61+
idxConstr{1,1} = 1:length(params.testSuite);
62+
m_constr = idxConstr{1,1}; % Initialize root node using all test cases
63+
options.cs.idzOutlier = [];
64+
num_nodes(1) = 1; % Root level has a single node
65+
66+
% Initial run with all constraints active to determine baseline performance
67+
options.cs.determineActive = true;
68+
options.cs.numOutlier = length(options.cs.idzOutlier);
69+
Timer = tic;
70+
[params_out,fval,p_opt,union_y_a,idzActive] = priv_conform_white(sys,params,options);
71+
results.sys = sys;
72+
results.unifiedOutputs = union_y_a;
73+
results.idzActive = idzActive;
74+
results.fval = fval;
75+
results.p = p_opt;
76+
T_all(1) = toc(Timer);
77+
resOut{1,1} = struct('params', params_out, 'results', results);
78+
fval_all(1,1) = results.fval; % Store initial cost function value
79+
fval_best(1,1) = results.fval;
80+
params_best{1} = params_out;
81+
correctPred = false;
82+
83+
if fval == 0
84+
% all predictions are already correct, removal of additional data
85+
% points does not imrpove the cost
86+
return
87+
end
88+
89+
% Begin search tree traversal to remove constraints
90+
for i=1:n_out % Iterate over each level of the tree (corresponding to outlier removals)
91+
num_nodes(i+1) = 0; % Initialize node count for the next level
92+
best_fval = Inf; % Track best fval for greedy selection
93+
best_idx_constr = [];
94+
best_results = [];
95+
if correctPred
96+
% copy results of previous level if all predictions are already correct
97+
num_nodes(i+1) = num_nodes(i+1) + 1; % Increment node count for next level
98+
idxConstr{i+1, num_nodes(i+1)} = m_constr; % Store the updated test suite
99+
fval_all(i+1, num_nodes(i+1)) = fval_l; % Store the updated cost function value
100+
resOut{i+1,num_nodes(i+1)} = struct('params', paramsRes_l, 'results', results_l);
101+
continue
102+
end
103+
104+
for h=1:num_nodes(i) % Iterate over each node at the current level
105+
if correctPred
106+
% cancel loop if all predictions are correct
107+
break
108+
end
109+
active_constraints = resOut{i, h}.results.idzActive; % Identify active constraints
110+
for l=1:length(active_constraints) % Iterate through active constraints
111+
c_index = active_constraints(l);
112+
% Remove the selected constraint from the parent node's test suite
113+
m_constr = idxConstr{i,h};
114+
m_constr(c_index) = [];
115+
116+
% Evaluate the modified parameter set
117+
options.cs.idzOutlier = 1:length(params.testSuite);
118+
options.cs.idzOutlier(m_constr) = [];
119+
options.cs.numOutlier = length(options.cs.idzOutlier);
120+
121+
[paramsRes_l,fval_l,p_opt,union_y_a,idzActive] = priv_conform_white(sys,params,options);
122+
results_l.sys = sys;
123+
results_l.unifiedOutputs = union_y_a;
124+
results_l.idzActive = idzActive;
125+
results_l.fval = fval_l;
126+
results_l.p = p_opt;
127+
128+
if greedy % Greedy approach: only keep track of the best option at each level
129+
if fval_l <= best_fval
130+
best_fval = fval_l;
131+
best_idx_constr = m_constr;
132+
best_results = results_l;
133+
best_params = paramsRes_l;
134+
135+
end
136+
else % Non-greedy approach: explore all possible nodes
137+
if (fval_l <= fval_all(i,h)) % If the new result improves the cost function
138+
% Consider the removed constraint as a support constraint and create a new node
139+
num_nodes(i+1) = num_nodes(i+1) + 1; % Increment node count for next level
140+
idxConstr{i+1, num_nodes(i+1)} = m_constr; % Store the updated test suite
141+
fval_all(i+1, num_nodes(i+1)) = fval_l; % Store the updated cost function value
142+
resOut{i+1,num_nodes(i+1)} = struct('params', paramsRes_l, 'results', results_l);
143+
end
144+
end
145+
if fval_l == 0
146+
correctPred = true;
147+
% no uncertainty in the point predictions
148+
break
149+
end
150+
end
151+
end
152+
153+
if greedy % If using greedy search, only keep the best node and discard others
154+
if ~isempty(best_idx_constr)
155+
num_nodes(i+1) = 1;
156+
idxConstr{i+1,1} = best_idx_constr;
157+
fval_all(i+1,1) = best_fval;
158+
resOut{i+1,1} = struct('params', best_params, 'results', best_results);
159+
end
160+
end
161+
T_all(i+1) = toc(Timer);
162+
[fval_best(i+1), idx] = min(fval_all(i+1,:),[],2);
163+
params_best{i+1} = resOut{i+1,idx}.params;
164+
end
165+
166+
% Identify the best result in the last level of the search tree
167+
for j=1:num_nodes(n_out+1)
168+
if resOut{n_out+1,j}.results.fval < results.fval % Check if a better result is found
169+
results = resOut{n_out+1,j}.results;
170+
params_out = resOut{n_out+1,j}.params;
171+
end
172+
end
173+
fval = results.fval;
174+
p_opt = results.p;
175+
union_y_a = results.unifiedOutputs;
176+
idzActive = results.idzActive;
177+
178+
end
179+
180+
% ------------------------------ END OF CODE ------------------------------

0 commit comments

Comments
 (0)