report_type keyword for Profile Likelihood

[kantundpeterpan]

Hello!

I performed the calculations for obtaining the profile likelihood as in the docs. Parsing the output gives the report for multi_parameter_estimation. I looked into the code, found the report_type keyword, but i couldn't figure out where to use it during the Parameter Estimation configuration.

[CiaranWelsh]

Hey!

It’s been some time since I’ve worked on this project but if you post some code examples and describe what you’re trying to do I may be able to help.

[kantundpeterpan]

Basically, I am following the tutorial to obtain the profile likelihoods here. The scans over the parameters work nicely, but the reports do not contain the values of the parameter that was scanned over, but only the values for the remaining optimized parameters. So I wondered if there's a way to get the "right" reports.

Example:

My variables are _k_neo _k_dig _k_trans

from pycotools3 import tasks, viz, model

mod = model.loada(antimony_str, 'test.cps')

with tasks.ParameterEstimation.Context(
    mod, experiment_filename,
    context='pl', parameters='a'
) as context:
    context.set('method', 'hooke_jeeves')
    context.set('pe_number', 5) # number of steps in each profile likelihood
    context.set('run_mode', 'parallel')
    context.set('separator', ',')
    context.set('number_of_iterations', 200)
    #context.set('report', 'profile_likelihood')
    context.set('prefix', '_')
    config = context.get_config()

pe = tasks.ParameterEstimation(config)

p = viz.Parse(pe)

print(p.data['_k_neo'])

Example output for _k_neo gives the values for the optimized values for the remaining parameters, but lacks the values for ´_k_neo` was scanned over.

         RSS     _k_dig  _k_trans
0    27.1121   0.000898  0.006912
1   822.0270   0.000564  0.016508
2   871.2460   0.000519  9.164050
3  1481.4700  12.288000  2.816190
4        inf   0.100000  0.100000

[CiaranWelsh]

Ah I see, if I recall correctly, the issue with storing the actual values of the scanned parameter is that they are calculated within COPASI directly, and the communication system between pycotools and COPASI (i.e. the copasiML) doesn't allow access to these values. Therefore (as you'll see in viz.PlotProfileLikelihood) when plotting these data I ended up recalculating them:

def compute_x(self):
    lb = self.pl.config.settings.pl_lower_bound
    ub = self.pl.config.settings.pl_upper_bound
    parameters = self.get_best_original_parameter_set()
    num_steps = self.pl.config.settings.pe_number
    dct = {}
    for i in parameters:
        val = parameters.loc[0, i]
        low = numpy.log10(val / lb)
        high = numpy.log10(val * ub)
        dct[i] = pandas.Series(numpy.logspace(low, high, num_steps).flatten())
    df = pandas.concat(dct, axis=1)

    return df

If you wanted to do stuff with the data, it might be an idea to either copy and paste the source code for viz.PlotProfileLikelihood and use this as a starting point or deriving a new class from it. I.e.

class MyComplicatedBespokeProfileLikelihoodAnalysis(viz.PlotProfileLikelihood):
    pass

Personally, I'd opt for the former. Hope this helps.

Ciaran