C++ GModelPar类代码示例

/***********************************************************************//**
 * @brief Read instrument scales from XML element
 *
 * @param[in] xml XML source element.
 *
 * Reads the instrument scale factors from a tag with the following format
 *
 *      <scaling>
 *         <instrument name="LAT" scale="1.0" min="0.1" max="10.0" value="1.0" free="0"/>
 *         <instrument name="CTA" scale="1.0" min="0.1" max="10.0" value="0.5" free="0"/>
 *      </scaling>
 *
 * The instrument name is case sensitive, but any leading and trailing white
 * space will be removed.
 *
 * If no scaling tag is found, all instrument scale factors will be cleared.
 ***************************************************************************/
void GModel::read_scales(const GXmlElement& xml)
{
    // Clear any existing scales
    m_scales.clear();

    // Continue only if there is a <scaling> tag
    if (xml.elements("scaling") != 0) {

        // Get pointer on first instrument scale factors
        const GXmlElement* scales = xml.element("scaling", 0);

        // Determine number of scale factors
        int nscales = scales->elements("instrument");

        // Read all scale factors
        for (int i = 0; i < nscales; ++i) {
            const GXmlElement* par = scales->element("instrument", i);
            GModelPar scale;
            scale.read(*par);
            scale.name(gammalib::strip_whitespace(par->attribute("name")));
            m_scales.push_back(scale);
        }

    }

    // Return
    return;
}

/***********************************************************************//**
 * @brief Test binned optimizer
 *
 * @param[in] datadir Directory of test data.
 * @param[in] irf Instrument response function.
 * @param[in] fit_results Expected fit result.
 *
 * Verifies the ability optimize binned Fermi/LAT data.
 ***************************************************************************/
void TestGLATOptimize::test_one_binned_optimizer(const std::string& datadir,
                                                 const std::string& irf,
                                                 const double*      fit_results)
{
    // Set filenames
    std::string lat_srcmap    = datadir+"/srcmap.fits";
    std::string lat_expmap    = datadir+"/binned_expmap.fits";
    std::string lat_ltcube    = datadir+"/ltcube.fits";
    std::string lat_model_xml = datadir+"/source_model.xml";

    // Setup GObservations for optimizing
    GObservations   obs;
    GLATObservation run;
    test_try("Setup for optimization");
    try {
        run.load_binned(lat_srcmap, lat_expmap, lat_ltcube);
        run.response(irf, lat_caldb);
        obs.append(run);
        test_try_success();
    }
    catch (std::exception &e) {
        test_try_failure(e);
    }

    // Load models from XML file
    obs.models(lat_model_xml);

    // Setup LM optimizer
    test_try("Perform LM optimization");
    try {
        GOptimizerLM opt;
        opt.max_iter(1000);
        obs.optimize(opt);
        obs.errors(opt);
        test_try_success();
        for (int i = 0, j = 0; i < obs.models().size(); ++i) {
            const GModel* model = obs.models()[i];
            for (int k = 0; k < model->size(); ++k) {
                GModelPar par  = (*model)[k];
                std::string msg = "Verify optimization result for " + par.print();
                test_value(par.value(), fit_results[j++], 5.0e-5, msg);
                test_value(par.error(), fit_results[j++], 5.0e-5, msg);
            }
        }
    }
    catch (std::exception &e) {
        test_try_failure(e);
    }

    // Exit test
    return;

}

/***********************************************************************//**
 * @brief Constructor
 *
 * @param[in] coeffs Vector of polynomial coefficients.
 ***************************************************************************/
GCTAModelRadialPolynom::GCTAModelRadialPolynom(const std::vector<double>& coeffs)
                                                          : GCTAModelRadial()
{
    // Initialise members
    init_members();

    // Assign coefficients
    for (int i = 0; i < coeffs.size(); ++i) {

        // Allocate parameter
        GModelPar par;

        // Set value
        par.real_value(coeffs[i]);

        // Set other attributes
        std::string name = "coeff"+str(i);
        par.name(name);
        par.unit("");
        par.free();
        par.scale(1.0);
        par.gradient(0.0);
        par.hasgrad(true);

        // Push coefficient on list
        m_coeffs.push_back(par);

    } // endfor: looped over coefficients

    // Update parameter mapping
    update_pars();

    // Return
    return;
}

/***********************************************************************//**
 * @brief Set model scale factor for a given instrument
 *
 * @param[in] par Model parameter for scaling.
 *
 * Sets the model parameter for a given instrument. The instrument name is
 * case sensitive, but any leading to trailing white space will be stripped.
 *
 * If the instrument is not yet defined it will be appended to the list of
 * instruments.
 ***************************************************************************/
void GModel::scale(const GModelPar& par)
{
    // String leading and trailing while space
    std::string instrument = gammalib::strip_whitespace(par.name());

    // Search for instrument and copy the model parameter if the instrument
    // has been found. Make sure that the instrument name is in upper case.
    bool found = false;
    for (int i = 0; i < m_scales.size(); ++i) {
        if (m_scales[i].name() == instrument) {
            found       = true;
            m_scales[i] = par;
            m_scales[i].name(instrument);
            break;
        }
    }

    // If instrument has not been found then append it now to the list
    // of instruments.
    if (!found) {
        m_scales.push_back(par);
        m_scales[m_scales.size()-1].name(instrument);
    }

    // Return
    return;
}

