New definition of complete_state for ProjectedDynamicalSystems

[awage]

When defining a ProjectedDynamicalSystem, the user must define a function complete_state that must returns a vector of the same dimension than the original dynamical system. Such that:

  ds = Systems.lorenz96(5)
  projection(u) = [sum(u), sqrt(u[1]^2 + u[2]^2)]
  complete_state(y) = repeat([y[1]/5], 5) 
  prods = ProjectedDynamicalSystem(ds, projection, complete_state)

The docstring says:

complete_state produces the state for the original system from
the projected state. complete_state can always be a function
that given the projected state returns a state in the original
space.

In the previous example this is not true, since the projected space is a 2D space and the complete_state function only takes the first dimension of the input vector y.

In fact the complete_state function can be a crazy input y, the only requisite is to return a valid vector on the original phase space.

I propose a minor change in the docstring and small modification of the source code:

function ProjectedDynamicalSystem(ds::DynamicalSystem, projection, complete_state)
    u0 = initial_state(ds)
    if projection isa AbstractVector{Int}
        all(1 .≤ projection .≤ dimension(ds)) || error("Dim. of projection must be in 1 ≤ np ≤ dimension(ds)")
        projection = SVector(projection...)
        y = u0[projection]
    else
        projection(u0) isa AbstractVector || error("Projected state must be an AbstracVector")
        y = projection(u0)
    end
    if complete_state isa AbstractVector
        projection isa AbstractVector{Int} || error("Projection vector must be of type Int")
        length(complete_state) + length(projection) == dimension(ds) ||
                                error("Wrong dimensions for complete_state and projection")
        remidxs = setdiff(1:dimension(ds), projection)
        !isempty(remidxs) || error("Error with the indices of the projection")
    else
        # Define a wrapper function to check the state since we do not know what will be the input y
        comp_state(y) = length(complete_state(y)) == dimension(ds) ? complete_state(y) :  error("The returned vector of complete_state must equal dimension(ds)")
        remidxs = nothing
    end
    u = zeros(dimension(ds))
    return ProjectedDynamicalSystem{
        typeof(projection), length(y), typeof(complete_state),
        typeof(remidxs), typeof(ds)}(projection, comp_state, u, remidxs, ds)
end

The change introduces a small overhead in the computation of the initial state. Maybe there is a better solution. Maybe a special error definition if the case arises that complete_state(y) != dimension(ds)

[Datseris]

where did you see this docstring? please provide versions.

The docstring i read says something diferent:

complete_state produces the state for the original system from the projected state. complete_state can always be a function that given the projected state returns a state in the original space. However, if projection isa AbstractVector{Int}, then complete_state can also be a vector that contains the values of the remaining variables of the system, i.e., those not contained in the projected space. In this case the projected space needs to be lower-dimensional than the original.

from https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystemsbase/stable/#DynamicalSystemsBase.ProjectedDynamicalSystem

[Datseris]

We should add a check in the source code that complete_state returns the same dimension as ds At the moment i don't check anything i think.

[awage]

This is the same docstring. I just quoted the first two sentences. We are on the same page.

What I think should be changed are these two sentences. If N = dimension(ds), the the function complete_state should be a function from R^M to R^N, where M can be any integer.

There is a check in the original source code: line 76 in projected_systems.jl:

 length(complete_state(y)) == dimension(ds) ||
                        error("The returned vector of complete_state must equal dimension(ds)")
        remidxs = nothing

The problem is that y is of the same dimension a the projected state. So the check fails if in the code the vector y is of dimension. lets says 2N for example. But the only true requisite is that the function returns a vector of dim N.

[Datseris]

the requirement is that the function takes a vector from the projected state and returns a vector in the original space. This is also what the docstring says. So M cannot be "any" integer, it has to be the dimension of the projected state. Feel free to make a PR that makes it clearer, but for me it is clear the way it is!

[Datseris]

to my knowledge there is no way to test, nor enforce, what is the input to a function. We can't check whether the user provided function complete_state "expects" an input of dimension M. We just hope it does because that what we document.

[awage]

The docstring is clear (I changed the name of the issue). My point is that we do not have to restrict the input the dimension of the projected space.

I have an example of mapper that tracks states on 2D space. The dynamical systems is N dimensional (coupled Kuramotos). I need to define the initial state of the mapper on the original space, not on the projected space. The attractors can be defiined on this 2D space so there is no problem.

Right now it is not possible because of the sanity check of ProjectedDynamicalSystem.

I will draft a small PR so can understand better what I mean.

[awage]

I have a system of n kuramoto oscillators described in DOI:10.1103/PhysRevLett.127.194101.

The stable states of this system are called twisted states: $\theta_i = 2\pi i q/n + C$ with $q$ the winding number of the state and $C$ a constant. You can notice that this is a linear equation $f(i) = a i + C$ with $a = 2\pi q/n$. Once you have converged to the stable state you can extract the winding number $q$ with a linear fit of the phases.

In order to track these states I set up a specific projection onto the $q$, $C$ values:

function proj_fun(y)
    @. model(x, p) = x*p[1]+p[2]
    N = length(y)     
    p0 = [0.5, 0.5]
    xdata = 1:N
    unwrap!(y)
   # this is the lsq fit:
    fit = curve_fit(model, xdata, y, p0)
    p = fit.param
    p[1] = p[1]/(2π/N)
    p[2] = 0.
    return p
end

So I have a projection function proj_fun from $\mathbb{R}^n \to \mathbb{R}^2$. The complete_state function is from $\mathbb{R}^2 \to \mathbb{R}^n$, but I have no idea what should be this function or even if it exists. I would need the property complete_state(proj_fun(u)) = u and I guess it is impossible in most cases.

Now, suppose I want to explore the fractions or basins of a subspace of $\mathbb{R}^n$, if I set up

u = proj_fun(u0)
l = mapper(u)

I am doomed because I don't now the if the dynamical system starts at the state u0 I am interested in.

This is an unusual case because I need to feed directly the initial conditions to the original dynamical system, not to the projection. This is why I proposed in the PR to be less restrictive about the complete_state function such that it can be any arbitrary function from $\mathbb{R}^p \to \mathbb{R}^n$.

Maybe I am missing something, but yet I don't know how to achieve this using the current API. Of course I can try to track the attractors directly in the original phase space $\mathbb{R}^n$ but it doesn't work well with the Kuramoto oscillators in general. This projection is very effective.

Nevertheless I could compute a basin on a two dimensional subspace using a trick, these basins are very cool:

[Datseris]

Why don't you use the featurizer approach...? Seems like exactly the type of case featurizing was developed for. The features will be C and q

function featurize(A, t)
u = A[end]
q, C = curve_fit(...)
return SVector(q, C)
end

[Datseris]

so you don't use a projected system at all.

[Datseris]

So I have a projection function proj_fun from $\mathbb{R}^n \to \mathbb{R}^2$. The complete_state function is from $\mathbb{R}^2 \to \mathbb{R}^n$, but I have no idea what should be this function or even if it exists. I would need the property complete_state(proj_fun(u)) = u and I guess it is impossible in most cases.

Well, I would say that this is a case where projected system is not something you can use.

Notice that for your particular example the complete state is actually simple: just make complete_state(y) = (q, C = y; [2*pi*q*i/n + C for i in 1:n]). I assume n is the number of oscillators?

[awage]

Why don't you use the featurizer approach...? Seems like exactly the type of case featurizing was developed for. The features will be C and q

When you have a hammer everything looks like a nail. I am always using the recurrence method even if it doesn't apply. I will try the featurizer and see how well it does for this system.

Thanks!