How to fix pbflux value for a sample of forward models?

[mwrona77]

Hi,
I want to compute a sample of forward models using distributions but I want to keep pbflux=2. Is there any way to do that using phoebe distributions and single run_compute(), assuming b.set_value('sample_num@compute', value=10)? I know I can do this by setting priors manually and running run_compute() 10 times, but it would be much more convenient if I can do this without a loop.

I can do this that way:

import phoebe
import numpy as np
from matplotlib import pyplot as plt

logger = phoebe.logger(clevel='WARNING')

b = phoebe.default_binary()


b.add_dataset(kind='lc', compute_phases=phoebe.linspace(-0.75, 0.25, 201), dataset='lc01')

b.set_value_all('ld_mode', 'manual')
b.set_value_all('atm', 'blackbody')
b.set_value_all('distortion_method', 'sphere')
b.set_value('irrad_method', 'none')


b.add_distribution({'ecc@component': phoebe.uniform(0.0,0.5),
                    'teff@secondary@component': phoebe.uniform(4000,6000)
                    }, distribution='my_dist')

e_sample, teff_sample = b['value@ecc@my_dist'].sample(10), b['value@teff@my_dist'].sample(10)

fig = plt.figure(figsize=(10,6))

for i in range(10):
    
    b.set_value('ecc@component', value=e_sample[i])
    b.set_value('teff@secondary@component', value=teff_sample[i])
    b.set_value('pblum@primary@lc01', value=4.*np.pi)
    
    pbflux = b.compute_pblums(pbflux=True)['pbflux@lc01'].value

    b.set_value('pblum@primary@lc01', value=4.*np.pi*2./pbflux)
    print(b.compute_pblums(pbflux=True))
    
    b.run_compute(model='model1', overwrite=True)
    plt.plot(b['value@times@model1'], b['value@fluxes@model1'])
    

plt.show()

But it would be great to do this in such a way:

import phoebe
import numpy as np
from matplotlib import pyplot as plt

logger = phoebe.logger(clevel='WARNING')

b = phoebe.default_binary()


b.add_dataset(kind='lc', compute_phases=phoebe.linspace(-0.75, 0.25, 201), dataset='lc01')

b.set_value_all('ld_mode', 'manual')
b.set_value_all('atm', 'blackbody')
b.set_value_all('distortion_method', 'sphere')
b.set_value('irrad_method', 'none')


b.add_distribution({'ecc@component': phoebe.uniform(0.0,0.5),
                    'teff@secondary@component': phoebe.uniform(4000,6000)
                    }, distribution='my_dist')

b.set_value('sample_from@compute', value='my_dist')
b.set_value('sample_mode@compute', value='all')
b.set_value('sample_num@compute', value=10)

# Option 1
b.set_value('pbflux@lc01', value=2.)
b.run_compute(model='model1', overwrite=True)
# Option 2
b.run_compute(model='model1', overwrite=True, pbflux_mode = 'const', pblux_value=2.)

b.plot(show=True)

Is there an option to this without using a loop for run_compute?

[kecnry]

Unfortunately pbflux is only supported as an output-parameter (since it requires the mesh at time t0 in order to compute). If you want to ignore distortion in the translation (which might not really be fair since you're sampling over eccentricity), you could write a custom constraint to set an approximate value for pblum within each iteration.

If your goal is just to have all the models "normalized" in some way, then my suggestion would be to just use a fixed value for pblum and handle the renormalization after-the-fact. You can use sample_mode='all' along with sample_mode=10 to expose all 10 models in the output and then either normalize each model based on its max-flux or to a specific phase/time-point. Alternatively you could keep your for-loop, but store each model separately by doing something like b.run_compute(model=f'model{i}', ...).

The only way I can think to do this within phoebe (in the future) is by allowing a custom callback function similar to the existing support for custom cost functions - do you think that would be useful or does that sound just as complicated?

[mwrona77]

You can use sample_mode='all' along with sample_mode=10 to expose all 10 models in the output and then either normalize each model based on its max-flux or to a specific phase/time-point.

That was my initial thought, but Andrej mentioned that it could be problematic for systems with ellipsoidal variability.

To provide context, I want to create a large set of light curves that will serve as the output for the neural network, with the system's parameters as the input. Aligning the light curves, at least on the baseline, will make it easier for the neural network to learn. Instead of using the original light curves:

I think it's better to use these (but maybe I'm completely wrong):

Alternatively you could keep your for-loop, but store each model separately by doing something like b.run_compute(model=f'model{i}', ...).

This is what I did, and it works, but it's much less convenient than having a single model with samples.

The only way I can think to do this within phoebe (in the future) is by allowing a custom callback function similar to the existing support for custom cost functions - do you think that would be useful or does that sound just as complicated?

I'm not sure how run_compute works under the hood, but I assume that when using samples, it must run some kind of loop. Therefore, if a user wants to use fixed pbflux, run_compute may just calculate pbflux at t0 and use this value to rescale pblum@primary, similar to what I did in the code at the beginning of this thread. This will increase computational time, but it would be more convenient for the user.

... or this is just only a one-time problem, and no one will use this feature in the future, so there's no need for implementing this.

[kecnry]

The complication is that the mapping from flux to luminosity would need computing the underlying mesh at time t0 and then iteratively finding the correct value of pblum to use internally (its not a simple mapping).

For your purposes here, if you want to anchor a specific time/phase to a flux value, you could add a fake observation (single time-flux pair to the dataset) and make use of pblum_mode = 'dataset-scaled' and then phoebe will do least-squares rescaling of the model to go through that one point... 🤔

[mwrona77]

The complication is that the mapping from flux to luminosity would need computing the underlying mesh at time t0 and then iteratively finding the correct value of pblum to use internally (its not a simple mapping).

Why isn't it a simple single mapping? First, I assume that pblum@primary is 4*pi. Since pblum_mode = 'component-coupled', pblum@secondary is rescaled automatically. Next, I obtain pbflux using b.compute_pblums(pbflux=True). Then, I simply rescale b.set_value('pblum@primary@lc01', value=4.*np.pi*2./pbflux) to obtain the desired pbflux = 2.0. Isn't the scaling factor between pblum@primary and pblum@secondary always the same for a given system at the specific time? Could it be a function of pblum@primary?

Sorry if it's a trivial question, I've started learning PHOEBE very recently.

For your purposes here, if you want to anchor a specific time/phase to a flux value, you could add a fake observation (single time-flux pair to the dataset) and make use of pblum_mode = 'dataset-scaled' and then phoebe will do least-squares rescaling of the model to go through that one point...

It works :) Unfortunately, it's much slower than:

    b.set_value('pblum@primary@lc01', value=4.*np.pi)
    
    pbflux = b.compute_pblums(pbflux=True)['pbflux@lc01'].value

    b.set_value('pblum@primary@lc01', value=4.*np.pi*2./pbflux)

For now I'll stay with this method, but I'll save all models using b.run_compute(model=f'model{i}', ...).

[kecnry]

Why isn't it a simple single mapping?

compute_pblums is creating a populating a mesh under-the-hood (so cannot be expressed as a simple equation) - although you're right, it likely wouldn't need iterating at least for the component-coupled case. I'll give this some thought and see if we can come up with a more convenient (and efficient) way to support this use-case, but its good to hear that the brute force way is doing the trick for now.