Readme update and trapezoidal integration of cost function #458

Ipuch · 2022-04-20T14:47:08Z

All Submissions:

Have you followed the guidelines in our Contributing document [docs/contribution.md]?
Have you checked to ensure there aren't other open [Pull Requests] for the same update/change?
Have you opened/linked the issue related to your pull request?
Have you used the tag [WIP] for on-going changes, and removed it when the pull request was ready?
When ready to merge, have you sent a comment pinging @pariterre in it?

New Feature Submissions:

Does your submission pass the tests (if not please explain why this is intended)?
Did you write a proper documentation (docstrings and ReadMe)
Have you linted your code locally prior to submission (using the command: black . -l120 --exclude "external/*")?

Changes to Core Features:

Have you added an explanation of what your changes do and why you'd like us to include them?
Have you written new examples for your core changes, as applicable?
Have you written new tests for your core changes, as applicable?

related to #416

This change is

first iteration

Ipuch · 2022-04-21T20:49:20Z

Still need to address some challenges:

Handle targets with trapezoidal integration
IPOPT seems to compute the objectives differently compared to sol.print_cost()
finish draft of the test.
if control are constant (ControlType.CONSTANT) then integration should be the same in trapezoidal and rectangle

Ipuch · 2022-04-25T21:59:03Z

still a weird bad file descriptor

… targets

Ipuch · 2022-04-26T22:21:09Z

@pariterre, I think it started to be something now. I have to finalize the tests for the targets for each use case. Let me know if I need to modify something in the core.

pariterre · 2022-04-27T12:48:42Z

@Ipuch please update your branch to master

pariterre

Reviewed 4 of 10 files at r1, 3 of 6 files at r2, 7 of 7 files at r3, 1 of 1 files at r4, all commit messages.
Reviewable status: all files reviewed, 21 unresolved discussions (waiting on @Ipuch)

README.md line 329 at r3 (raw file):

