User-defined model

[vdelapparent]

I am working on well resolved nearby galaxies and the availability of the priors is great!

However, other models would be better suited for spiral galaxies than Sersic (disk with inner/outer truncation), and a bar model would be helpful (more than 2/3 of galaxies have bars).

Actually, it would be wonderful to be able to add user-defined models...!

[louisquilley]

Hello, for the 0.16 release, it was mentioned that custom user-made models could now be included for the model-fitting. Unfortunately, I can't find any documentation about their implementation, and I am not familiar with machine-learning ONNX models, which is the format required, as said in 0.16 release. I am wondering if the problem mentioned here of having different profiles (truncated sérsic, bars), could be solved by this 0.16 add-on or if it is referring to something else, as ML seems overly complicated for the small change we are looking for. To be clear, we would like to change in the config file : add_model(measurement_group['band'], SersicModel(x, y, flux, rad, aspect, angle, nsersic)) to add_model(measurement_group['band'], CustomModel(x, y, flux, rad, aspect, angle, other parameters of the function)) for a simple analytical function that gives the flux as a function of radius and depends on a few parameters (that would be additional free parameters here).

Thank you in advance for your answer. Best regards, Louis

[marcschefer]

To avoid the performance issues of calling the Python interpreter for every pixel we have implemented custom models as ONNX compute graphs, this doesn't require using machine learning techniques but only the format. It's a bit convoluted but it allows the model function to be written in python and rendered fast enough for model fitting.

You can use any tool capable of creating ONNX format but here is an example using keras

x and y inputs are required, other inputs can be provided freely then it's just a matter of writing the python function

It's also recommended to save models with different size values and then SX++ can use the smallest one that fits the fitting region to avoid wasting performance.

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers

def exponential_model(inp):
    r = keras.backend.sqrt(inp[0]*inp[0]+inp[1]*inp[1])
    return inp[2] * keras.backend.exp(-inp[3] * r)

size = 200
k = layers.Input(shape=[1,], name='k')
i0 = layers.Input(shape=[1,], name='i0') # actually we don't need this input anymore as SourceXtractorPlusPlus will renormalize flux
x = layers.Input(shape=[size,size,1], name='x')
y = layers.Input(shape=[size,size,1], name='y')
out = layers.Lambda(exponential_model)([x, y, i0, k])

model = tf.keras.Model(inputs=[x, y, i0, k], outputs=[out])

model.summary()
model.save('exp_model')

# To convert the model model to ONNX format (tf2onnx is available as a pip package):
# python3 -m tf2onnx.convert --saved-model exp_model/ --output exp_model_200.onnx

This model can be used like this:

add_model(group, ComputeGraphModel(["exp25.onnx", "exp100.onnx", "exp200.onnx", "exp400.onnx"], x, y, flux[i], {"_scale": rad, "_aspect_ratio": ratio, "_angle": angle, "i0": 1.0, "k": 1.0}))

_scale, _aspect_ratio, _angle are special variables that handle those parameters automatically so you don't need to do that in your model.

It's a bit more complicated that writing the function directly in the configuration but I think it's still a lot simpler than having to write a custom C++ plugin while being close in performance.

Of course, we'd like to hear further feedback on whether that fits your needs and if the usability is adequate.

Cheers, Marc

[louisquilley]

Hello, thank you for your answer. I have tried using custom models following your recommendations and encountered a few problems.

First of all, as a test, I used your example of an exponential model to test the custom build method, and compared the results obtained in a 2-component model using a Sérsic and an exponential model. Replacing the built-in exp model by the one in onnx format retrieves extremely close values in all fields except for the rad_D value. I performed this check in a one band and in a multi-band model, and in both cases I have this discrepancy in the rad_D_band value. In all bands the rad_D_onnx/rad_D_default ratio has the same value (0.59). Could it be an effect of the size of the onnx model ? Because looking at the model images, they are restrained to a square inside the image when using onnx model. But this square is larger than the segmentation area so all relevant pixels are included.

This issue made me question something else regarding the effective radius of the exp (or Sérsic) profile in the case of an elliptical profile. Along which axis is computed the effective radius when using the built-in functions ? Major axis I guess, but it is never said clearly in the documentation.

Then I tried to build an exponential model with an inner truncation using the commands below.

def inner_truncated_exp_model(inp): r = keras.backend.sqrt(inp[0]*inp[0]+inp[1]*inp[1]) return inp[2] * keras.backend.exp( (r/inp[3]) + pow(inp[4]/r, inp[5]) )

size = 10000 r0 = layers.Input(shape=[1,], name='r0') rh = layers.Input(shape=[1,], name='rh') i0 = layers.Input(shape=[1,], name='i0') n = layers.Input(shape=[1,], name='n') x = layers.Input(shape=[size,size,1], name='x') y = layers.Input(shape=[size,size,1], name='y') out = layers.Lambda(inner_truncated_exp_model)([x, y, i0, r0, rh, n])

model = tf.keras.Model(inputs=[x, y, i0, r0, rh, n], outputs=[out])

model.summary() model.save('../Models/inn_truncated_exp_model_'+str(size))

Since there are 2 scale lengths here I was not sure of how to handle scale compared to the exponential model example, so I did the following in my config file : add_model(measurement_group['g'], ComputeGraphModel( ["inn_truncated_exp_model_800.onnx","inn_truncated_exp_model_3200.onnx", "inn_truncated_exp_model_6400.onnx"], x_D_g, y_D_g, flux_D_g, {"_scale": 1.0, "_aspect_ratio": aspect_D_g, "_angle": angle_D_g, "i0": 1.0, "r0": rad_D_g, "rh":rad_trunc_D_g, "n": 3.0})) and when launching a Sérsic + truncated exp with every parameter and prior copied from the "normal" exp case, except the addition of rad_trunc_D_g as a FreeParameter, it runs without any message error but returns fmf_stop_reason=7 in the catalogue, NaNs in the central region of the model and residual images, and all outputs as their initial value and 0 as associated error. I can send you the corresponding files if needed.

Thank you for your help. Best regards, Louis

[mkuemmel]

I am afraid we need @marcschefer to solve this, but he just started his vacation. He'll be back in August.

[mkuemmel]

@louisquilley and @vdelapparent,

I have applied for a focus demo at this years ADASS conference. If granted it would be great to present there also your user generated models. I have started a discussion on this at #499.

[marcschefer]

@louisquilley

Hi, I'm back sorry for the delay.

Yes, the factor 0.59 is to be expected for an exponential since we set k to 1, I did not include it in my example unfortunately but for the built-in models I have to convert the effective radius to k. For exponential it's k = 1.678 / Re and for a general Sersic profile with index n it's k = Bn / Re^(1/n) (we get Bn using the polynomial approximation from https://iopscience.iop.org/article/10.1086/344506/pdf )

We provide the scale parameter more as a convenience factor scale everything. I think setting it to 1 and using k as described would be the equivalent of the hard coded model. So I should have included that in my example.

As for your truncated exponential, it looks good overall but isn't there a missing minus sign inside the exponential? Could that be the issue?

Since it's been a few weeks, maybe you have made progress on that already, let me know if you need me to try to run it.

[louisquilley]

Hello, sorry for the delay, I just came back from vacation.

@mkuemmel I will keep working with these models and with @vdelapparent we will keep you updated in #499, but more generally I would be glad to participate and show to the community what I did with SE++.

@marcschefer Thank you for the clarification on scale length. And indeed the error came from the missing minus sign in the exponential. With it, model-fitting works in one band and multi-band so I can now explore further possibilities.