Line data Source code
1 : /**
2 : Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
3 : Initial version January 2010.
4 : Maintained by ESO since 2013.
5 :
6 : This file is part of LibAIR and is licensed under GNU Public
7 : License Version 2
8 :
9 : \file dtdlcoeffs.cpp
10 : Renamed dtdlcoeffs.cc 2023
11 :
12 : */
13 :
14 : #include <algorithm>
15 : #include <numeric>
16 :
17 : #include <cmath>
18 : #include <iostream>
19 : #include <stdio.h> /*** for sprintf(...) ***/
20 :
21 : #include "dtdlcoeffs.h"
22 : #include "almaresults.h"
23 :
24 : namespace LibAIR2 {
25 :
26 0 : dTdLCoeffsBase::dTdLCoeffsBase()
27 : {
28 0 : std::array<double, 4> _chmask={1,1,1,1};
29 0 : chmask=_chmask;
30 0 : }
31 :
32 0 : dTdLCoeffsBase::~dTdLCoeffsBase()
33 : {
34 0 : }
35 :
36 0 : dTdLCoeffsSingle::dTdLCoeffsSingle(const std::vector<double> &c,
37 0 : const std::vector<double> &e):
38 0 : c(c),
39 0 : c2(4, 0.0),
40 0 : e(e)
41 : {
42 0 : }
43 :
44 0 : dTdLCoeffsSingle::dTdLCoeffsSingle(const ALMAResBase &r):
45 0 : c(4),
46 0 : c2(4, 0.0),
47 0 : e(4)
48 : {
49 0 : for(size_t i=0; i<4; ++i)
50 : {
51 0 : c[i]=r.dTdL[i];
52 0 : e[i]=r.dTdL_err[i];
53 : }
54 0 : }
55 :
56 0 : void dTdLCoeffsSingle::get(size_t i,
57 : double time,
58 : double el,
59 : std::vector<double> &res,
60 : std::vector<double> &res2) const
61 : {
62 : (void) i; // use not yet implemented
63 : (void) time;
64 : (void) el;
65 0 : res.resize(c.size());
66 0 : std::copy(c.begin(),
67 : c.end(),
68 : res.begin());
69 :
70 0 : for(size_t j=0; j<4; ++j)
71 0 : res[j]*=chmask[j];
72 :
73 0 : res2.resize(c2.size());
74 0 : std::copy(c2.begin(),
75 : c2.end(),
76 : res2.begin());
77 :
78 0 : }
79 :
80 0 : void dTdLCoeffsSingle::print(std::ostream &os)
81 : {
82 0 : os<<"Coefficients: ";
83 0 : for(size_t i=0; i<4; ++i)
84 0 : os<<c[i]<<", ";
85 0 : os<<std::endl;
86 0 : }
87 :
88 0 : void dTdLCoeffsSingle::repr(std::vector<double> &res,
89 : std::vector<double> &err) const
90 : {
91 0 : res=c;
92 0 : err=e;
93 0 : }
94 :
95 :
96 0 : bool dTdLCoeffsSingle::isnan(void) const
97 : {
98 0 : return std::accumulate(c.begin(),c.end(),false,[](bool acc, double v) { return acc || std::isnan(v); });
99 : }
100 :
101 0 : dTdLCoeffsIndiv::dTdLCoeffsIndiv(const coeff_t &c):
102 0 : coeff(c)
103 : {
104 0 : }
105 :
106 0 : dTdLCoeffsIndiv::dTdLCoeffsIndiv(size_t nAnts):
107 0 : coeff(4,nAnts,3,0.)
108 : {
109 0 : }
110 :
111 0 : void dTdLCoeffsIndiv::set(size_t i,
112 : const std::vector<double> &c,
113 : const std::vector<double> &e)
114 : {
115 0 : for(size_t k=0; k<4; ++k)
116 : {
117 0 : coeff(k,i,0)=c[k];
118 0 : coeff(k,i,1)=e[k];
119 : }
120 0 : }
121 :
122 0 : void dTdLCoeffsIndiv::get(size_t i,
123 : double time,
124 : double el,
125 : std::vector<double> &res,
126 : std::vector<double> &res2) const
127 : {
128 : (void) time; // use not yet implemented
129 : (void) el;
130 :
131 0 : const size_t nc=4; //coeff.shape()[0];
132 0 : res.resize(nc);
133 0 : res2.resize(nc);
134 0 : for(size_t k=0; k<nc; ++k)
135 : {
136 0 : res[k]=coeff(k,i,0)*chmask[k];
137 0 : res2[k]=coeff(k,i,2);
138 : }
139 0 : }
140 :
141 0 : void dTdLCoeffsIndiv::print(std::ostream &os)
142 : {
143 0 : os<<"Coefficients: ";
144 0 : for(size_t i=0; i<4; ++i)
145 0 : os<<coeff(i,0,0)<<", ";
146 0 : os<<std::endl;
147 0 : }
148 :
149 0 : void dTdLCoeffsIndiv::repr(std::vector<double> &res,
150 : std::vector<double> &err) const
151 : {
152 0 : const size_t nc=4; //coeff.shape()[0];
153 0 : res.resize(nc);
154 0 : err.resize(nc);
155 0 : for(size_t k=0; k<nc; ++k)
156 : {
157 0 : res[k]=coeff(k,0,0);
158 0 : err[k]=coeff(k,0,1);
159 : }
160 0 : }
161 :
162 0 : bool dTdLCoeffsIndiv::isnan(void) const
163 : {
164 : // return True if one of the elements of coeff is a NaN
165 0 : if (coeff.contiguousStorage()) {
166 0 : return std::accumulate(coeff.cbegin(),
167 : coeff.cend(),
168 : false,
169 0 : [](bool acc, double v) { return acc || std::isnan(v); }
170 0 : );
171 : } else {
172 0 : return std::accumulate(coeff.begin(),
173 0 : coeff.end(),
174 : false,
175 0 : [](bool acc, double v) { return acc || std::isnan(v); }
176 0 : );
177 : }
178 :
179 : //return std::accumulate(coeff.origin(),
180 : // coeff.origin()+(coeff.shape()[0]*coeff.shape()[1]*coeff.shape()[2]),
181 : // false,
182 : // [](bool acc, double v) { return acc || std::isnan(v); }
183 : // );
184 : }
185 :
186 :
187 0 : dTdLCoeffsSingleInterpolated::dTdLCoeffsSingleInterpolated()
188 : {
189 0 : }
190 :
191 0 : void dTdLCoeffsSingleInterpolated::insert(double time,
192 : const std::array<double, 4> &coeffs,
193 : const std::array<double, 4> &err)
194 : {
195 : ret_t t;
196 0 : t.time=time;
197 0 : t.coeffs=coeffs;
198 0 : t.err=err;
199 0 : retrievals.insert(t);
200 0 : }
201 :
202 0 : void dTdLCoeffsSingleInterpolated::get(size_t i,
203 : double time,
204 : double el,
205 : std::vector<double> &res,
206 : std::vector<double> &c2) const
207 : {
208 : (void) i; // use not yet implemented
209 : (void) el;
210 :
211 0 : if (retrievals.size() == 0)
212 : {
213 0 : throw std::runtime_error("No retrievals have been entered.");
214 : }
215 :
216 0 : res.resize(4);
217 0 : c2.resize(4);
218 0 : if (time<retrievals.begin()->time)
219 : {
220 : // request is before the first retrieval
221 0 : std::copy(retrievals.begin()->coeffs.begin(),
222 0 : retrievals.begin()->coeffs.end(),
223 : res.begin());
224 0 : for(size_t j=0; j<4; ++j)
225 0 : res[j]*=chmask[j];
226 0 : std::copy(retrievals.begin()->c2.begin(),
227 0 : retrievals.begin()->c2.end(),
228 : c2.begin());
229 0 : return;
230 : }
231 0 : std::set<ret_t>::const_iterator prev;
232 0 : std::set<ret_t>::const_iterator next;
233 :
234 0 : for(next=retrievals.begin();
235 0 : next !=retrievals.end() and time >= next->time;
236 0 : ++next)
237 : {
238 : }
239 :
240 0 : if (next ==retrievals.end())
241 : {
242 : // requested time is after the last retrieval
243 0 : std::set<ret_t>::reverse_iterator last=retrievals.rbegin();
244 0 : std::copy(last->coeffs.begin(),
245 0 : last->coeffs.end(),
246 : res.begin());
247 0 : for(size_t j=0; j<4; ++j)
248 0 : res[j]*=chmask[j];
249 0 : std::copy(last->c2.begin(),
250 0 : last->c2.end(),
251 : c2.begin());
252 0 : return;
253 : }
254 :
255 0 : prev=next;
256 0 : prev--;
257 :
258 :
259 : //Finally if not before beginning or after end, interpolate
260 : //linearly
261 :
262 0 : double c_f=(time - prev->time)/(next->time - prev->time);
263 0 : double c_s=(next->time - time)/(next->time - prev->time);
264 :
265 0 : for(size_t k=0; k<4; ++k)
266 : {
267 0 : res[k]= c_s*(prev->coeffs[k]) + c_f*(next->coeffs[k]);
268 0 : res[k]*= chmask[k];
269 0 : c2[k]= c_s*(prev->c2[k]) + c_f*(next->c2[k]);
270 : }
271 : }
272 :
273 0 : void dTdLCoeffsSingleInterpolated::print(std::ostream &os)
274 : {
275 0 : for(std::set<ret_t>::const_iterator i=retrievals.begin();
276 0 : i!=retrievals.end();
277 0 : ++i)
278 : {
279 : char tstr[120];
280 0 : sprintf(tstr, "Retrieval at %f is: [", i->time);
281 0 : os << tstr;
282 0 : for(size_t j=0; j<4; ++j)
283 : {
284 0 : os << i->coeffs[j] <<",";
285 : }
286 0 : os << "]"
287 0 : << std::endl;
288 : }
289 0 : }
290 :
291 0 : void dTdLCoeffsSingleInterpolated::repr(std::vector<double> &res,
292 : std::vector<double> &err) const
293 : {
294 0 : std::set<ret_t>::const_iterator i=retrievals.begin();
295 : size_t j;
296 0 : for(j=0; j < retrievals.size()/2;++j)
297 : {
298 0 : ++i;
299 : };
300 0 : while(j < retrievals.size()-1) // make sure we have not hit an invalid coefficent
301 : {
302 0 : if(i->coeffs[0]>0.)
303 : {
304 0 : break;
305 : }
306 0 : ++j;
307 0 : ++i;
308 : };
309 :
310 0 : res.resize(4);
311 0 : err.resize(4);
312 0 : for(j=0; j<4; ++j)
313 : {
314 0 : res[j]=i->coeffs[j];
315 0 : err[j]=i->err[j];
316 : }
317 0 : }
318 :
319 0 : bool dTdLCoeffsSingleInterpolated::isnan() const
320 : {
321 0 : for(std::set<ret_t>::const_iterator i=retrievals.begin();
322 0 : i != retrievals.end();
323 0 : ++i)
324 : {
325 0 : for(size_t j=0; j<4; ++j)
326 : {
327 0 : if ( std::isnan(i->coeffs[j]) || (std::isnan(i->err[j])) || std::isnan(i->c2[j]) )
328 : {
329 0 : return true;
330 : }
331 : }
332 : }
333 0 : return false;
334 : }
335 :
336 : }
337 :
338 :
|
