Plot constraint violations

[Peter230655]

I tried to make the prints of the method plot_constraint violations a bit easier to read. In particular

• maximally 10 bars per plot. If there are more constraint violations, more plots will be made
• changed printing of the constraints, now in exponential writing, like e.g. $1.24 \cdot 10^{-3}$ instead of 0.000123456

[Peter230655]

I have a very basic question about testing such a thing:
When I started, I made a new branch, and in this branch, I opened direct_allocation.py
To test what I did, I simply copied a simulation below direct_allocation. and ran the wohle thing every time.

But this seems inefficient if I want to test it with different simulations.
In a normal simulation I wrtie:

• from opty.direct_allocation import Problem.
  But what do I write, if I want to use my modified direct_allocation ?

Thanks for any help!

[Peter230655]

The error Github found may be due to this:
If the user give the axes, he must give 'enough' axes, just two may not be the right amount.
It should be

• 2 if len(instant_constraints) $\leq$ 10
• len(instant_constraints) // 10 + 1 otherwise

[moorepants]

My approach to testing edits I make to a package, like opty, is to first create a conda environment with all the development dependencies. We include this file in the repo for such a thing: https://github.com/csu-hmc/opty/blob/master/opty-dev-env.yml. After that, I install the development version of opty in a way that any import of opty will always reflect the current state of the opty repo on my computer.

$ cd /path/to/opty/repo
$ conda env create -f opty-dev-env.yml
$ conda activate opty-dev
$ python -m pip install --no-deps -e .  # installs version in the current directory, i.e. the development version

Now as long as I have that conda environment activated it will use the opty version in /path/to/opty/repo.

So then, if I want to build the docs, just:

$ cd docs
$ make html

All of the constraint violations plots should generate using your edits to the direct_collocation.py file.

[moorepants]

Looks like your code doesn't work on the cycling or drone example:

Extension error:
Here is a summary of the problems encountered when running the examples:

