A (slightly more complex) PseudoPositioner problem case

[codedump]

Hello,

I'm looking for some pointers and hints as to how to solve a problem I have with Ophyd. If this is the wrong forum to ask, please advise as to where to turn to.

I'm implementing a simple Bluesky setup for an experiment which involves a signal generator. In a nutshell, every measurement point requires that a unique function is generated and uploaded to the signal generator. The "unique functions" are generated by parameters. For the sake of simplicity, let's assume that the function (signal) has two square peaks and it looks like this:

                  +-------------+
                  |             |
                  |             |
------------------+      A      +-----------+     B      +----------->
                                            |            |
                                            |            |
                                            +------------+

There are two peaks there, labeled "A" and "B".

Now I want the RunEngine in Bluesky to iterate over, say, different starting points, amplitudes, and/or widths of any of the peaks. Possibly only one. Let's say "A".

The underlying hardware has an EPICS-IOC which offers the basic capabilites one would expect (related to "management", e.g. triggering, start/stop, frequency and scaling of the signal in its entirety), as well as the ability to upload an arbitrary waveform shape for the signal generator to use.

There are obviously several ways of how one could be doing this, including modifying the IOC to actually generate the waveform on the fly and export parameters like FUNC:PEAKA:START, FUNC:PEAKA:WIDTH, ... FUNC:PEAKB:WIDTH... as EPICS variables for Ophyd to use.

But this would take away flexibility: needing to operate a different shape of function (say, with 3 peaks instead of 2, or with Gaussian peaks instead of squares) would require modification of the IOC, or at least restarts. Ugly.

My best idea so far was to implement all of this in Ophyd. The resulting API should look something like this:

# General function-generator control
class MyFunctionGenertor(Device):
    trigger = Component(...)
    freq = Component(...)
    ...

class Peak(MagicSubclassOfDevice): ...

# Specific shape of a square-peak
class SquarePeak(Peak):
    start = Component(...)
    duration = Component(...)
    amplitude = Component(...)

# ...or different type of peaks
class GaussPeak(Peak):
    fwhm = Component(...)
    ...

# Specific User-defined function
class MySpecificFunction(SubclassOfOphydDevice):
    peak_a = Component(SquarePeak, ...)
    peak_b = Component(SquarePeak, ...)

And then to use it like this:

detectors = [...]
dev = MyFunctionGenertor(...)
func = MySpecificFunction(...)

RE = RunEngine({})

RE( scan(
         [detectors],
         func.peak_a.amplitude, 0.0, 1.0, 10,
         func.peak_a.duration, 0.1, 0.5, 10, ...
) )

So this way I'd effectively use the peak parameters as motors. I don't care much for readback values at this point.

The essential part would be that the user can easily build their own function from building blocks like SquarePeak and GaussPeak:

class AnotherFunction(...):
    a = Component(SquarePeak, ...)
    b = Component(GaussPeak, ...)
    c = Component(SquarePeask, ...)

func2 = AnotherFunction(...)

RE( scan([...], func2.b.fwhm, 0, 1, 10, ... ))

My first attempt to do this was to make use of the PseudoPositioner class, along the lines of:

class SquarePeak(PseudoPositioner):
    start = Component(PseudoSingle)
    dur = Component(PseudoSingle)
    ampl = Component(PseudoSingle)

    single_peak_array = Component(SoftPositioner)

    @pseudo_position_argument
    def forward(self, pos):
        a = numpy.zeros(signal_length)[start:start+dur]=ampl
        return self.RealPosition(single_peak_array=a)

    @real_position_argument
    def inverse(self, pos):
        # This doesn't need to make sense. It not generally possible to reconstruct
        # original parameters from an array
        return self.PseudoPosition(start=0, duration=0, ampl=0)

The results are... mixed. It kind of does what it's supposed to in the sense that I can access fields:

In [1]: peak.start.set(3.14)
Out [1]: MoveStatus(done=False, pos=peak, elapsed=0.0, success=False, settle_time=0.0)

...and if I access other fields subsequently (start, dur, ampl) the results are as expected. The ugly side is that I'm forced to provide an inverse() even if it doesn't make sense. The inverse() can even be nonsensical, like here, and it still does the right ting, because apparently the forward() transformation receives all its parameters not from the current inverse(), but from previous setpoints. Still, this is a bit ugly.

But then it gets worse.

The actual data that gets written to single_peak_array is not an array at all, it's apparently only a 0 -- presumably the first array element. If I replace the definition by this, i.e. add the init_pos parameter:

