Analysis.cpp

/**
 * \file        Analysis.cpp
 * \date        Oct 10, 2014
 * \version     v0.7
 * \copyright   <2009-2016> Forschungszentrum Juelich GmbH. All rights reserved.
 *
 * \section License
 * This file is part of JuPedSim.
 *
 * JuPedSim is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * any later version.
 *
 * JuPedSim is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with JuPedSim. If not, see <http://www.gnu.org/licenses/>.
 *
 * \section Description
 * The Analysis class represents a process of analyzing groups of pedestrian
 * trajectories from experiment or simulation. Different measurement methods
 * can be used and are defined by various parameters and functions.
 *
 *
 **/

#include "Analysis.h"
#include "geometry/Room.h"
#include "geometry/SubRoom.h"

#include "methods/VoronoiDiagram.h"
#include "methods/Method_A.h"
#include "methods/Method_B.h"
#include "methods/Method_C.h"
#include "methods/Method_D.h"
#include "methods/Method_I.h"
#include "methods/PedData.h"

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cfloat>
#include <stdlib.h>
#include <algorithm>    // std::min_element, std::max_element

#ifdef __linux__
#include <sys/stat.h>
#include <dirent.h>
#elif   __APPLE__
#include <sys/stat.h>
#include <dirent.h>
#else
#include <direct.h>
#endif


using boost::geometry::dsv;
using namespace std;
OutputHandler* Log = new STDIOHandler();

/************************************************
 // Konstruktoren
 ************************************************/

Analysis::Analysis()
{
     _building = NULL;
     _projectRootDir="";
     _deltaF=5;   // half of the time interval that used to calculate instantaneous velocity of ped i. Here v_i = (X(t+deltaF) - X(t+deltaF))/(2*deltaF).   X is location.
     _DoesUseMethodA = false;                                           // Method A (Zhang2011a)
     _DoesUseMethodB = false;                   // Method B (Zhang2011a)
     _DoesUseMethodC = false;                                   // Method C //calculate and save results of classic in separate file
     _DoesUseMethodD = false;                                   // Method D--Voronoi method
     _cutByCircle = false;  //Adjust whether cut each original voronoi cell by a circle
     _getProfile = false;   // Whether make field analysis or not
     _outputGraph = false;   // Whether output the data for plot the fundamental diagram each frame
     _calcIndividualFD = false; //Adjust whether analyze the individual density and velocity of each pedestrian in stationary state (ALWAYS VORONOI-BASED)
     _vComponent = "B"; // to mark whether x, y or x and y coordinate are used when calculating the velocity
     _IgnoreBackwardMovement = false;
     _grid_size_X = 0.10;   // the size of the grid
     _grid_size_Y = 0.10;
     _lowVertexX = 0;// LOWest vertex of the geometry (x coordinate)
     _lowVertexY = 0; //  LOWest vertex of the geometry (y coordinate)
     _highVertexX = 10; // Highest vertex of the geometry
     _highVertexY = 10;

     _cutRadius=1.0;
     _circleEdges=6;
     _trajFormat=FileFormat::FORMAT_PLAIN;
     _isOneDimensional=false;
     _plotGraph=false;
}

Analysis::~Analysis()
{
     delete _building;
}

// file.txt ---> file
std::string Analysis::GetBasename (const std::string& str)
{
     std::cout << "Splitting: " << str << '\n';
     unsigned found = str.find_last_of(".");
     return  str.substr(0,found);
}
// c:\\windows\\winhelp.exe ---> winhelp.exe
std::string Analysis::GetFilename (const std::string& str)
{
     std::cout << "GetFilename: " << str << '\n';
     unsigned found = str.find_last_of("/\\");
     return   str.substr(found+1);
}