/***********************************************************************//**
 * @brief Test observations optimizer.
 *
 * @param[in] mode Testing mode.
 * 
 * This method supports two testing modes: 0 = unbinned and 1 = binned.
 ***************************************************************************/
void TestOpenMP::test_observations_optimizer(const int& mode)
{
    // Create Test Model
    GTestModelData model;

    // Create Models conteners
    GModels models;
    models.append(model);

    // Time iterval
    GTime tmin(0.0);
    GTime tmax(1800.0);

    // Rate : events/sec
    double rate = RATE;

    // Create observations
    GObservations obs;

    // Add some observation
    for (int i = 0; i < 6; ++i) {

        // Random Generator
        GRan ran;
        ran.seed(i);

        // Allocate events pointer
        GEvents *events;

        // Create either a event list or an event cube
        if (mode == UN_BINNED) {
            events = model.generateList(rate,tmin,tmax,ran);
        }
        else {
            events = model.generateCube(rate,tmin,tmax,ran);
        }

        // Create an observation
        GTestObservation ob;
        ob.id(gammalib::str(i));

        // Add events to the observation
        ob.events(*events);
        ob.ontime(tmax.secs()-tmin.secs());
        obs.append(ob);

        // Delete events pointer
        delete events;
        
    }

    // Add the model to the observation
    obs.models(models);

    // Create a GLog for show the interations of optimizer.
    GLog log;

    // Create an optimizer.
    GOptimizerLM opt(log);

    opt.max_stalls(50);

    // Optimize
    obs.optimize(opt);

    // Get the result
    GModelPar result = (*(obs.models()[0]))[0];

    // Check if converged
    test_assert(opt.status()==0, "Check if converged", 
                                 "Optimizer did not converge"); 

    // Check if value is correct
    test_value(result.factor_value(),RATE,result.factor_error()*3); 

    // Return
    return;
}

/***********************************************************************//**
 * @brief Return pivot energy
 *
 * @return Pivot energy.
 *
 * Returns the pivot energy.
 ***************************************************************************/
inline
GEnergy GModelSpectralPlaw::pivot(void) const
{
    GEnergy energy;
    energy.MeV(m_pivot.value());
    return energy;
}

/***********************************************************************//**
 * @brief Return exponential cut-off energy
 *
 * @return Exponential cut-off energy.
 *
 * Returns the exponential cut-off energy.
 ***************************************************************************/
inline
GEnergy GModelSpectralSuperExpPlaw::cutoff(void) const
{
    GEnergy energy;
    energy.MeV(m_ecut.value());
    return energy;
}

/***********************************************************************//**
 * @brief Return maximum energy
 *
 * @return Maximum energy.
 *
 * Returns the maximum energy.
 ***************************************************************************/
inline
GEnergy GModelSpectralPlaw2::emax(void) const
{
    GEnergy energy;
    energy.MeV(m_emax.value());
    return energy;
}

/***********************************************************************//**
 * @brief Return breakenergy energy
 *
 * @return breakenergy energy.
 *
 * Returns the breakenergy energy.
 ***************************************************************************/
inline
GEnergy GModelSpectralBrokenPlaw::breakenergy(void) const
{
    GEnergy energy;
    energy.MeV(m_breakenergy.value());
    return energy;
}

/***********************************************************************//**
 * @brief Set power law index2
 *
 * @param[in] index2 Power law index2.
 *
 * Sets the power law index2.
 ***************************************************************************/
inline
void GModelSpectralBrokenPlaw::index2(const double& index2)
{
    m_index2.value(index2);
    return;
}

/***********************************************************************//**
 * @brief Return power law index2
 *
 * @return Power law index2.
 *
 * Returns the power law index2.
 ***************************************************************************/
inline
double GModelSpectralBrokenPlaw::index2(void) const
{
    return (m_index2.value());
}

/***********************************************************************//**
 * @brief Set pre factor 
 *
 * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
 *
 * Sets the pre factor.
 ***************************************************************************/
inline
void GModelSpectralBrokenPlaw::prefactor(const double& prefactor)
{
    m_norm.value(prefactor);
    return;
}

/***********************************************************************//**
 * @brief Set normalization factor 
 *
 * @param[in] norm Normalization factor.
 *
 * Sets the normalization factor.
 ***************************************************************************/
inline
void GModelSpectralFunc::norm(const double& norm)
{
    m_norm.value(norm);
    return;
}

/***********************************************************************//**
 * @brief Set Declination of model centre
 *
 * @param[in] dec Declination (degrees).
 *
 * Sets the Declination of the model centre in degrees.
 ***************************************************************************/
inline
void GModelSpatialRadial::dec(const double& dec)
{
    m_dec.value(dec);
    return;
}

/***********************************************************************//**
 * @brief Set semi-major axis of ellipse
 *
 * @param[in] semimajor Semi-major axis of ellipse (degrees)
 *
 * Sets the semi-major axis of the ellipse in degrees.
 ***************************************************************************/
inline
void GModelSpatialEllipticalDisk::semimajor(const double& semimajor)
{
    m_semimajor.value(semimajor);
    return;
}