| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 << 26 // $Id: G4AdjointInterpolator.cc,v 1.4 2009/11/20 10:31:20 ldesorgh Exp $
>> 27 // GEANT4 tag $Name: geant4-09-03-patch-01 $
>> 28 //
>> 29 #include "G4AdjointCSMatrix.hh"
27 #include "G4AdjointInterpolator.hh" 30 #include "G4AdjointInterpolator.hh"
28 31
29 G4ThreadLocal G4AdjointInterpolator* G4Adjoint << 32 G4AdjointInterpolator* G4AdjointInterpolator::theInstance = 0;
30 <<
31 ////////////////////////////////////////////// 33 ///////////////////////////////////////////////////////
>> 34 //
32 G4AdjointInterpolator* G4AdjointInterpolator:: 35 G4AdjointInterpolator* G4AdjointInterpolator::GetAdjointInterpolator()
33 { << 36 { if(theInstance == 0) {
34 return GetInstance(); << 37 static G4AdjointInterpolator interpolator;
>> 38 theInstance = &interpolator;
>> 39 }
>> 40 return theInstance;
35 } 41 }
36 <<
37 ////////////////////////////////////////////// 42 ///////////////////////////////////////////////////////
>> 43 //
38 G4AdjointInterpolator* G4AdjointInterpolator:: 44 G4AdjointInterpolator* G4AdjointInterpolator::GetInstance()
39 { << 45 { if(theInstance == 0) {
40 if(!fInstance) << 46 static G4AdjointInterpolator interpolator;
41 { << 47 theInstance = &interpolator;
42 fInstance = new G4AdjointInterpolator; <<
43 } 48 }
44 return fInstance; << 49 return theInstance;
45 } 50 }
46 51
47 ////////////////////////////////////////////// 52 ///////////////////////////////////////////////////////
48 G4AdjointInterpolator::G4AdjointInterpolator() << 53 //
49 << 54 G4AdjointInterpolator::G4AdjointInterpolator()
>> 55 {;
>> 56 }
50 ////////////////////////////////////////////// 57 ///////////////////////////////////////////////////////
51 G4AdjointInterpolator::~G4AdjointInterpolator( << 58 //
52 << 59 G4AdjointInterpolator::~G4AdjointInterpolator()
>> 60 {;
>> 61 }
53 ////////////////////////////////////////////// 62 ///////////////////////////////////////////////////////
54 G4double G4AdjointInterpolator::LinearInterpol << 63 //
55 << 64 G4double G4AdjointInterpolator::LinearInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
56 << 65 { G4double res = y1+ (x-x1)*(y2-y1)/(x2-x1);
57 { << 66 //G4cout<<"Linear "<<res<<G4endl;
58 G4double res = y1 + (x - x1) * (y2 - y1) / ( << 67 return res;
59 return res; <<
60 } 68 }
61 <<
62 ////////////////////////////////////////////// 69 ///////////////////////////////////////////////////////
63 G4double G4AdjointInterpolator::LogarithmicInt << 70 //
64 G4double& x, G4double& x1, G4double& x2, G4d << 71 G4double G4AdjointInterpolator::LogarithmicInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
65 { << 72 { if (y1<=0 || y2<=0 || x1<=0) return LinearInterpolation(x,x1,x2,y1,y2);
66 if(y1 <= 0. || y2 <= 0. || x1 <= 0.) << 73 G4double B=std::log(y2/y1)/std::log(x2/x1);
67 return LinearInterpolation(x, x1, x2, y1, << 74 //G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<G4endl;
68 G4double B = std::log(y2 / y1) / std::log( << 75 G4double A=y1/std::pow(x1,B);
69 G4double A = y1 / std::pow(x1, B); << 76 G4double res=A*std::pow(x,B);
70 G4double res = A * std::pow(x, B); << 77 // G4cout<<"Log "<<res<<G4endl;
71 return res; 78 return res;
72 } 79 }
73 <<
74 ////////////////////////////////////////////// 80 ///////////////////////////////////////////////////////
75 G4double G4AdjointInterpolator::ExponentialInt << 81 //
76 G4double& x, G4double& x1, G4double& x2, G4d << 82 G4double G4AdjointInterpolator::ExponentialInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
77 { << 83 { G4double B=(std::log(y2)-std::log(y1));
78 G4double B = (std::log(y2) - std::log(y1)) << 84 B=B/(x2-x1);
79 G4double A = y1 * std::exp(-B * x1); << 85 G4double A=y1*std::exp(-B*x1);
80 G4double res = A * std::exp(B * x); << 86 G4double res=A*std::exp(B*x);
81 return res; 87 return res;
82 } 88 }
83 <<
84 ////////////////////////////////////////////// 89 ///////////////////////////////////////////////////////
85 G4double G4AdjointInterpolator::Interpolation( << 90 //
86 << 91 G4double G4AdjointInterpolator::Interpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2,G4String InterPolMethod)
87 <<
88 <<
89 { 92 {
90 if(InterPolMethod == "Log") << 93 if (InterPolMethod == "Log" ){
91 { << 94 return LogarithmicInterpolation(x,x1,x2,y1,y2);
92 return LogarithmicInterpolation(x, x1, x2, << 95 }
93 } << 96 else if (InterPolMethod == "Lin" ){
94 else if(InterPolMethod == "Lin") << 97 return LinearInterpolation(x,x1,x2,y1,y2);
95 { << 98 }
96 return LinearInterpolation(x, x1, x2, y1, << 99 else if (InterPolMethod == "Exp" ){
97 } << 100 return ExponentialInterpolation(x,x1,x2,y1,y2);
98 else if(InterPolMethod == "Exp") << 101 }
99 { << 102 else {
100 return ExponentialInterpolation(x, x1, x2, << 103 //G4cout<<"The interpolation method that you invoked does not exist!"<<G4endl;
101 } << 104 return -1111111111.;
102 else <<
103 { <<
104 G4ExceptionDescription ed; <<
105 ed << "The interpolation method that you i <<
106 G4Exception("G4AdjointInterpolator::Interp <<
107 FatalException, ed); <<
108 return 0.; <<
109 } 105 }
110 } 106 }
111 <<
112 ////////////////////////////////////////////// 107 ///////////////////////////////////////////////////////
113 // only valid if x_vec is monotically increasi << 108 //
114 std::size_t G4AdjointInterpolator::FindPositio << 109 size_t G4AdjointInterpolator::FindPosition(G4double& x,std::vector<G4double>& x_vec,size_t , size_t ) //only valid if x_vec is monotically increasing
115 std << 110 { //most rapid nethod could be used probably
116 std << 111 //It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
117 { << 112
118 // most rapid method could be used probably << 113
119 << 114 size_t ndim = x_vec.size();
120 std::size_t ndim = x_vec.size(); << 115 size_t ind1 = 0;
121 std::size_t ind1 = 0; << 116 size_t ind2 = ndim - 1;
122 std::size_t ind2 = ndim - 1; << 117 /* if (ind_max >= ind_min){
123 << 118 ind1=ind_min;
124 if(ndim > 1) << 119 ind2=ind_max;
125 { << 120
126 if(x_vec[0] < x_vec[1]) << 121
127 { // increasing << 122 }
128 do << 123 */
129 { << 124
130 std::size_t midBin = (ind1 + ind2) / 2 << 125
131 if(x < x_vec[midBin]) << 126 if (ndim >1) {
132 ind2 = midBin; << 127
133 else << 128 if (x_vec[0] < x_vec[1] ) { //increasing
134 ind1 = midBin; << 129 do {
135 } while(ind2 - ind1 > 1); << 130 size_t midBin = (ind1 + ind2)/2;
136 } << 131 if (x < x_vec[midBin])
137 else << 132 ind2 = midBin;
138 { << 133 else
139 do << 134 ind1 = midBin;
140 { << 135 } while (ind2 - ind1 > 1);
141 std::size_t midBin = (ind1 + ind2) / 2 << 136 }
142 if(x < x_vec[midBin]) << 137 else {
143 ind1 = midBin; << 138 do {
144 else << 139 size_t midBin = (ind1 + ind2)/2;
145 ind2 = midBin; << 140 if (x < x_vec[midBin])
146 } while(ind2 - ind1 > 1); << 141 ind1 = midBin;
147 } << 142 else
>> 143 ind2 = midBin;
>> 144 } while (ind2 - ind1 > 1);
>> 145 }
>> 146
148 } 147 }
149 << 148
150 return ind1; 149 return ind1;
151 } 150 }
152 151
153 ////////////////////////////////////////////// 152 ///////////////////////////////////////////////////////
154 // only valid if x_vec is monotically increasi << 153 //
155 std::size_t G4AdjointInterpolator::FindPositio << 154 size_t G4AdjointInterpolator::FindPositionForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec) //only valid if x_vec is monotically increasing
156 G4double& log_x, std::vector<G4double>& log_ << 155 { //most rapid nethod could be used probably
157 { << 156 //It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
158 // most rapid method could be used probably <<
159 return FindPosition(log_x, log_x_vec); 157 return FindPosition(log_x, log_x_vec);
160 } << 158 if (log_x_vec.size()>3){
161 << 159 size_t ind=0;
162 ////////////////////////////////////////////// << 160 G4double log_x1=log_x_vec[1];
163 G4double G4AdjointInterpolator::Interpolate(G4 << 161 G4double d_log =log_x_vec[2]-log_x1;
164 st << 162 G4double dind=(log_x-log_x1)/d_log +1.;
165 st << 163 if (dind <1.) ind=0;
166 co << 164 else if (dind >= double(log_x_vec.size())-2.) ind =log_x_vec.size()-2;
167 { << 165 else ind =size_t(dind);
168 std::size_t i = FindPosition(x, x_vec); << 166 return ind;
169 return Interpolation(x, x_vec[i], x_vec[i + << 167
170 InterPolMethod); << 168 }
171 } << 169 else return FindPosition(log_x, log_x_vec);
172 << 170
173 ////////////////////////////////////////////// << 171
174 G4double G4AdjointInterpolator::InterpolateWit << 172 }
175 G4double& x, std::vector<G4double>& x_vec, s << 173 ///////////////////////////////////////////////////////
176 std::vector<std::size_t>& index_vec, G4doubl << 174 //
177 G4double dx) // only linear interpolation p << 175 G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,G4String InterPolMethod)
178 { << 176 { size_t i=FindPosition(x,x_vec);
179 std::size_t ind = 0; << 177 //G4cout<<i<<G4endl;
180 if(x > x0) << 178 //G4cout<<x<<G4endl;
181 ind = G4int((x - x0) / dx); << 179 //G4cout<<x_vec[i]<<G4endl;
182 if(ind >= index_vec.size() - 1) << 180 return Interpolation( x,x_vec[i],x_vec[i+1],y_vec[i],y_vec[i+1],InterPolMethod);
183 ind = index_vec.size() - 2; << 181 }
184 std::size_t ind1 = index_vec[ind]; << 182
185 std::size_t ind2 = index_vec[ind + 1]; << 183 ///////////////////////////////////////////////////////
186 if(ind1 > ind2) << 184 //
187 { << 185 G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,
188 std::size_t ind11 = ind1; << 186 std::vector<size_t>& index_vec,G4double x0, G4double dx) //only linear interpolation possible
189 ind1 = ind2; << 187 { size_t ind=0;
190 ind2 = ind11; << 188 if (x>x0) ind=int((x-x0)/dx);
191 } << 189 if (ind >= index_vec.size()-1) ind= index_vec.size()-2;
192 ind = FindPosition(x, x_vec, ind1, ind2); << 190 size_t ind1 = index_vec[ind];
193 return Interpolation(x, x_vec[ind], x_vec[in << 191 size_t ind2 = index_vec[ind+1];
194 y_vec[ind + 1], "Lin"); << 192 if (ind1 >ind2) {
195 } << 193 size_t ind11=ind1;
196 << 194 ind1=ind2;
197 ////////////////////////////////////////////// << 195 ind2=ind11;
198 G4double G4AdjointInterpolator::InterpolateFor << 196
199 G4double& log_x, std::vector<G4double>& log_ << 197 }
200 std::vector<G4double>& log_y_vec) << 198
201 { << 199 ind=FindPosition(x,x_vec,ind1,ind2);
202 std::size_t i = FindPositionForLogVector(log << 200 return Interpolation( x,x_vec[ind],x_vec[ind+1],y_vec[ind],y_vec[ind+1],"Lin");
203 << 201
204 G4double log_y = LinearInterpolation(log_x, << 202 }
205 log_y_v << 203
206 return log_y; << 204
207 } << 205 ///////////////////////////////////////////////////////
>> 206 //
>> 207 G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec,std::vector<G4double>& log_y_vec)
>> 208 { //size_t i=0;
>> 209 size_t i=FindPositionForLogVector(log_x,log_x_vec);
>> 210 /*G4cout<<"In interpolate "<<G4endl;
>> 211 G4cout<<i<<G4endl;
>> 212 G4cout<<log_x<<G4endl;
>> 213 G4cout<<log_x_vec[i]<<G4endl;
>> 214 G4cout<<log_x_vec[i+1]<<G4endl;
>> 215 G4cout<<log_y_vec[i]<<G4endl;
>> 216 G4cout<<log_y_vec[i+1]<<G4endl;*/
>> 217
>> 218 G4double log_y=LinearInterpolation(log_x,log_x_vec[i],log_x_vec[i+1],log_y_vec[i],log_y_vec[i+1]);
>> 219 return log_y;
>> 220
>> 221 }
208 222