void Analysis::InitArgs(ArgumentParser* args)
{
     string s = "Parameter:\n";
     _building = new Building();
     _building->LoadGeometry(args->GetGeometryFilename());
     // create the polygons
     _building->InitGeometry();
     // _building->AddSurroundingRoom();

     if(args->GetIsMethodA()) {
          _DoesUseMethodA = true;
          vector<int> Measurement_Area_IDs = args->GetAreaIDforMethodA();
          for(unsigned int i=0; i<Measurement_Area_IDs.size(); i++)
          {
               _areaForMethod_A.push_back(dynamic_cast<MeasurementArea_L*>( args->GetMeasurementArea(Measurement_Area_IDs[i])));
          }
          _deltaT = args->GetTimeIntervalA();
          _plotTimeseriesA=args->GetIsPlotTimeSeriesA();
     }

     if(args->GetIsMethodB()) {
          _DoesUseMethodB = true;
          vector<int> Measurement_Area_IDs = args->GetAreaIDforMethodB();
          for(unsigned int i=0; i<Measurement_Area_IDs.size(); i++)
          {
               _areaForMethod_B.push_back(dynamic_cast<MeasurementArea_B*>( args->GetMeasurementArea(Measurement_Area_IDs[i])));
          }
     }

     if(args->GetIsMethodC()) {
          _DoesUseMethodC = true;
          vector<int> Measurement_Area_IDs = args->GetAreaIDforMethodC();
          for(unsigned int i=0; i<Measurement_Area_IDs.size(); i++)
          {
               _areaForMethod_C.push_back(dynamic_cast<MeasurementArea_B*>( args->GetMeasurementArea(Measurement_Area_IDs[i])));
          }
          _plotTimeseriesC=args->GetIsPlotTimeSeriesC();
     }

     if(args->GetIsMethodD()) {
          _DoesUseMethodD = true;
          vector<int> Measurement_Area_IDs = args->GetAreaIDforMethodD();
          for(unsigned int i=0; i<Measurement_Area_IDs.size(); i++)
          {
               _areaForMethod_D.push_back(dynamic_cast<MeasurementArea_B*>( args->GetMeasurementArea(Measurement_Area_IDs[i])));
          }
          _StartFramesMethodD = args->GetStartFramesMethodD();
          _StopFramesMethodD = args->GetStopFramesMethodD();
          _IndividualFDFlags = args->GetIndividualFDFlags();
          _plotTimeseriesD=args->GetIsPlotTimeSeriesD();
          _geoPolyMethodD = ReadGeometry(args->GetGeometryFilename(), _areaForMethod_D);

     }
     if(args->GetIsMethodI()) {
          _DoesUseMethodI = true;
          vector<int> Measurement_Area_IDs = args->GetAreaIDforMethodI();
          for(unsigned int i=0; i<Measurement_Area_IDs.size(); i++)
          {
               _areaForMethod_I.push_back(dynamic_cast<MeasurementArea_B*>( args->GetMeasurementArea(Measurement_Area_IDs[i])));
          }
          _StartFramesMethodI = args->GetStartFramesMethodI();
          _StopFramesMethodI = args->GetStopFramesMethodI();
          _IndividualFDFlags = args->GetIndividualFDFlags();
          _plotTimeseriesI=args->GetIsPlotTimeSeriesI();
          _geoPolyMethodI = ReadGeometry(args->GetGeometryFilename(), _areaForMethod_I);
     }

// ToDo: obsolete ?

//     if( _DoesUseMethodD &&  _DoesUseMethodI)
//     {
//          Log->Write("Warning:\t Using both method D and I is not safe!");
//          // because ReadGeometry() may be called twice (line 169 and 182)
//          // overwrite _geoPoly -> new value for each Method
//          Log->Write("Info:\t Using both method D and I is safe! :-)");
//     }

     _deltaF = args->GetDelatT_Vins();
     _cutByCircle = args->GetIsCutByCircle();
     _getProfile = args->GetIsGetProfile();
     _outputGraph = args->GetIsOutputGraph();
     _plotGraph = args->GetIsPlotGraph();
     _plotIndex = args->GetIsPlotIndex();
     _isOneDimensional=args->GetIsOneDimensional();
     _vComponent = args->GetVComponent();
     _IgnoreBackwardMovement =args->GetIgnoreBackwardMovement();
     _grid_size_X = int(args->GetGridSizeX());
     _grid_size_Y = int(args->GetGridSizeY());
     _geometryFileName=args->GetGeometryFilename();
     _projectRootDir=args->GetProjectRootDir();
     _trajFormat=args->GetFileFormat();
     _cutRadius=args->GetCutRadius();
     _circleEdges=args->GetCircleEdges();
     _scriptsLocation=args->GetScriptsLocation();
     _outputLocation=args->GetOutputLocation();
}