Unexpected failing examples (2):

    ../examples-gallery/plot_one_legged_time_trial.py failed leaving traceback:

    Traceback (most recent call last):
      File "/home/runner/work/opty/opty/examples-gallery/plot_one_legged_time_trial.py", line 748, in <module>
        problem.plot_constraint_violations(solution, axes=axes)
      File "/home/runner/work/opty/opty/opty/utils.py", line 218, in wrapper
        func(*args, **kwargs)
      File "/home/runner/work/opty/opty/opty/direct_collocation.py", line 524, in plot_constraint_violations
        axes[i+1].bar(
        ~~~~^^^^^
    IndexError: index 2 is out of bounds for axis 0 with size 2

    ../examples-gallery/plot_drone.py failed leaving traceback:

    Traceback (most recent call last):
      File "/home/runner/work/opty/opty/examples-gallery/plot_drone.py", line 235, in <module>
        prob.plot_constraint_violations(solution, axes=axes)
      File "/home/runner/work/opty/opty/opty/utils.py", line 218, in wrapper
        func(*args, **kwargs)
      File "/home/runner/work/opty/opty/opty/direct_collocation.py", line 524, in plot_constraint_violations
        axes[i+1].bar(
        ~~~~^^^^^
    IndexError: index 2 is out of bounds for axis 0 with size 2

-------------------------------------------------------------------------------
make: *** [Makefile:25: html] Error 2
Error: Process completed with exit code 2.

[moorepants]

If the user give the axes, he must give 'enough' axes, just two may not be the right amount.

As long as you document this, I think it is reasonable to ask the user to pass in the right number and combination of axes. They already have to do that if you create the axes first.

[moorepants]

A development setup of opty is explained in the README.

[Peter230655]

If the user give the axes, he must give 'enough' axes, just two may not be the right amount.

As long as you document this, I think it is reasonable to ask the user to pass in the right number and combination of axes. They already have to do that if you create the axes first.

This seems to be the problem with the drone and the cycling simulations. Both more than 10 constraints and give axes,
How should I document it? In the method itself?

[tjstienstra]

In addition to what Jason has already mentioned. You can easily run personal scripts, like some fun example, by creating a file like temp_playground.py in the root directory:

.
├── .git
│   ├── ...
│   ├── info
│   │    └── exclude
│   └── ...
├── ...
├── opty
│   ├── tests
│   ├── __init__.py
│   └── ...
├── ...
├── README.md
└── temp_playground.py

This file can just use from opty.direct_collocation import Problem and you can run it using python -m temp_playground.py. Just make sure not to add this file to git. To make my life easier I like to update the .git/info/exclude to also ignore all temp_*.py files. There are of course tons of other options, but this one I personally like.

[Peter230655]

In addition to what Jason has already mentioned. You can easily run personal scripts, like some fun example, by creating a file like temp_playground.py in the root directory:

.
├── .git
│   ├── ...
│   ├── info
│   │    └── exclude
│   └── ...
├── ...
├── opty
│   ├── tests
│   ├── __init__.py
│   └── ...
├── ...
├── README.md
└── temp_playground.py

This file can just use from opty.direct_collocation import Problem and you can run it using python -m temp_playground.py. Just make sure not to add this file to git. To make my life easier I like to update the .git/info/exclude to also ignore all temp_*.py files. There are of course tons of other options, but this one I personally like.

Thanks to both of you!
I have to think about these suggestions, I am a bit better with Github / git - but they are still far away from being my friends! :-)

[Peter230655]

The errors with cycling and with drone


were as they should be, I believe. I added a check if len(axes) is not correct.

[Peter230655]

My approach to testing edits I make to a package, like opty, is to first create a conda environment with all the development dependencies. We include this file in the repo for such a thing: https://github.com/csu-hmc/opty/blob/master/opty-dev-env.yml. After that, I install the development version of opty in a way that any import of opty will always reflect the current state of the opty repo on my computer.

$ cd /path/to/opty/repo
$ conda env create -f opty-dev-env.yml
$ conda activate opty-dev
$ python -m pip install --no-deps -e .  # installs version in the current directory, i.e. the development version

Now as long as I have that conda environment activated it will use the opty version in /path/to/opty/repo.

So then, if I want to build the docs, just:

$ cd docs
$ make html

All of the constraint violations plots should generate using your edits to the direct_collocation.py file.

I think I finally understood! Had to think about it for a good while.

[Peter230655]

Prints a warning if axes are given, with sharex = True (which makes no sense here).
If sharex = True it will generater plots, but they are nonsense.

[Peter230655]

I rotate the labels -90°. This way, they are always lined up with the ticks and never smear into each other.
If the state variables are defined as functions, their argument is also displayed as mantissa * 10^(exponent). If the argument is exactly 10, the display looks odd sometimes line q(--1.0 * 10^1) instead of q(1.0 * 10^1). I do not know why. (you can run the non-variable drone simulation with duration = 10.0 to see what I mean)
@moorepants

• which rotation do you prefer?
• should I keep the numbers as mantissa * 10^(exponent), or better like rounded numbers? Of course the mantissa is also rounded. I have a slight preference for rounded numbers, but this is your program!
  What I mean ist this: 50.12345 displayed as 5.01*10^1 or as 50.1

[Peter230655]

The failed check with cycling and drone is that they do not give the right number of axes, as before.

[moorepants]

We typically don't have print statements in library code because their display can't be controlled. Best to avoid warnings and just raise errors.

[moorepants]

I also don't think we need to warn about things like this. The best approach here is to explain in the docstring what the user must do if they want to pass in their own axes.

[Peter230655]

I also don't think we need to warn about things like this. The best approach here is to explain in the docstring what the user must do if they want to pass in their own axes.

If sharex = True, the program will not interrupt, but the plots are meaningless. Do you mean, I just let this happen, and explain in the method, that it should be set to False or removed? Should this be mentioned in the opty Documentation?
If YES, I have to find out how to do this.

[moorepants]

Do you mean, I just let this happen, and explain in the method, that it should be set to False or removed?

Yes, but I don't think you need to explain anything about sharex. plt.subplots() has a very long list of arguments and it is not our job to explain them all here. You simply have to state to the user the minimal information they need to know to pass in a correct axes object.

[moorepants]

which rotation do you prefer?

I don't know because we need to see how it looks on all the examples and any other test cases you are using.

[moorepants]

The failed check with cycling and drone is that they do not give the right number of axes, as before.

You have to correct that.

[Peter230655]

The failed check with cycling and drone is that they do not give the right number of axes, as before.

You have to correct that.

Do you mean I correct the drone example? If Yes, same would be true with the cycling example (?)
If yes, would I give the corrrect number of axes, or remove the keyword ax when the function is called?
Would I change drone and cycling in the same PR as this one, so the check will run, or change with separate PRs?

[moorepants]

Do you mean I correct the drone example?

Yes.

If Yes, same would be true with the cycling example (?)

Yes.

If yes, would I give the corrrect number of axes, or remove the keyword ax when the function is called?

Give correct number.

Would I change drone and cycling in the same PR as this one, so the check will run, or change with separate PRs?

All should happen in this PR.

[Peter230655]

1. I do not know, how my plot_two_body_skateboard got into this branch. This was intended for examples-gallery, but no PR, yet. git still not my friend! :-)
2. Shoud I make a PR with normal numbers, not exponential, so you can see the difference and decide which one you prefer?

