Skip to content

Commit

Permalink
fixed Linear structure function in Multiple structure
Browse files Browse the repository at this point in the history
  • Loading branch information
cristianlussana committed Sep 27, 2024
1 parent a73950a commit 10d9e33
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 83 deletions.
24 changes: 8 additions & 16 deletions include/gridpp.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1717,8 +1717,9 @@ namespace gridpp {
float test_vec2_argout(vec2& distances);

void test_not_implemented_exception();
/* vec2 test_args_for_R(const Points& bpoints, const StructureFunction& structure, const vec2& background);*/
vec2 test_args_for_R(const Points& bpoints, const vec2& background);
vec2 test_args_for_R(const Points& bpoints, const StructureFunction& structure, const vec2& background);
void test_args_for_R_1( const StructureFunction& structure);
/* vec2 test_args_for_R(const Points& bpoints, const vec2& background); */

/** Default value used to fill array in SWIG testing functions. Not useful for any other purpose. */
static const float swig_default_value = -1;
Expand Down Expand Up @@ -2145,10 +2146,10 @@ namespace gridpp {

/** Linear correlation function
* @param normdist Normalized distance between points. Must be in the range -1 to 1.
* @param min Minimum allowed value for the correlation (if less than 0, the return 1)
* @param length Minimum allowed value for the correlation (if less than 0, the return 1)
* @returns linear rho
*/
float linear_rho(float dist, float length, float min) const;
float linear_rho(float dist, float length) const;
float m_localization_distance;
};
class MultipleStructure: public StructureFunction {
Expand Down Expand Up @@ -2303,37 +2304,28 @@ namespace gridpp {
* @param h: Horizontal decorrelation length >=0 [m]
* @param v: Vertical decorrelation length >=0 [m]. If 0, disable decorrelation.
* @param w: Land/sea decorrelation length >=0 [1]. If 0, disable decorrelation.
* @param hmin: Horizontal decorrelation length >=0 [m]
* @param vmin: Vertical decorrelation length >=0 [m]. If 0, disable decorrelation.
* @param wmin: Land/sea decorrelation length >=0 [1]. If 0, disable decorrelation.
* @param hmax: Truncate horizontal correlation beyond this length [m]. If undefined, 3.64 * h.
*/
LinearStructure(float h, float v=0, float w=0, float hmin=0, float vmin=0, float wmin=0, float hmax=MV);
LinearStructure(float h, float v=0, float w=0, float hmax=MV);

/** Linear structure function where decorrelation varyies spatially
* @param grid: Grid of decorrelation field
* @param h: 2D vector of horizontal decorrelation lengths >=0, same size as grid [m]
* @param v: 2D vector of Vertical decorrelation lengths >=0 [m]. Set all to 0 to disable decorrelation.
* @param w: 2D vector of land/sea decorrelation lengths >=0 [1]. Set all to 0 to disable decorrelation.
* @param hmin: 2D vector of horizontal decorrelation lengths >=0, same size as grid [m]
* @param vmin: 2D vector of Vertical decorrelation lengths >=0 [m]. Set all to 0 to disable decorrelation.
* @param wmin: 2D vector of land/sea decorrelation lengths >=0 [1]. Set all to 0 to disable decorrelation.
* @param min_rho: Truncate horizontal correlation when rho is less than this value [m].
*/
LinearStructure(Grid grid, vec2 h, vec2 v, vec2 w, vec2 hmin, vec2 vmin, vec2 wmin, float min_rho=StructureFunction::default_min_rho);
LinearStructure(Grid grid, vec2 h, vec2 v, vec2 w, float min_rho=StructureFunction::default_min_rho);
float corr(const Point& p1, const Point& p2) const;
vec corr(const Point& p1, const std::vector<Point>& p2) const;
StructureFunctionPtr clone() const;
float localization_distance(const Point& p) const;
private:
float localization_distance(float h, float hmin) const;
float localization_distance(float h) const;
Grid m_grid;
vec2 mH;
vec2 mV;
vec2 mW;
vec2 mHmin;
vec2 mVmin;
vec2 mWmin;
float m_min_rho;
bool m_is_spatial;
};
Expand Down
96 changes: 36 additions & 60 deletions src/api/structure.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,16 +73,16 @@ float gridpp::StructureFunction::powerlaw_rho(float dist, float length) const {
return 1 / (1 + 0.5 * v * v);
}
/* linear function returns correlation between min and 1. normdist between -1 and 1 */
float gridpp::StructureFunction::linear_rho(float dist, float length, float min) const {
if(!gridpp::is_valid(min) || min < 0 || !gridpp::is_valid(length) || length < 0)
float gridpp::StructureFunction::linear_rho(float dist, float length) const {
if(!gridpp::is_valid(length) || length < 0)
// Disabled
return 1;
if(!gridpp::is_valid(dist))
return 0;
float normdist = abs(dist) / length;
float normdist = abs(dist);
if(normdist > 1)
normdist = 1;
return (1 - (1-min) * normdist);
return (1 - (1-length) * normdist);
}
float gridpp::StructureFunction::localization_distance(const Point& p) const {
return m_localization_distance;
Expand Down Expand Up @@ -382,7 +382,7 @@ float gridpp::MixAStructure::corr(const Point& p1, const Point& p2) const {
}
if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) {
float lafdist = p1.laf - p2.laf;
rho *= gridpp::StructureFunction::linear_rho(lafdist, 1., w);
rho *= gridpp::StructureFunction::linear_rho(lafdist, w);
}
}
else {
Expand All @@ -396,7 +396,7 @@ float gridpp::MixAStructure::corr(const Point& p1, const Point& p2) const {
}
if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) {
float lafdist = p1.laf - p2.laf;
rho *= gridpp::StructureFunction::linear_rho(lafdist, 1., mW[0][0]);
rho *= gridpp::StructureFunction::linear_rho(lafdist, mW[0][0]);
}
}
return rho;
Expand Down Expand Up @@ -425,7 +425,7 @@ vec gridpp::MixAStructure::corr(const Point& p1, const std::vector<Point>& p2) c
}
if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2[i].laf)) {
float lafdist = p1.laf - p2[i].laf;
rho *= gridpp::StructureFunction::linear_rho(lafdist, 1., w);
rho *= gridpp::StructureFunction::linear_rho(lafdist, w);
}
}
output[i] = rho;
Expand Down Expand Up @@ -596,11 +596,12 @@ float gridpp::SoarStructure::localization_distance(const Point& p) const {
}
float gridpp::SoarStructure::localization_distance(float h) const {
// To find an approximate analytical solution for the inverse of the function y=(1+x)exp(-x) when y is between 0 and 1, we can follow a similar approach used in approximations for transcendental functions.
// return sqrt(-2*log(m_min_rho)) * h;
if(m_min_rho > 0.5) {
return sqrt( 2 * (1-m_min_rho)) * h;
return sqrt( 1-m_min_rho) * h;
}
else {
return ( -log(m_min_rho) - log(log(m_min_rho))) * h;
return ( -log(m_min_rho) - log( -log(m_min_rho))) * h;
}
}
//-----------------------------------------------------------------------------
Expand Down Expand Up @@ -750,10 +751,10 @@ float gridpp::ToarStructure::localization_distance(const Point& p) const {
float gridpp::ToarStructure::localization_distance(float h) const {
// To find an approximate analytical solution for the inverse of the function y=(1+x+1/3*x^2)exp(-x) when y is between 0 and 1, we can follow a similar approach used in approximations for transcendental functions.
if(m_min_rho > 0.5) {
return sqrt( 6 * (1-m_min_rho)) * h;
return sqrt( 6/5 * (1-m_min_rho)) * h;
}
else {
return ( -log(m_min_rho)) * h;
return ( -log(m_min_rho) - log( -log(m_min_rho))) * h;
}
}
//-----------------------------------------------------------------------------
Expand Down Expand Up @@ -910,30 +911,19 @@ float gridpp::PowerlawStructure::localization_distance(float h) const {
//=============================================================================
/** Linear */
//=============================================================================
gridpp::LinearStructure::LinearStructure(float h, float v, float w, float hmin, float vmin, float wmin, float hmax) :
gridpp::LinearStructure::LinearStructure(float h, float v, float w, float hmax) :
m_is_spatial(false) {
if(gridpp::is_valid(hmax) && hmax < 0)
throw std::invalid_argument("min must be >= 0");
throw std::invalid_argument("hmax must be >= 0");
if(!gridpp::is_valid(h) || h < 0)
throw std::invalid_argument("h must be >= 0");
if(!gridpp::is_valid(v) || v < 0)
throw std::invalid_argument("v must be >= 0");
if(!gridpp::is_valid(w) || w < 0)
throw std::invalid_argument("w must be >= 0");
if(!gridpp::is_valid(hmin) || hmin < 0)
throw std::invalid_argument("hmin must be >= 0");
if(!gridpp::is_valid(vmin) || vmin < 0)
throw std::invalid_argument("vmin must be >= 0");
if(!gridpp::is_valid(wmin) || wmin < 0)
throw std::invalid_argument("wmin must be >= 0");

if(gridpp::is_valid(hmax))
if(hmax > h) {
m_min_rho = hmin;
}
else {
m_min_rho = (1 - (1-hmin) * hmax / h);
}
m_min_rho = default_min_rho;
else
m_min_rho = default_min_rho;
vec2 h2(1);
Expand All @@ -946,23 +936,20 @@ gridpp::LinearStructure::LinearStructure(float h, float v, float w, float hmin,
mV = v2;
mW = w2;
}
gridpp::LinearStructure::LinearStructure(Grid grid, vec2 h, vec2 v, vec2 w, vec2 hmin, vec2 vmin, vec2 wmin, float min_rho) :
gridpp::LinearStructure::LinearStructure(Grid grid, vec2 h, vec2 v, vec2 w, float min_rho) :
m_grid(grid),
m_min_rho(min_rho),
mH(h),
mV(v),
mW(w),
mHmin(h),
mVmin(v),
mWmin(w) {
if(mH.size() == 1 && mH[0].size() == 1 && mV.size() == 1 && mV[0].size() == 1 && mW.size() == 1 && mW[0].size() == 1 && mHmin.size() == 1 && mHmin[0].size() == 1 && mVmin.size() == 1 && mVmin[0].size() == 1 && mWmin.size() == 1 && mWmin[0].size() == 1) {
mW(w) {
if(mH.size() == 1 && mH[0].size() == 1 && mV.size() == 1 && mV[0].size() == 1 && mW.size() == 1 && mW[0].size() == 1) {
m_is_spatial = false;
}
else {
m_is_spatial = true;
if(grid.size()[0] != h.size() || grid.size()[0] != v.size() || grid.size()[0] != w.size() || grid.size()[0] != hmin.size() || grid.size()[0] != vmin.size() || grid.size()[0] != wmin.size())
if(grid.size()[0] != h.size() || grid.size()[0] != v.size() || grid.size()[0] != w.size())
throw std::invalid_argument("Grid size not the same as scale size");
if(grid.size()[1] != h[0].size() || grid.size()[1] != v[0].size() || grid.size()[1] != w[0].size() || grid.size()[1] != hmin[0].size() || grid.size()[1] != vmin[0].size() || grid.size()[1] != wmin[0].size())
if(grid.size()[1] != h[0].size() || grid.size()[1] != v[0].size() || grid.size()[1] != w[0].size())
throw std::invalid_argument("Grid size not the same as scale size");
}
}
Expand All @@ -980,37 +967,34 @@ float gridpp::LinearStructure::corr(const Point& p1, const Point& p2) const {
float h = mH[I[0]][I[1]];
float v = mV[I[0]][I[1]];
float w = mW[I[0]][I[1]];
float hmin = mHmin[I[0]][I[1]];
float vmin = mVmin[I[0]][I[1]];
float wmin = mWmin[I[0]][I[1]];

// Use h to compute this, so we don't have to get nearest neighbour twice
float loc_dist = localization_distance(h, hmin);
float loc_dist = localization_distance(h);
if(hdist > loc_dist)
return 0;

rho = gridpp::StructureFunction::linear_rho(hdist, h, hmin);
rho = gridpp::StructureFunction::linear_rho(hdist, h);
if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) {
float vdist = p1.elev - p2.elev;
rho *= gridpp::StructureFunction::linear_rho(vdist, v, vmin);
rho *= gridpp::StructureFunction::linear_rho(vdist, v);
}
if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) {
float lafdist = p1.laf - p2.laf;
rho *= gridpp::StructureFunction::linear_rho(lafdist, w, wmin);
rho *= gridpp::StructureFunction::linear_rho(lafdist, w);
}
}
else {
if(hdist > localization_distance(p1))
return 0;

rho = gridpp::StructureFunction::linear_rho(hdist, mH[0][0], mHmin[0][0]);
rho = gridpp::StructureFunction::linear_rho(hdist, mH[0][0]);
if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) {
float vdist = p1.elev - p2.elev;
rho *= gridpp::StructureFunction::linear_rho(vdist, mV[0][0], mVmin[0][0]);
rho *= gridpp::StructureFunction::linear_rho(vdist, mV[0][0]);
}
if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) {
float lafdist = p1.laf - p2.laf;
rho *= gridpp::StructureFunction::linear_rho(lafdist, mW[0][0], mWmin[0][0]);
rho *= gridpp::StructureFunction::linear_rho(lafdist, mW[0][0]);
}
}
return rho;
Expand All @@ -1027,22 +1011,19 @@ vec gridpp::LinearStructure::corr(const Point& p1, const std::vector<Point>& p2)
float h = mH[I[0]][I[1]];
float v = mV[I[0]][I[1]];
float w = mW[I[0]][I[1]];
float hmin = mHmin[I[0]][I[1]];
float vmin = mVmin[I[0]][I[1]];
float wmin = mWmin[I[0]][I[1]];
float loc_dist = localization_distance(h, hmin);
float loc_dist = localization_distance(h);
for(int i = 0; i < p2.size(); i++) {
float hdist = gridpp::KDTree::calc_straight_distance(p1, p2[i]);
float rho = 0;
if(hdist <= loc_dist) {
rho = gridpp::StructureFunction::linear_rho(hdist, h, hmin);
rho = gridpp::StructureFunction::linear_rho(hdist, h);
if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2[i].elev)) {
float vdist = p1.elev - p2[i].elev;
rho *= gridpp::StructureFunction::linear_rho(vdist, v, vmin);
rho *= gridpp::StructureFunction::linear_rho(vdist, v);
}
if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2[i].laf)) {
float lafdist = p1.laf - p2[i].laf;
rho *= gridpp::StructureFunction::linear_rho(lafdist, w, wmin);
rho *= gridpp::StructureFunction::linear_rho(lafdist, w);
}
}
output[i] = rho;
Expand All @@ -1056,24 +1037,19 @@ vec gridpp::LinearStructure::corr(const Point& p1, const std::vector<Point>& p2)
return output;
}
gridpp::StructureFunctionPtr gridpp::LinearStructure::clone() const {
return std::make_shared<gridpp::LinearStructure>(m_grid, mH, mV, mW, mHmin, mVmin, mWmin, m_min_rho);
return std::make_shared<gridpp::LinearStructure>(m_grid, mH, mV, mW, m_min_rho);
}
float gridpp::LinearStructure::localization_distance(const Point& p) const {
if(m_is_spatial) {
ivec I = m_grid.get_nearest_neighbour(p.lat, p.lon);
return localization_distance(mH[I[0]][I[1]], mHmin[I[0]][I[1]]);
return localization_distance(mH[I[0]][I[1]]);
}
else {
return localization_distance(mH[0][0], mHmin[0][0]);
return localization_distance(mH[0][0]);
}
}
float gridpp::LinearStructure::localization_distance(float h, float hmin) const {
if( m_min_rho < hmin) {
return h;
}
else {
return (1 - m_min_rho) / (1 - hmin) * h;
}
float gridpp::LinearStructure::localization_distance(float h) const {
return 0;
}
//-----------------------------------------------------------------------------
/** End Linear */
Expand Down
20 changes: 13 additions & 7 deletions src/api/swig.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,14 +99,14 @@ void gridpp::test_not_implemented_exception() {
throw gridpp::not_implemented_exception();
}
vec2 gridpp::test_args_for_R(const gridpp::Points& bpoints,
/* const gridpp::StructureFunction& structure,*/
const gridpp::StructureFunction& structure,
const vec2& background) {
int nS = bpoints.size();
int nEns = background[0].size();
int nY = background.size();
float missing_value = -99999.999;
/* gridpp::StructureFunction structure = gridpp::BarnesStructure(10000, 0, 0, 0);*/
BarnesStructure structure = BarnesStructure(10000,1000000,0.5,100000);
/* gridpp::StructureFunction structure = gridpp::BarnesStructure(10000, 0, 0, 0);
BarnesStructure structure = BarnesStructure(10000,1000000,0.5,100000); */
vec2 output = gridpp::init_vec2(nY, nEns, missing_value);
vec blats = bpoints.get_lats();
vec blons = bpoints.get_lons();
Expand All @@ -115,11 +115,11 @@ vec2 gridpp::test_args_for_R(const gridpp::Points& bpoints,
Point p2 = bpoints.get_point(0);
for(int i = 0; i < nS; i++) {
Point p1 = bpoints.get_point(i);
float localizationRadius = structure.localization_distance(p1);
float rhos = structure.corr_background(p1, p2);
/* float localizationRadius = structure.localization_distance(p1);
float rhos = structure.corr_background(p1, p2); */
std::cout << "----------------------------" << std::endl;
std::cout << i << " localizationRadius " << localizationRadius << std::endl;
std::cout << i << " rhos " << rhos << std::endl;
/* std::cout << i << " localizationRadius " << localizationRadius << std::endl;
std::cout << i << " rhos " << rhos << std::endl; */
std::cout << i << " blats " << blats[i] << std::endl;
std::cout << i << " blons " << blons[i] << std::endl;
std::cout << i << " belevs " << belevs[i] << std::endl;
Expand All @@ -131,3 +131,9 @@ vec2 gridpp::test_args_for_R(const gridpp::Points& bpoints,
}
return output;
}

void gridpp::test_args_for_R_1( const gridpp::StructureFunction& structure) {
float missing_value = -99999.999;
/* gridpp::StructureFunction structure = gridpp::BarnesStructure(10000, 0, 0, 0);
BarnesStructure structure = BarnesStructure(10000,1000000,0.5,100000); */
}

0 comments on commit 10d9e33

Please sign in to comment.