JuliaControl / ControlSystems.jl

A Control Systems Toolbox for Julia
https://juliacontrol.github.io/ControlSystems.jl/stable/
Other
508 stars 85 forks source link

Charpoly computation improved #804

Closed stepanoslejsek closed 1 year ago

stepanoslejsek commented 1 year ago

I took a task for improving computation of characteristic polynomial as it was mentioned in the #TODO note using the provided article. Testing withing runtest.jl was successful.

baggepinnen commented 1 year ago

Hello and thanks for your contribution :)

Would you mind submitting a version of the PR without any automatic formatting changes applied?

stepanoslejsek commented 1 year ago

Yeah for sure, I apologize for that formatting stuff.

codecov[bot] commented 1 year ago

Codecov Report

Merging #804 (d0e0e58) into master (0d06b70) will not change coverage. The diff coverage is n/a.

:exclamation: Current head d0e0e58 differs from pull request most recent head 21e35b2. Consider uploading reports for the commit 21e35b2 to get more accurate results

@@           Coverage Diff           @@
##           master     #804   +/-   ##
=======================================
  Coverage   97.14%   97.14%           
=======================================
  Files           4        4           
  Lines         315      315           
=======================================
  Hits          306      306           
  Misses          9        9           

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

baggepinnen commented 1 year ago

No worries, thanks for updating! It seems like github has been thrown off a bit by some of your changes, it shows the entire file as affected by the diff. I think the PR needs a squash, are you familiar with how to perform one?

stepanoslejsek commented 1 year ago

the PR needs a

No, I'm not familiar with PR squashing, but I'm looking forward to learn about it.

baggepinnen commented 1 year ago

hehe, it's quite simple, but it can be come messy at times. It boils down to combining several different commits so that they appear as a single commit. It can be done from the command-line interface, but I usually use a tool for the purpose, in my case, I use "git lens" in vscode.

stepanoslejsek commented 1 year ago

Thank you for the explanation, I'm going to squash the PR then. Hopefully it won't mess up.

baggepinnen commented 1 year ago

Hmm, your squash worked (nice ;) but there's still a whole-file diff. Base on this article https://stackoverflow.com/questions/19593909/git-diff-sees-whole-file-as-changed-when-its-not it seems like your editor has changed the whole file, are you by any chance using windows?

baggepinnen commented 1 year ago

Yeah, running this command git diff --ignore-space-at-eol master shows that your editor has inserted a bunch of invisible characters at the end of each line (all the ^Ms).

image

I'm not a windows person so I'm not sure how to fix that :/ I suggest googling how to turn that behavior off in your editor

stepanoslejsek commented 1 year ago

Oh I see now, I'm using linux, let me fix those invisible characters and squash it once again.

stepanoslejsek commented 1 year ago

It's still there, not exactly sure why, I'll try something else to deal with it.

baggepinnen commented 1 year ago

I think that your editor might be trying too hard to force its settings upon the file by auto formatting. I saw that you had one commit that did away with the line-ending edits, but it had instead inserted tabs instead of spaces for indentation 😃 Generally, these features are best left turned off when working on open-source packages, unless the package happens to follow exactly the same convention that your editor happens to be configured with.

Please don't let this discourage you, I do appreciate your PR and I'm sure you'll find the proper workflow soon!

mzaffalon commented 1 year ago

@stepanoslejsek What editor are you using? I found that Emacs + julia-mode + magit usually work, although I have never tried on Linux.

stepanoslejsek commented 1 year ago

I am using neovim, but I found the solution. It'll be without those invisible characters. I checked it with cat -A command.

stepanoslejsek commented 1 year ago

I did a benchmarking on BigFloat matrix using the old and new method and here is the result of benchmarking using @btime image

baggepinnen commented 1 year ago

Wonderful, that's a solid 10x improvement :) I think you might have misunderstood slightly what I meant with using BigFloat though, BigFloats are for comparing accuracy, while performance benchmarks are likely better run using the standard Float64. When you compute using BigFloats, you get an extremely accurate answer which can be considered "the truth", the old and new algorithms are then computed in Float64 and compared to the BigFloat result to see which one is closer to the truth.

stepanoslejsek commented 1 year ago

I see, sorry for misunderstanding, I took a random 10x10 Float64 matrix and use the old and new algorithm and pass trough a BigFloat representation and Float64 representation and then took a L2 norm on each coefficient and here is the result.

image

baggepinnen commented 1 year ago

That looks amazing :) Thanks for your nice contribution!

baggepinnen commented 1 year ago

May I ask how you performed the BigFloat computations? When I try to use any floating-point type that is not recognized by BLAS I get

julia> charpoly(big.(A))
ERROR: MethodError: no method matching hessenberg!(::Matrix{BigFloat})
stepanoslejsek commented 1 year ago

May I ask how you performed the BigFloat computations? When I try to use any floating-point type that is not recognized by BLAS I get

julia> charpoly(big.(A))
ERROR: MethodError: no method matching hessenberg!(::Matrix{BigFloat})

I ran into same problem with only using ControlSystemsBase, however I found out that in order to work, it is necessary to write using OrdinaryDiffEq, which is somewhat weird. I tested it with only using ControlSystems and then call ControlSystemsBase.charpoly(big.(A)).

baggepinnen commented 1 year ago

Thanks for the hint, OrdinaryDiffEq loads GenericSchur.jl which contains an implementation of hessenberg that works for BigFloat. I tried with GenericLinearAlgebra which has much poorer compatibility than GenericSchur, but now I know what to use instead :)