std::map<int, polygon_2d> Analysis::ReadGeometry(const fs::path& geometryFile, const std::vector<MeasurementArea_B*>& areas)
{
     Log->Write("INFO:\tReadGeometry with %s", geometryFile.string().c_str());
     double geo_minX  = FLT_MAX;
     double geo_minY  = FLT_MAX;
     double geo_maxX  = -FLT_MAX;
     double geo_maxY  = -FLT_MAX;

     std::map<int, polygon_2d> geoPoly;

     //loop over all areas
     for(auto&& area: areas)
     {
          //search for the subroom that contains that area
          for (auto&& it_room : _building->GetAllRooms())
          {
               for (auto&& it_sub : it_room.second->GetAllSubRooms())
               {
                    SubRoom* subroom = it_sub.second.get();
                    point_2d point(0,0);
                    boost::geometry::centroid(area->_poly,point);
                    //check if the area is contained in the obstacle
                    if(subroom->IsInSubRoom(Point(point.x()/M2CM,point.y()/M2CM)))
                    {
                         for (auto&& tmp_point : subroom->GetPolygon())
                         {
                              append(geoPoly[area->_id], make<point_2d>(tmp_point._x*M2CM, tmp_point._y*M2CM));
                              geo_minX = (tmp_point._x*M2CM<=geo_minX) ? (tmp_point._x*M2CM) : geo_minX;
                              geo_minY = (tmp_point._y*M2CM<=geo_minY) ? (tmp_point._y*M2CM) : geo_minY;
                              geo_maxX = (tmp_point._x*M2CM>=geo_maxX) ? (tmp_point._x*M2CM) : geo_maxX;
                              geo_maxY = (tmp_point._y*M2CM>=geo_maxY) ? (tmp_point._y*M2CM) : geo_maxY;
                         }
                         correct(geoPoly[area->_id]);
                         //cout<<"this is:\t"<<subroom->GetAllObstacles().size()<<endl;
                         //appen the holes/obstacles if any
                         int k=1;
                         for(auto&& obst: subroom->GetAllObstacles())
                         {
                              geoPoly[area->_id].inners().resize(k++);
                              geoPoly[area->_id].inners().back();
                              model::ring<point_2d>& inner = geoPoly[area->_id].inners().back();
                              for(auto&& tmp_point:obst->GetPolygon())
                              {
                                   append(inner, make<point_2d>(tmp_point._x*M2CM, tmp_point._y*M2CM));
                              }
                              correct(geoPoly[area->_id]);
                         }
                    }
               }
          }//room

          if(geoPoly.count(area->_id)==0)
          {
               Log->Write("ERROR: \t No polygon containing the measurement id [%d]", area->_id);
               geoPoly[area->_id] = area->_poly;
          }
     }

     _highVertexX = geo_maxX;
     _highVertexY = geo_maxY;
     _lowVertexX = geo_minX;
     _lowVertexY = geo_minY;
     return geoPoly;
}


