Method D: velocity calculation

[chraibi]

In Gitlab by @chraibi on Jun 21, 2018, 10:54 [origin]

In this line and this line

the velocity is calculated according to this definition

\quad v_{xy}={v_i(t)}\qquad {\rm{if}} (x,y) \in A_i,

\langle v \rangle_v=\frac{\iint{v_{xy}dxdy}}{b_\text{cor}\cdot\Delta x}.   

Is this correct?

Can we have a simple unit test example to check the correctness of this calculation?

[chraibi]

In Gitlab by @gjaeger on Jun 21, 2018, 10:58

@JuleAdrian FYI

[chraibi]

In Gitlab by @gjaeger on Jun 28, 2018, 11:19

The calculation of the individual velocity and density seems to depend on the measurement area. The output also seems to be limited to the measurement area.

methods/Method_D.cpp line 482 ff:

void Method_D::GetIndividualFD(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const vector<int>& Id, const polygon_2d& measureArea, const string& frid)
{
     double uniquedensity=0;
     double uniquevelocity=0;
     int uniqueId=0;
     int temp=0;
     for (const auto & polygon_iterator:polygon)
     {
          polygon_list v;
          intersection(measureArea, polygon_iterator, v);
          if(!v.empty()) {
               uniquedensity=1.0/(area(polygon_iterator)*CMtoM*CMtoM);
               uniquevelocity=Velocity[temp];
               uniqueId=Id[temp];
               fprintf(_fIndividualFD,"%s\t%d\t%.3f\t%.3f\n",frid.c_str(), uniqueId, uniquedensity,uniquevelocity);
          }
          temp++;
     }
}

line 491 - intersection(measureArea, polygon_iterator, v);: Why is the intersection calculated?

line 493 - uniquedensity=1.0/(area(polygon_iterator)*CMtoM*CMtoM);: Is this just the individual intersection?

[chraibi]

In Gitlab by @chraibi on Jul 25, 2018, 07:43

The calculation of the individual velocity and density seems to depend on the measurement area.

Of course, it does.

Why is the intersection calculated?

Because the calculations are specific to a geometrical polygon. This area can be e.g. a room.

Is this just the individual intersection?

$\rho = 1/A_i$ is the individual density.

[chraibi]

In Gitlab by @gjaeger on Jul 27, 2018, 11:47

The calculation of the velocity also includes agents outside the measurement area, for example (UNI_CORR_500_09, frame 88):

This leads to different values at the beginning:

In summary: The calculation of the velocity is not correct.

[chraibi]

In Gitlab by @gjaeger on Jul 27, 2018, 11:53

Is this just the individual intersection?

$\rho = 1/A_i$ is the individual density.

area(polygon_iterator) is the original polygon and thus independent of the measurement area. Is that correct?

[chraibi]

In Gitlab by @gjaeger on Jul 27, 2018, 12:23

The modifications can be found in the branch issue84.

[chraibi]

In Gitlab by @chraibi on Jul 27, 2018, 12:34

if(!v.empty())

This means that the polygon of agent overlaps with the measurement area. It does not assume that the agent itself is in the area.

I don't know how you calculate the blue data in your figure, but I think that the second condition (ped within area) is stronger than the first one (ped overlaps with the area).

I'm not sure which one is correct.

area(polygon_iterator) is the original polygon and thus independent of the measurement area. Is that correct?

Yes.

[chraibi]

In Gitlab by @gjaeger on Jul 27, 2018, 13:14

I calculate the blue data for all the agents who overlaps with the measurement area.

VEL_IFD = []
VEL = []
RHO = []

for f in range(df['Frame'].min(),df['Frame'].max()):

    meanV01 = []
    meanV02 = []
    density = []

    for i in df[df['Frame'] == f]['PersID'].iteritems():

        polygons_mA = []
        polygons_op = []
        polygons_ip = []

        # measurement area
        poly_ma = df['measurementArea'][i[0]].strip()
        # original polygon
        poly_op = df['polygon_iterator'][i[0]].strip()
        # intersected polygon
        poly_ip = df['v'][i[0]].strip()

        exec("p_mA = %s"%poly_ma)
        exec("p_op = %s"%poly_op)
        exec("p_ip = %s"%poly_ip)

        pmA = locals()['p_mA']
        pop = locals()['p_op']
        pip = locals()['p_ip']

        polygons_mA.append(pmA)
        polygons_op.append(pop)
        polygons_ip.append(pip)

        Poly_mA = np.array(polygons_mA[0]) / 10000
        Poly_op = np.array(polygons_op[0]) / 10000
        Poly_ip = np.array(polygons_ip[0]) / 10000

        # individual calculation
        # velocity
        # jupedsim: meanV+=Velocity[temp]*area(v[0]);
        vel = df['vel'][i[0]] * Polygon(polygons_ip[0]).area()

        meanV01.append(vel)
        meanV02.append(df['vel'][i[0]])

        # density
        # jupedsim: density+=area(v[0])/area(polygon_iterator);
        density_ratio = Polygon(polygons_ip[0]).area() / Polygon(polygons_op[0]).area()
        density.append(density_ratio)

    # calculation
    # velocity
    mV1 = np.array(meanV01)
    mV2 = np.array(meanV02)
    # jpsreport: meanV/area(measureArea);
    meanV1 = mV1.sum() / Polygon(polygons_mA[0]).area()
    meanV2 = mV2.mean()

    # jpsreport: density = density/(area(measureArea)*CMtoM*CMtoM);    
    rho = np.array(density)
    rho_mA = rho.sum() / (Polygon(polygons_mA[0]).area() * CMtoM * CMtoM)

    RHO.append(rho_mA)
    VEL.append(meanV2)
    VEL_IFD.append(meanV1)

    print('mean velocity IFD: {0:.3f} m/s'.format(meanV1))
    print('mean velocity: {0:.3f} m/s'.format(meanV2))
    print('mean density:  {0:.3f} P/m2'.format(rho_mA))

