CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RKIntegrator.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:
6 #include <climits>
7 #include <stdexcept>
8 namespace Genfun {
9 FUNCTION_OBJECT_IMP(RKIntegrator::RKFunction)
10 
11 RKIntegrator::RKFunction::RKFunction(RKData *data, unsigned int index)
12  :_data(data),
13  _index(index)
14 {
15  _data->ref();
16 }
17 
19 {
20  _data->unref();
21 }
22 
24  :AbsFunction(right),
25  _data(right._data),
26  _index(right._index)
27 {
28  _data->ref();
29 }
30 
31 
32 double RKIntegrator::RKFunction::operator() (double t) const {
33  if (t<0) return 0;
34  if (!_data->_locked) _data->lock();
35 
36  // Do this first, thereafter, just read the cache
37  _data->recache();
38 
39 
40  // If the cache is empty, make an entry for t=0;
41  size_t nvar = _data->_startingValParameter.size();
42  if (_data->_fx.empty()) {
43  RKData::Data d(nvar);
44  d.time=0;
45  Argument arg(nvar);
46  for (size_t f=0;f<nvar;f++) {
48  arg[f]=d.variable[f];
49  }
50  _data->_fx.insert(d);
51  }
52 
53  if (t==0) return (*_data->_fx.begin()).variable[_index];
54 
55  RKData::Data dt(nvar);
56  dt.time=t;
57  std::set<RKData::Data>::iterator l =_data->_fx.lower_bound(dt);
58 
59  // We did find an exact value (measure 0), just return it.
60  if (l!=_data->_fx.end() && (*l).time==t) {
61  return (*l).variable[_index];
62  }
63 
64  else {
65  std::set<RKData::Data>::iterator u =_data->_fx.upper_bound(dt);
66 
67  while (u==_data->_fx.end()) {
68  u--;
69  RKData::Data newData(nvar);;
70  _data->_stepper->step(_data,*u, newData, 0);
71  _data->_fx.insert(l,newData);
72  if (newData.time==t) return newData.variable[_index];
73  u = _data->_fx.upper_bound(dt);
74  }
75  u--;
76  _data->_stepper->step(_data,*u, dt, t);
77  return dt.variable[_index];
78  }
79 }
80 
81 
82 RKIntegrator::RKData::RKData():_locked(false) {
83 }
84 
85 RKIntegrator::RKData::~RKData() {
86  for (size_t i=0;i<_startingValParameter.size();i++) delete _startingValParameter[i];
87  for (size_t i=0;i<_controlParameter.size();i++) delete _controlParameter[i];
88  for (size_t i=0;i<_diffEqn.size(); i++) delete _diffEqn[i];
89  delete _stepper;
90 }
91 
93  _data(new RKData())
94 {
95  if (stepper) _data->_stepper=stepper->clone();
96  else _data->_stepper= new AdaptiveRKStepper();
97  _data->ref();
98 }
99 
101  _data->unref();
102  for (size_t i=0;i<_fcn.size();i++) delete _fcn[i];
103 }
104 
106  const std::string &variableName,
107  double defStartingValue,
108  double defValueMin,
109  double defValueMax) {
110  Parameter *par = new Parameter(variableName, defStartingValue, defValueMin, defValueMax);
111  _data->_startingValParameter.push_back(par);
112  _data->_diffEqn.push_back(diffEquation->clone());
113  _data->_startingValParameterCache.push_back(defStartingValue);
114  _fcn.push_back(new RKFunction(_data,_fcn.size()));
115  return par;
116 }
117 
118 
119 
120 
121 
122 Parameter * RKIntegrator::createControlParameter (const std::string & variableName,
123  double defStartingValue,
124  double startingValueMin,
125  double startingValueMax) {
126 
127  Parameter *par = new Parameter(variableName, defStartingValue, startingValueMin, startingValueMax);
128  _data->_controlParameter.push_back(par);
129  _data->_controlParameterCache.push_back(defStartingValue);
130  return par;
131 
132 }
133 
134 
135 
137  return _fcn[i];
138 }
139 
140 
141 
143  if (!_locked) {
144  unsigned int size = _diffEqn.size();
145  for (size_t i=0;i<size;i++) {
146  if (!(_diffEqn[i]->dimensionality()==size)) throw std::runtime_error("Runtime error in RKIntegrator");
147  }
148  _locked=true;
149  }
150 }
151 
152 
153 
155 
156  bool stale=false;
157  if (!stale) {
158  for (size_t p=0;p<_startingValParameter.size();p++) {
159  if (_startingValParameter[p]->getValue()!=_startingValParameterCache[p]) {
160  _startingValParameterCache[p]=_startingValParameter[p]->getValue();
161  stale=true;
162  break;
163  }
164  }
165  }
166 
167  if (!stale) {
168  for (size_t p=0;p<_controlParameter.size();p++) {
169  if (_controlParameter[p]->getValue()!=_controlParameterCache[p]) {
170  _controlParameterCache[p]=_controlParameter[p]->getValue();
171  stale=true;
172  break;
173  }
174  }
175  }
176 
177  if (stale) {
178  _fx.erase(_fx.begin(),_fx.end());
179  }
180 
181 }
182 
183 
184 
186 
187 
188 } // namespace Genfun
Genfun::RKIntegrator::RKData::Data::time
double time
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:112
l
long l
Definition: JamesRandomSeeding.txt:30
Genfun::RKIntegrator::RKData::RKData
RKData()
Definition: RKIntegrator.cc:82
AdaptiveRKStepper.hh
Genfun::RKIntegrator::RKData::_startingValParameter
std::vector< Parameter * > _startingValParameter
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:123
Genfun::RKIntegrator::RKFunction::RKFunction
RKFunction(RKData *data, unsigned int index)
Definition: RKIntegrator.cc:11
Genfun::RKIntegrator::RKData::lock
void lock()
Definition: RKIntegrator.cc:142
Genfun::AbsFunction
Definition: CLHEP/GenericFunctions/AbsFunction.hh:48
Genfun::RKIntegrator::RKFunction::~RKFunction
virtual ~RKFunction()
Definition: RKIntegrator.cc:18
RKIntegrator.hh
Genfun::RKIntegrator::RKFunction::operator()
virtual double operator()(double argument) const
Definition: RKIntegrator.cc:32
Genfun::RKIntegrator::RKData::_stepper
const RKStepper * _stepper
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:132
Genfun::RKIntegrator::RKData::_controlParameterCache
std::vector< double > _controlParameterCache
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:127
Genfun::RKIntegrator::RKData::_startingValParameterCache
std::vector< double > _startingValParameterCache
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:124
Genfun::RKIntegrator::RKData::_locked
bool _locked
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:131
Genfun::RCBase::unref
void unref() const
Definition: RCBase.cc:20
Genfun::RKIntegrator::addDiffEquation
Parameter * addDiffEquation(const AbsFunction *diffEquation, const std::string &variableName="anon", double defStartingValue=0.0, double startingValueMin=0.0, double startingValueMax=0.0)
Definition: RKIntegrator.cc:105
Genfun::RKIntegrator::RKStepper
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:172
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
Genfun::RKIntegrator::RKData::Data
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:108
Genfun::RKIntegrator::RKData::recache
void recache()
Definition: RKIntegrator.cc:154
Genfun::RKIntegrator::RKData::Data::variable
std::vector< double > variable
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:110
f
void f(void g())
Definition: excDblThrow.cc:38
Variable.hh
Genfun::AbsFunction::clone
virtual AbsFunction * clone() const =0
Genfun::RKIntegrator::RKFunction
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:140
Genfun::RKIntegrator::RKStepper::~RKStepper
virtual ~RKStepper()
Definition: RKIntegrator.cc:185
Genfun::RKIntegrator::RKData
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:102
Genfun::RKIntegrator::RKIntegrator
RKIntegrator(const RKStepper *stepper=NULL)
Definition: RKIntegrator.cc:92
Genfun::Argument
Definition: CLHEP/GenericFunctions/Argument.hh:17
Genfun::RKIntegrator::~RKIntegrator
virtual ~RKIntegrator()
Definition: RKIntegrator.cc:100
Genfun::RKIntegrator::RKStepper::step
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, double timeLimit=0) const =0
Genfun::RKIntegrator::RKStepper::clone
virtual RKStepper * clone() const =0
Genfun::RKIntegrator::RKData::_diffEqn
std::vector< const AbsFunction * > _diffEqn
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:129
Genfun::RKIntegrator::getFunction
const RKFunction * getFunction(unsigned int i) const
Definition: RKIntegrator.cc:136
Genfun::RKIntegrator::createControlParameter
Parameter * createControlParameter(const std::string &variableName="anon", double defStartingValue=0.0, double startingValueMin=0.0, double startingValueMax=0.0)
Definition: RKIntegrator.cc:122
Genfun::AdaptiveRKStepper
Definition: CLHEP/GenericFunctions/AdaptiveRKStepper.hh:10
i
long i
Definition: JamesRandomSeeding.txt:27
Genfun::RCBase::ref
void ref() const
Definition: RCBase.cc:15
Genfun::RKIntegrator::RKData::_controlParameter
std::vector< Parameter * > _controlParameter
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:126
Genfun::Parameter
Definition: CLHEP/GenericFunctions/Parameter.hh:35
FUNCTION_OBJECT_IMP
#define FUNCTION_OBJECT_IMP(classname)
Definition: CLHEP/GenericFunctions/AbsFunction.hh:156
Genfun::RKIntegrator
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:43
Genfun::RKIntegrator::RKData::_fx
std::set< Data > _fx
Definition: CLHEP/GenericFunctions/RKIntegrator.hh:130
Genfun
Definition: CLHEP/GenericFunctions/Abs.hh:14