Add ampl_repn option to Incidence Analysis

[Robbybp]

Fixes #3068

Summary/Motivation:

The main motivation is to add an option that greedily removes 0*x terms from nonlinear subexpressions. Secondarily, AMPLRepnVisitor likely has performance benefits compared to generate_standard_repn and is consistent with the problem structure sent to nl files.

Changes proposed in this PR:

• IncidenceMethod.ampl_repn option in pyomo.contrib.incidence_analysis.config

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

1. I agree my contributions are submitted under the BSD license.
2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

[Robbybp]

As written, this PR introduces a large performance bottleneck. For example:

import os
import examples.dae.distill_DAE as distill_module
import pyomo.environ as pyo 
from pyomo.contrib.incidence_analysis import IncidenceGraphInterface
from pyomo.contrib.incidence_analysis.config import IncidenceMethod
from pyomo.common.timing import TicTocTimer

datfile = os.path.join(
    os.path.dirname(distill_module.__file__),
    "distill.dat",
)

instance = distill_module.model.create_instance(datfile)
pyo.TransformationFactory("dae.finite_difference").apply_to(
    instance, wrt=instance.t, nfe=50,
)

timer = TicTocTimer()

timer.tic()
igraph = IncidenceGraphInterface(instance)
timer.toc("default (standard_repn)")
igraph = IncidenceGraphInterface(instance, method=IncidenceMethod.ampl_repn)
timer.toc("ampl_repn")

gives the output:

[    0.00] Resetting the tic/toc delta timer
[+   0.56] default (standard_repn)
[+  12.97] ampl_repn

[Robbybp]

As written, this PR introduces a large performance bottleneck.

We can make this a little better by re-using the visitor (and var_map, subexpression_cache, etc.) for all constraints:

[    0.00] Resetting the tic/toc delta timer
[+   0.57] default (standard_repn)
[+   6.68] ampl_repn

[Robbybp]

Changing to WIP until this performance hit is figured out.

[jsiirola]

Can you add an option to use / time pyomo.repn.linear.LinearRepnVisitor?

Some notes on the NL writer:

