From 66bfd22b217468004d49dbe3486f163dd89cc7c1 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Wed, 19 Mar 2025 12:00:17 +0100 Subject: [PATCH] fix find nonscalar when have params --- src/+mpecopt/Solver.m | 12 ++++++------ src/find_nonscalar.m | 14 +++++++++----- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/+mpecopt/Solver.m b/src/+mpecopt/Solver.m index 72f41d1..4b09e4c 100644 --- a/src/+mpecopt/Solver.m +++ b/src/+mpecopt/Solver.m @@ -955,8 +955,8 @@ function create_mpec_functions(obj) G_curr = mpec.G.(name)(); H_curr = mpec.H.(name)(); - [ind_scalar_G,ind_nonscalar_G, ind_map_G] = find_nonscalar(G_curr, mpec.w.sym); - [ind_scalar_H,ind_nonscalar_H, ind_map_H] = find_nonscalar(H_curr, mpec.w.sym); + [ind_scalar_G,ind_nonscalar_G, ind_map_G] = find_nonscalar(G_curr, mpec.w.sym, mpec.p.sym); + [ind_scalar_H,ind_nonscalar_H, ind_map_H] = find_nonscalar(H_curr, mpec.w.sym, mpec.p.sym); mpec.w.([name '_G_lift']) = {{'G', length(ind_nonscalar_G)}, 0, inf}; G_lift = G_curr(ind_nonscalar_G); @@ -1000,8 +1000,8 @@ function create_mpec_functions(obj) % TODO(@anton) this is in an if for performance reasons not to call find_nonscalar many times % however it is likely that this can be done in 1 shot? - [ind_scalar_G,ind_nonscalar_G, ind_map_G] = find_nonscalar(G_curr, mpec.w.sym); - [ind_scalar_H,ind_nonscalar_H, ind_map_H] = find_nonscalar(H_curr, mpec.w.sym); + [ind_scalar_G,ind_nonscalar_G, ind_map_G] = find_nonscalar(G_curr, mpec.w.sym, mpec.p.sym); + [ind_scalar_H,ind_nonscalar_H, ind_map_H] = find_nonscalar(H_curr, mpec.w.sym, mpec.p.sym); mpec.w.([name '_G_lift'])(curr{:}) = {{'G', length(ind_nonscalar_G)}, 0, inf}; G_lift = G_curr(ind_nonscalar_G); @@ -1134,7 +1134,7 @@ function create_mpec_functions(obj) else % lifting with only those that are not scaler % define lift vairables - [ind_scalar,ind_nonscalar_x1, ind_map] = find_nonscalar(G,x); + [ind_scalar,ind_nonscalar_x1, ind_map] = find_nonscalar(G,x,p); n_lift_x1 = length(ind_nonscalar_x1); if n_lift_x1 == 0 % TODO(@anton) Figure out what this does. @@ -1170,7 +1170,7 @@ function create_mpec_functions(obj) g = [g;x2-H]; else % lifting with only those that are not scaler - [ind_scalar,ind_nonscalar_x2, ind_map] = find_nonscalar(H,x); + [ind_scalar,ind_nonscalar_x2, ind_map] = find_nonscalar(H,x,p); n_lift_x2 = length(ind_nonscalar_x2); if n_lift_x2 == 0 % TODO(@anton) Figure out what this does. diff --git a/src/find_nonscalar.m b/src/find_nonscalar.m index 966f329..7e46e14 100644 --- a/src/find_nonscalar.m +++ b/src/find_nonscalar.m @@ -1,19 +1,23 @@ -function [ind_scalar,ind_nonscalar, ind_map] = find_nonscalar(g,w) +function [ind_scalar,ind_nonscalar, ind_map] = find_nonscalar(g,w,p) % Implemented by Anton Pozharskiy within nosnoc: https://github.com/nosnoc/nosnoc % This function returns the indices of g which contain elements of the symbolic vector w, along with its - % nonscalar indicies and the ind_map which are the indices in w which correspond to scalar elements of g. +% nonscalar indicies and the ind_map which are the indices in w which correspond to scalar elements of g. + if ~exist('p') + p = SX([]); + end import casadi.* % Take jacobian of g with respect to w. - ind_g_fun = Function('ind_G', {w}, {g.jacobian(w)}); + ind_g_fun = Function('ind_G', {w,p}, {g.jacobian(w)}); + p_val = SX(ones(length(p),1)); % HERE BE DRAGONS: % `ind_g_fun(w) == 1` creates a symbolic array where literal `1` are ones, other constants are zeros and % symbolics are nans when converted to a casadi.DM. % % We want to find rows which contain only exactly one nonstructural 0, 1, or nan. % We get this from using find on the sparsity patern of the DM. - sp = DM(ind_g_fun(w) == 1).sparsity; + sp = DM(ind_g_fun(w,p) == 1).sparsity; sub = sp.find; - [ind_g1, ind_g2] = ind2sub(size(ind_g_fun(w)),sub); + [ind_g1, ind_g2] = ind2sub(size(ind_g_fun(w,p)),sub); % transpose because groupcounts expects column vector c=groupcounts(ind_g1'); ind_scalar=ind_g1(find(c==1));