Method D (Voronoi) | Calculation of Voronoi velocity

[gjaeger]

JuPedSim - JPSreport

• Version : 0.8.4
• Commit hash : 74db190

Problem description When calculating the Voronoi density and Voronoi speed with the current version I can't get the results of J. Zhang 2012, p. 66, figure 4.6 for UO-080-300-300 (run of experiment of two-dimensional undirectional movement with open boundary conditions (UO), see Pedestrian Dynamics Data Archive).

The comparison of the time series shows that the density differs only minimally:

When comparing the time series for speed, differences are clearly visible:

I had expected that the time series for the speed should be similar.

Expected behavior I had expected that the time series for velocity and density would be similar to the Method C (Classical):

The differences could be numerical.

The two methods are currently implemented as follows:

Method C (Classical) (GetClassicVelocity in Method_C.cpp, line 115 et seqq.)

For the classical method, the spatial mean velocity is the average of the instantaneous velocities for all pedestrians in the measurement area at time:

This is implemented in the program code as follows

double Method_C::GetClassicalVelocity(const vector<double>& xs, const vector<double>& ys, const vector<double>& VInFrame, int pednum) const
{
     int pedsinMeasureArea=0;
     double velocity = 0;
     for(int i=0; i<pednum; i++)
     {
          if(within(make<point_2d>(xs[i], ys[i]), _areaForMethod_C->_poly))
          {
               velocity+=VInFrame[i];
               pedsinMeasureArea++;
          }
     }
     if(pedsinMeasureArea!=0)
     {
          velocity /= (1.0*pedsinMeasureArea);
     } else {
          velocity = 0;
     }
     return velocity;

}

This function was not edited in the last release.

Method D (Voronoi) (GetVoronoiDensityVelocity in Method_D.cpp, line 326 et seqq.)

For the Voronoi method, the speed (Voronoi velocity) for the measurement area is defined as:

std::tuple<double,double> Method_D::GetVoronoiDensityVelocity(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const polygon_2d & measureArea)
{
     double meanV=0;
     double density=0;
     int temp=0;
     for (auto && polygon_iterator:polygon)
     {
          polygon_list v;
          intersection(measureArea, polygon_iterator, v);
          if(!v.empty())
          {
               meanV+=Velocity[temp]*area(v[0]);
               density+=area(v[0])/area(polygon_iterator);
               if((area(v[0]) - area(polygon_iterator))>J_EPS)
               {
                    std::cout<<"----------------------Now calculating density-velocity!!!-----------------\n ";
                    std::cout<<"measure area: \t"<<std::setprecision(16)<<dsv(measureArea)<<"\n";
                    std::cout<<"Original polygon:\t"<<std::setprecision(16)<<dsv(polygon_iterator)<<"\n";
                    std::cout<<"intersected polygon: \t"<<std::setprecision(16)<<dsv(v[0])<<"\n";
                    std::cout<<"this is a wrong result in density calculation\t "<<area(v[0])<<'\t'<<area(polygon_iterator)<<  "  (diff=" << (area(v[0]) - area(polygon_iterator)) << ")" << "\n";
               }
          }
          temp++;
     }
     meanV = meanV/area(measureArea);
     density = density/(area(measureArea)*CMtoM*CMtoM);
     return std::make_tuple(density, meanV);
}

This function was edited in the last release (commits 694e497743ae5ae0f3656100d763e684858f7ae6 and ec35c62d3ccc2687b5f8cabe4c992f956065bd14).

Conclusion

We had agreed in issue #127 (gitlab issue 84) that

• Method D implements the abovementioned formula correctly. So the code is not buggy.
• When to use method D and when not is an important question. I hope it can be discussed in a publication. However, a short notice to warn users can be written in the documentation of method D.
• Other enhancements or changes of these formulae can be implemented. I would suggest a different method for that.

I just can't explain why I can't get the results from Zhang anymore.

@chraibi Were we wrong?

[gjaeger]

Before the modification (commit 694e497), the velocity was calculated as follows:

double Method_D::GetVoronoiVelocity(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const polygon_2d & measureArea)
{
     double meanV=0;
     int temp=0;
     for (auto && polygon_iterator:polygon)
     {
          polygon_list v;
          intersection(measureArea, polygon_iterator, v);
          if(!v.empty())
          {
               meanV+=(Velocity[temp]*area(v[0])/area(measureArea));
               if((area(v[0]) - area(polygon_iterator))>J_EPS)
               {
                    std::cout<<"this is a wrong result in calculating velocity\t"<<area(v[0])<<'\t'<<area(polygon_iterator)<< "  (diff=" << area(v[0]) - area(polygon_iterator) << ")"<< std::endl;
               }
          }
          temp++;
     }
     return meanV;
}

If I modify the function GetVoronoiDensityVelocity (current version, line 326 et seqq.) as follows:

• meanV+=(Velocity[temp]*area(v[0])/area(measureArea)); instead of meanV+=Velocity[temp]*area(v[0]);
• commented out: meanV = meanV/area(measureArea);

std::tuple<double,double> Method_D::GetVoronoiDensityVelocity(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const polygon_2d & measureArea)
{
     double meanV=0;
     double density=0;
     int temp=0;
     for (auto && polygon_iterator:polygon)
     {
          polygon_list v;
          intersection(measureArea, polygon_iterator, v);
          if(!v.empty())
          {
               meanV+=(Velocity[temp]*area(v[0])/area(measureArea));
               density+=area(v[0])/area(polygon_iterator);
               if((area(v[0]) - area(polygon_iterator))>J_EPS)
               {
                    std::cout<<"----------------------Now calculating density-velocity!!!-----------------\n ";
                    std::cout<<"measure area: \t"<<std::setprecision(16)<<dsv(measureArea)<<"\n";
                    std::cout<<"Original polygon:\t"<<std::setprecision(16)<<dsv(polygon_iterator)<<"\n";
                    std::cout<<"intersected polygon: \t"<<std::setprecision(16)<<dsv(v[0])<<"\n";
                    std::cout<<"this is a wrong result in density calculation\t "<<area(v[0])<<'\t'<<area(polygon_iterator)<<  "  (diff=" << (area(v[0]) - area(polygon_iterator)) << ")" << "\n";
               }
          }
          temp++;
     }
     // meanV = meanV/area(measureArea);
     density = density/(area(measureArea)*CMtoM*CMtoM);
     return std::make_tuple(density, meanV);
}

The result does not change:

[gjaeger]

When comparing another case (UO-300-300-120), the differences are not so strong.

and for completeness:

[gjaeger]

the files for analysis:

• UO-080-300-300
  • ini, geo and traj-file UO_080-300-300V2.zip
  • datafile from Zhang 2012 rho_v_uo-080-300-300v2-dat.zip
• UO-300-300-120
  • ini, geo and traj-file UO_300-300-120V2.zip
  • datafile from Zhang 2012 rho_v_uo-300-300-120-dat.zip

Note: The structure of datafiles from Zhang 2012:

• columns Frame | rho_voronoi | rho_classical | v_voronoi | v_classical

[gjaeger]

The result (Method Voronoi) with the proposed pull request #171 corresponds to my expectations for UO-080-300-300:

[gjaeger]

With the modification (PR #171) the result of Zhang can also be reproduced.