[moorepants]

You have mistakenly committed a skateboard example:

[Peter230655]

Those both look good. I'm fine with the normal numbers as long as the number of figures after the decimal are reduced (like you have it).

I will switch to normal numbers, and play around with exponentials (eliminate the double minus) and variable h next.
I will make a PR later today or tomorrow morning.

[moorepants]

These were ratio 3:4 aspect, I suspect on purpose. Is there a reason to change that?

[Peter230655]

These were ratio 3:4 aspect, I suspect on purpose. Is there a reason to change that?

No, I did not give any thought. I just changed to make it run with my proposal

[moorepants]

There is no way for the user to know what the value of bars_per_plot is in the docstring. Just state the number. This could be more simply stated in general. Something like, "six instance constraints will be plotted per axes, so you'll need to pass in enough rows to accomodate the number of instance constraints".

[Peter230655]

There is no way for the user to know what the value of bars_per_plot is in the docstring. Just state the number. This could be more simply stated in general. Something like, "six instance constraints will be plotted per axes, so you'll need to pass in enough rows to accomodate the number of instance constraints".

My approach was this: I select between 6 to 10 bars per plot to prevent situations like with the double pendulum., where there is only on bar on the last plot. If the wrong number of axes is given, the method will tell the user how many are needed - unless he only gives one, which is always wrong.

[moorepants]

It is better to leave the fig size for the user to create manually. Setting figsize inside a function like this can result in giant figures.

[moorepants]

Python has printing controls that work like this:

In [7]: '{:1.3f}'.format(1.2394359593939)
Out[7]: '1.239'

is there a reason to not use something like that instead of literally rounding the numerical values?

[Peter230655]

Python has printing controls that work like this:

In [7]: '{:1.3f}'.format(1.2394359593939)
Out[7]: '1.239'

is there a reason to not use something like that instead of literally rounding the numerical values?

Timo suggested this to me, so I used it, since it worked. I had no idea how to do such things.

[moorepants]

Ok, I see now there may not be a straight forward way to use Python printing methods on SymPy Floats in an expression. I thought there should be a setting on the StrPrinter that would control display of decimals, but there doesn't seem to be one.

[moorepants]

Maybe Timo already showed this (I struggle to find his comment), but this is a way to control printing of floats:

In [47]: from mpmath.libmp import prec_to_dps, to_str as mlib_to_str

In [48]: class MyPrinter(sm.StrPrinter):
    ...:     _default_settings = {
    ...:         "order": None,
    ...:         "full_prec": "auto",
    ...:         "sympy_integers": False,
    ...:         "abbrev": False,
    ...:         "perm_cyclic": True,
    ...:         "min": None,
    ...:         "max": None,
    ...:         "dps": None,
    ...:     }
    ...: 
    ...:     def _print_Float(self, expr):
    ...:         prec = expr._prec
    ...:         if prec < 5:
    ...:             dps = 0
    ...:         else:
    ...:             dps = prec_to_dps(expr._prec)
    ...:         if self._settings['dps']:
    ...:             dps = self._settings['dps']
    ...:         if self._settings["full_prec"] is True:
    ...:             strip = False
    ...:         elif self._settings["full_prec"] is False:
    ...:             strip = True
    ...:         elif self._settings["full_prec"] == "auto":
    ...:             strip = self._print_level > 1
    ...:         low = self._settings["min"] if "min" in self._settings else None
    ...:         high = self._settings["max"] if "max" in self._settings else None
    ...:         rv = mlib_to_str(expr._mpf_, dps, strip_zeros=strip, min_fixed=low, max_fixed=high)
    ...:         if rv.startswith('-.0'):
    ...:             rv = '-0.' + rv[3:]
    ...:         elif rv.startswith('.0'):
    ...:             rv = '0.' + rv[2:]
    ...:         if rv.startswith('+'):
    ...:             # e.g., +inf -> inf
    ...:             rv = rv[1:]
    ...:         return rv
    ...: 

