GitLab

Regression in 6303bd57. Forgot to include this particular call to
__res_vals in that commit.

A method in OptimizationProblem was returning a "LookupTable" type, with
the class not defined anywhere in scope. One of the CI/CD style checkers
became a bit more aggressive recently and started to catch this as an
error.

We now properly defined a base LookupTable class, from which the one in
csv_lookup_table_mixin.py inherits.

One of the CI/CD style checkers became a bit more aggressive recently,
and started to catch this "l" variable name as an error.

A call to math's isfinite() is about 20 times faster than calling
NumPy's isfinite() (~66 ns vs 1.2 us). We also save some attribute
lookup overhead (20 ns) by importing isfinite() directly, instead of
calling `math.isfinite()` (the latter takes 85 ns).

Instead of calling `canonical_state` every time via the use of the
AliasDict, we can just figure out all parameters and their aliases once.
That way we can rely on the (very) fast containment check that a native
Python set() offers.

The IO timeseries do not change in the update loop. Instead of
repeatedly requesting the timeseries, we can just store the values in an
internal cache.

Going to a NumPy array first is a faster way of achieving the same
result.

We do not know the size of the run a priori. Appending to NumPy arrays
is very slow (~5 us per call), whereas it is very fast for lists (~90
ns).

Most importantly is that AliasDicts need to be able to apply the
negative operator. A list does not support this operator itself, so we
have to work around it ourselves.

We can avoid the somewhat expensive `canonical_signed` call. Note that
this is typically only expensive because we tend to call get_var/set_var
millions of times.

Vertcat calls are relatively expensive to do. The vertcat call the
prepends `dt` was therefore taking a substantial amount of time, equal
or more than the actual solving of the equations with rootfinder for
linear models.

Instead of vertcatting, we can just specify the `dt` as additional
parameter to the ca.Function that we pass to rootfinder. We can then
also just pass it separately when calling `__do_step` to do a rootfinder
step.

The variables in the problem cannot change from update call to update
call. It therefore makes no sense to make the rather expensive call to
get all variables time and again. We can save some time by caching the
return values of these functions.

Calling `self.io.times_sec` would make a copy, to be sure that the user
could not overwrite the internal array of times. Especially in
simulation, where we can do many time steps, this copy() call was rather
expensive. The same intent can be accomplished by disallowing writes to
the array we return.

Note that the times in seconds are immutable as soon as they are
generated (which is usually done when a request is done like
`get_timeseries_sec`). At that point, they are no longer allowed to
change.

The method `get_timeseries_values` does not exist. Furthermore, the
check is supposed to check only collisions between the data specified
with CSV files. When the user is using more than one IO mixin,
overwriting initial state may actually be intended.

`io.times_sec` is a property, not a method.

`NetCDF4` depends on `cftime`, the latter released a version 1.1 which
changed the behavior of the `num2date` function. This function now no
longer returns Python `datetime` like objects. It also added a new
convenience method `num2pydate` which does.

We update our code to handle both old and new versions of `cftime`.

See e.g. https://github.com/Unidata/cftime/issues/134 for details on why
this breaking change was made.

`NetCDF4` depends on `cftime`. For some reason, older versions of cftime
had some rounding errors when handling arrays. For the mentioned test,
when comparing the NetCDF file to the PI-XML file, we inadvertently
relied on this behavior. In actuality, both files have the same
datetimes.

This commit makes sure that the three input files (CSV, PI-XML, and
NetCDF) have different datetimes. The CSV file has offsets of 0-5
minutes, the NetCDF file 0-5 hours, and the PI-XML now has 0-5 seconds.

See also the following issues/PRs on `cftime` which exposed this:
https://github.com/Unidata/cftime/issues/144
https://github.com/Unidata/cftime/pull/146

This reverts commit 30716146.

The second positional argument to integral() is the time in seconds from
which to start integrating. It is not the ensemble member, which is
easier passed as a keyword argument if we want to take the integral over
the entire optimization horizon.