int Analysis::RunAnalysis(const fs::path& filename, const fs::path& path)
{
     PedData data;
     if(data.ReadData(_projectRootDir, _outputLocation, path, filename, _trajFormat, _deltaF, _vComponent, _IgnoreBackwardMovement)==false)
     {
          Log->Write("ERROR:\tCould not parse the file %d", filename.c_str());
          return -1;
     }

     //-----------------------------check whether there is pedestrian outside the whole geometry--------------------------------------------
     std::map<int , std::vector<int> > _peds_t=data.GetPedsFrame();
     for(int frameNr = 0; frameNr < data.GetNumFrames(); frameNr++ )
     {
          vector<int> ids=_peds_t[frameNr];
          vector<int> IdInFrame = data.GetIdInFrame(frameNr, ids);
          vector<double> XInFrame = data.GetXInFrame(frameNr, ids);
          vector<double> YInFrame = data.GetYInFrame(frameNr, ids);
          for( unsigned int i=0;i<IdInFrame.size();i++)
          {
               bool IsInBuilding=false;
               for (auto&& it_room : _building->GetAllRooms())
               {
                    for (auto&& it_sub : it_room.second->GetAllSubRooms())
                    {
                         SubRoom* subroom = it_sub.second.get();
                         if(subroom->IsInSubRoom(Point(XInFrame[i]*CMtoM,YInFrame[i]*CMtoM)))
                         {
                              IsInBuilding=true;
                              break;
                         }
                    }
                    if(IsInBuilding)
                    {
                         break;
                    }
               }
               if(false==IsInBuilding)
               {
                    Log->Write("Warning:\tAt %dth frame pedestrian at <x=%.4f, y=%.4f> is not in geometry!", frameNr+data.GetMinFrame(), XInFrame[i]*CMtoM, YInFrame[i]*CMtoM );
               }
          }
     }
     //-----------------------------------------------------------------------------------------------------------------------------------------------------------------

     if(_DoesUseMethodA) //Method A
     {
          if(_areaForMethod_A.empty())
          {
               Log->Write("ERROR: Method A selected with no measurement area!");
               exit(EXIT_FAILURE);
          }
#pragma omp parallel for
          for(long unsigned int i=0; i < _areaForMethod_A.size(); i++)
          {
               Method_A method_A ;
               method_A.SetMeasurementArea(_areaForMethod_A[i]);
               method_A.SetTimeInterval(_deltaT[i]);
               method_A.SetPlotTimeSeries(_plotTimeseriesA[i]);
               bool result_A=method_A.Process(data,_scriptsLocation,_areaForMethod_A[i]->_zPos);
               if(result_A)
               {
                    Log->Write("INFO:\tSuccess with Method A using measurement area id %d!\n",_areaForMethod_A[i]->_id);
                    std::cout << "INFO:\tSuccess with Method A using measurement area id "<< _areaForMethod_A[i]->_id << "\n";
               }
               else
               {
                    Log->Write("INFO:\tFailed with Method A using measurement area id %d!\n",_areaForMethod_A[i]->_id);
               }
          }
     }

     if(_DoesUseMethodB) //Method_B
     {
          if(_areaForMethod_B.empty())
          {
               Log->Write("ERROR: Method B selected with no measurement area!");
               exit(EXIT_FAILURE);
          }

#pragma omp parallel for
          for(long unsigned int i=0; i < _areaForMethod_B.size(); i++)
          {
               Method_B method_B;
               method_B.SetMeasurementArea(_areaForMethod_B[i]);
               bool result_B = method_B.Process(data);
               if(result_B)
               {
                    Log->Write("INFO:\tSuccess with Method B using measurement area id %d!\n",_areaForMethod_B[i]->_id);
                    std::cout << "INFO:\tSuccess with Method B using measurement area id "<< _areaForMethod_B[i]->_id << "\n";
               }
               else
               {
                    Log->Write("INFO:\tFailed with Method B using measurement area id %d!\n",_areaForMethod_B[i]->_id);
               }
          }
     }

     if(_DoesUseMethodC) //Method C
     {
          if(_areaForMethod_C.empty())
          {
               Log->Write("ERROR: Method C selected with no measurement area!");
               exit(EXIT_FAILURE);
          }
#pragma omp parallel for
          for(long unsigned int i=0; i < _areaForMethod_C.size(); i++)
          {
               Method_C method_C;
               method_C.SetMeasurementArea(_areaForMethod_C[i]);
               bool result_C =method_C.Process(data,_areaForMethod_C[i]->_zPos);
               if(result_C)
               {
                    Log->Write("INFO:\tSuccess with Method C using measurement area id %d!\n",_areaForMethod_C[i]->_id);
                    std::cout << "INFO:\tSuccess with Method C using measurement area id "<< _areaForMethod_C[i]->_id << "\n";
                    if(_plotTimeseriesC[i])
                    {
                         string parameters_Timeseries=" " + _scriptsLocation.string()+
                              "/_Plot_timeseries_rho_v.py -p "+ _projectRootDir.string()+VORO_LOCATION + " -n "+filename.string()+
                              " -f "+boost::lexical_cast<std::string>(data.GetFps());
                         parameters_Timeseries = PYTHON + parameters_Timeseries;
                         int res=system(parameters_Timeseries.c_str());
                         Log->Write("INFO:\t time series result: %d ",res);
                    }
               }
               else
               {
                    Log->Write("INFO:\tFailed with Method C using measurement area id %d!\n",_areaForMethod_C[i]->_id);
               }
          }
     }

     if(_DoesUseMethodD) //method_D
     {
          if(_areaForMethod_D.empty())
          {
               Log->Write("ERROR: Method D selected with no measurement area!");
               exit(EXIT_FAILURE);
          }

#pragma omp parallel for
          for(long unsigned int i=0; i<_areaForMethod_D.size(); i++)
          {
               Method_D method_D;
               method_D.SetStartFrame(_StartFramesMethodD[i]);
               method_D.SetStopFrame(_StopFramesMethodD[i]);
               method_D.SetCalculateIndividualFD(_IndividualFDFlags[i]);
               method_D.SetGeometryPolygon(_geoPolyMethodD[_areaForMethod_D[i]->_id]);
               method_D.SetGeometryFileName(_geometryFileName);
               method_D.SetGeometryBoundaries(_lowVertexX, _lowVertexY, _highVertexX, _highVertexY);
               method_D.SetGridSize(_grid_size_X, _grid_size_Y);
               method_D.SetOutputVoronoiCellData(_outputGraph);
               method_D.SetPlotVoronoiGraph(_plotGraph);
               method_D.SetPlotVoronoiIndex(_plotIndex);
               method_D.SetDimensional(_isOneDimensional);
               method_D.SetCalculateProfiles(_getProfile);
               method_D.SetTrajectoriesLocation(path);
               if(_cutByCircle)
               {
                    method_D.Setcutbycircle(_cutRadius, _circleEdges);
               }
               method_D.SetMeasurementArea(_areaForMethod_D[i]);
               bool result_D = method_D.Process(data,_scriptsLocation, _areaForMethod_D[i]->_zPos);
               if(result_D)
               {
                    Log->Write("INFO:\tSuccess with Method D using measurement area id %d!\n",_areaForMethod_D[i]->_id);
                    std::cout << "INFO:\tSuccess with Method D using measurement area id "<< _areaForMethod_D[i]->_id << "\n";
                    if(_plotTimeseriesD[i])
                    {
                         string parameters_Timeseries= " " +_scriptsLocation.string()+"/_Plot_timeseries_rho_v.py -p "+ _projectRootDir.string()+VORO_LOCATION + " -n "+filename.string()+
                              " -f "+boost::lexical_cast<std::string>(data.GetFps());
                         parameters_Timeseries = PYTHON + parameters_Timeseries;
                         std::cout << parameters_Timeseries << "\n;";

                         int res=system(parameters_Timeseries.c_str());
                         Log->Write("INFO:\t time series result: %d ",res);
                    }
               }
               else
               {
                    Log->Write("INFO:\tFailed with Method D using measurement area id %d!\n",_areaForMethod_D[i]->_id);
               }
          }
     }

     if(_DoesUseMethodI) //method_I
     {
          if(_areaForMethod_I.empty())
          {
               Log->Write("ERROR: Method I selected with no measurement area!");
               exit(EXIT_FAILURE);
          }

#pragma omp parallel for
          for(long unsigned int i=0; i<_areaForMethod_I.size(); i++)
          {
               Method_I method_I;
               method_I.SetStartFrame(_StartFramesMethodI[i]);
               method_I.SetStopFrame(_StopFramesMethodI[i]);
               method_I.SetCalculateIndividualFD(_IndividualFDFlags[i]);
               method_I.SetGeometryPolygon(_geoPolyMethodI[_areaForMethod_I[i]->_id]);
               method_I.SetGeometryFileName(_geometryFileName);
               method_I.SetGeometryBoundaries(_lowVertexX, _lowVertexY, _highVertexX, _highVertexY);
               method_I.SetGridSize(_grid_size_X, _grid_size_Y);
               method_I.SetOutputVoronoiCellData(_outputGraph);
               // method_I.SetPlotVoronoiGraph(_plotGraph);
               method_I.SetPlotVoronoiIndex(_plotIndex);
               method_I.SetDimensional(_isOneDimensional);
               // method_I.SetCalculateProfiles(_getProfile);
               method_I.SetTrajectoriesLocation(path);
               if(_cutByCircle)
               {
                    method_I.Setcutbycircle(_cutRadius, _circleEdges);
               }
               method_I.SetMeasurementArea(_areaForMethod_I[i]);
               bool result_I = method_I.Process(data,_scriptsLocation, _areaForMethod_I[i]->_zPos);
               if(result_I)
               {
                    Log->Write("INFO:\tSuccess with Method I using measurement area id %d!\n",_areaForMethod_I[i]->_id);
                    std::cout << "INFO:\tSuccess with Method I using measurement area id "<< _areaForMethod_I[i]->_id << "\n";
                    if(_plotTimeseriesI[i])
                    {
                         string parameters_Timeseries= " " +_scriptsLocation.string()+"/_Plot_timeseries_rho_v.py -p "+ _projectRootDir.string()+VORO_LOCATION + " -n "+filename.string()+
                              " -f "+boost::lexical_cast<std::string>(data.GetFps());
                         parameters_Timeseries = PYTHON + parameters_Timeseries;
                         std::cout << parameters_Timeseries << "\n;";

                         int res=system(parameters_Timeseries.c_str());
                         Log->Write("INFO:\t time series result: %d ",res);
                    }
               }
               else
               {
                    Log->Write("INFO:\tFailed with Method I using measurement area id %d!\n",_areaForMethod_I[i]->_id);
               }
          }
     }


     return 0;
}

