diff --git a/covid-sim.vcxproj b/covid-sim.vcxproj
index b24b4149f..b107ec981 100644
--- a/covid-sim.vcxproj
+++ b/covid-sim.vcxproj
@@ -107,6 +107,7 @@
+
@@ -134,6 +135,7 @@
+
diff --git a/docs/Doxyfile.in b/docs/Doxyfile.in
index 35e33c9d6..693e7e92b 100644
--- a/docs/Doxyfile.in
+++ b/docs/Doxyfile.in
@@ -791,6 +791,7 @@ WARN_LOGFILE = warnings.txt
# Note: If this tag is empty the current directory is searched.
INPUT = @CMAKE_CURRENT_SOURCE_DIR@/src \
+ @CMAKE_CURRENT_SOURCE_DIR@/include \
@CMAKE_CURRENT_SOURCE_DIR@/data \
@CMAKE_CURRENT_SOURCE_DIR@/docs \
@CMAKE_CURRENT_SOURCE_DIR@/Rscripts \
diff --git a/include/geometry/CMakeLists.txt b/include/geometry/CMakeLists.txt
index 266b10e37..1f3b8086c 100644
--- a/include/geometry/CMakeLists.txt
+++ b/include/geometry/CMakeLists.txt
@@ -2,4 +2,5 @@ set(GEOMETRY_HDR_FILES
${CMAKE_CURRENT_SOURCE_DIR}/Size.h
${CMAKE_CURRENT_SOURCE_DIR}/Vector2.h
${CMAKE_CURRENT_SOURCE_DIR}/BoundingBox.h
+ ${CMAKE_CURRENT_SOURCE_DIR}/Distance.h
CACHE INTERNAL "Headers" FORCE)
diff --git a/include/geometry/Distance.h b/include/geometry/Distance.h
new file mode 100644
index 000000000..e3afce9a2
--- /dev/null
+++ b/include/geometry/Distance.h
@@ -0,0 +1,58 @@
+#ifndef COVIDSIM_GEOMETRY_DISTANCE_H_INCLUDED_
+#define COVIDSIM_GEOMETRY_DISTANCE_H_INCLUDED_
+
+#include "BoundingBox.h"
+#include "Vector2.h"
+#include "Size.h"
+
+#include
+
+namespace CovidSim
+{
+ namespace Geometry
+ {
+ /// ABC for distance calculations
+ class Distance
+ {
+ public:
+ /// Calculate the distance squared between two rectangular grid indices
+ /// \param l index
+ /// \param m index
+ /// \param stride number of indices in a row
+ /// \param scales conversion to global co-ordinates
+ /// \param minimum whether to calculate the minimum between cells
+ /// \return distance squared
+ virtual double distance_squared(int l, int m, int stride, const Size& scales, bool minimum = false) = 0;
+
+ /// Calculate the distance squared between two points
+ /// \param a point in global co-ordinates
+ /// \param b point in global co-ordinates
+ /// \return distance squared
+ virtual double distance_squared(const Vector2d& a, const Vector2d& b) = 0;
+
+ /// Calculate the distance squared between two points
+ /// \param a point in global co-ordinates
+ /// \param b point in global co-ordinates
+ /// \return distance squared
+ double distance_squared(const Vector2f& a, const Vector2f& b)
+ {
+ return distance_squared(Vector2d(a.x, a.y), Vector2d(b.x, b.y));
+ }
+ };
+
+ /// Provide a factory method for creating a concrete Distance class
+ class DistanceFactory
+ {
+ public:
+ /// Create a concrete distance
+ /// \param utm whether to perform UTM calculations
+ /// \param periodic_boundaries if not UTM whether the geometry is periodic or Eulerian
+ /// \param spatial_bounding_box for UTM the bounds to which any lat long co-ordinates are relative
+ /// \param global for periodic the horizontal and vertical periods in global units
+ /// \return instance of a concrete distance
+ static std::shared_ptr create(bool utm, bool periodic_boundaries, const BoundingBox2d& spatial_bounding_box, const Size& global);
+ };
+ }
+}
+
+#endif
diff --git a/src/Constants.h b/src/Constants.h
index 17a40fbb8..a07a49027 100644
--- a/src/Constants.h
+++ b/src/Constants.h
@@ -13,56 +13,6 @@
*/
constexpr double PI = 3.1415926535; // full double precision: 3.14159265358979323846
-/**
- * An arc minute of latitude along any line of longitude in meters.
- *
- * Also known as the International Nautical Mile.
- *
- * @see https://en.wikipedia.org/wiki/Nautical_mile
- */
-constexpr int NMI = 1852;
-
-/**
- * The number of arc minutes in one degree.
- *
- * @see https://en.wikipedia.org/wiki/Minute_and_second_of_arc
- */
-constexpr int ARCMINUTES_PER_DEGREE = 60;
-
-/**
- * The number of degrees in a complete rotation.
- *
- * @see https://en.wikipedia.org/wiki/Turn_(angle)
- */
-constexpr int DEGREES_PER_TURN = 360;
-
-/**
- * The earth's circumference in meters.
- *
- * The units of cancellation:
- * meters/minute * minutes/degree * degrees = meters
- */
-constexpr int EARTH_CIRCUMFERENCE = NMI * ARCMINUTES_PER_DEGREE * DEGREES_PER_TURN;
-
-/**
- * The earth's diameter in meters.
- */
-constexpr double EARTH_DIAMETER = EARTH_CIRCUMFERENCE / PI;
-
-/**
- * The Earth's radius in meters.
- *
- * The previous hardcoded value used 6366707 which was derived from the
- * following formula:
- *
- * Earth's radius (m) = Earth's circumference / 2 * Pi
- *
- * where Earth's circumference can be derived with the following formula:
- *
- * Earth's circumference (m) = NMI * ARCMINUTES_PER_DEGREE * DEGREES_PER_TURN
- */
-constexpr double EARTHRADIUS = EARTH_DIAMETER / 2;
-
const int OUTPUT_DIST_SCALE = 1000;
const int MAX_PLACE_SIZE = 20000;
const int MAX_NUM_SEED_LOCATIONS = 10000;
diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp
index 34755eb1c..ed8519539 100644
--- a/src/CovidSim.cpp
+++ b/src/CovidSim.cpp
@@ -2074,20 +2074,9 @@ void ReadParams(std::string const& ParamFile, std::string const& PreParamFile, s
P.usCaseIsolationDuration = ((unsigned short int) (P.CaseIsolationDuration * P.TimeStepsPerDay));
P.usCaseAbsenteeismDuration = ((unsigned short int) (P.CaseAbsenteeismDuration * P.TimeStepsPerDay));
P.usCaseAbsenteeismDelay = ((unsigned short int) (P.CaseAbsenteeismDelay * P.TimeStepsPerDay));
- if (P.DoUTM_coords)
- {
- for (i = 0; i <= 1000; i++)
- {
- asin2sqx[i] = asin(sqrt(i / 1000.0));
- asin2sqx[i] = asin2sqx[i] * asin2sqx[i];
- }
- for (i = 0; i <= DEGREES_PER_TURN; i++)
- {
- t = PI * i / 180;
- sinx[i] = sin(t);
- cosx[i] = cos(t);
- }
- }
+
+ P.distance_ = CovidSim::Geometry::DistanceFactory::create(P.DoUTM_coords, P.DoPeriodicBoundaries, P.SpatialBoundingBox, P.in_degrees_);
+
if (P.R0scale != 1.0)
{
P.HouseholdTrans *= P.R0scale;
@@ -2475,7 +2464,7 @@ void ReadAirTravel(std::string const& air_travel_file, std::string const& output
for (j = 0; j < Airports[i].num_connected; j++)
{
k = (int)Airports[i].conn_airports[j];
- traf = floor(sqrt(dist2_raw(Airports[i].loc.x, Airports[i].loc.y, Airports[k].loc.x, Airports[k].loc.y)) / OUTPUT_DIST_SCALE);
+ traf = floor(sqrt(P.distance_->distance_squared(Airports[i].loc, Airports[k].loc)) / OUTPUT_DIST_SCALE);
l = (int)traf;
//fprintf(stderr,"%(%i) ",l);
if (l < MAX_DIST)
@@ -3210,7 +3199,7 @@ void SaveDistribs(std::string const& output_file_base)
((AdUnits[Mcells[Hosts[i].mcell].adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor == (P.OutputPlaceDistAdunit % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor))
{
k = Hosts[i].PlaceLinks[j];
- s = sqrt(dist2_raw(Households[Hosts[i].hh].loc.x, Households[Hosts[i].hh].loc.y, Places[j][k].loc.x, Places[j][k].loc.y)) / OUTPUT_DIST_SCALE;
+ s = sqrt(P.distance_->distance_squared(Households[Hosts[i].hh].loc, Places[j][k].loc)) / OUTPUT_DIST_SCALE;
k = (int)s;
if (k < MAX_DIST) PlaceDistDistrib[j][k]++;
}
diff --git a/src/Dist.cpp b/src/Dist.cpp
index 781c3ae26..d1db7a559 100644
--- a/src/Dist.cpp
+++ b/src/Dist.cpp
@@ -1,177 +1,38 @@
-#include
-#include
-
-#include "Constants.h"
#include "Dist.h"
#include "Param.h"
#include "Model.h"
-double sinx[DEGREES_PER_TURN + 1], cosx[DEGREES_PER_TURN + 1], asin2sqx[1001];
-
-//// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// ****
-//// **** DISTANCE FUNCTIONS (return distance-squared, which is input for every Kernel function)
-//// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// ****
-
-double periodic_xy(double x, double y) {
- if (P.DoPeriodicBoundaries)
- {
- if (x > P.in_degrees_.width * 0.5) x = P.in_degrees_.width - x;
- if (y > P.in_degrees_.height * 0.5) y = P.in_degrees_.height - y;
- }
- return x * x + y * y;
-}
-
-double dist2UTM(double x1, double y1, double x2, double y2)
-{
- double x, y, cy1, cy2, yt, xi, yi;
-
- x = fabs(x1 - x2) / 2;
- y = fabs(y1 - y2) / 2;
- xi = floor(x);
- yi = floor(y);
- x -= xi;
- y -= yi;
- x = (1 - x) * sinx[(int)xi] + x * sinx[((int)xi) + 1];
- y = (1 - y) * sinx[(int)yi] + y * sinx[((int)yi) + 1];
- yt = fabs(y1 + P.SpatialBoundingBox.bottom_left().y);
- yi = floor(yt);
- cy1 = yt - yi;
- cy1 = (1 - cy1) * cosx[((int)yi)] + cy1 * cosx[((int)yi) + 1];
- yt = fabs(y2 + P.SpatialBoundingBox.bottom_left().y);
- yi = floor(yt);
- cy2 = yt - yi;
- cy2 = (1 - cy2) * cosx[((int)yi)] + cy2 * cosx[((int)yi) + 1];
- x = fabs(1000 * (y * y + x * x * cy1 * cy2));
- xi = floor(x);
- x -= xi;
- y = (1 - x) * asin2sqx[((int)xi)] + x * asin2sqx[((int)xi) + 1];
- return 4 * EARTHRADIUS * EARTHRADIUS * y;
-}
double dist2(Person* a, Person* b)
{
- double x, y;
-
- if (P.DoUTM_coords)
- return dist2UTM(Households[a->hh].loc.x, Households[a->hh].loc.y, Households[b->hh].loc.x, Households[b->hh].loc.y);
- else
- {
- x = fabs(Households[a->hh].loc.x - Households[b->hh].loc.x);
- y = fabs(Households[a->hh].loc.y - Households[b->hh].loc.y);
- return periodic_xy(x, y);
- }
+ return P.distance_->distance_squared(Households[a->hh].loc, Households[b->hh].loc);
}
+
double dist2_cc(Cell* a, Cell* b)
{
- double x, y;
- int l, m;
+ int l = (int)(a - Cells);
+ int m = (int)(b - Cells);
- l = (int)(a - Cells);
- m = (int)(b - Cells);
- if (P.DoUTM_coords)
- return dist2UTM(P.in_cells_.width * fabs((double)(l / P.nch)), P.in_cells_.height * fabs((double)(l % P.nch)),
- P.in_cells_.width * fabs((double)(m / P.nch)), P.in_cells_.height * fabs((double)(m % P.nch)));
- else
- {
- x = P.in_cells_.width * fabs((double)(l / P.nch - m / P.nch));
- y = P.in_cells_.height * fabs((double)(l % P.nch - m % P.nch));
- return periodic_xy(x, y);
- }
+ return P.distance_->distance_squared(l, m, P.nch, P.in_cells_);
}
+
double dist2_cc_min(Cell* a, Cell* b)
{
- double x, y;
- int l, m, i, j;
+ int l = (int)(a - Cells);
+ int m = (int)(b - Cells);
- l = (int)(a - Cells);
- m = (int)(b - Cells);
- i = l; j = m;
- if (P.DoUTM_coords)
- {
- if (P.in_cells_.width * ((double)abs(m / P.nch - l / P.nch)) > PI)
- {
- if (m / P.nch > l / P.nch)
- j += P.nch;
- else if (m / P.nch < l / P.nch)
- i += P.nch;
- }
- else
- {
- if (m / P.nch > l / P.nch)
- i += P.nch;
- else if (m / P.nch < l / P.nch)
- j += P.nch;
- }
- if (m % P.nch > l % P.nch)
- i++;
- else if (m % P.nch < l % P.nch)
- j++;
- return dist2UTM(P.in_cells_.width * fabs((double)(i / P.nch)), P.in_cells_.height * fabs((double)(i % P.nch)),
- P.in_cells_.width * fabs((double)(j / P.nch)), P.in_cells_.height * fabs((double)(j % P.nch)));
- }
- else
- {
- if ((P.DoPeriodicBoundaries) && (P.in_cells_.width * ((double)abs(m / P.nch - l / P.nch)) > P.in_degrees_.width * 0.5))
- {
- if (m / P.nch > l / P.nch)
- j += P.nch;
- else if (m / P.nch < l / P.nch)
- i += P.nch;
- }
- else
- {
- if (m / P.nch > l / P.nch)
- i += P.nch;
- else if (m / P.nch < l / P.nch)
- j += P.nch;
- }
- if ((P.DoPeriodicBoundaries) && (P.in_degrees_.height * ((double)abs(m % P.nch - l % P.nch)) > P.in_degrees_.height * 0.5))
- {
- if (m % P.nch > l % P.nch)
- j++;
- else if (m % P.nch < l % P.nch)
- i++;
- }
- else
- {
- if (m % P.nch > l % P.nch)
- i++;
- else if (m % P.nch < l % P.nch)
- j++;
- }
- x = P.in_cells_.width * fabs((double)(i / P.nch - j / P.nch));
- y = P.in_cells_.height * fabs((double)(i % P.nch - j % P.nch));
- return periodic_xy(x, y);
- }
+ return P.distance_->distance_squared(l, m, P.nch, P.in_cells_, true);
}
+
double dist2_mm(Microcell* a, Microcell* b)
{
- double x, y;
- int l, m;
+ int l = (int)(a - Mcells);
+ int m = (int)(b - Mcells);
- l = (int)(a - Mcells);
- m = (int)(b - Mcells);
- if (P.DoUTM_coords)
- return dist2UTM(P.in_microcells_.width * fabs((double)(l / P.total_microcells_high_)), P.in_microcells_.height * fabs((double)(l % P.total_microcells_high_)),
- P.in_microcells_.width * fabs((double)(m / P.total_microcells_high_)), P.in_microcells_.height * fabs((double)(m % P.total_microcells_high_)));
- else
- {
- x = P.in_microcells_.width * fabs((double)(l / P.total_microcells_high_ - m / P.total_microcells_high_));
- y = P.in_microcells_.height * fabs((double)(l % P.total_microcells_high_ - m % P.total_microcells_high_));
- return periodic_xy(x, y);
- }
+ return P.distance_->distance_squared(l, m, P.total_microcells_high_, P.in_microcells_);
}
double dist2_raw(double ax, double ay, double bx, double by)
{
- double x, y;
-
- if (P.DoUTM_coords)
- return dist2UTM(ax, ay, bx, by);
- else
- {
- x = fabs(ax - bx);
- y = fabs(ay - by);
- return periodic_xy(x, y);
- }
+ return P.distance_->distance_squared(CovidSim::Geometry::Vector2d(ax, ay), CovidSim::Geometry::Vector2d(bx, by));
}
diff --git a/src/Dist.h b/src/Dist.h
index 930cdc0d0..96dde86c2 100644
--- a/src/Dist.h
+++ b/src/Dist.h
@@ -6,7 +6,6 @@
#include "Models/Microcell.h"
#include "Constants.h"
-extern double sinx[DEGREES_PER_TURN + 1], cosx[DEGREES_PER_TURN + 1], asin2sqx[1001];
double dist2UTM(double, double, double, double);
double dist2(Person*, Person*);
double dist2_cc(Cell*, Cell*);
diff --git a/src/Geometry/CMakeLists.txt b/src/Geometry/CMakeLists.txt
index 428e5c0a6..8aaca8d39 100644
--- a/src/Geometry/CMakeLists.txt
+++ b/src/Geometry/CMakeLists.txt
@@ -1,4 +1,4 @@
-set(GEOMETRY_SRC_FILES Vector2.cpp)
+set(GEOMETRY_SRC_FILES Distance.cpp Vector2.cpp)
source_group(covidsim\\geometry FILES ${GEOMETRY_SRC_FILES} ${GEOMETRY_HDR_FILES})
add_library(geometrylib ${GEOMETRY_SRC_FILES} ${GEOMETRY_HDR_FILES})
diff --git a/src/Geometry/Distance.cpp b/src/Geometry/Distance.cpp
new file mode 100644
index 000000000..74e5fdcd5
--- /dev/null
+++ b/src/Geometry/Distance.cpp
@@ -0,0 +1,301 @@
+#include "geometry/Distance.h"
+
+namespace CovidSim
+{
+ namespace Geometry
+ {
+ /// Concrete class for performing UTM distance calculations
+ class UTMDistance : public Distance
+ {
+ /// Math constant defined as the ratio of a circle's circumference to its diameter.
+ ///
+ /// TODO: since all calculations using this constant are being automatically
+ /// type-casted to double, should the precision be extended for more accuracy in
+ /// the simulations?
+ ///
+ /// Eventually could be replaced with C++20's std::numbers::pi.
+ /// https://en.cppreference.com/w/cpp/header/numbers
+ static constexpr double pi_ = 3.1415926535; // full double precision: 3.14159265358979323846
+
+ /// An arc minute of latitude along any line of longitude in meters.
+ ///
+ /// Also known as the International Nautical Mile.
+ ///
+ /// @see https://en.wikipedia.org/wiki/Nautical_mile
+ static constexpr int nmi_ = 1852;
+
+ /// The number of arc minutes in one degree.
+ ///
+ /// @see https://en.wikipedia.org/wiki/Minute_and_second_of_arc
+ static constexpr int arcminutes_per_degree_ = 60;
+
+ /// The number of degrees in a complete rotation.
+ ///
+ /// @see https://en.wikipedia.org/wiki/Turn_(angle)
+ static constexpr int degrees_per_turn_ = 360;
+
+ /// The earth's circumference in meters.
+ ///
+ /// The units of cancellation:
+ /// meters/minute * minutes/degree * degrees = meters
+ static constexpr int earth_circumference_ = nmi_ * arcminutes_per_degree_ * degrees_per_turn_;
+
+ /// The earth's diameter in meters.
+ static constexpr double earth_diameter_ = earth_circumference_ / pi_;
+
+ /// The Earth's radius in meters.
+ static constexpr double earth_radius_ = earth_diameter_ / 2;
+
+ /// The bounds to which any lat long co-ordinates are relative
+ const BoundingBox2d& spatial_bounding_box_;
+
+ /// The lookup table for sin(x)
+ double sin_x_[degrees_per_turn_ + 1];
+
+ /// The lookup table for cos(x)
+ double cos_x_[degrees_per_turn_ + 1];
+
+ /// The lookup table for asin(sqrt(x)) * asin(sqrt(x))
+ double asin2_sq_x_[1001];
+
+ public:
+ /// Construct a class for UTM distance calculations
+ /// \param spatial_bounding_box the bounds to which any lat long co-ordinates are relative
+ UTMDistance(const BoundingBox2d& spatial_bounding_box)
+ : spatial_bounding_box_(spatial_bounding_box)
+ {
+ }
+
+ /// Create the lookup tables
+ void initialise()
+ {
+ for (int i = 0; i <= 1000; i++)
+ {
+ asin2_sq_x_[i] = asin(sqrt(i / 1000.0));
+ asin2_sq_x_[i] = asin2_sq_x_[i] * asin2_sq_x_[i];
+ }
+
+ for (int i = 0; i <= degrees_per_turn_; i++)
+ {
+ double t = pi_ * i / 180;
+ sin_x_[i] = sin(t);
+ cos_x_[i] = cos(t);
+ }
+ }
+
+ /// Calculate the distance squared between two rectangular grid indices
+ /// \param l index
+ /// \param m index
+ /// \param stride number of indices in a row
+ /// \param scales conversion to global co-ordinates
+ /// \param minimum whether to calculate the minimum between cells
+ /// \return distance squared
+ double distance_squared(int l, int m, int stride, const Size& scales, bool minimum) override
+ {
+ int i = l;
+ int j = m;
+ if (minimum)
+ {
+ if (scales.width * ((double)abs(m / stride - l / stride)) > pi_)
+ {
+ if (m / stride > l / stride)
+ j += stride;
+ else if (m / stride < l / stride)
+ i += stride;
+ }
+ else
+ {
+ if (m / stride > l / stride)
+ i += stride;
+ else if (m / stride < l / stride)
+ j += stride;
+ }
+
+ if (m % stride > l % stride)
+ i++;
+ else if (m % stride < l % stride)
+ j++;
+ }
+
+ CovidSim::Geometry::Vector2d u(scales.width * fabs((double)(i / stride)), scales.height * fabs((double)(i % stride)));
+ CovidSim::Geometry::Vector2d v(scales.width * fabs((double)(j / stride)), scales.height * fabs((double)(j % stride)));
+ return distance_squared(u, v);
+ }
+
+ /// Calculate the distance squared between two points
+ /// \param a lat long point in degrees relative to the spatial bounding box
+ /// \param b lat long point in degrees relative to the spatial bounding box
+ /// \return distance squared
+ double distance_squared(const Vector2d& a, const Vector2d& b) override
+ {
+ double x1 = a.x;
+ double y1 = a.y;
+ double x2 = b.x;
+ double y2 = b.y;
+
+ double x, y, cy1, cy2, yt, xi, yi;
+
+ x = fabs(x1 - x2) / 2;
+ y = fabs(y1 - y2) / 2;
+ xi = floor(x);
+ yi = floor(y);
+ x -= xi;
+ y -= yi;
+ x = (1 - x) * sin_x_[(int)xi] + x * sin_x_[((int)xi) + 1];
+ y = (1 - y) * sin_x_[(int)yi] + y * sin_x_[((int)yi) + 1];
+ yt = fabs(y1 + spatial_bounding_box_.bottom_left().y);
+ yi = floor(yt);
+ cy1 = yt - yi;
+ cy1 = (1 - cy1) * cos_x_[((int)yi)] + cy1 * cos_x_[((int)yi) + 1];
+ yt = fabs(y2 + spatial_bounding_box_.bottom_left().y);
+ yi = floor(yt);
+ cy2 = yt - yi;
+ cy2 = (1 - cy2) * cos_x_[((int)yi)] + cy2 * cos_x_[((int)yi) + 1];
+ x = fabs(1000 * (y * y + x * x * cy1 * cy2));
+ xi = floor(x);
+ x -= xi;
+ y = (1 - x) * asin2_sq_x_[((int)xi)] + x * asin2_sq_x_[((int)xi) + 1];
+ return 4 * earth_radius_ * earth_radius_ * y;
+ }
+ };
+
+ /// Concrete class for performing periodic distance calculations
+ class PeriodicDistance : public Distance
+ {
+ /// The horizontal and vertical periods in global units
+ const Size& global_;
+
+ public:
+ /// Construct a class for periodic distance calculations
+ /// \param global the horizontal and vertical periods in global units
+ PeriodicDistance(const Size& global)
+ : global_(global)
+ {
+ }
+
+ /// Calculate the distance squared between two rectangular grid indices
+ /// \param l index
+ /// \param m index
+ /// \param stride number of indices in a row
+ /// \param scales conversion to global co-ordinates
+ /// \param minimum whether to calculate the minimum between cells
+ /// \return distance squared
+ double distance_squared(int l, int m, int stride, const Size& scales, bool minimum) override
+ {
+ int i = l;
+ int j = m;
+ if (minimum)
+ {
+ if (scales.width * ((double)abs(m / stride - l / stride)) > global_.width * 0.5)
+ {
+ if (m / stride > l / stride)
+ j += stride;
+ else if (m / stride < l / stride)
+ i += stride;
+ }
+ else
+ {
+ if (m / stride > l / stride)
+ i += stride;
+ else if (m / stride < l / stride)
+ j += stride;
+ }
+
+ if (scales.height * ((double)abs(m % stride - l % stride)) > global_.height * 0.5)
+ {
+ if (m % stride > l % stride)
+ j++;
+ else if (m % stride < l % stride)
+ i++;
+ }
+ else
+ {
+ if (m % stride > l % stride)
+ i++;
+ else if (m % stride < l % stride)
+ j++;
+ }
+ }
+
+ CovidSim::Geometry::Vector2d u(scales.width * fabs((double)(i / stride)), scales.height * fabs((double)(i % stride)));
+ CovidSim::Geometry::Vector2d v(scales.width * fabs((double)(j / stride)), scales.height * fabs((double)(j % stride)));
+ return distance_squared(u, v);
+ }
+
+ /// Calculate the distance squared between two points
+ /// \param a point in global co-ordinates
+ /// \param b point in global co-ordinates
+ /// \return distance squared
+ double distance_squared(const Vector2d& a, const Vector2d& b) override
+ {
+ auto diff = a - b;
+ if (diff.x > global_.width * 0.5)
+ diff.x = global_.width - diff.x;
+
+ if (diff.y > global_.height * 0.5)
+ diff.y = global_.height - diff.y;
+
+ return diff.length_squared();
+ }
+ };
+
+ /// Concrete class for performing Eulerian distance calculations
+ class EulerDistance : public Distance
+ {
+ public:
+ /// Calculate the distance squared between two rectangular grid indices
+ /// \param l index
+ /// \param m index
+ /// \param stride number of indices in a row
+ /// \param scales conversion to global co-ordinates
+ /// \param minimum whether to calculate the minimum between cells
+ /// \return distance squared
+ double distance_squared(int l, int m, int stride, const Size& scales, bool minimum) override
+ {
+ int i = l;
+ int j = m;
+ if (minimum)
+ {
+ if (m / stride > l / stride)
+ i += stride;
+ else if (m / stride < l / stride)
+ j += stride;
+
+ if (m % stride > l % stride)
+ i++;
+ else if (m % stride < l % stride)
+ j++;
+ }
+
+ CovidSim::Geometry::Vector2d u(scales.width * fabs((double)(i / stride)), scales.height * fabs((double)(i % stride)));
+ CovidSim::Geometry::Vector2d v(scales.width * fabs((double)(j / stride)), scales.height * fabs((double)(j % stride)));
+ return distance_squared(u, v);
+ }
+
+ /// Calculate the distance squared between two points
+ /// \param a point in global co-ordinates
+ /// \param b point in global co-ordinates
+ /// \return distance squared
+ double distance_squared(const Vector2d& a, const Vector2d& b) override
+ {
+ auto diff = a - b;
+ return diff.length_squared();
+ }
+ };
+
+ std::shared_ptr DistanceFactory::create(bool utm, bool periodic_boundaries, const BoundingBox2d& spatial_bounding_box, const Size& global)
+ {
+ if (utm)
+ {
+ auto dist = std::make_shared(spatial_bounding_box);
+ dist->initialise();
+ return dist;
+ }
+
+ if (periodic_boundaries)
+ return std::make_shared(global);
+
+ return std::make_shared();
+ }
+ }
+}
diff --git a/src/Param.h b/src/Param.h
index d5e7499af..d51389579 100644
--- a/src/Param.h
+++ b/src/Param.h
@@ -10,6 +10,7 @@
#include "MicroCellPosition.hpp"
#include "geometry/BoundingBox.h"
+#include "geometry/Distance.h"
#include "geometry/Size.h"
/** @brief Enumeration of bitmap formats. */
@@ -50,6 +51,8 @@ struct Param
int NMCP, ncw, nch, DoUTM_coords, nsp, DoSeasonality, DoCorrectAgeDist, DoPartialImmunity;
int total_microcells_wide_, total_microcells_high_;
+ std::shared_ptr distance_;
+
MicroCellPosition get_micro_cell_position_from_cell_index(int cell_index) const;
int get_micro_cell_index_from_position(MicroCellPosition const& position) const;
bool is_in_bounds(MicroCellPosition const& position) const;
diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp
index d7f0d2ad1..3930c1ba6 100644
--- a/src/SetupModel.cpp
+++ b/src/SetupModel.cpp
@@ -1604,7 +1604,7 @@ void SetupAirports(void)
for (int j = 0; j < P.Nairports; j++)
if (Airports[j].total_traffic > 0)
{
- t = P.KernelLookup.num(dist2_raw(x, y, Airports[j].loc.x, Airports[j].loc.y)) * Airports[j].total_traffic;
+ t = P.KernelLookup.num(P.distance_->distance_squared(CovidSim::Geometry::Vector2d(x, y), CovidSim::Geometry::Vector2d(Airports[j].loc))) * Airports[j].total_traffic;
if (k < NNA)
{
Mcells[i].AirportList[k].id = j;
@@ -1703,16 +1703,14 @@ void SetupAirports(void)
tmin += 0.25 * MAX_DIST_AIRPORT_TO_HOTEL * MAX_DIST_AIRPORT_TO_HOTEL;
Airports[i].num_place = 0;
for (int j = 0; j < P.Nplace[P.HotelPlaceType]; j++)
- if (dist2_raw(Airports[i].loc.x, Airports[i].loc.y,
- Places[P.HotelPlaceType][j].loc.x, Places[P.HotelPlaceType][j].loc.y) < tmin)
+ if (P.distance_->distance_squared(Airports[i].loc, Places[P.HotelPlaceType][j].loc) < tmin)
Airports[i].num_place++;
} while (Airports[i].num_place < m);
if (tmin > MAX_DIST_AIRPORT_TO_HOTEL * MAX_DIST_AIRPORT_TO_HOTEL) fprintf(stderr_shared, "*** %i : %lg %i ***\n", i, sqrt(tmin), Airports[i].num_place);
Airports[i].DestPlaces = (IndexList*)Memory::xcalloc(Airports[i].num_place, sizeof(IndexList));
Airports[i].num_place = 0;
for (int j = 0; j < P.Nplace[P.HotelPlaceType]; j++)
- if ((t = dist2_raw(Airports[i].loc.x, Airports[i].loc.y,
- Places[P.HotelPlaceType][j].loc.x, Places[P.HotelPlaceType][j].loc.y)) < tmin)
+ if ((t = P.distance_->distance_squared(Airports[i].loc, Places[P.HotelPlaceType][j].loc)) < tmin)
{
Airports[i].DestPlaces[Airports[i].num_place].prob = (float)P.KernelLookup.num(t);
Airports[i].DestPlaces[Airports[i].num_place].id = j;
@@ -2203,8 +2201,7 @@ void AssignPeopleToPlaces()
auto const place_idx = cur_cell.places[tp][cnt];
if (place_idx >= P.Nplace[tp]) fprintf(stderr, "#%i %i %i ", tp, ic, cnt);
auto const& cur_place = Places[tp][place_idx];
- t = dist2_raw(Households[Hosts[i].hh].loc.x, Households[Hosts[i].hh].loc.y,
- cur_place.loc.x, cur_place.loc.y);
+ t = P.distance_->distance_squared(Households[Hosts[i].hh].loc, cur_place.loc);
s = P.KernelLookup.num(t);
if (tp < P.nsp)
{
@@ -2343,7 +2340,7 @@ void AssignPeopleToPlaces()
fprintf(stderr, "*%i %i: %i %i\n", k, tp, j, P.Nplace[tp]);
ERR_CRITICAL("Out of bounds place link\n");
}
- t = dist2_raw(Households[Hosts[k].hh].loc.x, Households[Hosts[k].hh].loc.y, Places[tp][j].loc.x, Places[tp][j].loc.y);
+ t = P.distance_->distance_squared(Households[Hosts[k].hh].loc, Places[tp][j].loc);
s = ((double)ct->S) / ((double)ct->S0) * P.KernelLookup.num(t) / Cells[i].max_trans[l];
if ((P.DoAdUnits) && (P.InhibitInterAdunitPlaceAssignment[tp] > 0))
{
diff --git a/src/Sweep.cpp b/src/Sweep.cpp
index d344a5d5b..efb2ec962 100644
--- a/src/Sweep.cpp
+++ b/src/Sweep.cpp
@@ -130,7 +130,7 @@ void TravelDepartSweep(double t)
k = Airports[i].conn_airports[l];
if (bm)
{
- if (dist2_raw(Airports[i].loc.x, Airports[i].loc.y, Airports[k].loc.x, Airports[k].loc.y) > P.MoveRestrRadius2)
+ if (P.distance_->distance_squared(Airports[i].loc, Airports[k].loc) > P.MoveRestrRadius2)
{
if (ranf_mt(tn) > P.MoveRestrEffect)
{
@@ -227,7 +227,7 @@ void TravelDepartSweep(double t)
}
if (f3)
{
- double s2 = dist2_raw(Households[Hosts[i2].hh].loc.x, Households[Hosts[i2].hh].loc.y, Places[P.HotelPlaceType][i].loc.x, Places[P.HotelPlaceType][i].loc.y);
+ double s2 = P.distance_->distance_squared(Households[Hosts[i2].hh].loc, Places[P.HotelPlaceType][i].loc);
int f2 = 1;
if ((bm) && (s2 > P.MoveRestrRadius2))
{
@@ -440,8 +440,7 @@ void InfectSweep(double t, int run) //added run number as argument in order to r
{
// if distance between si's household and linked place
// is greater than movement restriction radius
- if ((dist2_raw(Households[si->hh].loc.x, Households[si->hh].loc.y,
- Places[k][l].loc.x, Places[k][l].loc.y) > P.MoveRestrRadius2))
+ if ((P.distance_->distance_squared(Households[si->hh].loc, Places[k][l].loc) > P.MoveRestrRadius2))
{
// multiply infectiousness of place by movement restriction effect
s3 *= P.MoveRestrEffect;
@@ -559,8 +558,7 @@ void InfectSweep(double t, int run) //added run number as argument in order to r
if (bm)
{
// if potential infectee i3's household is further from selected place
- if ((dist2_raw(Households[Hosts[i3].hh].loc.x, Households[Hosts[i3].hh].loc.y,
- Places[k][l].loc.x, Places[k][l].loc.y) > P.MoveRestrRadius2))
+ if ((P.distance_->distance_squared(Households[Hosts[i3].hh].loc, Places[k][l].loc) > P.MoveRestrRadius2))
{
// multiply susceptibility by movement restriction effect
s *= P.MoveRestrEffect;
@@ -661,8 +659,7 @@ void InfectSweep(double t, int run) //added run number as argument in order to r
if (bm)
{
// if potential infectees household is farther away from hotel than restriction radius
- if ((dist2_raw(Households[Hosts[i3].hh].loc.x, Households[Hosts[i3].hh].loc.y,
- Places[k][l].loc.x, Places[k][l].loc.y) > P.MoveRestrRadius2))
+ if ((P.distance_->distance_squared(Households[Hosts[i3].hh].loc, Places[k][l].loc) > P.MoveRestrRadius2))
{
// multiply susceptibility by movement restriction effect
s *= P.MoveRestrEffect;
@@ -889,8 +886,7 @@ void InfectSweep(double t, int run) //added run number as argument in order to r
s *= CalcPersonSusc(i3, ts, ci);
if (bm)
{
- if ((dist2_raw(Households[si->hh].loc.x, Households[si->hh].loc.y,
- Households[Hosts[i3].hh].loc.x, Households[Hosts[i3].hh].loc.y) > P.MoveRestrRadius2))
+ if ((P.distance_->distance_squared(Households[si->hh].loc, Households[Hosts[i3].hh].loc) > P.MoveRestrRadius2))
s *= P.MoveRestrEffect;
}
else if ((mt->moverest != mi->moverest) && ((mt->moverest == 2) || (mi->moverest == 2)))