Skip to content

Commit

Permalink
Separate out the distance calculations
Browse files Browse the repository at this point in the history
  • Loading branch information
zebmason committed Aug 6, 2020
1 parent 87127e7 commit efd5056
Show file tree
Hide file tree
Showing 12 changed files with 323 additions and 240 deletions.
2 changes: 2 additions & 0 deletions covid-sim.vcxproj
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,7 @@
<ClCompile Include="src\Update.cpp" />
<ClCompile Include="src\Param.cpp" />
<ClCompile Include="src\Direction.cpp" />
<ClCompile Include="src\Geometry\Distance.cpp" />
<ClCompile Include="src\Geometry\Vector2.cpp" />
<ClCompile Include="src\InverseCdf.cpp" />
</ItemGroup>
Expand Down Expand Up @@ -134,6 +135,7 @@
<ClInclude Include="src\MicroCellPosition.hpp" />
<ClInclude Include="src\Direction.hpp" />
<ClInclude Include="include\geometry\BoundingBox.h" />
<ClInclude Include="include\geometry\Distance.h" />
<ClInclude Include="include\geometry\Size.h" />
<ClInclude Include="include\geometry\Vector2.h" />
<ClInclude Include="src\Models\Cell.h" />
Expand Down
2 changes: 1 addition & 1 deletion include/geometry/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1 +1 @@
set(GEOMETRY_HDR_FILES Size.h Vector2.h BoundingBox.h)
set(GEOMETRY_HDR_FILES Size.h Vector2.h BoundingBox.h Distance.h)
35 changes: 35 additions & 0 deletions include/geometry/Distance.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
#ifndef COVIDSIM_GEOMETRY_DISTANCE_H_INCLUDED_
#define COVIDSIM_GEOMETRY_DISTANCE_H_INCLUDED_

#include "BoundingBox.h"
#include "Vector2.h"
#include "Size.h"

#include <memory>

namespace CovidSim
{
namespace Geometry
{
class Distance
{
public:
virtual double distance_squared(int l, int m, int stride, const Size<double>& scales, bool minimum = false) = 0;

virtual double distance_squared(const Vector2d& a, const Vector2d& b) = 0;

double distance_squared(const Vector2f& a, const Vector2f& b)
{
return distance_squared(Vector2d(a.x, a.y), Vector2d(b.x, b.y));
}
};

class DistanceFactory
{
public:
static std::shared_ptr<Distance> create(bool utm, bool periodic_boundaries, const BoundingBox2d& spatial_bounding_box, const Size<double>& in_degrees);
};
}
}

#endif
50 changes: 0 additions & 50 deletions src/Constants.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
21 changes: 5 additions & 16 deletions src/CovidSim.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -2473,7 +2462,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)
Expand Down Expand Up @@ -3208,7 +3197,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]++;
}
Expand Down
167 changes: 14 additions & 153 deletions src/Dist.cpp
Original file line number Diff line number Diff line change
@@ -1,177 +1,38 @@
#include <cstdlib>
#include <cmath>

#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));
}
1 change: 0 additions & 1 deletion src/Dist.h
Original file line number Diff line number Diff line change
Expand Up @@ -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*);
Expand Down
2 changes: 1 addition & 1 deletion src/Geometry/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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})
Loading

0 comments on commit efd5056

Please sign in to comment.