`NetCDF4` depends on `cftime`, the latter released a version 1.1 which
changed the behavior of the `num2date` function. It also added a new
convenience method `num2pydate` which seems to sometimes lose the
millisecond information.

We only implicitly depend on `cftime` via `NetCDF4`, but with this
breaking change on their end we need to make sure that the correct
version is installed. We therefore add `cftime` to our list of
(optional) dependencies, and rely on the version we know is correct
(which excludes 1.1.0).

See e.g. https://github.com/Unidata/cftime/issues/134 for details on why
this breaking change was made.

See https://github.com/Unidata/cftime/issues/144 for the issue on
num2pydate losing millisecond information.

It is often useful to make sure that the constraints are linear, such
that convexity is guaranteed. Depending in the options of for example
GoalProgrammingMixin, it can be possible that a quadratic objective in
one priority leads to a quadratic constraint in the next priority, which
a solver may not be able to solve (or may not even be convex anymore).

Eventually this check might be rewritten such that it (also) logs
whether a problem is (MI)LP/QP/QCQP/NLP, and also say something about
the convexity.

Interior point solvers need the linear independence constraint
qualification (LICQ) to be satisfied. Because a solver like IPOPT is
callback based, it cannot detect even the simplest of dependencies that
violate this qualification.

This commit adds some checks on typical linear dependencies that arise.
Currently these checks consist of:

1. Duplicate constraints that set a variable to a certain value.
   For example the constraint `x = 1` twice.

2. Setting a equality constraint on a variable that has equal bounds.
   For example having a constraint `x = 1` when x has bounds [1, 1].

Many solvers (IPOPT, CPLEX, etc) would prefer the coefficients in the
Jacobian to be in the order of [1E-2, 1E2], and also for the range of
coefficients for a certain variable or in a certain constraint to be
rather limited (say 1E3). The check added in this commit gives feedback
to the user on what constraints/variables do not meet these targets.

Note that not satisfying these limits/ranges does not mean that the
problem will not solve or even be slower, because solvers tend to be
quite smart in how they themselves rescale equations or evade numerical
instabilities. That said, the larger (more variables, more constraints)
a problem gets, the more sensitive solvers tend to become to bad
scaling.

When setting up or debugging a model it can be useful to do typical
checks to see if you made a mistake. These checks can be rather
computationally expensive, which is why it should be easy to disable
them when running operationally. A (for now private) flag
`_debug_check_level` is added to OptimizationProblem/SimulationProblem
that can be used to enable/disable certain checks.

Foreseeing quite a few of these checks, it should be easy to filter
checks based on verbosity and how expensive they are to perform. A
typical way of doing this is with an enum, just like logging levels
(DEBUG, INFO, WARNING, etc). To this end the value of
`_debug_check_level` can be set to the value of a `DebugLevel` enum.
This enum has with 5 levels, with room for future refinement. The
default check level is in the middle (MEDIUM), with two available below
(NONE, LOW) and two above (HIGH, VERYHIGH).

With more of these checks, and considering how expensive and useful they
might be for a certain type of issue, it is also possible to filter more
precise by setting the value of `_debug_check_level` to a function. That
way one can filter a specific test on its level and fully qualified
name.

For some checks it might make sense to accept some arguments. For
example, a certain tolerance value above which messages should be
logged. For this purpose the `_debug_check_options` class variable of
OptimizationProblem/SimulationProblem can be used. This is a dictionary
that maps a fully specified check name to (a dictionary of) additional
keyword arguments.

Note that it is likely that some details of these debug checks will
change, which is why it is only part of the private API for now.

Closes #1042

Allow users to specify the name of the input files used to extract that
initial states, parameters and ensemble information. Note that this was
already possible when using PIMixin.

We use np.genfromtxt with dtype=None, which means that it will guess. If
a column does not have any value specified, it will guess the dtype
boolean with values `False`. This would eventually be turned into a
Timeseries of `0.0`. We however want a Timeseries full of `np.nan`.

