Script to Apply Pose as Rest Pose, without corrupting animation action keyframes

I was watching a recent Ian Hubert video where he has to use the “retarget” operation in the GYAZ Animation Tools Plugin to reset the animation, and when Ian can’t figure out an elegant way to solve this, you know the situation is dire.

For the uninformed, when you apply a pose as a rest pose on a rig with animation data, you get this warning:
“Warning: Actions on this armature will be destroyed by this new rest pose as the transforms stored are relative to the old rest pose”
And it’s no lie. None of the animation data gets updated, which usually means the animation is WRECT.
Here’s my first shot at a script to rescue the rig animation:

import bpy

a = bpy.context.object.animation_data.action
curframe = bpy.context.scene.frame_current

opstrs = ("[1]","pnt.handle_left[1]","pnt.handle_right[1]")
def noop(pnt,val): pass
def subop(pnt,val):
    for op in opstrs: exec(f"{op} -= val")
def divop(pnt,val):
    if val == 0: return
    for op in opstrs: exec(f"{op} /= val")
op_dict = {"location":subop, "rotation_euler":subop, "scale":divop}
#op_dict = {}

    obop = bpy.ops.object
    ed_mode =
    if ed_mode: obop.editmode_toggle()
    else:    obop.posemode_toggle()
if a is not None:
    print("No no!\nIgnore that noise, we've got magic that saves your animation data!")
    old_quats = {}
    old_euler = {}
    old_scale = {}
    old_loc = {}
    old_xfrms = (old_quats, old_euler, old_scale, old_loc)

    for fc in a.fcurves:
        val = fc.evaluate(curframe)
        dp = fc.data_path
        dps = dp.split(".")
        if dps[0] != "pose":
            #this is not pose animation, so leave it alone
        dptype = dps[-1]
        #process them now if possible
        # this only works for position and scale
        # when you add rotation everything breaks
        if dptype in op_dict:
            for pnt in fc.keyframe_points:
        # if not, store the transforms for later
        elif dptype == "rotation_quaternion":
            if dp not in old_quats: old_quats[dp] = [1,0,0,0]
            old_quats[dp][fc.array_index] = val
        elif dptype == "rotation_euler":
            if dp not in old_euler: old_euler[dp] = [0,0,0]
            old_euler[dp][fc.array_index] = val
        elif dptype == "location":
            if dp not in old_loc: old_loc[dp] = [0,0,0]
            old_loc[dp][fc.array_index] = val
        elif dptype == "scale":
            if dp not in old_scale: old_scale[dp] = [0,0,0]
            old_scale[dp][fc.array_index] = val

It works great!
Well, that is, it works when all you’re doing is animating location and scale. Once you throw rotation in there, all bets are off. Since rotation is most of what rigs are used for last time I checked, this is something of a profound failure.

I think the way forward is going to be to capture a full transform matrix for each keyframe location, and then invert the transform at the current frame. After applying the pose, I think we can just apply that invert transform matrix to each keyframe transform before decomposing the resulting matrix and plugging the new values into the appropriate keyframes.
Yes, I realize non-uniform scaling is going to result in lost child shear, but there’s no getting around that.

I’ve gotten as far as capturing the current frame transform for animated parameters, but haven’t pulled the non-animated transform data (do I need to?) or done anything with it (do I need to?). Also, I don’t handle quaternion transformations, though handling that should be a fairly straightforward sub-set of the transform matrix solution.
If there are any matrix wizards out there who can give me a few pointers (or outright fix my code), I’d appreciate it.

Otherwise, it’s pretty fun to scrub through the animation, run the script while in edit mode, and see the rest pose update right in front of you! Enjoy!