FILE* Analysis::CreateFile(const string& filename)
{
     //first try to create the file
     FILE* fHandle= fopen(filename.c_str(),"w");
     if(fHandle) return fHandle;

     unsigned int found=filename.find_last_of("/\\");
     string dir = filename.substr(0,found)+"/";
     //string file= filename.substr(found+1);

     // the directories are probably missing, create it
     if (mkpath((char*)dir.c_str())==-1) {
          Log->Write("ERROR:\tcannot create the directory <%s>",dir.c_str());
          return NULL;
     }
     //second and last attempt
     return fopen(filename.c_str(),"w");
}

#if defined(_WIN32)
// @todo: rewrite using boost
int Analysis::mkpath(char* file_path)
{
     assert(file_path && *file_path);
     char* p;
     for (p=strchr(file_path+1, '/'); p; p=strchr(p+1, '/')) {
          *p='\0';

          if (_mkdir(file_path)==-1) {

               if (errno!=EEXIST) {
                    *p='/';
                    return -1;
               }
          }
          *p='/';
     }
     return 0;
}

#else

// @todo: rewrite using boost
int Analysis::mkpath(char* file_path, mode_t mode)
{
     assert(file_path && *file_path);
     char* p;
     for (p=strchr(file_path+1, '/'); p; p=strchr(p+1, '/')) {
          *p='\0';

          if (mkdir(file_path, mode)==-1) {

               if (errno!=EEXIST) {
                    *p='/';
                    return -1;
               }
          }
          *p='/';
     }
     return 0;
}
// delete
#endif