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