Mantid
Loading...
Searching...
No Matches
GetTimeSeriesLogInformation.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8#include "MantidAPI/Run.h"
15#include "MantidHistogramData/Histogram.h"
18#include <algorithm>
19#include <fstream>
20
21using namespace Mantid::Kernel;
22using namespace Mantid::API;
23using namespace Mantid::DataObjects;
24using namespace Mantid::HistogramData;
25using Mantid::Types::Core::DateAndTime;
26
27using namespace std;
28
29namespace Mantid::Algorithms {
30
31DECLARE_ALGORITHM(GetTimeSeriesLogInformation)
32
33
36 : API::Algorithm(), m_dataWS(), mRunStartTime(), mFilterT0(), mFilterTf(), m_intInfoMap(), m_dblInfoMap(),
37 m_log(nullptr), m_timeVec(), m_valueVec(), m_starttime(), m_endtime(), m_ignoreNegativeTime(false) {}
38
43 std::make_unique<API::WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "Anonymous", Direction::InOut),
44 "Input EventWorkspace. Each spectrum corresponds to 1 pixel");
45
46 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "Dummy", Direction::Output),
47 "Name of the workspace of log delta T distribution. ");
48
49 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("InformationWorkspace", "", Direction::Output),
50 "Name of optional log statistic information output Tableworkspace.");
51
52 declareProperty("LogName", "", "Log's name to filter events.");
53
54 std::vector<std::string> timeoptions{"Absolute Time (nano second)", "Relative Time (second)"};
55 declareProperty("TimeRangeOption", "Relative Time (second)", std::make_shared<StringListValidator>(timeoptions),
56 "User defined time range (T0, Tf) is of absolute time (second). ");
57
58 declareProperty("FilterStartTime", EMPTY_DBL(),
59 "Earliest time of the events to be selected. "
60 "It can be absolute time (ns), relative time (second) or percentage.");
61 declareProperty("FilterStopTime", EMPTY_DBL(),
62 "Latest time of the events to be selected. "
63 "It can be absolute time (ns), relative time (second) or percentage.");
64
65 declareProperty("TimeStepBinResolution", 0.0001, "Time resolution (second) for time stamp delta T disibution. ");
66
67 declareProperty("IgnoreNegativeTimeInterval", false,
68 "If true, then the time interval with negative number will "
69 "be neglected in doing statistic.");
70}
71
75 // 1. Get wrokspace, log property and etc.
76 m_dataWS = this->getProperty("InputWorkspace");
77 if (!m_dataWS) {
78 throw runtime_error("Inputworkspace cannot be parsed to a MatrixWorkspace.");
79 }
80
81 string logname = getProperty("LogName");
82 if (logname.empty())
83 throw runtime_error("Input log value cannot be an empty string. ");
84
85 Kernel::Property *log = m_dataWS->run().getProperty(logname);
86 if (!log) {
87 stringstream errmsg;
88 errmsg << "Property " << logname << " does not exit in sample of workspace " << m_dataWS->getName() << ".";
89 g_log.error(errmsg.str());
90 throw std::invalid_argument(errmsg.str());
91 } else {
92 m_log = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(log);
93 if (!m_log) {
94 stringstream errmsg;
95 errmsg << "Log " << logname << " is found, but is not a double type time series log";
96 g_log.error(errmsg.str());
97 throw std::invalid_argument(errmsg.str());
98 }
99 }
100
103
104 m_starttime = m_dataWS->run().startTime();
105 m_endtime = m_dataWS->run().endTime();
106
107 m_ignoreNegativeTime = getProperty("IgnoreNegativeTimeInterval");
108
109 // 2. Process start time and end time
111
112 // 3. Check time stamps
114
115 // 4. Calculate distribution of delta T
116 double resolution = getProperty("TimeStepBinResolution");
117 Workspace2D_sptr timestatws = calDistributions(m_timeVec, resolution);
118
119 // 5. Check whether the log is alternating
121
122 // 6. Export error log
123 if (false) {
124 double userinputdt = 1 / 240.1;
125 exportErrorLog(m_dataWS, m_timeVec, userinputdt);
126 }
127
128 // -1. Finish set output properties.
130 setProperty("InformationWorkspace", infows);
131
132 this->setProperty("OutputWorkspace", std::dynamic_pointer_cast<MatrixWorkspace>(timestatws));
133
134 // 4. Do more staticts (examine)
135 // std::string outputdir = this->getProperty("OutputDirectory");
136 // examLog(logname, outputdir);
137}
138
143 // Orignal
144 m_intInfoMap.emplace("Items", m_log->size());
145
146 // Input time
147 double t0r = this->getProperty("FilterStartTime");
148 double tfr = this->getProperty("FilterStopTime");
149
150 // Time unit option
151 string timeoption = this->getProperty("TimeRangeOption");
152 int timecase = 0;
153 if (timeoption == "Absolute Time (nano second)")
154 timecase = 1;
155 else if (timeoption == "Relative Time (second)")
156 timecase = 0;
157 else
158 timecase = -1;
159
160 double duration = static_cast<double>(m_timeVec.back().totalNanoseconds() - m_timeVec[0].totalNanoseconds()) * 1.0E-9;
161
162 // Process start time
163 if (t0r == EMPTY_DBL()) {
164 // Default
165 mFilterT0 = m_timeVec[0];
166 } else {
167 switch (timecase) {
168 case 0:
169 // Relative time (second)
171 break;
172 case 1:
173 // Absolute time (nano second)
175 break;
176 case -1:
177 mFilterT0 = calculateRelativeTime(t0r * duration);
178 break;
179 default:
180 throw runtime_error("Coding error!");
181 break;
182 }
183 } // Filter start time
184
185 // Process filter end time
186 if (tfr == EMPTY_DBL()) {
187 // Default
189 } else {
190 switch (timecase) {
191 // Set with input
192 case 0:
193 // Relative time (second)
195 break;
196 case 1:
197 // Absolute time (nano second)
199 break;
200 case -1:
201 mFilterTf = calculateRelativeTime(tfr * duration);
202 break;
203 default:
204 throw runtime_error("Coding error!");
205 break;
206 }
207 }
208
209 // Check validity
210 if (mFilterTf.totalNanoseconds() <= mFilterT0.totalNanoseconds()) {
211 stringstream errmsg;
212 errmsg << "User defined filter starting time @ " << mFilterT0 << " (T = " << t0r
213 << ") is later than or equal to filer ending time @ " << mFilterTf << " (T = " << tfr << ").";
214 g_log.error(errmsg.str());
215 throw std::invalid_argument(errmsg.str());
216 }
217}
218
223Types::Core::DateAndTime GetTimeSeriesLogInformation::getAbsoluteTime(double abstimens) {
224 DateAndTime temptime(static_cast<int64_t>(abstimens));
225
226 return temptime;
227}
228
232Types::Core::DateAndTime GetTimeSeriesLogInformation::calculateRelativeTime(double deltatime) {
233 int64_t totaltime = m_starttime.totalNanoseconds() + static_cast<int64_t>(deltatime * 1.0E9);
234 DateAndTime abstime(totaltime);
235
236 return abstime;
237}
238
242 auto tablews = std::make_shared<TableWorkspace>();
243
244 tablews->addColumn("str", "Name");
245 tablews->addColumn("double", "Value");
246
247 // 1. Integer part
248 for (auto &intmapiter : m_intInfoMap) {
249 string name = intmapiter.first;
250 size_t value = intmapiter.second;
251
252 TableRow newrow = tablews->appendRow();
253 newrow << name << static_cast<double>(value);
254 }
255
256 // 2. Double part
257 map<string, double>::iterator dblmapiter;
258 for (dblmapiter = m_dblInfoMap.begin(); dblmapiter != m_dblInfoMap.end(); ++dblmapiter) {
259 string name = dblmapiter->first;
260 double value = dblmapiter->second;
261
262 TableRow newrow = tablews->appendRow();
263 newrow << name << value;
264 }
265
266 return tablews;
267}
268
277void GetTimeSeriesLogInformation::exportErrorLog(const MatrixWorkspace_sptr &ws, vector<DateAndTime> abstimevec,
278 double dts) {
279 std::string outputdir = getProperty("OutputDirectory");
280 if (!outputdir.empty() && outputdir.back() != '/')
281 outputdir += "/";
282
283 std::string ofilename = outputdir + "errordeltatime.txt";
284 g_log.notice() << ofilename << '\n';
285 std::ofstream ofs;
286 ofs.open(ofilename.c_str(), std::ios::out);
287
288 Types::Core::DateAndTime t0(ws->run().getProperty("run_start")->value());
289
290 for (size_t i = 1; i < abstimevec.size(); i++) {
291 double tempdts =
292 static_cast<double>(abstimevec[i].totalNanoseconds() - abstimevec[i - 1].totalNanoseconds()) * 1.0E-9;
293 double dev = (tempdts - dts) / dts;
294
295 if (fabs(dev) > 0.5) {
296 double deltapulsetimeSec1 =
297 static_cast<double>(abstimevec[i - 1].totalNanoseconds() - t0.totalNanoseconds()) * 1.0E-9;
298 double deltapulsetimeSec2 =
299 static_cast<double>(abstimevec[i].totalNanoseconds() - t0.totalNanoseconds()) * 1.0E-9;
300 auto index1 = static_cast<int>(deltapulsetimeSec1 * 60);
301 auto index2 = static_cast<int>(deltapulsetimeSec2 * 60);
302
303 ofs << "Error d(T) = " << tempdts << " vs Correct d(T) = " << dts << '\n';
304 ofs << index1 << "\t\t" << abstimevec[i - 1].totalNanoseconds() << "\t\t" << index2 << "\t\t"
305 << abstimevec[i].totalNanoseconds() << '\n';
306 }
307 }
308
309 ofs.close();
310}
311
318Workspace2D_sptr GetTimeSeriesLogInformation::calDistributions(std::vector<Types::Core::DateAndTime> timevec,
319 double stepsize) {
320 // 1. Get a vector of delta T (in unit of seconds)
321 double dtmin = static_cast<double>(timevec.back().totalNanoseconds() - timevec[0].totalNanoseconds()) * 1.0E-9;
322 double dtmax = 0.0;
323
324 vector<double> vecdt(timevec.size() - 1, 0.0);
325 for (size_t i = 1; i < timevec.size(); ++i) {
326 vecdt[i - 1] = static_cast<double>(timevec[i].totalNanoseconds() - timevec[i - 1].totalNanoseconds()) * 1.0E-9;
327 if (vecdt[i - 1] < dtmin)
328 dtmin = vecdt[i - 1];
329 else if (vecdt[i - 1] > dtmax)
330 dtmax = vecdt[i - 1];
331 }
332
333 // 2. Create a vector of counts
334 size_t numbins;
335 if (m_ignoreNegativeTime && dtmin < 0) {
336 numbins = static_cast<size_t>(ceil((dtmax) / stepsize)) + 2;
337 } else {
338 numbins = static_cast<size_t>(ceil((dtmax - dtmin) / stepsize)) + 2;
339 }
340
341 g_log.notice() << "Distribution has " << numbins << " bins. Delta T = (" << dtmin << ", " << dtmax << ")\n";
342
343 Workspace2D_sptr distws = create<Workspace2D>(1, Points(numbins));
344 auto &vecDeltaT = distws->mutableX(0);
345 auto &vecCount = distws->mutableY(0);
346
347 double countmin = dtmin;
348 if (m_ignoreNegativeTime && dtmin < 0)
349 countmin = 0;
350
351 for (size_t i = 0; i < numbins; ++i)
352 vecDeltaT[i] = countmin + (static_cast<double>(i) - 1) * stepsize;
353 for (size_t i = 0; i < numbins; ++i)
354 vecCount[i] = 0;
355
356 // 3. Count
357 for (double dt : vecdt) {
358 int index;
359 if (dt < 0 && m_ignoreNegativeTime) {
360 index = 0;
361 } else {
362 auto viter = lower_bound(vecDeltaT.begin(), vecDeltaT.end(), dt);
363 index = static_cast<int>(viter - vecDeltaT.begin());
364 if (index >= static_cast<int>(vecDeltaT.size())) {
365 // Out of upper boundary
366 g_log.error() << "Find index = " << index << " > vecX.size = " << vecDeltaT.size() << ".\n";
367 } else if (dt < vecDeltaT[index]) {
368 --index;
369 }
370
371 if (index < 0)
372 throw runtime_error("How can this happen.");
373 }
374 vecCount[index] += 1;
375 }
376
377 return distws;
378}
379
383 // 1. Time correctness: same time, disordered time
384 size_t countsame = 0;
385 size_t countinverse = 0;
386 for (size_t i = 1; i < m_timeVec.size(); i++) {
387 Types::Core::DateAndTime tprev = m_timeVec[i - 1];
388 Types::Core::DateAndTime tpres = m_timeVec[i];
389 if (tprev == tpres)
390 countsame++;
391 else if (tprev > tpres)
392 countinverse++;
393 }
394
395 // Written to summary map
396 /*
397 Types::Core::time_duration dts = m_timeVec[0]-m_starttime;
398 Types::Core::time_duration dtf = m_timeVec.back() - m_timeVec[0];
399 size_t f = m_timeVec.size()-1;
400 */
401
402 m_intInfoMap.emplace("Number of Time Stamps", m_timeVec.size());
403 m_intInfoMap.emplace("Number of Equal Time Stamps", countsame);
404 m_intInfoMap.emplace("Number of Reversed Time Stamps", countinverse);
405
406 // 2. Average and standard deviation (delta t)
407 double runduration_sec = static_cast<double>(m_endtime.totalNanoseconds() - m_starttime.totalNanoseconds()) * 1.0E-9;
408 if (runduration_sec < 0.0) {
409 g_log.warning() << "It shows that the run duration is not right. "
410 << "Run start = " << m_starttime.toFormattedString() << "; "
411 << "Run End = " << m_endtime.toFormattedString() << ".\n";
412 g_log.notice() << "Log start time = " << m_timeVec[0].toFormattedString() << "; "
413 << "Log end time = " << m_timeVec.back().toFormattedString() << ".\n";
414 }
415
416 double sum_deltaT1 = 0.0; // sum(dt ) in second
417 double sum_deltaT2 = 0.0; // sum(dt^2) in second^2
418 double max_dt = 0;
419 double min_dt = 0;
420 if (runduration_sec > 0)
421 min_dt = runduration_sec;
422
423 size_t numpts = m_timeVec.size();
424 for (size_t i = 0; i < numpts - 1; ++i) {
425 int64_t dtns = m_timeVec[i + 1].totalNanoseconds() - m_timeVec[i].totalNanoseconds();
426 double dt = static_cast<double>(dtns) * 1.0E-9; // in second
427
428 if (dt < 0) {
429 g_log.warning() << "Reversed dT: dt = " << dt << " between " << m_timeVec[i].toFormattedString() << " and "
430 << m_timeVec[i + 1].toFormattedString() << " @ index = " << i << ".\n";
431 }
432
433 sum_deltaT1 += dt;
434 sum_deltaT2 += dt * dt;
435
436 if (dt > max_dt)
437 max_dt = dt;
438 if (dt < min_dt)
439 min_dt = dt;
440 }
441
442 double avg_dt = sum_deltaT1 / static_cast<double>(numpts - 1);
443 double std_dt = sqrt(sum_deltaT2 / static_cast<double>(numpts - 1) - avg_dt * avg_dt);
444
445 m_dblInfoMap.emplace("Average(dT)", avg_dt);
446 m_dblInfoMap.emplace("Sigma(dt)", std_dt);
447 m_dblInfoMap.emplace("Min(dT)", min_dt);
448 m_dblInfoMap.emplace("Max(dT)", max_dt);
449
450 // 3. Count number of time intervals beyond 10% of deviation
451 /* Temporarily disabled
452 for (size_t i = 1; ; i ++)
453 {
454 int64_t dtns =
455 m_timeVec[i].totalNanoseconds()-m_timeVec[i-1].totalNanoseconds();
456 double dtms = static_cast<double>(dtns)*1.0E-3;
457
458 if (dtms > dtmsA10p)
459 numdtabove10p ++;
460 else if (dtms < dtmsB10p)
461 numdtbelow10p ++;
462
463 } // ENDFOR
464 */
465
466 // 4. Output
467 /* Temporily disabled
468 g_log.notice() << "Run Start = " << t0.totalNanoseconds() << '\n';
469 g_log.notice() << "First Log: " << "Absolute Time = " <<
470 m_timeVec[0].totalNanoseconds() << "(ns), "
471 << "Relative Time = " <<
472 DateAndTime::nanosecondsFromDuration(dts) << "(ns) \n";
473 g_log.notice() << "Last Log: " << "Absolute Time = " <<
474 m_timeVec[f].totalNanoseconds() << "(ns), "
475 << "Relative Time = " <<
476 DateAndTime::nanosecondsFromDuration(dtf) << "(ns) \n";
477 g_log.notice() << "Normal dt = " << numnormal << '\n';
478 g_log.notice() << "Zero dt = " << numsame << '\n';
479 g_log.notice() << "Negative dt = " << numinvert << '\n';
480 g_log.notice() << "Avg d(T) = " << dt << " seconds +/- " << stddt << ",
481 Frequency = " << 1.0/dt << '\n';
482 g_log.notice() << "d(T) (unit ms) is in range [" << mindtms << ", " << maxdtms
483 << "]"<< '\n';
484 g_log.notice() << "Number of d(T) 10% larger than average = " <<
485 numdtabove10p << '\n';
486 g_log.notice() << "Number of d(T) 10% smaller than average = " <<
487 numdtbelow10p << '\n';
488
489 g_log.notice() << "Size of timevec = " << m_timeVec.size() << '\n';
490 */
491}
492
501void GetTimeSeriesLogInformation::checkLogValueChanging(vector<DateAndTime> timevec, vector<double> values,
502 double delta) {
503 std::stringstream ss;
504 ss << "Alternating Threashold = " << delta << '\n';
505
506 size_t numchange = 0;
507 for (size_t i = 1; i < values.size(); i++) {
508 double tempdelta = values[i] - values[i - 1];
509 if (fabs(tempdelta) < delta) {
510 // Value are 'same'
511 numchange++;
512
513 // An error message
514 ss << "@ " << i << "\tDelta = " << tempdelta << "\t\tTime From " << timevec[i - 1].totalNanoseconds() << " to "
515 << timevec[i].totalNanoseconds() << '\n';
516 }
517 }
518
519 m_intInfoMap.insert(make_pair("Number of adjacent time stamp w/o value change", numchange));
520 g_log.debug() << ss.str();
521}
522
523} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
GetTimeSeriesLogInformation : Read a TimeSeries log and return some information required by users.
void processTimeRange()
Do statistic on user proposed range and examine the log inside the given time range.
Types::Core::DateAndTime calculateRelativeTime(double deltatime)
Calculate the time from a given relative time from run start.
std::vector< Types::Core::DateAndTime > m_timeVec
Types::Core::DateAndTime getAbsoluteTime(double abstimens)
Convert a value in nanosecond to DateAndTime.
void exportErrorLog(const API::MatrixWorkspace_sptr &ws, std::vector< Types::Core::DateAndTime > abstimevec, double dts)
Export time stamps looking erroreous.
const std::string name() const override
function to return a name of the algorithm, must be overridden in all algorithms
void checkLogValueChanging(std::vector< Types::Core::DateAndTime > timevec, std::vector< double > values, double delta)
Check whether log values are changing from 2 adjacent time stamps.
DataObjects::Workspace2D_sptr calDistributions(std::vector< Types::Core::DateAndTime > timevec, double stepsize)
Calcualte the distribution of delta T in time stamps.
DataObjects::TableWorkspace_sptr generateStatisticTable()
Generate statistic information table workspace.
void checkLogBasicInforamtion()
Check log in workspace including ... ...
void init() override
Definition of all input arguments.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
Base class for properties.
Definition: Property.h:94
A specialised Property class for holding a series of time-value pairs.
int size() const override
Returns the number of values at UNIQUE time intervals in the time series.
std::vector< TYPE > valuesAsVector() const
Return the time series's values (unfiltered) as a vector<TYPE>
std::vector< Types::Core::DateAndTime > timesAsVector() const override
Return the time series's times as a vector<DateAndTime>
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Output
An output workspace.
Definition: Property.h:54