See also snippet.

[chraibi]

In Gitlab by @chraibi on Jul 28, 2018, 04:52

In this line you calculate a velocity, which eventually would not be calculated in jpsreport. (because v.empty())

Does this explain why your blue line does not start from 0, as it normally should be?

PS: For readability,​ I edited your code as a snippet. Hope its :cool:

[chraibi]

In Gitlab by @chraibi on Jul 28, 2018, 04:53

changed time estimate to 1w

[chraibi]

In Gitlab by @gjaeger on Jul 30, 2018, 12:22

PS: For readability,​ I edited your code as a snippet. Hope its :cool:

It's not :cool:. You have to open the script in a new window. This disturbs the reading.

The calculation in line 43 and line 58

       vel = df['vel'][i[0]] * Polygon(polygons_ip[0]).area()

meanV1 = mV1.sum() / Polygon(polygons_mA[0]).area()

based on Method_D.cpp, line 323:

meanV+=Velocity[temp]*area(v[0])

This is not the explanation why the mean velocity starts at 0.

[chraibi]

In Gitlab by @chraibi on Jul 31, 2018, 04:37

I don't understand your meaning.

Your reimplementation of this part of method D is not accurate in my opinion. method D considers only polygons that overlap with the measurement area, while your calculation does not.

Is that correct?

PS: your code is back in the comment as plain text. ;-)

[chraibi]

In Gitlab by @gjaeger on Jul 31, 2018, 08:37

Your reimplementation of this part of method D is not accurate in my opinion. method D considers only polygons that overlap with the measurement area, while your calculation does not.

My implementation also considers only polygons that overlap with the measurement area.

Would you like to see the individual plots?

[chraibi]

In Gitlab by @chraibi on Jul 31, 2018, 09:03

show me please the line in the code that considers these..

[chraibi]

In Gitlab by @gjaeger on Jul 31, 2018, 11:20

There is not one line. In the first step I merge the individual data with the trajectories. This dataframe contains only data from agents that come in contact with the measurement area. Then I compare the calculation of the mean velocity of these agents/people:

    # calculation
    # velocity
    mV1 = np.array(meanV01)
    mV2 = np.array(meanV02)
    # jpsreport: meanV/area(measureArea);
    meanV1 = mV1.sum() / Polygon(polygons_mA[0]).area()
    meanV2 = mV2.mean()

[chraibi]

In Gitlab by @gjaeger on Jul 31, 2018, 11:44

My reimplementation (voronoi velocity (green line) based on individual data from the folder _/Output/Fundamental_Diagram/IndividualFD) matches with the classical calculation (classical voronoi velocity (red line) based on the data from the folder _/Output/Fundamental_Diagram/ClassicalVoronoi).

I don't understand why the voronoi velocity is calculated depending on the intersection.

[chraibi]

In Gitlab by @gjaeger on Aug 3, 2018, 07:57

After our phonecall: The velocity definition

\quad v_{xy}={v_i(t)}\qquad {\rm{if}} (x,y) \in A_i,

\langle v \rangle_v=\frac{\iint{v_{xy}dxdy}}{b_\text{cor}\cdot\Delta x}.   

is implemented correctly in method D (line 312 ff.).

[chraibi]

In Gitlab by @gjaeger on Aug 3, 2018, 08:42

The question remains whether the definition is correct.

The method calculates the mean velocity as a function of the intersection. Should not the mean velocity of the agents associated with the measurement area be determined?

 \bar{v}  =\frac{1}{n} \sum_{i = 1}^{n} v_i.   

This question is illustrated by a simple example:

An agent passes through a measurement area at a constant speed (1 m/s). The mean velocity in the measurement area should be 1 m/s if the agent moves over it.

