bpp-core  2.1.0
 All Classes Namespaces Files Functions Variables Typedefs Friends
ThreePointsNumericalDerivative.cpp
Go to the documentation of this file.
1 //
2 // File: ThreePointsNumericalDerivative.cpp
3 // Created by: Julien Dutheil
4 // Created on: Thu Aug 17 15:00 2006
5 //
6 
7 /*
8  Copyright or © or Copr. Bio++ Development Team, (November 17, 2004)
9 
10  This software is a computer program whose purpose is to provide classes
11  for numerical calculus.
12 
13  This software is governed by the CeCILL license under French law and
14  abiding by the rules of distribution of free software. You can use,
15  modify and/ or redistribute the software under the terms of the CeCILL
16  license as circulated by CEA, CNRS and INRIA at the following URL
17  "http://www.cecill.info".
18 
19  As a counterpart to the access to the source code and rights to copy,
20  modify and redistribute granted by the license, users are provided only
21  with a limited warranty and the software's author, the holder of the
22  economic rights, and the successive licensors have only limited
23  liability.
24 
25  In this respect, the user's attention is drawn to the risks associated
26  with loading, using, modifying and/or developing or reproducing the
27  software by the user in light of its specific status of free software,
28  that may mean that it is complicated to manipulate, and that also
29  therefore means that it is reserved for developers and experienced
30  professionals having in-depth computer knowledge. Users are therefore
31  encouraged to load and test the software's suitability as regards their
32  requirements in conditions enabling the security of their systems and/or
33  data to be ensured and, more generally, to use and operate it in the
34  same conditions as regards security.
35 
36  The fact that you are presently reading this means that you have had
37  knowledge of the CeCILL license and that you accept its terms.
38  */
39 
40 #include "ThreePointsNumericalDerivative.h"
41 
42 using namespace bpp;
43 using namespace std;
44 
45 void ThreePointsNumericalDerivative::updateDerivatives(const ParameterList parameters)
46 throw (ParameterNotFoundException, ConstraintException)
47 {
48  if (computeD1_ && variables_.size() > 0)
49  {
50  if (function1_) function1_->enableFirstOrderDerivatives(false);
51  if (function2_) function2_->enableSecondOrderDerivatives(false);
52  function_->setParameters(parameters);
53  f2_ = function_->getValue();
54  if ((abs(f2_) >= NumConstants::VERY_BIG()) || isnan(f2_)){
55  for (unsigned int i = 0; i < variables_.size(); i++){
56  der1_[i]=log(-1);
57  der2_[i]=log(-1);
58  }
59  return;
60  }
61 
62  string lastVar;
63  bool functionChanged = false;
64  ParameterList p;
65  bool start = true;
66  for (unsigned int i = 0; i < variables_.size(); i++)
67  {
68  string var = variables_[i];
69  if (!parameters.hasParameter(var)) continue;
70  if (!start)
71  {
72  vector<string> vars(2);
73  vars[0] = var;
74  vars[1] = lastVar;
75  p = parameters.subList(vars);
76  }
77  else
78  {
79  p = parameters.subList(var);
80  start = false;
81  }
82  lastVar = var;
83  functionChanged = true;
84  double value = function_->getParameterValue(var);
85  double h = -(1. + std::abs(value)) * h_;
86  if (abs(h)<p[0].getPrecision())
87  h=h<0?-p[0].getPrecision():p[0].getPrecision();
88  double hf1(0), hf3(0);
89  unsigned int nbtry=0;
90 
91  //Compute f1_
92  while (hf1==0){
93  try
94  {
95  p[0].setValue(value + h);
96  function_->setParameters(p); //also reset previous parameter...
97 
98  p = p.subList(0);
99  f1_ = function_->getValue();
100  if ((abs(f1_) >= NumConstants::VERY_BIG()) || isnan(f1_))
101  throw ConstraintException("f1_ too large", &p[0], f1_);
102  else
103  hf1=h;
104  }
105  catch (ConstraintException& ce)
106  {
107  if (++nbtry==10) // no possibility to compute derivatives
108  break;
109  else
110  if (h<0)
111  h=-h; // try on the right
112  else
113  h/=-2; // try again on the left with smaller interval
114  }
115  }
116 
117  if (hf1!=0){
118  //Compute f3_
119  if (h<0)
120  h=-h; // on the right
121  else
122  h/=2; // on the left with smaller interval
123 
124  nbtry=0;
125  while (hf3==0){
126  try
127  {
128  p[0].setValue(value + h);
129  function_->setParameters(p); //also reset previous parameter...
130 
131  p = p.subList(0);
132  f3_ = function_->getValue();
133  if ((abs(f3_) >= NumConstants::VERY_BIG()) || isnan(f3_))
134  throw ConstraintException("f3_ too large", &p[0], f3_);
135  else
136  hf3=h;
137  }
138  catch (ConstraintException& ce)
139  {
140  if (++nbtry==10) // no possibility to compute derivatives
141  break;
142  else
143  if (h<0)
144  h=-h; // try on the right
145  else
146  h/=-2; // try again on the left with smaller interval
147  }
148  }
149  }
150 
151  if (hf3==0){
152  der1_[i]=log(-1);
153  der2_[i]=log(-1);
154  }
155  else {
156  der1_[i] = (f1_ - f3_) / (hf1-hf3);
157  der2_[i] = ((f1_ - f2_)/hf1 - (f3_ - f2_)/hf3)*2/(hf1-hf3);
158  }
159  }
160 
161 
162 
163 
164  if (computeCrossD2_)
165  {
166  string lastVar1, lastVar2;
167  for (unsigned int i = 0; i < variables_.size(); i++)
168  {
169  string var1 = variables_[i];
170  if (!parameters.hasParameter(var1)) continue;
171  for (unsigned int j = 0; j < variables_.size(); j++)
172  {
173  if (j == i)
174  {
175  crossDer2_(i, j) = der2_[i];
176  continue;
177  }
178  string var2 = variables_[j];
179  if (!parameters.hasParameter(var2)) continue;
180 
181  vector<string> vars(2);
182  vars[0] = var1;
183  vars[1] = var2;
184  if (i > 0 && j > 0)
185  {
186  if (lastVar1 != var1 && lastVar1 != var2) vars.push_back(lastVar1);
187  if (lastVar2 != var1 && lastVar2 != var2) vars.push_back(lastVar2);
188  }
189  p = parameters.subList(vars);
190 
191  double value1 = function_->getParameterValue(var1);
192  double value2 = function_->getParameterValue(var2);
193  double h1 = (1. + std::abs(value1)) * h_;
194  double h2 = (1. + std::abs(value2)) * h_;
195 
196  //Compute 4 additional points:
197  try
198  {
199  p[0].setValue(value1 - h1);
200  p[1].setValue(value2 - h2);
201  function_->setParameters(p); //also reset previous parameter...
202  vector<size_t> tmp(2);
203  tmp[0] = 0;
204  tmp[1] = 1;
205  p = p.subList(tmp); //removed the previous parameters.
206  f11_ = function_->getValue();
207 
208  p[1].setValue(value2 + h2);
209  function_->setParameters(p.subList(1));
210  f12_ = function_->getValue();
211 
212  p[0].setValue(value1 + h1);
213  function_->setParameters(p.subList(0));
214  f22_ = function_->getValue();
215 
216  p[1].setValue(value2 - h2);
217  function_->setParameters(p.subList(1));
218  f21_ = function_->getValue();
219 
220  crossDer2_(i, j) = ((f22_ - f21_) - (f12_ - f11_)) / (4 * h1 * h2);
221  }
222  catch (ConstraintException& ce)
223  {
224  throw Exception("ThreePointsNumericalDerivative::setParameters. Could not compute cross derivatives at limit.");
225  }
226 
227  lastVar1 = var1;
228  lastVar2 = var2;
229  }
230  }
231  }
232 
233  //Reset last parameter and compute analytical derivatives if any.
234  if (function1_) function1_->enableFirstOrderDerivatives(computeD1_);
235  if (function2_) function2_->enableSecondOrderDerivatives(computeD2_);
236  if (functionChanged)
237  function_->setParameters(parameters.subList(lastVar));
238  }
239  else
240  {
241  //Reset initial value and compute analytical derivatives if any.
242  if (function1_) function1_->enableFirstOrderDerivatives(computeD1_);
243  if (function2_) function2_->enableSecondOrderDerivatives(computeD2_);
244  function_->setParameters(parameters);
245  //Just in case derivatives are not computed:
246  f2_ = function_->getValue();
247  }
248 }
249