• the list of nonlinear variable IDs can contain duplicates (it is the list of placeholders in the NL expression, so duplicates are expected)
• the list of nonlinear variable IDs does NOT any nonlinear variables from named expressions (but would include IDs from named subexpressions appearing nonlinearly (or containing nonlinear components).

[Robbybp]

• the list of nonlinear variable IDs does NOT any nonlinear variables from named expressions (but would include IDs from named subexpressions appearing nonlinearly (or containing nonlinear components).

This is only the case if export_defined_variables=True, right? I am running with export_defined_variables=False for simplicity. The distillation model has no named expression, so this should not be the cause of the performance hit in this case.

[Robbybp]

Can you add an option to use / time pyomo.repn.linear.LinearRepnVisitor?

I will try this. In the meantime, these tests

timer = TicTocTimer()
timer.tic()
igraph = IncidenceGraphInterface(instance)
timer.toc("default (standard_repn)")
igraph = IncidenceGraphInterface(instance, method=IncidenceMethod.standard_repn_compute_va
timer.toc("standard_repn_compute_values")
igraph = IncidenceGraphInterface(instance, method=IncidenceMethod.ampl_repn)
timer.toc("ampl_repn")
pyo.SolverFactory("ipopt").solve(instance)
timer.toc("solve with Ipopt")

give these results

[    0.00] Resetting the tic/toc delta timer
[+   0.57] default (standard_repn)
[+   0.56] standard_repn_compute_values
[+   6.65] ampl_repn
[+   0.21] solve with Ipopt

We can write the nl file and solve with Ipopt in a fraction of a second, so the extremely long compilation time in IncidenceGraphInterface is likely due to some error using AMPLRepnVisitor on my part.

[Robbybp]

We can make this a little better by re-using the visitor (and var_map, subexpression_cache, etc.) for all constraints

I had missed passing the visitor to another spot where get_incident_variables was called. With this fix, we get:

[    0.00] Resetting the tic/toc delta timer
[+   0.58] default (standard_repn)
[+   0.58] standard_repn_compute_values
[+   0.60] ampl_repn
[+   0.21] solve with Ipopt

This couples IncidenceGraphInterface and get_incident_variables way more tightly than I initially intended, and will require a bit more thought before I'm ready to merge this PR. In particular, every time we parse a group of constraints for the incidence graph, I need to make sure that if AMPLRepnVisitor is used, we re-use the visitor among all constraints. Configuration of the visitor (e.g. with an initialized var_map or with export_defined_variables set) should happen in one location so we don't have multiple locations that get out of sync.

Given that generate_standard_repn(expr, compute_values=True) seems to achieve my goals as well, AMPLRepnVisitor may not be worth the additional complexity at this point.

[jsiirola]

It's surprising to me that the visitor's ramp-up is so significant, but I could see it being an issue. You could consider making get_incidence_variables be a functor that you create instead of a function (then it can be a different class depending on which walker you want to use): it can create / manage it's own context, and might reduce the amount of "if method is ...".

Separately, if you are going to walk / collect a bunch of expressions, it might make sense to to it all within a PauseGC() context (that speeds thing up for the writers a significant amount).

[Robbybp]

Separately, if you are going to walk / collect a bunch of expressions, it might make sense to to it all within a PauseGC() context (that speeds thing up for the writers a significant amount).

I did not notice any difference when using PauseGC here.

[Robbybp]

I just pushed a change that moves construction of the AMPLRepnVisitor to the config system. The motivation for this is that I don't want to have code in the main interface class and functions handle a particular IncidenceMethod specially. However, this involves some hacking of ConfigDict, so I may revert this change. The alternative would be to call a method like kwds = _construct_visitor_if_necessary(kwds) before passing kwds to IncidenceConfig.

[Robbybp]

At some point I needed a singleton to trigger re-construction of IncidenceConfig's visitor. I'm not sure if this is still necessary, and will re-visit.

[Robbybp]

By default, _ReconstructVisitor is passed to the validator, and we construct a new visitor. However, as written, cloning this ConfigValue (via __call__) will create a new value pointing to the same visitor. This effectively means that, unless a new visitor is explicitly received, the same default visitor is used for all new ConfigValues. I.e. we re-use the visitor for multiple models, or even worse, for the same model before and after modification. The code below implements a hack around ConfigDict to re-construct the visitor instead of using the default.

[Robbybp]

If we are using ampl_repn and did not provide a visitor, explicitly trigger reconstruction of the visitor.

[Robbybp]

The most recent changes still give the timing results that I expect:

[    0.00] Resetting the tic/toc delta timer
[+   0.64] default (standard_repn)
[+   0.66] standard_repn_compute_values
[+   0.60] ampl_repn
[+   0.97] solve with Ipopt

[Robbybp]

I've updated this to use a get_config_from_kwds function that just constructs a new AMPLRepnVisitor if method=ampl_repn and a visitor is not provided. This removes the subclass ConfigDict with a hacky re-implementation of __call__.

I think this PR is ready for review.

[codecov]

Codecov Report

Attention: 3 lines in your changes are missing coverage. Please review.

Comparison is base (835a2df) 88.29% compared to head (5a5fa33) 88.30%.
Report is 15 commits behind head on main.

Files	Patch %	Lines
pyomo/contrib/incidence_analysis/incidence.py	94.44%	2 Missing ⚠️
pyomo/contrib/incidence_analysis/config.py	96.15%	1 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #3069   +/-   ##
=======================================
  Coverage   88.29%   88.30%           
=======================================
  Files         832      832           
  Lines       92234    92293   +59     
=======================================
+ Hits        81438    81495   +57     
- Misses      10796    10798    +2     

Flag	Coverage Δ
linux	86.05% <95.65%> (+<0.01%)	⬆️
osx	75.56% <95.65%> (+0.01%)	⬆️
other	86.25% <95.65%> (+<0.01%)	⬆️
win	83.46% <95.65%> (+<0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Robbybp]

Following our conversation at the dev meeting today, I tested a linear_repn method that uses LinearRepnVisitor, and got the following results on a larger version of the distillation model that I benchmarked above (nfe=300):

[    0.00] Resetting the tic/toc delta timer
[+   4.19] default (standard_repn)
[+   4.24] standard_repn_compute_values
[+   4.02] ampl_repn
[+   4.29] linear_repn
[+   6.25] solve with Ipopt

It seems LinearRepnVisitor is slightly slower than AmplRepnVisitor for this nonlinear model, which surprised me. I think it should be possible to create a visitor that creates and reduces the linear part of the expression, but stores the nonlinear part just as a set of variable ids. I think I tried this last year, but was tripped up by (a) trying to support None-valued parameters and (b) something about it being hard to test. Maybe I'll give this another shot at some point, although I do like the idea of incidence graphs being consistent with the structure exported to NL files.

[blnicho]

I noted a few imports that can be cleaned up but otherwise I think this looks good!