It is difficult to make np.genfromtxt read in floats, when the first
column can sometimes be a string. Therefore, we detect the case where
the dtype of a column is boolean, and force a conversion to float.
Basically, we therefore only support reading in floats (and a optionally
string for the first column).

Note that CSVMixin (both in optimization and simulation) forces a
conversion to float already anyway, so there is no freedom lost.

Closes #1128

One of the CI/CD style checkers became a bit more aggressive recently,
and started to catch this redundancy as an error.

This function can be very costly to regenerate. For NZV, it can easily
take 30 seconds. Because SXFunctions can be cached, we can also just
cache this function to save a lot of time.

It is necessary to overwrite the seed within the HomotopyMixin class
whenever the results of the latest optimization run are stored here.
Note that a necessary requirement for this is that theta must be greater
than zero. Moreover, there are only two situations in which we need to
overwrite the seed: when the class GoalProgrammingMixin is not used, or
whenever it is the first run of the goal programming loop. To detect
this, and to make it clear that we internally rely on this variable, we
rename __first_run to _gp_first_run.

If the start attribute is the default of zero, we assume that it has not
been explicitly set. The type of the start attribute can however be more
than just a float; it can also be a bool, int (depending on the Python
type of the variable). We reformulate the check such that it catches all
these cases.

Note that the default seed value is zero anyway in CollInt, so in
practice the result will often be the same. The only difference is that,
previously, failing check on zero (e.g. because the start attribute was
in int) would _explicitly_ set the seed to zero. This could trigger code
that checks if a seed for a certain variable has been specified.

In the cascading channel example, one of the water level initial
conditions was overdetermined, which lead to a degenerate Jacobian
matrix. We remove one of these initial conditions.

Closes #1099

The generation of plots requires matplotlib, which is not an explicit
dependency of our test environment. Instead of adding matplotlib to our
test dependencies, we can just disable the generation of plots, as it
does not add anything to the tests.

Note that the generation of plots might make sense when debugging why
the test is failing, but that is a different use case for which the user
can easily just install matplotlib.

Pymoca requires (as a workaround for a CasADi shortcoming) that
`expand_mx` is set when using `eliminable_variable_expression`.

The behavior in GoalProgrammingMixin with keep_soft_constraints=True is
that the constraint_relaxation of priority 2 is applied to the objective
of priority 2. This is because we call `goal_programming_options` in
GoalProgrammingMixin at the end of priority 2 to make this constraint.

In SinglePassGoalProgrammingMixin, we would get the
constraint_relaxation by calling `goal_programming_options` at the start
of the next priority, so the constraint_relaxation of priority 3 would
be applied to the objective of priority 2. This is not necessarily
wrong, but does differ from GoalProgrammingMixin. We therefore fix the
discrepancy by making SinglePassGoalProgrammingMixin behave equal in
this sense to GoalProgrammingMixin; the constraint relaxation at
priority 2 is applied to the objective of priority 2.

The naive implementation until this commit was O(N), which means that
state lookup could take up a significant portion of the computational
time when there were many variables/parameters. We implement a O(1)
state index lookup for SimulationProblem using a dictionary, just like
fb6bbf08 did for optimization's CollocatedIntegratedOptimizationProblem.

Closes #1052

If there are many timeseries and/or a lot of reference times, we can
save quite some time by converting only the times of the current series
to a set. Imagine a case where a lot of timeseries are only specified at
t0 (for the initial state) or are empty, but the total length of all
times is a series of length 100,000+.

Although simulation mode only makes sense for a single ensemble member,
PI-XML files can contain multiple ensemble members. PIMixin in
simulation would always read the first ensemble member before this
commit. Instead of forcing the user to write a separate file for each
ensemble member, we now allow a new attribute `pi_ensemble_member` that
specifies what ensemble member should be read from the PI-XML file.

Setting the goal programming option 'mu_reinit' to False would trigger a
CasADi RuntimeError, as the `gather_stats` option no longer exists.
CasADi instead just always collects stats, so we can simply remove the
`gather_stats` option.

Note that the `gather_stats` option was deprecated by CasADi about 4
years ago at time of writing this commit.