[ExpandingMan]

This adds a more generic scaling operator that accepts arbitrary functions and fields as arguments. This obviates the need to materialize new field arrays for more complicated scaling operations.

It's hard to guarantee the properties of this operator because it's so generic, so for now I impose that it's Hermitian (though we must somewhat trust the user to obey this).

I kept the original ScalingOp both because it can be non-Hermitian and because the code of the new op is complicated enough to contain foot-guns and I think it's good to keep a super simple op.

Not at all attached to the awkward name GHScalingOp (Generic Hermitian).

[kylebeggs]

I think we should try to merge these types. I like what we have here and I think we should make this the new ScalingOp and then we can create aliases for special sub cases. maybe something like

struct ScalingOp{F,A<:Tuple} <: AbstractOperator
    func::F
    args::A
end
ScalingOp(c) = ScalingOp(identity, (c,))          # old one-arg constructor keeps working
const CoeffScaling{C} = ScalingOp{typeof(identity),Tuple{C}}

or whatever to name the simple scaling that ScalarOp is now

[ExpandingMan]

Alright, I had claude do this quick, it's a pretty simple thing it's basically just the same code with the CoeffScaling alias.

The one thing that's a little dubious is that now if ScalingOp is not CoeffScaling it automatically gets declared hermitian and it expects the user to obey that. The reason is that whether it's hermitian depends on the function body and that's not realistically something we can inspect. This is noted in the doc strings.

[kylebeggs]

The reason is that whether it's hermitian depends on the function body and that's not realistically something we can inspect. This is noted in the doc strings.

right, i wonder if instead we should declare it non-hermitian by default as thats the safer assumption (at least for correctness reasons) - I think there should be a way to do this via dispatch? We could also just actually check f is hermitian but maybe thats overkill and could be a significant computational cost if the grid is big enough.

I could be convinced to just trust the user and clearly document it. your thoughts?

[ExpandingMan]

I don't see how we reliably check if f is hermitian that seems complicated way beyond the scope of this package. I think realistically the only options are to either leave it as it is or add an option for the user to specify whether it's hermitian or not (i.e. so that the struct carries a boolean).

I suspect that hermitian is the more common use case, I have no clue what optimizations, if any, are in this package for hermitian operators so I thought it safer at least for my own immediate purposes to enforce that it is rather than that it is not.

[kylebeggs]

after conversing with claude some, I think since this is a diagonal op the adjoint is just conj of the diagonal, so we don't need to assume anything. I think we just need to check if f returns real or complex for isselfadjoint? Something like

isselfadjoint(S::ScalingOp) = Base.promote_op(S.func, map(_argtype, S.args)...) <: Real

adjoint_operator(S::ScalingOp) = isselfadjoint(S) ? S : ScalingOp(conj ∘ S.func, S.args)

function apply_adjoint!(y::Field, S::ScalingOp, x::Field, g::AbstractGrid, α, β)
    return apply!(y, adjoint_operator(S), x, g, α, β)
end

[ExpandingMan]

I don't think this is a good idea. For one, promote_op is internal and the docs explicitly warn that its use should be avoided. For another, the result in general depends on the field value. For example
$\phi_1 + i \phi_2$
with $\phi_1, \phi_2 \in \mathbb{R}$ is hermitian iff $\phi_2 = 0$ (github is doing a pretty shitty job typesetting today) which may be true in some regions but not others (granted the return type could not be Real but in any case our check would be a bit dishonest because it doesn't check Complex cases). And if you care about whether the operator is hermitian, there's a good chance you are returning Complex.

Oh, and adjoint_operator as you've written above is then type unstable... and to make matters worse it's type unstable depending on what promote_op and constant propagation decides to do, so it's very unpredictable.

[kylebeggs]

Right about promote_op and type instability.

Regarding the hermitian trait my main issue is hard coding hermitian to be true and putting it in the docs because all it takes is a user (or LLM) to not digest the docs properly then they would get silently wrong results. I'm generally OK with leaning on just documenting things but I think this could really bite a user and after they find out their output is silently wrong it will leave a really bad taste in their mouth even if it was their fault.

If you are to assume anything, I think hermitian = false is the safer choice, then we can add an opt in to override. We can make it a type param like this:

struct ScalingOp{F,A,SA}
    func::F
    args::A
    ScalingOp{SA}(func::F, args::A) where {F,A,SA} = new{F,A,SA::Bool}(func, args)
end

thoughts?

This whole conversation has made it painfully obvious we need a higher level refactor of the API because this is not as straightforward as I would have hoped...

[kylebeggs]

then we have

adjoint_operator(S::ScalingOp{<:Any,<:Any,false}) = ScalingOp(conj ∘ S.func, S.args)
adjoint_operator(S::ScalingOp{<:Any,<:Any,true}) = S
isselfadjoint(::ScalingOp{F,A,SA}) where {F,A,SA} = SA

function apply_adjoint!(y::Field, S::ScalingOp, x::Field, g::AbstractGrid, α, β)
    return apply!(y, adjoint_operator(S), x, g, α, β)
end

[kylebeggs]

maybe this is getting too convoluted, out of scope, but that goes back to my "we need to reconsider the higher level API here"