In [49]: MyPrinter(settings={'dps': 3}).doprint(f(1.329294))
Out[49]: 'f(1.33)'

[moorepants]

Hah, I found Timo's solution which was the same.

[Peter230655]

Maybe Timo already showed this (I struggle to find his comment), but this is a way to control printing of floats:

In [47]: from mpmath.libmp import prec_to_dps, to_str as mlib_to_str

In [48]: class MyPrinter(sm.StrPrinter):
    ...:     _default_settings = {
    ...:         "order": None,
    ...:         "full_prec": "auto",
    ...:         "sympy_integers": False,
    ...:         "abbrev": False,
    ...:         "perm_cyclic": True,
    ...:         "min": None,
    ...:         "max": None,
    ...:         "dps": None,
    ...:     }
    ...: 
    ...:     def _print_Float(self, expr):
    ...:         prec = expr._prec
    ...:         if prec < 5:
    ...:             dps = 0
    ...:         else:
    ...:             dps = prec_to_dps(expr._prec)
    ...:         if self._settings['dps']:
    ...:             dps = self._settings['dps']
    ...:         if self._settings["full_prec"] is True:
    ...:             strip = False
    ...:         elif self._settings["full_prec"] is False:
    ...:             strip = True
    ...:         elif self._settings["full_prec"] == "auto":
    ...:             strip = self._print_level > 1
    ...:         low = self._settings["min"] if "min" in self._settings else None
    ...:         high = self._settings["max"] if "max" in self._settings else None
    ...:         rv = mlib_to_str(expr._mpf_, dps, strip_zeros=strip, min_fixed=low, max_fixed=high)
    ...:         if rv.startswith('-.0'):
    ...:             rv = '-0.' + rv[3:]
    ...:         elif rv.startswith('.0'):
    ...:             rv = '0.' + rv[2:]
    ...:         if rv.startswith('+'):
    ...:             # e.g., +inf -> inf
    ...:             rv = rv[1:]
    ...:         return rv
    ...: 

In [49]: MyPrinter(settings={'dps': 3}).doprint(f(1.329294))
Out[49]: 'f(1.33)'

He did show this or something similar, but I liked the simple round(a, 3) much better! :-)

[moorepants]

I'll merge if you remove the figsize kwarg inside the function and try to simplify the docstring explanation to how many rows they'll need in the axes.

[Peter230655]

I'll merge if you remove the figsize kwarg inside the function and try to simplify the docstring explanation to how many rows they'll need in the axes.

Great!
Removing the kwarg is easy, frankly I do not recall why I did it.
As regards the number of axes, would this be o.k.:
"Select at least 2 axes. If more are needed, the method will tell the use how many are needed." ?

I just removed the kwarg. Then the output of the drone looks like this.
This is probably why I changed it.

My idea, based on you correct complaint about the double pendulum was that I select the number of bars per plot so all plots have approximately the same number of bars. It seems to me that in this case I have to give a figsize depending on the number of plots?

I can find no way to avoid giving figsize=(....., .....).
I tried to copy what was done in plot_trajectories, but it did not work.

[Peter230655]

My idea, based on your (very justified!) complaint about the double pendulum is this:
bars_per_plot is between 6 and 10.
Say, there are 11 instance constraints, as in the double pendulum. Then there will be two plots for the instance constraints, the first has 6 bars, the second one has 5 bars
Say, there are 16 instance constraints, then both plots have 8 bats
That is, I select the bars_per_plot so, that (bars_per_plot - num_instance_constraints % bars_per_plot) is minimal.
(Of course, I have to separately consider when (bars_per_plot - num_instance_constraints % bars_per_plot) = 0, etc, hence my program looks maybe a bit confusing. Also my execution is likely not the most elegant one)

• I found no way to get the plots right unless I give figsize = (..., ...). The method plot_trajectories does not give this and the plots are not good unless the user provides the axes. Since the philosphy is that the two methods should give reasonable (not perfect) plots, I think I should add the figsize=(..., ...) there, too

• If a user wants to give his own axes, he should give at least two axes, the program will tell him how many are needed.

[Peter230655]