simulation:

measurement area:

v-t-diagram:

• red line: based on method D output data (for plot_time_series v-t-diagram)
• blue line: based on new definition (see above)

[chraibi]

In Gitlab by @gjaeger on Sep 13, 2018, 17:43

The files are available here (expires on 30.09.2018)

[chraibi]

In Gitlab by @gjaeger on Sep 14, 2018, 08:23

comment to: https://gitlab.version.fz-juelich.de/jupedsim/jpsreport/issues/84#note_17069

• FDS file: UNI_CORR_500_09.fds
• trajectories: traj_UNI_CORR_500_09_0001.txt

[chraibi]

In Gitlab by @chraibi on Sep 15, 2018, 09:19

I think method D, as published in Jun's articles, is indeed different from the "classical" definition you are suggesting here.

However, there is a logic in Jun's definition. In analogy to the density calculation, the speed is a continuous function and not a step function.

It considers pedestrian's "contribution" to the overall speed (analogy to Voronoi density).

For one pedestrian moving with 1 m/s the calculated speed looks like.

The speed is maximal when the pedestrian is fully contained in the area.

Please note, the velocity is not exactly 1 m/s, because the overlap area is slightly smaller than the measurement area. See also

By the way, the step function you are showing is given in the file Individual_FD/*.dat, but that's not important.

Now, what happens if we make the measurement area bigger? (probably that's what you did in your test)

Here you can see if the area is enlarged in the y-axis we get different results for the same test scenario.

To summarise, I think Jun's definition of a "Voronoi speed" is good, however, it has limitations when the measurement area is bigger than the pedestrian's size (see case two).

Therefore, in my opinion,​ it needs to be changed in a way to be less sensitive to measurement areas.

For example, what about this:

\quad v_{xy}={v_i(t)}\qquad {\rm{if}} (x,y) \in A_i,

\langle v \rangle_v=\frac{\iint{v_{xy}dxdy}}{\min(P_i, b_\text{cor}\cdot\Delta x)},

where $P_i$ is the area of the polygon of pedestrian $i$.

@PaulGeoerg, @JuleAdrian, @gjaeger: Any suggestions?

---

Here are the files for my test.

Note in the simulation we need at least 4 pedestrians, otherwise jpsreport won't calculate Voronoi densities and velocities. See this line

• geometry file

• jpscore inifile

• jpsreport inifile

• trajectory file

[chraibi]

In Gitlab by @gjaeger on Sep 15, 2018, 11:17

The course of the voronoi velocity relative to the measurement area suggests an acceleration.

Back to the original question:

I think, Method D is suitable for considering a steady state. The manual should be supplemented with a note. Also, the limitations should be clearly described in the manual.

The velocity calculation as described above is not suitable for the analysis of a time series.

For a time series analysis, jpsreport need a new method for data output.

[chraibi]

In Gitlab by @chraibi on Sep 18, 2018, 08:55

We can document that this method, although correct, does not fit all cases, especially when the measurement area is big.

Apart from that, any suggestions on​ how we can generalize it?

[chraibi]

In Gitlab by @gjaeger on Sep 18, 2018, 09:35

We can document that this method, although correct, does not fit all cases, especially when the measurement area is big.

I would not say that the method is correct. It has its application limits.

Apart from that, any suggestions on​ how we can generalize it?

Yes, that's another issue (follows soon).

[chraibi]

In Gitlab by @gjaeger on Oct 9, 2018, 18:00

Do we want to take the knowledge from issue #87? if yes, how?

[chraibi]

In Gitlab by @chraibi on Oct 9, 2018, 18:15

If I may summarise (please correct me, if I forgot something):

• 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.

Do you think this branch should be merged in develop? It means that the output of individual_FD/IndividualFDtaj*.dat contains the different polygons.

[chraibi]

In Gitlab by @gjaeger on Oct 9, 2018, 18:34

For your summary, I have no supplement.

[chraibi]

In Gitlab by @gjaeger on Oct 9, 2018, 18:35

Do you think this branch should be merged in develop? It means that the output of individual_FD/IndividualFDtaj*.dat contains the different polygons.

I used the different polygons to improve the method D. For the develop branch, I would reduce the output to the original polygon and the intersected polygon.

[chraibi]

In Gitlab by @chraibi on Oct 9, 2018, 18:47

mentioned in commit 4627d35ab1c502a421b1faf744f50797cdf77c08

[chraibi]

In Gitlab by @chraibi on Oct 9, 2018, 18:49

closed via merge request !12

[chraibi]

In Gitlab by @chraibi on Oct 9, 2018, 18:49

closed via commit 4627d35ab1c502a421b1faf744f50797cdf77c08

[chraibi]

In Gitlab by @chraibi on Oct 9, 2018, 18:49

mentioned in commit 3dc868aea0916d7185247da310b20e0c5fac3da6