| 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 69844 2013-05-16 09:19:33Z gcosmo $
>> 27 //
>> 28 #include "G4AdjointCSMatrix.hh"
27 #include "G4AdjointInterpolator.hh" 29 #include "G4AdjointInterpolator.hh"
28 30
29 G4ThreadLocal G4AdjointInterpolator* G4Adjoint << 31 G4AdjointInterpolator* G4AdjointInterpolator::theInstance = 0;
30 <<
31 ////////////////////////////////////////////// 32 ///////////////////////////////////////////////////////
>> 33 //
32 G4AdjointInterpolator* G4AdjointInterpolator:: 34 G4AdjointInterpolator* G4AdjointInterpolator::GetAdjointInterpolator()
33 { << 35 { if(theInstance == 0) {
34 return GetInstance(); << 36 static G4AdjointInterpolator interpolator;
>> 37 theInstance = &interpolator;
>> 38 }
>> 39 return theInstance;
35 } 40 }
36 <<
37 ////////////////////////////////////////////// 41 ///////////////////////////////////////////////////////
>> 42 //
38 G4AdjointInterpolator* G4AdjointInterpolator:: 43 G4AdjointInterpolator* G4AdjointInterpolator::GetInstance()
39 { << 44 { if(theInstance == 0) {
40 if(!fInstance) << 45 static G4AdjointInterpolator interpolator;
41 { << 46 theInstance = &interpolator;
42 fInstance = new G4AdjointInterpolator; <<
43 } 47 }
44 return fInstance; << 48 return theInstance;
45 } 49 }
46 50
47 ////////////////////////////////////////////// 51 ///////////////////////////////////////////////////////
48 G4AdjointInterpolator::G4AdjointInterpolator() << 52 //
49 << 53 G4AdjointInterpolator::G4AdjointInterpolator()
>> 54 {;
>> 55 }
50 ////////////////////////////////////////////// 56 ///////////////////////////////////////////////////////
51 G4AdjointInterpolator::~G4AdjointInterpolator( << 57 //
52 << 58 G4AdjointInterpolator::~G4AdjointInterpolator()
>> 59 {;
>> 60 }
53 ////////////////////////////////////////////// 61 ///////////////////////////////////////////////////////
54 G4double G4AdjointInterpolator::LinearInterpol << 62 //
55 << 63 G4double G4AdjointInterpolator::LinearInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
56 << 64 { G4double res = y1+ (x-x1)*(y2-y1)/(x2-x1);
57 { << 65 //G4cout<<"Linear "<<res<<G4endl;
58 G4double res = y1 + (x - x1) * (y2 - y1) / ( << 66 return res;
59 return res; <<
60 } 67 }
61 <<
62 ////////////////////////////////////////////// 68 ///////////////////////////////////////////////////////
63 G4double G4AdjointInterpolator::LogarithmicInt << 69 //
64 G4double& x, G4double& x1, G4double& x2, G4d << 70 G4double G4AdjointInterpolator::LogarithmicInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
65 { << 71 { if (y1<=0 || y2<=0 || x1<=0) return LinearInterpolation(x,x1,x2,y1,y2);
66 if(y1 <= 0. || y2 <= 0. || x1 <= 0.) << 72 G4double B=std::log(y2/y1)/std::log(x2/x1);
67 return LinearInterpolation(x, x1, x2, y1, << 73 //G4cout<<"x1,x2,y1,y2 "<<x1<<'\t'<<x2<<'\t'<<y1<<'\t'<<y2<<'\t'<<G4endl;
68 G4double B = std::log(y2 / y1) / std::log( << 74 G4double A=y1/std::pow(x1,B);
69 G4double A = y1 / std::pow(x1, B); << 75 G4double res=A*std::pow(x,B);
70 G4double res = A * std::pow(x, B); << 76 // G4cout<<"Log "<<res<<G4endl;
71 return res; 77 return res;
72 } 78 }
73 <<
74 ////////////////////////////////////////////// 79 ///////////////////////////////////////////////////////
75 G4double G4AdjointInterpolator::ExponentialInt << 80 //
76 G4double& x, G4double& x1, G4double& x2, G4d << 81 G4double G4AdjointInterpolator::ExponentialInterpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2)
77 { << 82 { G4double B=(std::log(y2)-std::log(y1));
78 G4double B = (std::log(y2) - std::log(y1)) << 83 B=B/(x2-x1);
79 G4double A = y1 * std::exp(-B * x1); << 84 G4double A=y1*std::exp(-B*x1);
80 G4double res = A * std::exp(B * x); << 85 G4double res=A*std::exp(B*x);
81 return res; 86 return res;
82 } 87 }
83 <<
84 ////////////////////////////////////////////// 88 ///////////////////////////////////////////////////////
85 G4double G4AdjointInterpolator::Interpolation( << 89 //
86 << 90 G4double G4AdjointInterpolator::Interpolation(G4double& x,G4double& x1,G4double& x2,G4double& y1,G4double& y2,G4String InterPolMethod)
87 <<
88 <<
89 { 91 {
90 if(InterPolMethod == "Log") << 92 if (InterPolMethod == "Log" ){
91 { << 93 return LogarithmicInterpolation(x,x1,x2,y1,y2);
92 return LogarithmicInterpolation(x, x1, x2, << 94 }
93 } << 95 else if (InterPolMethod == "Lin" ){
94 else if(InterPolMethod == "Lin") << 96 return LinearInterpolation(x,x1,x2,y1,y2);
95 { << 97 }
96 return LinearInterpolation(x, x1, x2, y1, << 98 else if (InterPolMethod == "Exp" ){
97 } << 99 return ExponentialInterpolation(x,x1,x2,y1,y2);
98 else if(InterPolMethod == "Exp") << 100 }
99 { << 101 else {
100 return ExponentialInterpolation(x, x1, x2, << 102 //G4cout<<"The interpolation method that you invoked does not exist!"<<G4endl;
101 } << 103 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 } 104 }
110 } 105 }
111 <<
112 ////////////////////////////////////////////// 106 ///////////////////////////////////////////////////////
113 // only valid if x_vec is monotically increasi << 107 //
114 std::size_t G4AdjointInterpolator::FindPositio << 108 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 << 109 { //most rapid nethod could be used probably
116 std << 110 //It is important to put std::vector<G4double>& such that the vector itself is used and not a copy
117 { << 111
118 // most rapid method could be used probably << 112
119 << 113 size_t ndim = x_vec.size();
120 std::size_t ndim = x_vec.size(); << 114 size_t ind1 = 0;
121 std::size_t ind1 = 0; << 115 size_t ind2 = ndim - 1;
122 std::size_t ind2 = ndim - 1; << 116 /* if (ind_max >= ind_min){
123 << 117 ind1=ind_min;
124 if(ndim > 1) << 118 ind2=ind_max;
125 { << 119
126 if(x_vec[0] < x_vec[1]) << 120
127 { // increasing << 121 }
128 do << 122 */
129 { << 123
130 std::size_t midBin = (ind1 + ind2) / 2 << 124
131 if(x < x_vec[midBin]) << 125 if (ndim >1) {
132 ind2 = midBin; << 126
133 else << 127 if (x_vec[0] < x_vec[1] ) { //increasing
134 ind1 = midBin; << 128 do {
135 } while(ind2 - ind1 > 1); << 129 size_t midBin = (ind1 + ind2)/2;
136 } << 130 if (x < x_vec[midBin])
137 else << 131 ind2 = midBin;
138 { << 132 else
139 do << 133 ind1 = midBin;
140 { << 134 } while (ind2 - ind1 > 1);
141 std::size_t midBin = (ind1 + ind2) / 2 << 135 }
142 if(x < x_vec[midBin]) << 136 else {
143 ind1 = midBin; << 137 do {
144 else << 138 size_t midBin = (ind1 + ind2)/2;
145 ind2 = midBin; << 139 if (x < x_vec[midBin])
146 } while(ind2 - ind1 > 1); << 140 ind1 = midBin;
147 } << 141 else
>> 142 ind2 = midBin;
>> 143 } while (ind2 - ind1 > 1);
>> 144 }
>> 145
148 } 146 }
149 << 147
150 return ind1; 148 return ind1;
151 } 149 }
152 150
153 ////////////////////////////////////////////// 151 ///////////////////////////////////////////////////////
154 // only valid if x_vec is monotically increasi << 152 //
155 std::size_t G4AdjointInterpolator::FindPositio << 153 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_ << 154 { //most rapid nethod could be used probably
157 { << 155 //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); 156 return FindPosition(log_x, log_x_vec);
160 } << 157 /*
161 << 158 if (log_x_vec.size()>3){
162 ////////////////////////////////////////////// << 159 size_t ind=0;
163 G4double G4AdjointInterpolator::Interpolate(G4 << 160 G4double log_x1=log_x_vec[1];
164 st << 161 G4double d_log =log_x_vec[2]-log_x1;
165 st << 162 G4double dind=(log_x-log_x1)/d_log +1.;
166 co << 163 if (dind <1.) ind=0;
167 { << 164 else if (dind >= double(log_x_vec.size())-2.) ind =log_x_vec.size()-2;
168 std::size_t i = FindPosition(x, x_vec); << 165 else ind =size_t(dind);
169 return Interpolation(x, x_vec[i], x_vec[i + << 166 return ind;
170 InterPolMethod); << 167
171 } << 168 }
172 << 169 else return FindPosition(log_x, log_x_vec);
173 ////////////////////////////////////////////// << 170 */
174 G4double G4AdjointInterpolator::InterpolateWit << 171
175 G4double& x, std::vector<G4double>& x_vec, s << 172
176 std::vector<std::size_t>& index_vec, G4doubl << 173 }
177 G4double dx) // only linear interpolation p << 174 ///////////////////////////////////////////////////////
178 { << 175 //
179 std::size_t ind = 0; << 176 G4double G4AdjointInterpolator::Interpolate(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,G4String InterPolMethod)
180 if(x > x0) << 177 { size_t i=FindPosition(x,x_vec);
181 ind = G4int((x - x0) / dx); << 178 //G4cout<<i<<G4endl;
182 if(ind >= index_vec.size() - 1) << 179 //G4cout<<x<<G4endl;
183 ind = index_vec.size() - 2; << 180 //G4cout<<x_vec[i]<<G4endl;
184 std::size_t ind1 = index_vec[ind]; << 181 return Interpolation( x,x_vec[i],x_vec[i+1],y_vec[i],y_vec[i+1],InterPolMethod);
185 std::size_t ind2 = index_vec[ind + 1]; << 182 }
186 if(ind1 > ind2) << 183
187 { << 184 ///////////////////////////////////////////////////////
188 std::size_t ind11 = ind1; << 185 //
189 ind1 = ind2; << 186 G4double G4AdjointInterpolator::InterpolateWithIndexVector(G4double& x,std::vector<G4double>& x_vec,std::vector<G4double>& y_vec,
190 ind2 = ind11; << 187 std::vector<size_t>& index_vec,G4double x0, G4double dx) //only linear interpolation possible
191 } << 188 { size_t ind=0;
192 ind = FindPosition(x, x_vec, ind1, ind2); << 189 if (x>x0) ind=int((x-x0)/dx);
193 return Interpolation(x, x_vec[ind], x_vec[in << 190 if (ind >= index_vec.size()-1) ind= index_vec.size()-2;
194 y_vec[ind + 1], "Lin"); << 191 size_t ind1 = index_vec[ind];
195 } << 192 size_t ind2 = index_vec[ind+1];
196 << 193 if (ind1 >ind2) {
197 ////////////////////////////////////////////// << 194 size_t ind11=ind1;
198 G4double G4AdjointInterpolator::InterpolateFor << 195 ind1=ind2;
199 G4double& log_x, std::vector<G4double>& log_ << 196 ind2=ind11;
200 std::vector<G4double>& log_y_vec) << 197
201 { << 198 }
202 std::size_t i = FindPositionForLogVector(log << 199
203 << 200 ind=FindPosition(x,x_vec,ind1,ind2);
204 G4double log_y = LinearInterpolation(log_x, << 201 return Interpolation( x,x_vec[ind],x_vec[ind+1],y_vec[ind],y_vec[ind+1],"Lin");
205 log_y_v << 202
206 return log_y; << 203 }
207 } << 204
>> 205
>> 206 ///////////////////////////////////////////////////////
>> 207 //
>> 208 G4double G4AdjointInterpolator::InterpolateForLogVector(G4double& log_x,std::vector<G4double>& log_x_vec,std::vector<G4double>& log_y_vec)
>> 209 { //size_t i=0;
>> 210 size_t i=FindPositionForLogVector(log_x,log_x_vec);
>> 211 /*G4cout<<"In interpolate "<<G4endl;
>> 212 G4cout<<i<<G4endl;
>> 213 G4cout<<log_x<<G4endl;
>> 214 G4cout<<log_x_vec[i]<<G4endl;
>> 215 G4cout<<log_x_vec[i+1]<<G4endl;
>> 216 G4cout<<log_y_vec[i]<<G4endl;
>> 217 G4cout<<log_y_vec[i+1]<<G4endl;*/
>> 218
>> 219 G4double log_y=LinearInterpolation(log_x,log_x_vec[i],log_x_vec[i+1],log_y_vec[i],log_y_vec[i+1]);
>> 220 return log_y;
>> 221
>> 222 }
208 223