WSL/SLF GitLab Repository

WindowedFilter.cc 9.56 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
18
#include <meteoio/meteofilters/WindowedFilter.h>
19
20
21
22
23

using namespace std;

namespace mio {

24
WindowedFilter::WindowedFilter(const std::string& name)
25
	: FilterBlock(name), min_time_span(0.0, 0.), centering(WindowedFilter::center),
26
	  last_start(0), last_end(0), min_data_points(1), vec_window(), is_soft(false)
27
28
{}

29
30
unsigned int WindowedFilter::get_centering(std::vector<std::string>& vec_args)
{
31
	if (!vec_args.empty()){
32
		if (vec_args.front() == "left"){
33
34
			vec_args.erase(vec_args.begin());
			return WindowedFilter::left;
35
		} else if (vec_args.front() == "right"){
36
37
			vec_args.erase(vec_args.begin());
			return WindowedFilter::right;
38
		} else if (vec_args.front() == "center"){
39
40
41
42
			vec_args.erase(vec_args.begin());
			return WindowedFilter::center;
		}
	}
43

44
45
46
47
48
49
50
51
52
53
54
	return WindowedFilter::center; //the default
}

/**
 * @brief A function that cuts out the desired window for the 'index' element within
 *        ivec, the window elements are stored into vec_window
 *        Calls to this function have to start with index 0, then 1, 2, 3, ...
 *        vec_window is not allowed to be changed between two calls
 * @param index The index of the element in ivec that requires a window
 * @param ivec The original sequence of data points
 */
55
const std::vector<const MeteoData*>& WindowedFilter::get_window(const size_t& index,
56
                                                                const std::vector<MeteoData>& ivec)
57
{
58
	if(index==0) {
59
		vec_window.clear();
60
		vec_window.reserve(min_data_points*2); //to have enough margin for the time criteria
61
62
	}

63
64
65
66
67
68
	size_t start, end;
	if(!get_window_specs(index, ivec, start, end)) {
		vec_window.clear();
		vec_window.reserve(min_data_points*2); //to have enough margin for the time criteria
		vec_window.push_back( &ivec[index] );
	}
69

70
71
72
73
74
75
76
	if(index==0) {
		for(size_t ii=start; ii<=end; ii++) vec_window.push_back( &ivec[ii] );
	} else {
		if(last_start<start) {
			vec_window.erase( vec_window.begin(),  vec_window.begin()+(start-last_start));
		}
		if(last_start>start) {
77
			const std::vector<const MeteoData*>::iterator it = vec_window.begin();
78
79
80
81
82
83
84
			for(size_t ii=start; ii<=last_start; ii++) vec_window.insert(it, &ivec[ii]);
		}
		if(last_end<end) {
			for(size_t ii=last_end+1; ii<=end; ii++) vec_window.push_back( &ivec[ii] );
		}
		if(last_end>end) {
			vec_window.erase( vec_window.end()-(last_end-end), vec_window.end() );
85
86
87
		}
	}

88
89
90
	//save window metadata
	last_start = start;
	last_end = end;
91
92
93
	return vec_window;
}

94
95
96
97
98
99
100
101
/**
 * @brief A function that computes the start and end for a window for the 'index' element from ivec
 * @param index The index of the element in ivec that requires a window
 * @param ivec The original sequence of data points
 * @param start the start index of the window
 * @param end the end index of the window
 * @return true if success, false if a window could not be computed
 */
102
bool WindowedFilter::get_window_specs(const size_t& index, const std::vector<MeteoData>& ivec, size_t& start, size_t& end) const
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
{	/*
	The principle is too compute the first index that matches the minimum number of points criteria,
	the the one that matches the minimum time window,
	then combine them (with the equivalent of OR: we take the MIN index).
	Afterward, we compute the last index [...] for number of points
	and the last index [...] for the time window
	and combine them (with the equivalent of OR: we take the MIN index).
	(or vice versa for right centering)
	*/
	const Date date = ivec[index].date;
	start = end = index; //for proper initial value, so we can bail out without worries

	if(centering == WindowedFilter::left) {
		//get start of window
		size_t start_elements = min_data_points - 1; //start elements criteria
		if(start_elements>index) {
			if(!is_soft) return false;
			start_elements = index; //as many as possible
121
		}
122
123
		const Date start_date = date - min_time_span;
		size_t start_time_idx = IOUtils::seek(start_date, ivec, false); //start time criteria
124
		if(start_time_idx!=IOUtils::npos) start_time_idx = (start_time_idx>0)?start_time_idx-1:IOUtils::npos;
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
		if(start_time_idx==IOUtils::npos) {
			if(!is_soft) return false;
			start_time_idx=0; //first possible element
		}
		const size_t elements_left = MAX(index - start_time_idx, start_elements);
		start = index - elements_left;

		//get end of window
		if(!is_soft) return true; //with end=index
		size_t end_elements = (min_data_points>(elements_left+1))?min_data_points - (elements_left + 1):0;
		const Date end_date = ivec[start].date+min_time_span;
		size_t end_time_idx = (end_date>date)?IOUtils::seek(end_date, ivec, false):index; //end time criteria
		if(end_time_idx==IOUtils::npos) {
			if(!is_soft) return false;
			end_time_idx=ivec.size()-1; //last possible element
		}
		const size_t elements_right = MAX(end_time_idx - index, end_elements);
		end = index + elements_right;
143
144
	}

145
146
147
148
149
150
	if(centering == WindowedFilter::right) {
		//get end of window
		size_t end_elements = min_data_points - 1; //end elements criteria
		if(end_elements>(ivec.size()-1-index)) {
			if(!is_soft) return false;
			end_elements = (ivec.size()-1-index); //as many as possible
151
		}
152
153
154
155
156
		const Date end_date = date+min_time_span;
		size_t end_time_idx = IOUtils::seek(end_date, ivec, false); //end time criteria
		if(end_time_idx==IOUtils::npos) {
			if(!is_soft) return false;
			end_time_idx=ivec.size()-1; //last possible element
157
		}
158
159
160
161
162
163
164
165
		const size_t elements_right = MAX(end_time_idx - index, end_elements);
		end = index + elements_right;

		//get start of window
		if(!is_soft) return true; //with start=index
		size_t start_elements = (min_data_points>(elements_right+1))?min_data_points - (elements_right + 1):0;
		const Date start_date = ivec[end].date-min_time_span;
		size_t start_time_idx = (start_date<date)?IOUtils::seek(start_date, ivec, false):index; //start time criteria
166
		if(start_time_idx!=IOUtils::npos && start_time_idx!=index) start_time_idx = (start_time_idx>0)?start_time_idx-1:IOUtils::npos;
167
168
169
		if(start_time_idx==IOUtils::npos) {
			if(!is_soft) return false;
			end_time_idx=0; //first possible element
170
		}
171
172
		const size_t elements_left = MAX(index - start_time_idx, start_elements);
		start = index - elements_left;
173
174
	}

175
176
177
178
179
180
	if(centering == WindowedFilter::center) {
		//get start of ideal window
		size_t start_elements = min_data_points/2; //start elements criteria
		if(start_elements>index) {
			if(!is_soft) return false;
			start_elements = index; //as many as possible
181
		}
182
183
		const Date start_date = date - min_time_span/2;
		size_t start_time_idx = IOUtils::seek(start_date, ivec, false); //start time criteria
184
		if(start_time_idx!=IOUtils::npos) start_time_idx = (start_time_idx>0)?start_time_idx-1:IOUtils::npos;
185
186
187
		if(start_time_idx==IOUtils::npos) {
			if(!is_soft) return false;
			start_time_idx=0; //first possible element
188
		}
189
190
191
192
193
194
195
196
		const size_t elements_left = MAX(index - start_time_idx, start_elements);
		start = index - elements_left;

		//get end of ideal window
		size_t end_elements = min_data_points/2; //end elements criteria
		if(end_elements>(ivec.size()-1-index)) {
			if(!is_soft) return false;
			end_elements = (ivec.size()-1-index); //as many as possible
197
		}
198
199
200
201
202
203
204
205
206
207
208
209
210
211
		const Date end_date = date+min_time_span/2;
		size_t end_time_idx = IOUtils::seek(end_date, ivec, false); //end time criteria
		if(end_time_idx==IOUtils::npos) {
			if(!is_soft) return false;
			end_time_idx=ivec.size()-1; //last possible element
		}
		const size_t elements_right = MAX(end_time_idx - index, end_elements);
		end = index + elements_right;

		//now, check (and modify) if the window could not be centered
		if(elements_left==elements_right) return true;
		if(!is_soft) return false;
		if(elements_left<elements_right) { //we hit the left border
			//get again the end of window
212
213
214
215
216
			size_t end_elems = (min_data_points>(elements_left+1))?min_data_points - (elements_left + 1):0;
			const Date end_dt = ivec[start].date+min_time_span;
			size_t end_tm_idx = (end_dt>date)?IOUtils::seek(end_dt, ivec, false):index; //end time criteria
			if(end_tm_idx==IOUtils::npos) {
				end_tm_idx=ivec.size()-1; //last possible element
217
			}
218
219
			const size_t elems_right = MAX(end_tm_idx - index, end_elems);
			end = index + elems_right;
220
221
		} else { //we hit the right border
			//get again the start of window
222
223
224
225
226
			size_t start_elems = (min_data_points>(elements_right+1))?min_data_points - (elements_right + 1):0;
			const Date start_dt = ivec[end].date-min_time_span;
			size_t start_tm_idx = (start_dt<date)?IOUtils::seek(start_dt, ivec, false):index; //start time criteria
			if(start_tm_idx!=IOUtils::npos && start_tm_idx!=index) start_tm_idx = (start_tm_idx>0)?start_tm_idx-1:IOUtils::npos;
			if(start_tm_idx==IOUtils::npos) {
227
228
				end_time_idx=0; //first possible element
			}
229
230
			const size_t elems_left = MAX(index - start_tm_idx, start_elems);
			start = index - elems_left;
231
232
		}
	}
233
234

	return true;
235
236
}

237
} //namespace