ChainRule.h
1//-*-C++-*-
2/***************************************************************************
3 *
4 * Copyright (C) 2004 by Willem van Straten
5 * Licensed under the Academic Free License version 2.1
6 *
7 ***************************************************************************/
8
9// psrchive/More/MEAL/MEAL/ChainRule.h
10
11#ifndef __MEAL_ChainRule_H
12#define __MEAL_ChainRule_H
13
14#include "MEAL/ProjectGradient.h"
15#include "MEAL/Composite.h"
16#include "MEAL/Scalar.h"
17
18namespace MEAL {
19
21 class ConstrainedParameter {
22
23 public:
24 ConstrainedParameter (unsigned iparam, bool infit, Scalar* func)
25 : scalar (func)
26 {
27 previously_infit = infit;
28 parameter = iparam;
29 }
30
32 unsigned parameter;
33
36
39
41 std::vector<double> gradient;
42
44 unsigned index;
45 };
46
48
56
57 template<class T>
58 class ChainRule : public T
59 {
60
61 public:
62
63 typedef typename T::Result Result;
64
66 ChainRule () : composite (this) { }
67
69 ChainRule (const ChainRule& rule) : composite (this) { operator = (rule); }
70
73
75 void set_model (T* model);
76 T* get_model () { return model; }
77
79 void set_constraint (unsigned iparam, Scalar* scalar);
80 bool has_constraint (unsigned iparam);
81 Scalar* get_constraint (unsigned iparam);
82
83 // ///////////////////////////////////////////////////////////////////
84 //
85 // Function implementation
86 //
87 // ///////////////////////////////////////////////////////////////////
88
90 std::string get_name () const;
91
92 protected:
93
95 void calculate (Result& result, std::vector<Result>* gradient);
96
98 std::vector<ConstrainedParameter> constraints;
99
102
105
106 private:
107
109 Composite composite;
110
111 };
112
115
116 void covariance ( Scalar* function, unsigned index,
117 std::vector<unsigned>& imap,
118 std::vector< std::vector<double> >& covar);
119
120}
121
122template<class T>
124{
125 return "ChainRule<" + std::string(T::Name)+ ">";
126}
127
128template<class T>
131{
132 if (this == &rule)
133 return *this;
134
135 set_model (rule.model);
136
137 for (unsigned ic=0; ic < rule.constraints.size(); ic++)
138 set_constraint (rule.constraints[ic].parameter,
139 rule.constraints[ic].scalar);
140
141 return *this;
142}
143
144template<class T>
147{
148 for (unsigned ifunc=0; ifunc<constraints.size(); ifunc++)
149 if (constraints[ifunc].parameter == iparam)
150 {
151 constraints[ifunc].index = ifunc;
152 return &constraints[ifunc];
153 }
154
155 return 0;
156}
157
158template<class T>
159void MEAL::ChainRule<T>::set_constraint (unsigned iparam, Scalar* scalar)
160{
161 if (this->get_verbose())
162 std::cerr << get_name() + "::set_constraint iparam=" << iparam
163 << std::endl;
164
165
166 ConstrainedParameter* existing = find_constraint (iparam);
167
168 // only one scalar may control a parameter
169 if (existing)
170 {
171 if (this->get_verbose())
172 std::cerr << get_name() + "::set_constraint"
173 " remove param=" << iparam << std::endl;
174
175 composite.unmap (existing->scalar);
176
177 if (model)
178 model->set_infit (iparam, existing->previously_infit);
179
180 constraints.erase (constraints.begin() + existing->index);
181 }
182
183 if (!scalar)
184 return;
185
186 bool infit = false;
187
188 if (model)
189 infit = model->get_infit (iparam);
190
191 constraints.push_back (ConstrainedParameter (iparam, infit, scalar));
192
193 if (this->get_verbose())
194 std::cerr << get_name() + "::set_constraint"
195 " add param=" << iparam << std::endl;
196
197 composite.map (constraints.back().scalar);
198
199 if (model)
200 model->set_infit (iparam, false);
201}
202
203template<class T>
204MEAL::Scalar* MEAL::ChainRule<T>::get_constraint (unsigned iparam)
205{
206 ConstrainedParameter* existing = find_constraint (iparam);
207 if (!existing)
208 throw Error (InvalidParam, "MEAL::ChainRule<T>::get_constraint",
209 "iparam=%u is not constrained", iparam);
210
211 return existing->scalar;
212}
213
214template<class T>
215bool MEAL::ChainRule<T>::has_constraint (unsigned iparam)
216{
217 return find_constraint (iparam) != 0;
218}
219
220template<class T>
222{
223 if (!_model)
224 return;
225
226 if (model)
227 {
228 if (this->get_verbose())
229 std::cerr << get_name() + "::set_model"
230 " unmap old model" << std::endl;
231 composite.unmap (model);
232 }
233
234 model = _model;
235
236 if (this->get_verbose())
237 std::cerr << get_name() + "::set_model"
238 " map new model" << std::endl;
239
240 composite.map (model);
241
242 for (unsigned ifunc=0; ifunc<constraints.size(); ifunc++)
243 {
244 unsigned iparam = constraints[ifunc].parameter;
245 constraints[ifunc].previously_infit = model->get_infit (iparam);
246 model->set_infit (iparam, false);
247 }
248}
249
250template<class T>
251void MEAL::ChainRule<T>::calculate (Result& result,
252 std::vector<Result>* grad)
253{
254 if (!model)
255 throw Error (InvalidState, get_name() + "::calculate","no model");
256
257 if (this->get_verbose())
258 std::cerr << get_name() + "::calculate" << std::endl;
259
260 for (unsigned ifunc=0; ifunc<constraints.size(); ifunc++)
261 {
262 if (T::very_verbose)
263 std::cerr << get_name() + "::calculate constraint on iparam="
264 << constraints[ifunc].parameter << " by "
265 << constraints[ifunc].scalar->get_name () << std::endl;
266
267 std::vector<double>* fgrad = 0;
268
269 if (grad)
270 fgrad = &(constraints[ifunc].gradient);
271
272 model->set_param (constraints[ifunc].parameter,
273 constraints[ifunc].scalar->evaluate(fgrad));
274 }
275
276 std::vector<Result> model_grad;
277 std::vector<Result>* model_grad_ptr = 0;
278 if (grad)
279 model_grad_ptr = & model_grad;
280
281 result = model->evaluate (model_grad_ptr);
282
283 if (grad)
284 {
285 unsigned ngrad = this->get_nparam();
286 grad->resize (ngrad);
287
288 unsigned igrad;
289 for (igrad=0; igrad<ngrad; igrad++)
290 (*grad)[igrad] = 0.0;
291
292 // map the model gradient
293 ProjectGradient (model, model_grad, *(grad));
294
295 // map the scalar gradients
296 std::vector<Result> fgrad;
297
298 for (unsigned ifunc=0; ifunc<constraints.size(); ifunc++)
299 {
300 unsigned iparam = constraints[ifunc].parameter;
301
302 ngrad = constraints[ifunc].gradient.size();
303 fgrad.resize (ngrad);
304
305 // dM/dxi = dM/df * df/dxi, where M=model, f=func, xi=func.param[i]
306 for (igrad=0; igrad<ngrad; igrad++)
307 fgrad[igrad] = model_grad[iparam] * constraints[ifunc].gradient[igrad];
308
309 ProjectGradient (constraints[ifunc].scalar, fgrad, *(grad));
310 }
311 }
312
313 if (this->get_verbose())
314 {
315 std::cerr << get_name() + "::calculate result\n"
316 " " << result << std::endl;
317 if (grad)
318 {
319 std::cerr << get_name() + "::calculate gradient\n";
320 for (unsigned i=0; i<grad->size(); i++)
321 std::cerr << " " << i << ":" << this->get_infit(i) << "="
322 << (*grad)[i] << std::endl;
323 }
324 }
325
326}
327
328#endif
Parameterizes a Function by one or more Scalar ordinates.
Definition ChainRule.h:59
ChainRule & operator=(const ChainRule &rule)
Assignment operator.
Definition ChainRule.h:130
ChainRule()
Default constructor.
Definition ChainRule.h:66
std::vector< ConstrainedParameter > constraints
Scalars used to constrain the parameters.
Definition ChainRule.h:98
void set_constraint(unsigned iparam, Scalar *scalar)
Set the Scalar instance used to constrain the specified parameter.
Definition ChainRule.h:159
Project< MEAL::Scalar > model
Definition ChainRule.h:104
ConstrainedParameter * find_constraint(unsigned iparam)
Find the matching constraint.
Definition ChainRule.h:146
std::string get_name() const
Return the name of the class.
Definition ChainRule.h:123
ChainRule(const ChainRule &rule)
Copy constructor.
Definition ChainRule.h:69
void calculate(Result &result, std::vector< Result > *gradient)
Return the Result and its gradient.
Definition ChainRule.h:251
void set_model(T *model)
Set the Function to be constrained by Scalar ordinates.
Definition ChainRule.h:221
Parameter policy for composite functions.
Definition Composite.h:20
A parameter and its constraining Scalar instance.
Definition ChainRule.h:21
unsigned parameter
The index of the Function parameter constrained by scalar.
Definition ChainRule.h:32
Project< Scalar > scalar
The constraining Scalar instance.
Definition ChainRule.h:38
std::vector< double > gradient
The last calculated gradient of the Scalar.
Definition ChainRule.h:41
bool previously_infit
The value of the fit flag before the parameter was constrained.
Definition ChainRule.h:35
unsigned index
Index within ChainRule vector.
Definition ChainRule.h:44
Template combines a reference to a Component and its Projection.
Definition Projection.h:60
Pure virtual base class of scalar functions.
Definition Scalar.h:20
Namespace in which all modeling and calibration related code is declared.
Definition ExampleComplex2.h:16
void covariance(Scalar *function, unsigned index, std::vector< unsigned > &imap, std::vector< std::vector< double > > &covar)
Definition ChainRule.C:12

Generated using doxygen 1.14.0