more general `bymap` (RFC)

[aplavin]

I drafted a reasonably simple and more generic bymap version that works for arbitrary inputs that contain Particles inside.
It assumes (and asserts) all Particles have the same number of particles and iterates over them.
In the end, it gathers all results from applying f to arguments stripped of Particles, and transforms all numeric arguments there into Particles, keeping non-numbers from the first output. Assumes that all f results are consistent, have the same structure.

import MonteCarloMeasurements as MCM
using AccessorsExtra

function bymap₊(f, args...)
	Ns = @getall args |> RecursiveOfType(MCM.AbstractParticles) |> length(_.particles)
	N = uniqueonly(Ns)
	vals = map(1:N) do i
		curargs = @modify(args |> RecursiveOfType(MCM.AbstractParticles)) do p
			p.particles[i]
		end
		f(curargs...)
	end
	v = first(vals)
	numoptics = AccessorsExtra.flat_concatoptic(v, RecursiveOfType(Number))
	valps = map(o -> MCM.Particles(o.(vals)), AccessorsExtra._optics(numoptics))
	return setall(v, numoptics, valps)
end

x = (MCM.Particles(randn(1024)), (a=MCM.Particles(randn(1024)), b=123))
bymap₊(x) do x
	vals = (x[1], x[2].a)
	(;x[2].b, i=argmax(vals))
end
# returns:
# (b = 123.0 ± 0.0, i = 1.52 ± 0.5)

I'm not proposing to include it here, but wonder if you have some suggestions/recommendations on the specific semantics given your experience with bymap/Workspace. This bymap₊ is supposed to be much easier to use than the workspace thing, and applicable in many cases.

[baggepinnen]

The most difficult thing to support I've found is when particles are located in the fields of structs, the Workspace object is intended to handle this case, but does so in a somewhat non-robust fashion. Does this approach using Accessors handle this use case? If it does, we could perhaps make it the default way of handling difficult cases.

Where is uniqueonly defined? The example does not run for me

[baggepinnen]

It does handle this example, which is quite impressive :)

struct Bar{T1, T2}
    a::T1
    b::T2
end

x = Bar(MCM.Particles(randn(1024)), 0)

y = bymap₊(x) do x
    val = if x.a > 0
	    x.a + x.b
    else
        0
    end
    Bar(val, val)
end

using Test
@test all(>=(0), y.a.particles)
@test y.a == y.b

but it fails if the fields of Bar are untyped

[aplavin]

Where is uniqueonly defined? The example does not run for me

Sorry, I'm so used to DataManipulation.jl being loaded "almost always", and forgot to check if the code works in a fresh session. It's just the efficient version of uniqueonly(x) = only(unique(x)).

but it fails if the fields of Bar are untyped

Yes, it does require everything to be type-stable. With untyped fields, the loop/map would be slow anyway – f() uses and creates a Bar object N times.

Does this approach using Accessors handle this use case? If it does, we could perhaps make it the default way of handling difficult cases.
...
It does handle this example, which is quite impressive :)

The Accessors approach to various problems has proven to be general and widely applicable :) This case is no exception.
Accessors.jl itself is a very stable package with the intent to avoid any breaking changes and remain foundational. That's why more opinionated functions/those that requite Julia internals/where the details of the semantics is not 100% established go to AccessorsExtra.jl.
Earlier, we discussed upstreaming stuff like this recursive optic to Accessors, but that's quite challenging for the full functionality including setall and flattening. AccessorsExtra is also well-tested and doesn't make major breaking changes – but I admit it's quite a heavy dependency. I personally use it liberally, and recommend for endusers – but lightweight packages probably don't want to depend on it.

This raises the question where to put such a bymap₊ function. I'm glad that you like the semantics and it works for you as well :) Maybe, add it here as an experimental extension? Like, define the function itself in MonteCarloMeasurements.jl, but implement in MCM -> AccessorsExtraExt.jl. Then, users can just load AccessorsExtra and use bymap₊...

[baggepinnen]

Maybe, add it here as an experimental extension?

Yeah, that sounds like a good option! (preferably without the unusual ₊ character)

[BambOoxX]

Hi ! I have just stumbled on a case where bymap currently fails. Might be a good test case as well.

using LinearAlgebra, MonteCarloMeasurements
H = (10± 0.1) .*rand(ComplexF64,10,5)
bymap(x -> qr(x).R, H)

ERROR: ArgumentError: tuple must be non-empty
Stacktrace:
[1] first(::Tuple{})
@ Base .\tuple.jl:196
[2] bymap(f::var"#39#40", args::Matrix{Complex{Particles{Float64, 2000}}})
@ MonteCarloMeasurements D:\bamboo.julia\packages\MonteCarloMeasurements\hHq6h\src\bymap.jl:67
[3] top-level scope
@ REPL[111]:1

bymap(x-> qr(x).R, real(H))

works fine

[baggepinnen]

Yeah that would e a nice addition to the test cases. In the meantime, you can try

fun = x -> qr(x).R
MonteCarloMeasurements.ℂⁿ2ℂⁿ_function(fun, H)

[BambOoxX]

@baggepinnen Ha ! Just in time for #167 !

[aplavin]

As is, bymap₊ from this thread doesn't support this case – although it would probably work if that was a StaticMatrix.
It's possible to extend it to this case, especially if we don't care about efficiency and extraneous allocations. But to be reasonably efficient while staying general, the transformation from "list of matrices" to "matrix of particles with lists" should be carefully thought out.

[aplavin]

resolved by #173!