I added figsize(...) to plot_trajectories, so it gives nice pictures without the user giving axes. (I see no way to get by without this)
I changed the description on how the number of axes is determined in plot_constraint_violations a little bit. With the way I determine the number of axis, I see no simple way to tell the user how many are needed.

[Peter230655]

h = solution[-1], correct? I hope I can change this,

Yes, but only if the problem is variable time. This is not necessary in this PR, but could be done if you like (or in a later PR).

I may have found a solution to print the correct time with variable h.
I do it in a different branch. ( I am much more sure I found a solution for variable h than I am about this different branch...)

[moorepants]

One possible way around it is not to use the layout='compressed' flag. That flag breaks down if there are too many subplots in the default size of the figure.

Let's keep the automated figsize and the layout=compressed. I typically do not muck with figsize in library code due to figures becoming too large, but I think this is as good as we can get now.

[moorepants]

Thanks!

[Peter230655]

Thanks!

Thank YOU! I bet you spent more time coaching me compared to having it done yourself! :-)

A general question:
This is now part of opty 1.4.0dev0 ?
I tested it the best I could, but if users find problems I will learn about it so I can correct them?
Now that you merged it do I delete this branch?

[moorepants]

Thank YOU! I bet you spent more time coaching me compared to having it done yourself! :-)

I'm a teacher by trade, so that's just the way it always is :)

This is now part of opty 1.4.0dev0 ?

Yes

I tested it the best I could, but if users find problems I will learn about it so I can correct them?

They will open issues on this repository.

Now that you merged it do I delete this branch?

Yes, you can safely delete, but make sure to pull the master from this repository in your local master to get the updates.

[Peter230655]

Thank YOU! I bet you spent more time coaching me compared to having it done yourself! :-)

I'm a teacher by trade, so that's just the way it always is :)

This is now part of opty 1.4.0dev0 ?

Yes

I tested it the best I could, but if users find problems I will learn about it so I can correct them?

They will open issues on this repository.

Now that you merged it do I delete this branch?

Yes, you can safely delete, but make sure to pull the master from this repository in your local master to get the updates.

Yeah, part of your job! :-)
I may have found a way to give the time with variable h. Seems to work so far, but still not all exceptions caught. I will make a PR when I feel much more confident.

[Peter230655]

I think, I found a pretty good way to print the rounded number with variable h. I simply look around the symbols tree (I guess that is what it is called). It works for many ways of giving the duration, like
duration = (num_nodes - 1) * h
duration = h * (num_nodes - 1)
duration = (num_nodes - 1) * h / 2.0

It does NOT work for
duration = (num_nodes - 1) * h / 2 In this case it prints 125h/2 if num_nodes is 126

The variable inteval may have any name, not just h

I'd like to make a PR to see if it passes all the test (I think it would..) but I am afraid, I am swamping you with PRs!

[moorepants]

The variable inteval may have any name, not just h

h is stored in Problem, so you'll always have the symbol inside problem.

[Peter230655]

The variable inteval may have any name, not just h

h is stored in Problem, so you'll always have the symbol inside problem.

At the top of class Problem, I set, in line 155,
self.zeit = node_time_interval (I used zeit = German for time, to ensure I would not mess up anything, as presumably you not use German names )

[moorepants]

self.zeit = node_time_interval (I used zeit = German for time, to ensure I would not mess up anything, as presumably you not use German names )

presumably not, but you are teaching me new words every day :)

[Peter230655]

self.zeit = node_time_interval (I used zeit = German for time, to ensure I would not mess up anything, as presumably you not use German names )

presumably not, but you are teaching me new words every day :)

:-) :-) German is very similar to Dutch.
It seems to work. Should I make a PR so you can look at it at your convenience?

[moorepants]

Just a note for future, in the next sympy version you can do:

StrPrinter(settings={"dps": 3}).doprint(Function('f')(1.329294))

which returns

f(1.33)

[Peter230655]

Just a note for future, in the next sympy version you can do:

StrPrinter(settings={"dps": 3}).doprint(Function('f')(1.329294))

which returns

f(1.33)

Thanks!
I did not know this at all, never seen StrPrinter.
This works with the 'development version' of sympy?

[moorepants]

Yes, we just merged the change.

[moorepants]

I did not know this at all, never seen StrPrinter.

Both Timo and I showed you this above at least twice.

[Peter230655]

I did not know this at all, never seen StrPrinter.

Both Timo and I showed you this above at least twice.

Seems I am getting old! :-(