From 8f6ba2b6719900574d5734c548e190fb7c8b6e81 Mon Sep 17 00:00:00 2001 From: Anton Edvinovich Pozharskiy Date: Wed, 2 Apr 2025 11:44:44 +0200 Subject: [PATCH 1/2] Add logic to deduplicate complementarity vector efficiently --- src/+mpecopt/Solver.m | 40 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/src/+mpecopt/Solver.m b/src/+mpecopt/Solver.m index 9059274..22d89bb 100644 --- a/src/+mpecopt/Solver.m +++ b/src/+mpecopt/Solver.m @@ -1077,10 +1077,26 @@ function create_mpec_functions(obj) % TODO(@anton) do we actually need this here or can we calculate these "analytically" ind_x1 = []; ind_x2 = []; - ind_x1_fun = Function('ind_1',{x},{x.jacobian(x1)}); + % HERE BE DRAGONS: This is using some casadi internals to deduplicate x1, x2, expressions. + % Get element cell arrays. + x1_elem = x1.elements(); + x2_elem = x2.elements(); + % Get element hashes. + x1_hashes = cellfun(@(x) x.element_hash(), x1_elem); + x2_hashes = cellfun(@(x) x.element_hash(), x2_elem); + % Get indicies of unique elements by comparing the element hashes. + % WARNING: Technically these hashes may be nonunique but looking at how they are generated this should + % never practically be the case. If it happens though the symptoms should be obvious. + [~,x1_uniq,x1_uniq_backmap] = unique(x1_hashes, 'stable'); + [~,x2_uniq,x2_uniq_backmap] = unique(x2_hashes, 'stable'); + % Get indicies with jacobian trick + ind_x1_fun = Function('ind_1',{x},{x.jacobian(x1(x1_uniq))}); + ind_x2_fun = Function('ind_2',{x},{x.jacobian(x2(x2_uniq))}); [ind_x1,~] = find(sparse(ind_x1_fun(zeros(size(x))))==1); - ind_x2_fun = Function('ind_2',{x},{x.jacobian(x2)}); [ind_x2,~] = find(sparse(ind_x2_fun(zeros(size(x))))==1); + % Resize indices back into correct size by mapping back the nonunique elements + ind_x1 = ind_x1(x1_uniq_backmap); + ind_x2 = ind_x2(x2_uniq_backmap); opts.nlp_is_mpec = 1; % TODO: check is this still used? (its purpose: if not an mpec, just make single nlp call without mpec machinery); else opts.nlp_is_mpec = 0; @@ -1205,10 +1221,26 @@ function create_mpec_functions(obj) % TODO(@anton) do we actually need this here or can we calculate these "analytically" ind_x1 = []; ind_x2 = []; - ind_x1_fun = Function('ind_1',{x},{x.jacobian(x1)}); + % HERE BE DRAGONS: This is using some casadi internals to deduplicate x1, x2, expressions. + % Get element cell arrays. + x1_elem = x1.elements(); + x2_elem = x2.elements(); + % Get element hashes. + x1_hashes = cellfun(@(x) x.element_hash(), x1_elem); + x2_hashes = cellfun(@(x) x.element_hash(), x2_elem); + % Get indicies of unique elements by comparing the element hashes. + % WARNING: Technically these hashes may be nonunique but looking at how they are generated this should + % never practically be the case. If it happens though the symptoms should be obvious. + [~,x1_uniq,x1_uniq_backmap] = unique(x1_hashes, 'stable'); + [~,x2_uniq,x2_uniq_backmap] = unique(x2_hashes, 'stable'); + % Get indicies with jacobian trick + ind_x1_fun = Function('ind_1',{x},{x.jacobian(x1(x1_uniq))}); + ind_x2_fun = Function('ind_2',{x},{x.jacobian(x2(x2_uniq))}); [ind_x1,~] = find(sparse(ind_x1_fun(zeros(size(x))))==1); - ind_x2_fun = Function('ind_2',{x},{x.jacobian(x2)}); [ind_x2,~] = find(sparse(ind_x2_fun(zeros(size(x))))==1); + % Resize indices back into correct size by mapping back the nonunique elements + ind_x1 = ind_x1(x1_uniq_backmap); + ind_x2 = ind_x2(x2_uniq_backmap); opts.nlp_is_mpec = 1; % TODO: check is this still used? (its purpose: if not an mpec, just make single nlp call without mpec machinery); else opts.nlp_is_mpec = 0; From d17127b974cb438915ef6bbedb1c9c3c5c5c9d07 Mon Sep 17 00:00:00 2001 From: Anton Edvinovich Pozharskiy Date: Wed, 2 Apr 2025 12:38:52 +0200 Subject: [PATCH 2/2] Also fix mathworks breaking ind2sub backwards compatibility for seemingly no reason --- src/+mpecopt/Solver.m | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/+mpecopt/Solver.m b/src/+mpecopt/Solver.m index 22d89bb..e2dd49e 100644 --- a/src/+mpecopt/Solver.m +++ b/src/+mpecopt/Solver.m @@ -1901,7 +1901,11 @@ function create_lpec_functions(obj) X = zeros(prod(sizeThisSet),numSets); for i=1:size(X,1) ixVect = cell(length(sizeThisSet),1); - [ixVect{:}] = ind2sub(flip(sizeThisSet),i); + sz = flip(sizeThisSet); + if length(sz) == 1 + sz = [sz,1]; + end + [ixVect{:}] = ind2sub(sz,i); ixVect = flip([ixVect{:}]); vect = zeros(1, numSets); for jj=1:numSets