Justine Devin / gammalib

Commit 5853ecbb37c4bedc2bc36631c3fd1f76c2aea87c

Authored by knodlseder
1 parent f5c94a68

add parameter limits and NaN/Inf checking

Inline Side-by-side
Showing 1 changed file with 59 additions and 8 deletions
src/model/GModelSpectralExpPlaw.cpp
  View file @5853ecb
... ... @@ -4,10 +4,18 @@
4 4 * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7   - * This program is free software; you can redistribute it and/or modify *
8   - * it under the terms of the GNU General Public License as published by *
9   - * the Free Software Foundation; either version 2 of the License, or *
10   - * (at your option) any later version. *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
11 19 * *
12 20 ***************************************************************************/
13 21 /**
... ... @@ -241,6 +249,24 @@ double GModelSpectralExpPlaw::eval(const GEnergy&amp; srcEng) const
241 249 double power = std::pow(e_norm, index()) * std::exp(-e_cut);
242 250 double value = norm() * power;
243 251  
  252 + // Compile option: Check for NaN/Inf
  253 + #if defined(G_NAN_CHECK)
  254 + if (std::isnan(value) || std::isinf(value)) {
  255 + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval";
  256 + std::cout << "(srcEng=" << srcEng << "):";
  257 + std::cout << " NaN/Inf encountered";
  258 + std::cout << " (value=" << value;
  259 + std::cout << ", energy=" << energy;
  260 + std::cout << ", index=" << index();
  261 + std::cout << ", ecut=" << ecut();
  262 + std::cout << ", pivot=" << pivot();
  263 + std::cout << ", e_norm=" << e_norm;
  264 + std::cout << ", e_cut=" << e_cut;
  265 + std::cout << ", power=" << power;
  266 + std::cout << ")" << std::endl;
  267 + }
  268 + #endif
  269 +
244 270 // Return
245 271 return value;
246 272 }
... ... @@ -294,6 +320,28 @@ double GModelSpectralExpPlaw::eval_gradients(const GEnergy&amp; srcEng) const
294 320 ((GModelSpectralExpPlaw*)this)->m_ecut.gradient(g_ecut);
295 321 ((GModelSpectralExpPlaw*)this)->m_pivot.gradient(g_pivot);
296 322  
  323 + // Compile option: Check for NaN/Inf
  324 + #if defined(G_NAN_CHECK)
  325 + if (std::isnan(value) || std::isinf(value)) {
  326 + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients";
  327 + std::cout << "(srcEng=" << srcEng << "):";
  328 + std::cout << " NaN/Inf encountered";
  329 + std::cout << " (value=" << value;
  330 + std::cout << ", energy=" << energy;
  331 + std::cout << ", index=" << index();
  332 + std::cout << ", ecut=" << ecut();
  333 + std::cout << ", pivot=" << pivot();
  334 + std::cout << ", e_norm=" << e_norm;
  335 + std::cout << ", e_cut=" << e_cut;
  336 + std::cout << ", power=" << power;
  337 + std::cout << ", g_norm=" << g_norm;
  338 + std::cout << ", g_index=" << g_index;
  339 + std::cout << ", g_ecut=" << g_ecut;
  340 + std::cout << ", g_pivot=" << g_pivot;
  341 + std::cout << ")" << std::endl;
  342 + }
  343 + #endif
  344 +
297 345 // Return
298 346 return value;
299 347 }
... ... @@ -569,8 +617,9 @@ void GModelSpectralExpPlaw::init_members(void)
569 617 m_norm.clear();
570 618 m_norm.name("Prefactor");
571 619 m_norm.unit("ph/cm2/s/MeV");
572   - m_norm.value(1.0);
573 620 m_norm.scale(1.0);
  621 + m_norm.value(1.0); // default: 1.0
  622 + m_norm.min(0.0); // min: 0.0
574 623 m_norm.free();
575 624 m_norm.gradient(0.0);
576 625 m_norm.hasgrad(true);
... ... @@ -578,8 +627,9 @@ void GModelSpectralExpPlaw::init_members(void)
578 627 // Initialise powerlaw index
579 628 m_index.clear();
580 629 m_index.name("Index");
581   - m_index.value(-2.0);
582 630 m_index.scale(1.0);
  631 + m_index.value(-2.0); // default: -2.0
  632 + m_index.range(-10.0,+10.0); // range: [-10,+10]
583 633 m_index.free();
584 634 m_index.gradient(0.0);
585 635 m_index.hasgrad(true);
... ... @@ -588,8 +638,9 @@ void GModelSpectralExpPlaw::init_members(void)
588 638 m_ecut.clear();
589 639 m_ecut.name("Cutoff");
590 640 m_ecut.unit("MeV");
591   - m_ecut.value(1000.0);
592 641 m_ecut.scale(1.0);
  642 + m_ecut.value(1000.0); // default: 1000.0
  643 + m_norm.min(0.1); // min: 0.1
593 644 m_ecut.free();
594 645 m_ecut.gradient(0.0);
595 646 m_ecut.hasgrad(true);
... ... @@ -598,8 +649,8 @@ void GModelSpectralExpPlaw::init_members(void)
598 649 m_pivot.clear();
599 650 m_pivot.name("PivotEnergy");
600 651 m_pivot.unit("MeV");
601   - m_pivot.value(100.0);
602 652 m_pivot.scale(1.0);
  653 + m_pivot.value(100.0);
603 654 m_pivot.fix();
604 655 m_pivot.gradient(0.0);
605 656 m_pivot.hasgrad(true);
... ...