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