Coverage for wsimod/core/core.py: 18%
147 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-01-11 16:39 +0000
1# -*- coding: utf-8 -*-
2"""Created on Wed Apr 7 15:44:48 2021.
4@author: Barney
6Converted to totals on Thur Apr 21 2022
7"""
8from math import log10
10from wsimod.core import constants
13class WSIObj:
14 """"""
16 def __init__(self):
17 """WSIObj is the base object of everything in WSIMOD. It is used to perform VQIP
18 operations and mass balance checking behaviour.
20 RSE has suggested that it would make more sense to leave VQIP operations as
21 regular functions in a module or associated them with a VQIP class.
23 Predefining empty_vqip_predefined in a class object is sensible though because
24 it is the foundation of many operations, and copying a dict is many times
25 quicker than copying a class.
27 For now I will leave WSIObj as the base object, but this may change.
28 """
29 # Predefine empty concentrations because copying is quicker than defining
30 self.empty_vqip_predefined = dict.fromkeys(constants.POLLUTANTS + ["volume"], 0)
32 def empty_vqip(self):
33 """Return a copy of the predefined empty vqip. All pollutants and volume
34 initialised in a dict and set to 0.
36 Returns:
37 empty_vqip_predefined (dict): Copy of empty_vqip_predefined
39 Examples:
40 >>> obj = WSIObj()
41 >>> obj.empty_vqip()
42 """
43 return self.empty_vqip_predefined.copy()
45 def copy_vqip(self, t):
46 """Wrapper to copy VQIP.
48 Args:
49 t (dict): A VQIP
51 Returns:
52 (dict): A copy of t
53 """
54 return t.copy()
56 def blend_vqip(self, c1, c2):
57 """Blends together two VQIPs that are assumed to have pollutant entries set as
58 pollution concentrations, blending occurs with proportionate mixing.
60 NOTE: VQIPs in WSIMOD in general store pollution as a total rather than a
61 concentration. So you should only blend if you are doing it intentionally and
62 know what you're doing. This won't do anything on VQIPs with 0 volume.
64 Args:
65 c1 (dict): A VQIP where pollutant entries are concentrations c2 (dict): A
66 VQIP where pollutant entries are concentrations
68 Returns:
69 c (dict): A new VQIP where c1 and c2 have been proportionately blended
70 """
71 # Blend two vqips given as concentrations
72 c = self.empty_vqip()
74 c["volume"] = c1["volume"] + c2["volume"]
75 if c["volume"] > 0:
76 for pollutant in constants.POLLUTANTS:
77 c[pollutant] = (
78 c1[pollutant] * c1["volume"] + c2[pollutant] * c2["volume"]
79 ) / c["volume"]
81 return c
83 def sum_vqip(self, t1, t2):
84 """Combines two VQIPs where pollutant entries are assumed to be given as mass.
85 Volume and additive pollutants are summed while non additive pollutants are
86 proportionately blended.
88 Args:
89 t1 (dict): A VQIP where pollutant entries are mass totals t2 (dict): A VQIP
90 where pollutant entries are mass totals
92 Returns:
93 t (dict): A VQIP that is the sum of t1 and t2 (except for non-additive
94 pollutants)
96 Examples:
97 >>> t1 = {'phosphate' : 0.25, 'volume' : 100, 'temperature' : 10}
98 >>> t2 = {'phosphate' : 0.25, 'volume' : 10, 'temperature' : 15}
99 >>> t = sum_vqip(t1, t2)
100 >>> print(t)
101 {'phosphate' : 0.5, 'volume' : 110, 'temperature' : 10.45}
102 """
103 # Sum two vqips given as totals
104 t = self.copy_vqip(t1)
105 t["volume"] += t2["volume"]
106 for pollutant in constants.ADDITIVE_POLLUTANTS:
107 t[pollutant] += t2[pollutant]
109 if t["volume"] > 0:
110 # Assume proportional blending of non additive pollutants
111 for pollutant in constants.NON_ADDITIVE_POLLUTANTS:
112 t[pollutant] = (
113 t2[pollutant] * t2["volume"] + t1[pollutant] * t1["volume"]
114 ) / t["volume"]
116 return t
118 def concentration_to_total(self, c):
119 """Convert a VQIP that has pollutant entries as concentrations into mass totals.
121 Args:
122 c (dict): A VQIP where pollutant entries are concentrations
124 Returns:
125 c (dict): A VQIP where pollutant entries are mass totals
126 """
127 c = self.copy_vqip(c)
128 for pollutant in constants.ADDITIVE_POLLUTANTS:
129 # Multiply concentration by volume to get mass for additive pollutants
130 c[pollutant] *= c["volume"]
131 return c
133 def total_to_concentration(self, t):
134 """Converts a VQIP that has pollutant entries as mass totals into
135 concentrations. Note, that this won't work for VQIPs with 0 volume.
137 Args:
138 t (dict): A VQIP where pollutant entries are mass totals
140 Returns:
141 c (dict): A VQIP where pollutant entries are concentrations
142 """
143 c = self.copy_vqip(t)
144 for pollutant in constants.ADDITIVE_POLLUTANTS:
145 # Divide concentration by volume to get concentration for additive
146 # pollutants
147 c[pollutant] /= c["volume"]
148 return c
150 def extract_vqip(self, t1, t2):
151 """Extract one VQIP from another where both VQIPs have pollutants as mass
152 totals. Each volume and additive pollutant is directly subtracted.
154 Args:
155 t1 (dict): A VQIP where pollutant entries are mass totals to subtract
156 from
157 t2 (dict): A VQIP where pollutant entries are mass totals to subtract
159 Returns:
160 t (dict): A copy of t1 where each additive pollutant and volume has had
161 t2 subtracted from it
163 Examples:
164 >>> t1 = {'phosphate' : 0.25, 'volume' : 100, 'temperature' : 10}
165 >>> t2 = {'phosphate' : 0.25, 'volume' : 10, 'temperature' : 15}
167 >>> t = extract_vqip(t1, t2)
168 >>> print(t)
169 {'phosphate' : 0, 'volume' : 90, 'temperature' : 10}
170 """
171 # TODO should probably be called 'subtract_vqip' TODO need to analyse uses of
172 # this to see if it is sensible to do something for non additive
173 t = self.copy_vqip(t1)
174 # Directly subtract t2 from t1 for vol and additive pollutants
175 for pol in constants.ADDITIVE_POLLUTANTS + ["volume"]:
176 t[pol] -= t2[pol]
178 return t
180 def extract_vqip_c(self, c1, c2):
181 """Extract one VQIP from another where both VQIPs have pollutants as
182 concentrations. Operation performed for volume and additive pollutants.
184 NOTE: VQIPs in WSIMOD in general store pollution as a total rather than a
185 concentration. So you should only work with concentrations if you are doing it
186 intentionally and know what you're doing.
188 Args:
189 c1 (dict): A VQIP where pollutant entries are concentrations to subtract
190 from
191 c2 (dict): A VQIP where pollutant entries are concentrations to subtract
193 Returns:
194 c (dict): A copy of c1 where each additive pollutant and volume has had
195 c2 proportionately extracted from it
196 """
197 c = self.copy_vqip(c1)
199 c1 = self.concentration_to_total(c1)
200 c2 = self.concentration_to_total(c2)
201 c["volume"] = c1["volume"] - c2["volume"]
202 if c["volume"] > 0:
203 for pollutant in constants.ADDITIVE_POLLUTANTS:
204 # Subtract c2 from c1 for vol and additive pollutants
205 c[pollutant] = (c1[pollutant] - c2[pollutant]) / c["volume"]
207 return c
209 def v_distill_vqip(self, t, v):
210 """Directly remove a volume from a VQIP.
212 Args:
213 t (dict): A VQIP where pollutant entries are mass totals to remove
214 volume from
215 v (float): Volume to remove
217 Returns:
218 t (dict): Updated VQIP
219 """
220 # Distill v from t
221 t = self.copy_vqip(t)
222 t["volume"] -= v
223 return t
225 def v_distill_vqip_c(self, c, v):
226 """Directly remove a volume from a VQIP, where pollutant entries are
227 concentrations.
229 NOTE: VQIPs in WSIMOD in general store pollution as a total rather than a
230 concentration. So you should only work with concentrations if you are doing it
231 intentionally and know what you're doing.
233 Args:
234 c (dict): A VQIP where pollutant entries are concentrations to remove
235 volume from
236 v (float): Volume to remove
238 Returns:
239 c (dict): Updated VQIP
240 """
241 # Distill v from c
242 c = self.copy_vqip(c)
243 d = self.empty_vqip()
244 d["volume"] = -v
245 c_ = self.blend_vqip(c, d)
246 # Directly copy non additive pollutants
247 for pollutant in constants.NON_ADDITIVE_POLLUTANTS:
248 c_[pollutant] = c[pollutant]
249 return c_
251 def v_change_vqip(self, t, v):
252 """Change the volume of a VQIP, where pollutants are mass totals, and update
253 pollutant values in proportion to the change in volume.
255 Args:
256 t (dict): A VQIP where pollutant entries are mass totals to get
257 pollutant concentrations from
258 v (float): Volume from t to get proportionate pollutant values in
260 Returns:
261 (dict): A VQIP with v volume and pollutions in proportion to t
263 Examples:
264 You want to extract 10m3 from 100m3 of water (store), to do this you need to
265 understand how much phosphate to extract in addition to volume.
267 >>> store = {'volume' : 100, 'phosphate' : 0.25}
268 >>> to_extract = v_change_vqip(store, 10)
270 >>> print(to_extract)
271 {'volume': 10, 'phosphate': 0.025}
272 """
273 t = self.copy_vqip(t)
274 if t["volume"] > 0:
275 # change all values of t by volume v in proportion to volume of t
276 ratio = v / t["volume"]
277 t["volume"] *= ratio
278 for pol in constants.ADDITIVE_POLLUTANTS:
279 t[pol] *= ratio
281 else:
282 # Assign volume directly
283 t["volume"] = v
284 return t
286 def v_change_vqip_c(self, c, v):
287 """Change the volume of a VQIP, where pollutants are concentrations.
289 NOTE: VQIPs in WSIMOD in general store pollution as a total rather than a
290 concentration. So you should only work with concentrations if you are doing it
291 intentionally and know what you're doing.
293 Args:
294 c (dict): A VQIP where pollutant entries are concentrations v (float):
295 Volume to change c's volume to
297 Returns:
298 c (dict): A new VQIP with volume udpated
299 """
300 # Change volume of vqip
301 c = self.copy_vqip(c)
302 c["volume"] = v
303 return c
305 def ds_vqip(self, t, t_):
306 """Get difference between each additive pollutant and volume for VQIPs where
307 pollutants are given as mass totals.
309 Args:
310 t (dict): A VQIP where pollutant entries are mass totals to subtract
311 values from
312 t_ (_type_): A VQIP where pollutant entries are mass totals to subtract
314 Returns:
315 ds (dict): Difference between t and t_ in mass totals
317 Examples:
318 >>> t1 = {'phosphate' : 0.25, 'volume' : 100, 'temperature' : 10}
319 >>> t2 = {'phosphate' : 0.2, 'volume' : 90, 'temperature' : 9}
321 >>> t = ds_vqip(t1, t2)
322 >>> print(t)
323 {'phosphate' : 0.05, 'volume' : 10, 'temperature' : 0}
324 """
325 ds = self.empty_vqip()
326 for pol in constants.ADDITIVE_POLLUTANTS + ["volume"]:
327 ds[pol] = t[pol] - t_[pol]
328 return ds
330 def ds_vqip_c(self, c, c_):
331 """Get difference between each additive pollutant and volume for VQIPs where
332 pollutants are given as concentrations but difference is given as mass totals.
334 NOTE: VQIPs in WSIMOD in general store pollution as a total rather than a
335 concentration. So you should only work with concentrations if you are doing it
336 intentionally and know what you're doing.
338 Args:
339 c (dict): A VQIP where pollutant entries are concentrations to subtract
340 values from
341 c_ (_type_): A VQIP where pollutant entries are concentrations to
342 subtract
344 Returns:
345 ds (dict): Difference between c and c_ in mass totals
346 """
347 ds = self.empty_vqip()
348 ds["volume"] = c["volume"] - c_["volume"]
349 for pol in constants.ADDITIVE_POLLUTANTS:
350 ds[pol] = c["volume"] * c[pol] - c_["volume"] * c_[pol]
351 # TODO what about non-additive ...
352 return ds
354 def compare_vqip(self, t1, t2):
355 """Compare two VQIPs and check if the difference between each key is less ' than
356 constants.FLOAT_ACCURACY.
358 Args:
359 t1 (dict): A VQIP t2 (dict): A VQIP
361 Returns:
362 bool: True if the difference is less for each key, False otherwise
363 """
364 reply = True
365 for v in t1.keys():
366 if abs(t1[v] - t2[v]) > constants.FLOAT_ACCURACY:
367 reply = False
368 return reply
370 def mass_balance(self):
371 """Call all mass balance functions and compare to see if discrepancy (i.e., if
372 in_ != (out_ + ds_) for volume or for any additive pollutant).
374 Comparison is performed in the magnitude of the largest value of in_, ds_ or
375 out_. And so judgement should be exercised as to whether a mass balance has
376 actually occurred
378 Returns:
379 in_ (dict): A VQIP of the total from mass_balance_in functions ds_ (dict): A
380 VQIP of the total from mass_balance_ds functions out_ (dict): A VQIP of the
381 total from mass_balance_out functions
383 Raises:
384 Message if mass balance does not close to constants.FLOAT_ACCURACY
385 """
386 # Iterate over mass_balance_in functions, summing values in in_
387 in_ = self.empty_vqip()
388 for f in self.mass_balance_in:
389 in_ = self.sum_vqip(in_, f())
391 # Iterate over mass_balance_out functions, summing values in out_
392 out_ = self.empty_vqip()
393 for f in self.mass_balance_out:
394 out_ = self.sum_vqip(out_, f())
396 # Iterate over mass_balance_ds functions, summing values in ds_
397 ds_ = self.empty_vqip()
398 for f in self.mass_balance_ds:
399 ds_f = f()
400 for v in constants.ADDITIVE_POLLUTANTS + ["volume"]:
401 ds_[v] += ds_f[v]
403 # Iterate over volume and additive pollutants to perform comparison
404 for v in ["volume"] + constants.ADDITIVE_POLLUTANTS:
405 # Find the largest value of in_, out_, ds_
406 largest = max(in_[v], out_[v], ds_[v])
408 if largest > constants.FLOAT_ACCURACY:
409 # Convert perform comparison in a magnitude to match the largest value
410 magnitude = 10 ** int(log10(largest))
411 in_10 = in_[v] / magnitude
412 out_10 = out_[v] / magnitude
413 ds_10 = ds_[v] / magnitude
414 else:
415 in_10 = in_[v]
416 ds_10 = ds_[v]
417 out_10 = out_[v]
419 if abs(in_10 - ds_10 - out_10) > constants.FLOAT_ACCURACY:
420 # Print mass balance error Print actual difference rather than magnitude
421 # comparison to enable user judgement
423 print(
424 "mass balance error for {0} of {1} in {2}".format(
425 v, in_[v] - ds_[v] - out_[v], self.name
426 )
427 )
429 return in_, ds_, out_
432class DecayObj(WSIObj):
433 """"""
435 # TODO - internet says this is a bad idea (diamond will occur when a Node - a type
436 # of WSIObj inherits a DecayObj - also a type of WSIObj). The reason diamonds are
437 # problems is because there can be conflicts in functions. But I don't want anyone
438 # to overwrite WSIObj functions so I don't see an issue?
439 def __init__(self, decays):
440 """A WSIObj that has decay functions built in.
442 Args:
443 decays (dict): A dict of dicts containing a key for each pollutant that
444 decays
445 and, within that, a key for each parameter (a constant and exponent)
447 Examples:
448 The 'constant' parameter represents what proportion of an amount will
449 decrease each time make_decay is called. Lower value will reduce decay.
450 Bounded between 0 and 1. The 'exponent' parameter represents how temperature
451 sensitive the decay is. The higher the value, the more pollution occurs at
452 higher values. Values expected to vary between 1 (no temperature
453 sensitivity) and 1.1 (high temperature sensitivity).
455 >>> decays = {'phosphate' : {'constant' : 0.001, 'exponent' : 1.005}}
457 Raises:
458 Message if no access to temperature data
459 """
460 # Store decays
461 self.decays = decays
462 super().__init__()
464 # Identify parent object to read temperature data
465 if "parent" in dir(self):
466 self.data_input_object = self.parent
467 elif "in_port" in dir(self):
468 self.data_input_object = self.in_port
469 else:
470 print("warning: decay object cannot access temperature data")
472 self.total_decayed = self.empty_vqip()
474 def make_decay(self, vqip):
475 """Make decay, reading tempature and updating pollutant amounts. A wrapper for
476 generic_temperature_decay.
478 Args:
479 vqip (dict): A VQIP to decay where pollutants are given as mass totals
481 Returns:
482 vqip_ (dict): A VQIP with pollutant amounts updated
483 """
484 # Read temperature data
485 temperature = self.data_input_object.data_input_dict[
486 ("temperature", self.data_input_object.t)
487 ]
488 # Make decay
489 vqip_, diff = self.generic_temperature_decay(vqip, self.decays, temperature)
490 # Update total_decayed for mass balance checking
491 self.total_decayed = self.sum_vqip(self.total_decayed, diff)
492 return vqip_
494 def generic_temperature_decay(self, t, d, temperature):
495 """Performs temperature sensitive pollutant decay calculations for a VQIP where
496 pollutants are given as mass totals.
498 Args:
499 t (dict): A VQIP to decay where pollutants are given as mass totals d
500 (dict): decays in a DecayObj temperature (float): temperature
502 Returns:
503 t (dict): A VQIP with updated pollutant values diff (dict): A VQIP storing
504 the change in pollutant values (decreases
505 stored as positive numbers)
506 """
507 t = self.copy_vqip(t)
508 diff = self.empty_vqip()
509 # Iterate over pollutants in d (keys)
510 for pol, pars in d.items():
511 # Perform calculation
512 diff[pol] = t[pol] * min(
513 pars["constant"]
514 * pars["exponent"]
515 ** (temperature - constants.DECAY_REFERENCE_TEMPERATURE),
516 1,
517 )
518 # Update VQIP
519 t[pol] -= diff[pol]
521 return t, diff
523 def generic_temperature_decay_c(self, c, d, temperature):
524 """Performs temperature sensitive pollutant decay calculations for a VQIP where
525 pollutants are given as concentrations.
527 Args:
528 c (dict): A VQIP to decay where pollutants are given as concentrations. d
529 (dict): decays in a DecayObj temperature (float): temperature
531 Returns:
532 t (dict): A VQIP with updated pollutant values (pollutants as
533 concentrations)
534 diff (dict): A VQIP storing the change in pollutant values (decreases
535 stored as positive numbers). Pollutants as mass totals.
536 """
537 c = self.copy_vqip(c)
538 diff = self.empty_vqip()
539 for pol, pars in d.items():
540 diff[pol] = c[pol] * min(
541 pars["constant"]
542 * pars["exponent"]
543 ** (temperature - constants.DECAY_REFERENCE_TEMPERATURE),
544 1,
545 )
546 c[pol] -= diff[pol]
548 diff[pol] *= c["volume"]
549 return c, diff