class SquarePeak(PseudoPositioner):
    ...
    single_peak_array = Component(SoftPositioner, init_pos=np.zeros(25))
    ...

Then things fall apart. First, it's even uglier to need to know the size of the array at writing-time. If anything, I should specify it when defining something very specific, like MySpecificFunction, not when defining a general component like SquarePeak.

But even if we ignore that for a moment, what happens is this:

In [6]: s = SquarePeak(name="peak")                                                                                                                                                          
peak_partial_peak: < (array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),) {}                                                                                                       
peak_partial_peak: move (array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),) {'wait': False, 'moved_cb': None, 'timeout': None}                                                    
Subscription readback callback exception (PartialPeak(name='peak_partial_peak', parent='peak', settle_time=0.0, timeout=None, egu='', limits=(0, 0), source='computed'))                     
Traceback (most recent call last):                                                                                                                                                           
  File "/var/home/florin/.local/lib/python3.12/site-packages/ophyd/ophydobj.py", line 492, in inner                                                                                          
    cb(*args, **kwargs)                                                                                                                                                                      
  File "/var/home/florin/.local/lib/python3.12/site-packages/ophyd/pseudopos.py", line 734, in _real_pos_update                                                                              
    self._update_position()                                                                                                                                                                  
  File "/var/home/florin/.local/lib/python3.12/site-packages/ophyd/pseudopos.py", line 714, in _update_position                                                                              
    if None in real_cur_pos:                                                                                                                                                                 
       ^^^^^^^^^^^^^^^^^^^^                                                                                                                                                                  
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() 

So apparently in pseudopos.py somewhere there's an assumption that the position can never be an array. (Right?)

At this point I'm already wondering am I just "holding it wrong"?

And then, even if I overcome these obstacles, I don't know how to continue. Once I have the SquarePeak working, with its partial function / sampling array, I need to bundle an arbitrary number of those into an actual MyFunction object. The obvious PseudoPositioner approach will fail, because SquarePeak is not a PseudoSingle

class MyFunction(PseudoPositionerSubclass):
    peak_a = Component(SquarePeak) # remember that this is a PseudoPositioner itself, not a PseudoSingle!
    peak_b = Component(SquarePeak)

    # This is where we eventually want to end up...
    device_function_upload = Component(EpicsSignal, ...)

So.

I can't shake the feeling that I'm awfully off-course here, but I have no idea where to go.

Do I:

• just need to find "the right documentation"? (I've read pretty much everything Google turned up, but it's very limited "official" Ophyd docs only)
• scrap through the Ophyd sources and try to puzzle it together until it makes sense?
• forget PseudoPositioner, but reimplement what feels like 80-90% of parallel functionality, with the particular caveats in mind for my problem?
• forget implementing this as an Ophyd device? (What else? How to I make use of the Ophyd-Device-ness when feeding this to RunEngine?)

The core of the issue seems to be that I need to access sub-components that go more than 1level deep (2 or 3, e.g. func.peak.param ...), but which need to directly influence the top-level data (func), and the data is an array.

Any pointers?

Thanks & Cheers, F.

[codedump]

The core of the issue seems to be that I need to access sub-components that go more than 1level deep (2 or 3, e.g. func.peak.param ...), but which need to directly influence the top-level data (func), and the data is an array.

...oh, and one particular aspect is that I need to pass down specify information like array lengths, which is available to the "user", at "terminal time", down to intermediate components; and the current Ophyd API seems all about specifying everything at "library editing time", i.e. when defining those intermediate components (which I don't want my users to have to ever touch).

[codedump]

On further thought, is there a sanctioned way to programatically add positioners to Device (as opposed to declaratively)? Something along the lines of:

class Function(Device):
    def addSquare(self, start, duration, amplitude): ...

pulse = Function(...)
pulse.addSquare("a", ...)
pulse.addSquare("b", ...)

RE( scan([...], pulse.a_amplitude, ..., pulse.a_duration, ...) )

Because technically, my intermediate namespaces from above (i.e. the "peaks") aren't any kind of components or positioners that will ever be used independently. It's really just the individual peak parameters that I'm ever going to want to modify. But I'll want my users to be able to define a number of functions with wildly different number of peaks. Having them exist flatly, directly as children of the Function device is actually closer to the reality of my problem.

The preferred way of adding peaks directly to the definition of Function by use of Component(...) surely seems elegant, but is overengineered and actually hampering for my particular use case.