Optimized shooting method

[JonasKoziorek]

I implemented recent paper that searches for UPOs via optimization (link). The paper also includes python implementation in the appendix.

My implementation has a few potential problems, that should be addressed:

• I want to avoid creating ds in every iteration

  PeriodicOrbits.jl/src/algorithms/optimized_shooting.jl

  Lines 94 to 95 in 507234f

  	f = dynamic_rule(ds)
  	nds = CoupledODEs(new_rule(f, T), current_state(ds), current_parameters(ds); diffeq = ds.diffeq)

  This should be possible by simply not integrating from 0 to 1 but from 0 to T. However I couldn't make it to work for some reason.
• I would like to use other solvers than RKO65

  PeriodicOrbits.jl/src/algorithms/optimized_shooting.jl

  Lines 50 to 51 in 507234f

  	alg = OptimizedShooting(Δt=1/(2^6), p=3)
  	ds = CoupledODEs(dynamic_rule, state, params; diffeq = (alg=RKO65(), dt=alg.Δt))

  but all other solvers I tried simply weren't working. In the paper they use simple explicit 5-th order runge kutta. I think that choice of solver might be the cause of problems in Saiki too.
• I know that DynamicalSystemsBase.jl has OrdinaryDiffEq.jl as a dependency but the docs here say that to access the solvers, it has to be imported by the user. I wanted to change the default solver so I added the OrdinaryDiffEq.jl as a dependency again. I don't know if this is optimal. Is there a way to access the various solvers internally through DynamicalSystemsBase.jl?

[Datseris]

cc also @rveltz : if you know this algorithm, let's make sure that we aren't duplicating effort. If something like this is in BFKit, or parts of the code have been tried in BFKit, let's reuse the knowledge.

[JonasKoziorek]

Everything works as we discussed today. true in step was causing those problems.

[Datseris]

perfect!

[Datseris]

You can optimize the code. If you don't do it until our next meeting, then I will show you how to use profiling, and we can see the big difference in profiling before and after the optimizations.

[Datseris]

Suggested change

	to find periodic orbits of continuous systems.
	to find periodic orbits of continuous time systems.

[Datseris]

Hm, isn't this here too ineficient? Every time you call the cost function you are creating R1, R2, err. Does the optimize function run in parallel? If not, we can re-use existing vectors R1, R2, err instead of creating them all the time. For example, you can make costfunc be a type and use function-like objects. Something like

struct CostFunc
R1::Vector
R2::Vector
err:Vector
end

function (cf::CostFunc)(v, ds, alg)
   (; R1, R2, err) = cf
   # same code as in your function
end

[Datseris]

Same story. You need to make this compute_residual!(R, ...) that uses an existing array R instead of creating a new one every time.

[JonasKoziorek]

I tried to make it faster but with no success. Even though I reduced allocations almost in half, the runtime wasn't faster. Let's try profiling together next time then!

[JonasKoziorek]

The minimal_period was messing up one of the tests so I disabled it for now. I will allow passing kwargs to minimal_period in a different PR to be able to affect the precision.

[JonasKoziorek]

I suddenly can't build the docs. I attach the stacktrace.txt here. It originates in build_docs_with_style.jl line 30. Commenting this line solves the problem. Could you help me with this @Datseris?

[Datseris]

docbuild problem comes from Documenter.jl. We can't resolve it :( :( :( cc @Tortar , I searched throught the JuliaDynamics/doctheme repo, we don't have $todo anywhere. It has to be from Documenter.jl.