Hi Justin,
I’ll play with your script over the weekend!
Simple things doable in Nuke (and work if you are happy to manually nudge things) are along the lines of rotation to Luminance&Chrominance/Luma&Chroma space (CIE LCh, YCoCg + Polar Transformation), and then adjust the Chrominance/Chroma plane to wish.
set cut_paste_input [stack 0]
version 12.1 v1
ColorWheel {
inputs 0
format "512 512 0 0 512 512 1 square_512"
name ColorWheel
selected true
xpos 550
ypos 4
}
Colorspace {
colorspace_out CIE-LCH
name To_CIE_LCh_Colorspace
selected true
xpos 550
ypos 76
}
set N78849000 [stack 0]
Gamma {
value 0.5
name chrominance_masking_Gamma
selected true
xpos 660
ypos 100
}
push $N78849000
Multiply {
inputs 1+1
channels {-rgba.red rgba.green -rgba.blue none}
maskChannelMask rgba.green
name compression_Multiply
selected true
xpos 550
ypos 104
}
Colorspace {
colorspace_in CIE-LCH
name From_CIE_LCh_Colorspace
selected true
xpos 550
ypos 147
}
Colorspace {
colorspace_out CIE-Yxy
name To_CIE_xyY_Colorspace
selected true
xpos 550
ypos 171
}
Shuffle2 {
fromInput1 {{0} B}
fromInput2 {{0} B}
mappings "4 rgba.blue 0 2 rgba.blue 0 2 rgba.alpha 0 3 rgba.alpha 0 3 rgba.red 0 0 rgba.green 0 1 rgba.green 0 1 rgba.red 0 0"
name Axis_Reorder_Shuffle
selected true
xpos 550
ypos 209
}
PositionToPoints2 {
display textured
render_mode textured
P_channel rgb
detail 1
pointSize 4
name PositionToPoints
selected true
xpos 550
ypos 253
}
This will effectively compress and expand the gamut non-linearly, from there you can introduce some heuristics to find the compression factor and the ideal mask for the Chrominance plane. This is the hard part depending on the requirements/constraints.
This is not dependent on the linearity of your inputs but on the linearity of the transformation you apply on those inputs. For example, if you are using a non-linear mask, then your transformation is non-linear and thus most likely not exposure invariant. Whether it is a problem or not is another question 
We have a lot of code for that in Colour which will be easier to translate (or test directly): https://github.com/colour-science/colour/blob/develop/colour/characterisation/correction.py
What you will find is that it is not trivial to port the polynomial methods into Nuke unless somebody can share code to apply a nxn matrix to colour data efficiently. Nuke, out of the box, is limited to 3x3 via the ColorMatrix node.
You will also find that none of the linear or polynomial methods is able to bring the data where you would like it to go. The polynomial methods will get you closer, for example, the light blue dots in that image:
No cigar, unfortunately! That being said, it could be the first step to make the work easier.
Finally, polynomials are EXTREMELY dangerous, they behave quite well on the training data set, but for anything not within it, they are bound to Cosmical Explosions! One of my favorite shirt, so much that I have it 
Cheers,
Thomas