To solve the ocp, you simply have to call the `solve()` method of the `ocp` class
```python
solver=Solver.IPOPT(show_online_optim=True)

"r = S"

bioptim/__init__.py line 184 at r3 (raw file):

    VariableType,
    SolutionIntegrator,
    IntegralApproximation,

I feel that "Integral" is the value of an "Integration", but I may be wrong... Inuitively it makes more sense to me to call it "IntegrationApproximation" as it is the process of integration which is approximate

bioptim/examples/getting_started/pendulum.py line 111 at r3 (raw file):

    # Custom plots
    ocp.add_plot_penalty(CostType.OBJECTIVES)

Please do not change, apart for good reasons, the behavior of example. If you want to try something you can create a "sandbox" folder and it will automatically be ignore for your push

bioptim/examples/getting_started/pendulum.py line 117 at r3 (raw file):

    # --- Solve the ocp --- #
    sol = ocp.solve(Solver.IPOPT(show_online_optim=False))

Same here

bioptim/gui/plot.py line 62 at r3 (raw file):

        label: list = None,
        compute_derivative: bool = False,
        integration_rule: IntegralApproximation = IntegralApproximation.RECTANGLE,

Should trapezoidal be the default?

bioptim/interfaces/ipopt_interface.py line 214 at r3 (raw file):

    def __get_all_penalties(self, nlp, penalties):
        def format_target(target_in: list) -> np.array:

Thx for that :)

bioptim/interfaces/ipopt_interface.py line 294 at r3 (raw file):

                            target0 = format_target(penalty.target[0:1])
                            target1 = format_target(penalty.target[1:])
                            target = np.vstack((target0, target1)).T

This should be
target0 = penalty.target[0]
target1 = penalty.target[1]

Otherwise, what is the point of separating them in a list?

Code quote:

                            target0 = format_target(penalty.target[0:1])
                            target1 = format_target(penalty.target[1:])
                            target = np.vstack((target0, target1)).T

bioptim/limits/objective_functions.py line 260 at r3 (raw file):

            raise ValueError("'list_index' must be defined properly")

        ocp_or_nlp.J[list_index].target = [new_target]

What if the user already sent the second target (already a list)?

bioptim/limits/penalty.py line 85 at r3 (raw file):

            penalty.quadratic = True if penalty.quadratic is None else penalty.quadratic
            if penalty.integration_rule == IntegralApproximation.RECTANGLE:
                penalty.add_target_to_plot(all_pn=all_pn, combine_to=f"{key}_controls")

Please add a todo so we remember some day to implement when it is not rectangle

bioptim/limits/penalty_option.py line 169 at r3 (raw file):

            if len(self.target) == 1 and self.integration_rule == IntegralApproximation.TRAPEZOIDAL:
                if self.node == Node.ALL or self.node == Node.DEFAULT:
                    self.target = [self.target[0][:, :-1], self.target[0][:, 1:]]

See, you are doing it here already, that should not be done in IPOPT_INTERFACE

bioptim/limits/penalty_option.py line 173 at r3 (raw file):

                    raise NotImplementedError(
                        f"target with IntegralApproximation.TRAPEZOIDAL doesn't work with {self.node}, "
                        "it only works if node is Node.NODE_ALL and Node.NODE_DEFAULT"

That is not true, the "automatic" does not work, but one can provide the two columns and it will (should) work. Please rephrase ("Trap needs the beginning and end column values except for NODE_ALL which can be automatically generated" or something like that

Code quote:

                    raise NotImplementedError(
                        f"target with IntegralApproximation.TRAPEZOIDAL doesn't work with {self.node}, "
                        "it only works if node is Node.NODE_ALL and Node.NODE_DEFAULT"
                    )

bioptim/limits/penalty_option.py line 179 at r3 (raw file):

        self.target_to_plot = None
        # not implemented yet for trapezoidal integration
        self.plot_target = False if self.integration_rule == IntegralApproximation.TRAPEZOIDAL else True

Add a todo

bioptim/limits/penalty_option.py line 298 at r3 (raw file):

        elif self.integration_rule == IntegralApproximation.TRAPEZOIDAL:

            for i in range(2):

for target in self.target

As a rule of thumb, if you find yourself writting "for ... in range(...)" in Python, you are probably doing something wrong

bioptim/limits/penalty_option.py line 313 at r3 (raw file):

            )

            for i in range(2):

Same. At least len(target).. How do you know it is 2, maybe we implemented polynomial integration?

bioptim/limits/penalty_option.py line 325 at r3 (raw file):

                    and all_pn.nlp.ns in all_pn.t
                    and self.target[0].shape[-1] == all_pn.nlp.ns - 1
                    and self.target[1].shape[-1] == all_pn.nlp.ns - 1

Isn't it always true?
Again, how do you know there are exaclty 2 targets? If someone sends 3 or 4, what happens? I haven't see any control for that yet (maybe later)

Code quote:

                    and self.target[0].shape[-1] == all_pn.nlp.ns - 1
                    and self.target[1].shape[-1] == all_pn.nlp.ns - 1

bioptim/limits/penalty_option.py line 426 at r3 (raw file):

        dt_cx = nlp.cx.sym("dt", 1, 1)
        target_shape = list(sub_fcn.shape)

Why would you save that information?

bioptim/limits/penalty_option.py line 428 at r3 (raw file):

        target_shape = list(sub_fcn.shape)
        target_shape[1] += 1 if self.integration_rule == IntegralApproximation.TRAPEZOIDAL else 0
        target_cx = nlp.cx.sym("target", tuple(target_shape))

why not do
...("target", self.rows, self.cols + 1 if self.integration_rule == IntegralApproximation.TRAPEZOIDAL else 0)? which replaces the 5 lines before

Code quote:

        target_shape[1] += 1 if self.integration_rule == IntegralApproximation.TRAPEZOIDAL else 0
        target_cx = nlp.cx.sym("target", tuple(target_shape))

bioptim/limits/penalty_option.py line 430 at r3 (raw file):

        target_cx = nlp.cx.sym("target", tuple(target_shape))
        weight_cx = nlp.cx.sym("weight", 1, 1)
        n = 2 if self.quadratic and self.weight else 1

Please give better naming

bioptim/limits/penalty_option.py line 467 at r3 (raw file):

        self.weighted_function = Function(
            name, [state_cx, control_cx, param_cx, weight_cx, target_cx, dt_cx], [modified_fcn]
        ).expand()

please remove this.
expand only if the expand parameter is to True (see 10 lines later)

bioptim/misc/enums.py line 136 at r3 (raw file):

    """

    NONE = None

I am not sure None makes sense...

bioptim/optimization/optimal_control_program.py line 742 at r3 (raw file):

                    or penalty.transition
                ):
                    plot_params["plot_type"] = PlotType.POINT

Just to confirm, this was removed because it shows the wrong result?

Code quote:

                    isinstance(penalty.type, ObjectiveFcn.Mayer)
                    or isinstance(penalty.type, ConstraintFcn)
                    or penalty.transition
                ):
                    plot_params["plot_type"] = PlotType.POINT
                    plot_params["node_idx"] = penalty.node_idx
                else:
                    plot_params["plot_type"] = PlotType.INTEGRATED

Ipuch · 2022-04-27T14:32:31Z

README.md line 329 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

"r = S"

I didn't understand

Ipuch · 2022-04-27T14:34:19Z

bioptim/__init__.py line 184 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

I feel that "Integral" is the value of an "Integration", but I may be wrong... Inuitively it makes more sense to me to call it "IntegrationApproximation" as it is the process of integration which is approximate

IntegralApproximation is fine to me, What about NumericalIntegration ?

Ipuch · 2022-04-27T14:36:03Z

bioptim/examples/getting_started/pendulum.py line 117 at r3 (raw file):

Restored, Done.
Restored, Done.

Ipuch · 2022-04-27T18:40:51Z

bioptim/limits/penalty_option.py line 298 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

for target in self.target

As a rule of thumb, if you find yourself writting "for ... in range(...)" in Python, you are probably doing something wrong

yes, done.

Ipuch · 2022-04-27T18:50:08Z

bioptim/limits/penalty_option.py line 426 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

Why would you save that information?

It has been simplified to ease the readiness of the code.

Ipuch · 2022-04-27T18:51:15Z

bioptim/limits/penalty_option.py line 428 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

why not do
...("target", self.rows, self.cols + 1 if self.integration_rule == IntegralApproximation.TRAPEZOIDAL else 0)? which replaces the 5 lines before

Done, self.rows is an array of indexes. I needed the size.

Ipuch

ok done.

Reviewable status: 13 of 14 files reviewed, 21 unresolved discussions (waiting on @pariterre)

bioptim/interfaces/ipopt_interface.py line 294 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

This should be
target0 = penalty.target[0]
target1 = penalty.target[1]

Otherwise, what is the point of separating them in a list?

Done.

bioptim/limits/penalty_option.py line 169 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

See, you are doing it here already, that should not be done in IPOPT_INTERFACE

I think I resolved that too

bioptim/limits/penalty_option.py line 325 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

Isn't it always true?
Again, how do you know there are exaclty 2 targets? If someone sends 3 or 4, what happens? I haven't see any control for that yet (maybe later)

Now it is checked.

I guess i'm done with this PR.

pariterre

Reviewed 10 of 10 files at r6, all commit messages.
Reviewable status: all files reviewed, 7 unresolved discussions (waiting on @Ipuch)

bioptim/__init__.py line 184 at r3 (raw file):

Previously, Ipuch (Pierre Puchaud) wrote…

IntegralApproximation is fine to me, What about NumericalIntegration ?

Let's leave it as is :)

bioptim/gui/plot.py line 62 at r3 (raw file):

Previously, Ipuch (Pierre Puchaud) wrote…

If so it will change the value for every single test on an OCP and it will deprecate all the examples designed outside. But yes, let us try it as an option for now and when everybody agreed it is better, let's make the refactor.

To be done later

bioptim/limits/objective_functions.py line 260 at r3 (raw file):

Previously, Ipuch (Pierre Puchaud) wrote…

Idk. Can you explain more? I've never put my finger into MHE until yesterday.

This comment is not related to MHE, it is related to the user input

bioptim/limits/penalty_option.py line 298 at r3 (raw file):

Previously, Ipuch (Pierre Puchaud) wrote…

yes, done.

for target in self.target

Code quote:

range(len

bioptim/limits/penalty_option.py line 426 at r3 (raw file):

Previously, Ipuch (Pierre Puchaud) wrote…

It has been simplified to ease the readiness of the code.

Please use the "is_integrated" I talk later

bioptim/limits/penalty_option.py line 659 at r6 (raw file):

                all_pn.t[:-1]
                if (
                    self.integration_rule == IntegralApproximation.TRAPEZOIDAL

You find yourself doing this check over and over and it is confusing.
create a variable:
is_integrated = self.integration_rule == IntegralApproximation.TRAPEZOIDAL or self.integration_rule == IntegralApproximation.TRUE_TRAPEZOIDAL

and job done

Code quote:

                    self.integration_rule == IntegralApproximation.TRAPEZOIDAL
                    or self.integration_rule == IntegralApproximation.TRUE_TRAPEZOIDAL

tests/test_cost_function_integration.py line 148 at r6 (raw file):

    ocp = prepare_ocp(
        biorbd_model_path=bioptim_folder + "/models/pendulum.bioMod",
        n_shooting=30,

that seems a lot... 5? We don't care that the test actually converges... Ideally yes, but not at the expance of a test that is more than 5 to 10 seconds

tests/test_cost_function_integration.py line 381 at r6 (raw file):

    solver = Solver.IPOPT()
    solver.set_maximum_iterations(1)

I see, that was for to compensate the 30? Maybe 10 iterations so it leaves the beginning enough?

Ipuch

Reviewable status: all files reviewed, 7 unresolved discussions (waiting on @pariterre)

bioptim/__init__.py line 184 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

Let's leave it as is :)

ok

bioptim/gui/plot.py line 62 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

To be done later

ok

bioptim/limits/objective_functions.py line 260 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

This comment is not related to MHE, it is related to the user input

modified as,
[new_target] if not isinstance(new_target, list) else new_target

bioptim/limits/penalty_option.py line 298 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

for target in self.target

Thx, done.

bioptim/limits/penalty_option.py line 426 at r3 (raw file):

Previously, pariterre (Pariterre) wrote…

Please use the "is_integrated" I talk later

done, with is_trapezoidal

bioptim/limits/penalty_option.py line 659 at r6 (raw file):

Previously, pariterre (Pariterre) wrote…

You find yourself doing this check over and over and it is confusing.
create a variable:
is_integrated = self.integration_rule == IntegralApproximation.TRAPEZOIDAL or self.integration_rule == IntegralApproximation.TRUE_TRAPEZOIDAL

and job done

done.

tests/test_cost_function_integration.py line 148 at r6 (raw file):

Previously, pariterre (Pariterre) wrote…

that seems a lot... 5? We don't care that the test actually converges... Ideally yes, but not at the expance of a test that is more than 5 to 10 seconds

ok.

tests/test_cost_function_integration.py line 381 at r6 (raw file):

Previously, pariterre (Pariterre) wrote…

I see, that was for to compensate the 30? Maybe 10 iterations so it leaves the beginning enough?

ok. if you want, I picked 5 like the other one.

…e is used with targets

codecov · 2022-05-04T17:57:46Z

Codecov Report

Merging #458 (57dce8b) into master (17f40d7) will increase coverage by 0.07%.
The diff coverage is 87.41%.

@@            Coverage Diff             @@
##           master     #458      +/-   ##
==========================================
+ Coverage   80.48%   80.56%   +0.07%     
==========================================
  Files          85       85              
  Lines        8998     9080      +82     
==========================================
+ Hits         7242     7315      +73     
- Misses       1756     1765       +9

Impacted Files	Coverage Δ
bioptim/__init__.py	`100.00% <ø> (ø)`
bioptim/examples/getting_started/pendulum.py	`71.42% <ø> (ø)`
bioptim/gui/graph.py	`96.18% <ø> (+1.17%)`	⬆️
bioptim/optimization/optimal_control_program.py	`88.08% <72.72%> (-0.32%)`	⬇️
bioptim/limits/penalty_option.py	`89.51% <80.00%> (-2.68%)`	⬇️
bioptim/interfaces/ipopt_interface.py	`91.07% <95.83%> (+0.68%)`	⬆️
bioptim/gui/plot.py	`76.00% <100.00%> (ø)`
bioptim/interfaces/acados_interface.py	`93.44% <100.00%> (ø)`
bioptim/limits/objective_functions.py	`93.75% <100.00%> (+0.32%)`	⬆️
bioptim/limits/penalty.py	`93.41% <100.00%> (+0.05%)`	⬆️
... and 5 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 48e5b64...57dce8b. Read the comment docs.

pariterre

Reviewed 5 of 5 files at r7, all commit messages.
Reviewable status: all files reviewed, 4 unresolved discussions (waiting on @Ipuch and @pariterre)

bioptim/limits/objective_functions.py line 275 at r7 (raw file):

            raise ValueError("'list_index' must be defined properly")

        ocp_or_nlp.J[list_index].target = [new_target] if not isinstance(new_target, list) else new_target

Test for list and tuple

bioptim/limits/penalty_option.py line 316 at r7 (raw file):

                raise RuntimeError(f"targets with trapezoidal integration rule need to get a list of two elements.")

            for t in self.target:

please call "t" "target"

Ipuch

Reviewable status: 13 of 15 files reviewed, 4 unresolved discussions (waiting on @pariterre)

bioptim/limits/objective_functions.py line 275 at r7 (raw file):

Previously, pariterre (Pariterre) wrote…

Test for list and tuple

done

Code quote:

        ocp_or_nlp.J[list_index].target = [new_target] if not isinstance(new_target, list) else new_target

bioptim/limits/penalty_option.py line 316 at r7 (raw file):

Previously, pariterre (Pariterre) wrote…

please call "t" "target"

done

Code quote:

for t in self.target:

pariterre

Reviewed 2 of 2 files at r8, all commit messages.
Reviewable status: complete! all files reviewed, all discussions resolved (waiting on @Ipuch)

Ipuch added 8 commits April 19, 2022 11:24

fix a bug in linear continuous

a52c859

taking into account some other constraints to plot it

71c86a9

solver readme update

5291f49

not working yet

d8c01f4

Merge remote-tracking branch 'origin/master' into trapeze

2cd89a7

docstring and handling graph

d22244f

first iteration

plot objectives in plotstep is more rigorous

5037224

print of cost funtcion fixed

d7bc494

Ipuch changed the title ~~Readme update~~ Readme update and trapezoidal integration of cost function Apr 21, 2022

Ipuch mentioned this pull request Apr 21, 2022

Plot of contacts #363

Closed

Ipuch added 2 commits April 21, 2022 15:05

restored pendulum example

66040c8

blacked and found ipopt objectives values to be different ...

007cda7

Ipuch added 4 commits April 25, 2022 09:34

dot and line is more rigorous for objective functions

90b375c

piecewise constant control compatible with trapz integration

7d44e9f

linear continuous also works

db6266f

blacked

a7b19cf

Ipuch added 3 commits April 26, 2022 15:20

Handling targets in list instead of numpy array to prepare trapz with…

a9fa4f2

… targets

handling trapz with target in lagrange

d7e0133

fixing tests

60e4e14

Merge remote-tracking branch 'pyomeca/master'

e7f40ab

pariterre requested changes Apr 27, 2022

View reviewed changes

small fix on expanding penalty functions

3c09023

restored pendulum

815da29

Ipuch commented Apr 27, 2022

View reviewed changes

Ipuch added 6 commits April 27, 2022 15:17

Answering all the comments in reviewable.

b86e725

new tests for targets and blacked

5e1c882

new tests

72febc0

fix a print in a test

365943c

True trapezoidal integration

4109221

I guess i'm done with this PR.

Readme to go with it

4c8487e

pariterre changed the base branch from master to CrazyDynamics April 29, 2022 18:57

pariterre requested changes Apr 29, 2022

View reviewed changes

answering a first wave of comments.

b772444

Ipuch commented Apr 29, 2022

View reviewed changes

Ipuch changed the base branch from CrazyDynamics to master April 29, 2022 20:19

Ipuch added 7 commits May 3, 2022 17:52

small commit to handle bugs. tests still doesn't work.

9814f14

shall the test pass in windows ?

096b8c4

blacked

7b443cc

extra bug fixed which is the consequence of using list when trapz rul…

629727a

…e is used with targets

The test going with it

2bd6c11

to pass windows tests

4543818

blacked

e969c14

pariterre requested changes May 4, 2022

View reviewed changes

final trapz push

57dce8b

Ipuch commented May 4, 2022

View reviewed changes

pariterre approved these changes May 5, 2022

View reviewed changes

pariterre merged commit 6aff2c1 into pyomeca:master May 5, 2022

pariterre mentioned this pull request May 5, 2022

test the tests #462

Closed

Ipuch mentioned this pull request Apr 27, 2023

Lagrange objective should be calculated at each integration steps #249

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Readme update and trapezoidal integration of cost function #458

Readme update and trapezoidal integration of cost function #458

Ipuch commented Apr 20, 2022 •

edited

Ipuch commented Apr 21, 2022 •

edited

Ipuch commented Apr 25, 2022

Ipuch commented Apr 26, 2022

pariterre commented Apr 27, 2022

pariterre left a comment

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch left a comment

pariterre left a comment

Ipuch left a comment

codecov bot commented May 4, 2022 •

edited

pariterre left a comment

Ipuch left a comment

pariterre left a comment

Readme update and trapezoidal integration of cost function #458

Readme update and trapezoidal integration of cost function #458

Conversation

Ipuch commented Apr 20, 2022 • edited

All Submissions:

New Feature Submissions:

Changes to Core Features:

Ipuch commented Apr 21, 2022 • edited

Ipuch commented Apr 25, 2022

Ipuch commented Apr 26, 2022

pariterre commented Apr 27, 2022

pariterre left a comment

Choose a reason for hiding this comment

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch commented Apr 27, 2022

Ipuch left a comment

Choose a reason for hiding this comment

pariterre left a comment

Choose a reason for hiding this comment

Ipuch left a comment

Choose a reason for hiding this comment

codecov bot commented May 4, 2022 • edited

Codecov Report

pariterre left a comment

Choose a reason for hiding this comment

Ipuch left a comment

Choose a reason for hiding this comment

pariterre left a comment

Choose a reason for hiding this comment

Ipuch commented Apr 20, 2022 •

edited

Ipuch commented Apr 21, 2022 •

edited

codecov bot commented May 4, 2022 •

edited