Miind
BinEstimator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2005 - 2012 Marc de Kamps
2 // All rights reserved.
3 //
4 // Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
5 //
6 // * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7 // * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation
8 // and/or other materials provided with the distribution.
9 // * Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
10 // without specific prior written permission.
11 //
12 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
13 // WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY
14 // DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
15 // USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
16 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
17 //
18 // If you use this software in work leading to a scientific publication, you should include a reference there to
19 // the 'currently valid reference', which can be found at http://miind.sourceforge.net
20 #include <gsl/gsl_errno.h>
21 #include <gsl/gsl_spline.h>
22 #include <GeomLib/BinEstimator.hpp>
23 #include "../NumtoolsLib/NumtoolsLib.h"
24 #include "GeomLibException.hpp"
25 
26 using namespace GeomLib;
27 using namespace NumtoolsLib;
28 
30 (
31  const vector<double>& vec_interpretation,
32  const OdeParameter& par_ode
33 ):
34 _par_ode(par_ode),
35 _vec_interpretation(vec_interpretation)
36 {
37 }
38 
40 }
41 
43  if( v > _par_ode._par_pop._theta )
44  throw GeomLibException("BinEstimator: v larger than V_max");
45  if( v < _par_ode._V_min )
46  throw GeomLibException("BinEstimator: v smaller than V_min");
47  Number n = _vec_interpretation.size();
48  if ( v > _vec_interpretation.back() )
49  return n-1;
50  // changed search: now from one onwards
51  for (MPILib::Index i = 1; i < n; i++)
52  if ( v < _vec_interpretation[i] )
53  return i-1;
54 
55  throw GeomLibException("No sensible return for SearchBin");
56 }
57 
59  int sg = (dv > 0) - (dv < 0);
60  int n = static_cast<int>(_vec_interpretation.size());
61  for (int i = 0; i < n; i++){
62  MPILib::Index j = modulo(ind + sg*i, n);
63  if (_vec_interpretation[j] <= v && (j == static_cast<MPILib::Index>(n-1) || _vec_interpretation[j+1] > v) )
64  return j;
65  }
66  assert(false);
67  throw GeomLibException("BinEstimator search failed");
68 }
69 
71 
72  int i = ind;
73  if ( v < _par_ode._V_min)
74  return -1;
75  int n = static_cast<int>(_vec_interpretation.size());
76  if ( v > _par_ode._par_pop._theta)
77  return n;
78 
79  if (dv == 0)
80  return i;
81 
82  int i_ret = Search(i,v, dv);
83  return i_ret;
84 }
85 
87 {
88  Potential diff = v + delta_v;
89  return diff;
90 }
91 
92 double BinEstimator::BinLowFraction(Potential v, int i_tr_low) const {
93  int n = _vec_interpretation.size();
94  if (i_tr_low < 0 || i_tr_low > n-1)
95  return 0.0;
96 
97  Potential low_bin = _vec_interpretation[i_tr_low];
98  Potential high_bin = (i_tr_low != static_cast<int>(_vec_interpretation.size()) - 1)? _vec_interpretation[i_tr_low + 1] : _par_ode._par_pop._theta;
99  assert( low_bin <= v && high_bin >= v);
100  return (high_bin -v )/(high_bin - low_bin);
101 }
102 
103 double BinEstimator::BinHighFraction(Potential v, int i_tr_high) const {
104  int n = _vec_interpretation.size();
105  if (i_tr_high < 0 || i_tr_high > n-1)
106  return 0.0;
107 
108  Potential low_bin = _vec_interpretation[i_tr_high];
109  Potential high_bin = (i_tr_high != static_cast<int>(_vec_interpretation.size()) - 1)? _vec_interpretation[i_tr_high + 1] : _par_ode._par_pop._theta;
110  assert( low_bin <= v && high_bin >= v);
111  return (v - low_bin)/(high_bin - low_bin);
112 }
113 
115  assert(i < _vec_interpretation.size());
116  CoverPair pair_ret;
117  Potential low = _vec_interpretation[i];
118  Potential high = (i != _vec_interpretation.size() - 1) ? _vec_interpretation[i+1] : _par_ode._par_pop._theta;
119 
120  Potential trans_low = Translate(low, delta_v);
121  Potential trans_high = Translate(high, delta_v);
122 
123  // delta_v is being passed on to determine the sign of the translation
124  int i_tr_low = this->SearchBin(i,trans_low, delta_v);
125  int i_tr_high = this->SearchBin(i,trans_high,delta_v);
126 
127  double alpha_low = BinLowFraction (trans_low, i_tr_low);
128  double alpha_high = BinHighFraction(trans_high, i_tr_high);
129 
130  pair_ret.first._alpha = alpha_low;
131  pair_ret.first._index = i_tr_low;
132  pair_ret.second._alpha = alpha_high;
133  pair_ret.second._index = i_tr_high;
134 
135  return pair_ret;
136 }
BinEstimator(const vector< Potential > &, const OdeParameter &)
Constructor.
double Potential
Contains the parameters necessary to configure a concrete OdeSystem instance. See AbstractOdeSystem a...
double BinHighFraction(Potential, int) const
CoverPair CalculateBinCover(MPILib::Index i, Potential h) const
Cover pair looks at the boundaries of the current bins and calculates where they would end up if a po...
unsigned int Number
unsigned int Index
Base class for all exceptions thrown in GeomLib.
MPILib::Index SearchBin(Potential) const
Use for testing purposes, for any normal usage CalculateBinCover should be used. This function return...
int Search(MPILib::Index, Potential, Potential) const
std::pair< BinCover, BinCover > CoverPair
Type representing begin an end bins of a translated bins and the respective covering fraction...
Potential Translate(Potential, Potential) const
The objective is find a numerical solution for this equation This requires a numerical representation of the density We will work in the state space of a two dimensional and define a mesh there We first give two examples and then define the general procedure and given in Table for given fixed Delta g see Fig and we will denote coordinates in this dimension by a small letter $v The second dimension can be used to represent parameters as varied as and will represented by $w A strip is constructed by choosing two neighbouring points in state e g and integrating the vector field for a time $T that is assumed to be an integer multiple of a period of time Delta which we assume to be a defining characteristic of the grid Let then the set of points the set of points which is quadrilateral in shape The quadrilateral should be but not necessarily as long as they are but it is convenient to number them in order of creation In the we will assume that strip numbers created by the integration procedure start and are so that the numbers i in each identify a unique strip Strip no is reserved for stationary points There may or more cells in strip The number of cells in strip $i denoted by with $i the strip number and $j the cell as the i
Definition: 2D.hpp:145
double BinLowFraction(Potential, int) const