#include "main.h" #include "mesh.h" #include "filter.h" #include "Settings.h" #include "meshPlotter.h" #include "Plotty.h" #include "Plotta.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "FtmKalman.h" #include "misc.h" #include using namespace std::chrono_literals; enum class AggregateMethod { None, Median, MovingMedian }; struct MovingMedianTS2 { private: struct TimeValue { Timestamp timestamp; double value; }; int timeWindow; // ms std::vector values; public: MovingMedianTS2() : timeWindow(0) {} MovingMedianTS2(const Timestamp window) : timeWindow(window.ms()) {} void add(Timestamp ts, double value) { values.push_back(TimeValue{ ts, value }); } bool tryGet(const Timestamp ts, double& value) { int wnd = timeWindow; values.erase(std::remove_if(values.begin(), values.end(), [wnd, ts](TimeValue tv) { return std::abs((ts - tv.timestamp).ms()) >= wnd; }), values.end()); if (values.size() == 0) return false; Stats::Median median; for (auto tv : values) { median.add(tv.value); } value = median.get(); return true; } }; std::vector> getFtmValues(Offline::FileReader& fr, Interpolator& gtInterpolator, const MACAddress nuc) { std::vector> result; for (const Offline::Entry& e : fr.getEntries()) { if (e.type == Offline::Sensor::WIFI_FTM) { const Timestamp ts = Timestamp::fromMS(e.ts); Point3 gtPos = gtInterpolator.get(static_cast(ts.ms())) + Point3(0, 0, 1.3); auto wifi = fr.getWifiFtm()[e.idx].data; if (wifi.getAP().getMAC() == nuc) { Point3 apPos = Settings::CurrentPath.NUCs.find(wifi.getAP().getMAC())->second.position; float apDist = gtPos.getDistance(apPos); float ftmDist = wifi.getFtmDist(); float rssi = wifi.getRSSI(); result.push_back({ apDist, ftmDist, rssi }); } } } return result; } std::pair optimizeFtm(std::vector>& values) { std::vector> error; for (float offset = 0; offset < 10.0f; offset += 0.25) { Stats::Statistics diffs; for (const auto& data : values) { float apDist = std::get<0>(data); float ftmDist = std::get<1>(data); ftmDist += offset; float diff = (apDist - ftmDist); diffs.add(diff); } error.push_back({ offset, diffs.getSquaredSumAvg() }); } auto minElement = std::min_element(error.begin(), error.end(), [](std::pair a, std::pair b) { return a.second < b.second; }); std::cout << "Min ftm offset \t" << minElement->first << "\t" << minElement->second << "\n"; return *minElement; } std::pair optimizeRssi(std::vector>& values) { std::vector> error; for (float pathLoss = 2.0f; pathLoss < 4.0f; pathLoss += 0.125) { Stats::Statistics diffs; for (const auto& data : values) { float apDist = std::get<0>(data); float rssi = std::get<2>(data); float rssiDist = LogDistanceModel::rssiToDistance(-40, pathLoss, rssi); float diff = (apDist - rssiDist); diffs.add(diff); } error.push_back({ pathLoss, diffs.getSquaredSumAvg() }); } auto minElement = std::min_element(error.begin(), error.end(), [](std::pair a, std::pair b) { return a.second < b.second; }); std::cout << "Min path loss \t" << minElement->first << "\t" << minElement->second << "\n"; return *minElement; } void optimizeWifiParameters(Offline::FileReader& fr, Interpolator& gtInterpolator) { int i = 1; for (auto nuc : { Settings::NUC1, Settings::NUC2, Settings::NUC3, Settings::NUC4 }) { auto values = getFtmValues(fr, gtInterpolator, nuc); std::cout << "NUC" << i++ << "\n"; optimizeFtm(values); optimizeRssi(values); } } void exportFtmValues(Offline::FileReader& fr, Interpolator& gtInterpolator) { std::fstream fs; fs.open("test.txt", std::fstream::out); fs << "timestamp;nucid;dist;rssiDist;ftmDist;ftmStdDev;numMeas;numSuccesMeas" << "\n"; for (const Offline::Entry& e : fr.getEntries()) { if (e.type == Offline::Sensor::WIFI_FTM) { const Timestamp ts = Timestamp::fromMS(e.ts); Point3 gtPos = gtInterpolator.get(static_cast(ts.ms())) + Point3(0, 0, 1.3); auto wifi = fr.getWifiFtm()[e.idx].data; int nucid = Settings::CurrentPath.NUCs.at(wifi.getAP().getMAC()).ID; float ftm_offset = Settings::CurrentPath.NUCs.at(wifi.getAP().getMAC()).ftm_offset; float rssi_pathloss = Settings::CurrentPath.NUCs.at(wifi.getAP().getMAC()).rssi_pathloss; float rssiDist = LogDistanceModel::rssiToDistance(-40, rssi_pathloss, wifi.getRSSI()); float ftmDist = wifi.getFtmDist() + ftm_offset; //in m; plus static offset float ftmStdDev = wifi.getFtmDistStd(); int numMeas = wifi.getNumAttemptedMeasurements(); int numSuccessMeas = wifi.getNumSuccessfulMeasurements(); Point3 apPos = Settings::CurrentPath.NUCs.find(wifi.getAP().getMAC())->second.position; float apDist = gtPos.getDistance(apPos); fs << ts.ms() << ";" << nucid << ";" << apDist << ";" << rssiDist << ";" << ftmDist << ";" << ftmStdDev << ";" << numMeas << ";" << numSuccessMeas << "\n"; } } fs.close(); } template struct TimeSeries { std::vector t; std::vector values; void add(const Timestamp ts, const T value) { t.push_back(ts); values.push_back(value); } }; static float kalman_procNoiseDistStdDev = 1.2f; // standard deviation of distance for process noise static float kalman_procNoiseVelStdDev = 0.1f; // standard deviation of velocity for process noise static CombinedStats run(Settings::DataSetup setup, int walkIdx, std::string folder) { // reading file std::string currDir = std::filesystem::current_path().string(); Floorplan::IndoorMap* map = Floorplan::Reader::readFromFile(setup.map); Offline::FileReader fr(setup.training[walkIdx], setup.HasNanoSecondTimestamps); // ground truth std::vector gtPath = setup.gtPath; Interpolator gtInterpolator = fr.getGroundTruthPath(map, gtPath); CombinedStats errorStats; //calculate distance of path std::vector::InterpolatorEntry> gtEntries = gtInterpolator.getEntries(); double gtTotalDistance = 0; Stats::Statistics gtWalkingSpeed; for (int i = 1; i < gtEntries.size(); ++i) { double distance = gtEntries[i].value.getDistance(gtEntries[i - 1].value); double timeDiff = static_cast(gtEntries[i].key - gtEntries[i - 1].key); double walkingSpeed = distance / (timeDiff / 1000.0f); // m / s gtWalkingSpeed.add(walkingSpeed); gtTotalDistance += distance; } std::cout << "Distance of Path: " << gtTotalDistance << std::endl; std::cout << "GT walking speed: " << gtWalkingSpeed.getAvg() << "m/s (" << gtWalkingSpeed.getAvg()*3.6f << "km/h)" << std::endl; // error file const long int t = static_cast(time(NULL)); auto evalDir = std::filesystem::path(Settings::errorDir); evalDir.append(folder); if (!std::filesystem::exists(evalDir)) { std::filesystem::create_directory(evalDir); } std::ofstream errorFile; errorFile.open(evalDir.string() + "/" + std::to_string(walkIdx) + "_" + std::to_string(t) + ".csv"); // Output dir auto outputDir = std::filesystem::path(Settings::outputDir); outputDir.append(Settings::CurrentPath.name + "_" + std::to_string(walkIdx)); if (!std::filesystem::exists(outputDir)) { std::filesystem::create_directories(outputDir); } // wifi auto kalmanMap = std::make_shared>(); for (auto& nucConfig : setup.NUCs) { kalmanMap->insert({ nucConfig.first, Kalman(nucConfig.second.ID, nucConfig.second.kalman_measStdDev, kalman_procNoiseDistStdDev, kalman_procNoiseVelStdDev) }); } std::cout << "Optimal wifi parameters for " << setup.training[walkIdx] << "\n"; optimizeWifiParameters(fr, gtInterpolator); // mesh NM::NavMeshSettings set; set.maxQuality_m = 0.10; // because of narrow hallways and small rooms reduce min. triangle size (default is 0.2) MyNavMesh mesh; MyNavMeshFactory fac(&mesh, set); fac.build(map); // debug show //MeshPlotter dbg; //dbg.addFloors(map); //dbg.addOutline(map); //dbg.addMesh(mesh); ////dbg.addDijkstra(mesh); //dbg.draw(); Plotty plot(map); plot.buildFloorplan(); plot.setGroundTruth(gtPath); plot.setView(30, 0); // APs Positions for (auto& nucConfig : setup.NUCs) { plot.addCircle(10000 + nucConfig.second.ID, nucConfig.second.position.xy(), 0.1); } plot.plot(); // particle-filter const int numParticles = 5000; //auto init = std::make_unique(&mesh, srcPath0); // known position auto init = std::make_unique(&mesh); // uniform distribution auto eval = std::make_unique(); eval->ftmKalmanFilters = kalmanMap; auto trans = std::make_unique(); //auto trans = std::make_unique(); auto resample = std::make_unique>(); auto estimate = std::make_unique>(); //auto estimate = std::make_unique>(); // setup MyFilter pf(numParticles, std::move(init)); pf.setEvaluation(std::move(eval)); pf.setTransition(std::move(trans)); pf.setResampling(std::move(resample)); pf.setEstimation(std::move(estimate)); pf.setNEffThreshold(0.85); // sensors MyControl ctrl; MyObservation obs; Timestamp lastTimestamp = Timestamp::fromMS(0); std::vector errorValuesFtm, errorValuesRssi; std::vector timestamps; std::vector timestampsDist; std::vector> gtDistances, rssiDistances; // distance per AP std::array>, 4> ftmDistances; TimeSeries> ftmPredictions; Plotta::Plotta errorPlot("errorPlot", outputDir.string() + "/errorData.py"); Plotta::Plotta distsPlot("distsPlot", outputDir.string() + "/distances.py"); std::unordered_map movMedianPerAP; for (auto& nucConfig : setup.NUCs) { movMedianPerAP[nucConfig.first] = MovingMedianTS2(Timestamp::fromMS(500)); } for (const Offline::Entry& e : fr.getEntries()) { if (e.type != Offline::Sensor::WIFI_FTM) { continue; } // TIME const Timestamp ts = Timestamp::fromMS(e.ts); auto wifiFtm = fr.getWifiFtm()[e.idx].data; obs.ftm.push_back(wifiFtm); if (ts - lastTimestamp >= Timestamp::fromMS(500)) { // Do update step Point2 gtPos = gtInterpolator.get(static_cast(ts.ms())).xy(); plot.setGroundTruth(Point3(gtPos.x, gtPos.y, 0.1)); // TODO //gtDistances.push_back({ gtPos.getDistance(Settings::CurrentPath.nucInfo(0).position.xy()), // gtPos.getDistance(Settings::CurrentPath.nucInfo(1).position.xy()), // gtPos.getDistance(Settings::CurrentPath.nucInfo(2).position.xy()), // gtPos.getDistance(Settings::CurrentPath.nucInfo(3).position.xy()) }); Point3 estPos; float distErrorFtm = 0; float distErrorRssi = 0; const AggregateMethod aggrMethod = AggregateMethod::None; if (aggrMethod == AggregateMethod::Median) { // Compute median of observations std::unordered_map> apMeas; for (WiFiMeasurement wifi : obs.ftm) { apMeas[wifi.getAP().getMAC()].add(wifi.getFtmDist()); } obs.ftm.clear(); for (auto& pair : apMeas) { double median = pair.second.get(); obs.ftm.push_back(WiFiMeasurement(AccessPoint(pair.first), NAN, ts, median, NAN, 3, 3)); } } else if (aggrMethod == AggregateMethod::MovingMedian) { for (WiFiMeasurement wifi : obs.ftm) { movMedianPerAP[wifi.getAP().getMAC()].add(wifi.getTimestamp(), wifi.getFtmDist()); } obs.ftm.clear(); for (auto& pair : movMedianPerAP) { double median = 0; if (pair.second.tryGet(ts, median)) { obs.ftm.push_back(WiFiMeasurement(AccessPoint(pair.first), NAN, ts, median, NAN, 3, 3)); } } } // Store measurements //for (WiFiMeasurement wifi : obs.ftm) //{ // if (wifi.getNumSuccessfulMeasurements() < 3) // { // continue; // } // Point2 gtPos2 = gtInterpolator.get(static_cast(wifi.getTimestamp().ms())).xy(); // Point2 apPos2 = Settings::CurrentPath.NUCs[wifi.getAP().getMAC()].position.xy(); // float gtDist2 = gtPos2.getDistance(apPos2); // // store distances // const int nucIdx = Settings::nucIndex(wifi.getAP().getMAC()); // ftmDistances[nucIdx].add(wifi.getTimestamp(), { wifi.getFtmDist(), gtDist2, wifi.getRSSI() }); // TODO //} // Kalman debugging (can't be used with active PF) //{ // // Kalman predict & update for available measurments // for (WiFiMeasurement wifi : obs.ftm) // { // kalmanMap->at(wifi.getAP().getMAC()).predictAndUpdate(wifi.getTimestamp(), wifi.getFtmDist()); // } // // Kalman prediction only for current timestamp // std::array pred; // for (size_t i = 0; i < 4; i++) // { // KalmanPrediction prediction = kalmanMap->at(Settings::nucFromIndex(i)).predict(ts); // prediction.P[0] = kalmanMap->at(Settings::nucFromIndex(i)).P[0]; // prediction.P[1] = kalmanMap->at(Settings::nucFromIndex(i)).P[1]; // prediction.P[2] = kalmanMap->at(Settings::nucFromIndex(i)).P[2]; // prediction.P[3] = kalmanMap->at(Settings::nucFromIndex(i)).P[3]; // pred[i] = prediction; // } // ftmPredictions.add(ts, pred); //} // Run PF obs.currentTime = ts; ctrl.currentTime = ts; MyState est = pf.update(&ctrl, obs); ctrl.afterEval(); lastTimestamp = ts; estPos = est.pos.pos; ctrl.lastEstimate = estPos; plot.setCurEst(Point3(estPos.x, estPos.y, 0.1)); plot.addEstimationNode(Point3(estPos.x, estPos.y, 0.1)); // Error distErrorFtm = gtPos.getDistance(estPos.xy()); errorStats.ftm.add(distErrorFtm); // draw wifi ranges plot.clearDistanceCircles(); for (size_t i = 0; i < obs.ftm.size(); i++) { WiFiMeasurement wifi2 = obs.ftm[i]; Point3 apPos = Settings::CurrentPath.nuc(wifi2.getAP().getMAC()).position; K::GnuplotColor color; switch (Settings::CurrentPath.nuc(wifi2.getAP().getMAC()).ID) { case 1: color = K::GnuplotColor::fromRGB(0, 255, 0); break; case 2: color = K::GnuplotColor::fromRGB(0, 0, 255); break; case 3: color = K::GnuplotColor::fromRGB(255, 255, 0); break; default: color = K::GnuplotColor::fromRGB(255, 0, 0); break; } //float rssiDist = LogDistanceModel::rssiToDistance(-40, 2.5, wifi2.getRSSI()); //plot.addDistanceCircle(apPos.xy(), rssiDist, color); plot.addDistanceCircle(apPos.xy(), wifi2.getFtmDist(), color); } obs.wifi.clear(); obs.ftm.clear(); errorValuesFtm.push_back(distErrorFtm); errorValuesRssi.push_back(distErrorRssi); timestamps.push_back(ts.ms()); // Error plot errorPlot.add("t", timestamps); errorPlot.add("errorFtm", errorValuesFtm); errorPlot.add("errorRssi", errorValuesRssi); errorPlot.frame(); // Distances plot //distsPlot.add("t", timestamps); //distsPlot.add("gtDists", gtDistances); //distsPlot.add("ftmDists", ftmDistances); //distsPlot.frame(); // Plotting plot.showParticles(pf.getParticles()); plot.setCurEst(estPos); plot.setGroundTruth(Point3(gtPos.x, gtPos.y, 0.1)); plot.addEstimationNode(estPos); //plot.setActivity((int)act.get()); //plot.splot.getView().setEnabled(false); //plot.splot.getView().setCamera(0, 0); //plot.splot.getView().setEqualXY(true); plot.plot(); std::this_thread::sleep_for(100ms); } } printErrorStats(errorStats); std::ofstream plot_out; plot_out.open(outputDir.string() + "/plot.gp"); plot.clearDistanceCircles(); plot.saveToFile(plot_out); std::ofstream errorStats_out; errorStats_out.open(outputDir.string() + "/error_stats.txt"); printErrorStats(errorStats_out, errorStats); errorPlot.frame(); // MATLAB output //{ // std::ofstream matlab_error_out; // matlab_error_out.open(outputDir.string() + "/error.csv"); // matlab_error_out << "t;ftmError" << "\n"; // for (size_t i = 0; i < timestamps.size(); i++) // { // matlab_error_out << timestamps[i] << ";" << errorValuesFtm[i] << "\n"; // } //} //{ // std::ofstream matlab_gt_out; // matlab_gt_out.open(outputDir.string() + "/distance_gt.csv"); // matlab_gt_out << "t;distGT1;distGT2;distGT3;distGT4" << "\n"; // for (size_t i = 0; i < gtDistances.size(); i++) // { // matlab_gt_out << timestamps[i] << ";" << gtDistances[i][0] << ";" << gtDistances[i][1] << ";" << gtDistances[i][2] << ";" << gtDistances[i][3] << "\n"; // } //} //for (size_t i = 0; i < 4; i++) //{ // std::ofstream matlab_out; // matlab_out.open(outputDir.string() + "/distance_ap" + std::to_string(i+1) + ".csv"); // matlab_out << "t;distAp;distGT;rssi" << "\n"; // for (size_t j = 0; j < ftmDistances[i].t.size(); j++) // { // matlab_out << ftmDistances[i].t[j].ms() // << ";" << ftmDistances[i].values[j][0] // << ";" << ftmDistances[i].values[j][1] // << ";" << ftmDistances[i].values[j][2] // << "\n"; // } //} //{ // std::ofstream matlab_prediction_out; // matlab_prediction_out.open(outputDir.string() + "/predictions.csv"); // matlab_prediction_out << "t;pAP1d;pAP1dDev;pAP1s;pAP1sDev;pAP2d;pAP2dDev;pAP2s;pAP2sDev;pAP3d;pAP3dDev;pAP3s;pAP3sDev;pAP4d;pAP4dDev;pAP4s;pAP4sDev" << "\n"; // for (size_t i = 0; i < ftmPredictions.values.size(); i++) // { // matlab_prediction_out << ftmPredictions.t[i].ms(); // for (size_t j = 0; j < 4; j++) // { // const KalmanPrediction v = ftmPredictions.values[i][j]; // if (isnan(v.distance)) // matlab_prediction_out << ";nan"; // else // matlab_prediction_out << ";" << v.distance; // if (isnan(v.P[0])) // matlab_prediction_out << ";nan"; // else // matlab_prediction_out << ";" << std::sqrt(v.P[0]); // if (isnan(v.speed)) // matlab_prediction_out << ";nan"; // else // matlab_prediction_out << ";" << v.speed; // if (isnan(v.P[2])) // matlab_prediction_out << ";nan"; // else // matlab_prediction_out << ";" << std::sqrt(v.P[2]); // } // matlab_prediction_out << "\n"; // } //} return errorStats; } int main(int argc, char** argv) { CmdArguments args(argc, argv); if (args.hasFlag("prob")) { std::cout << "Probabilistic" << "\n"; return mainProp(argc, argv); } else if (args.hasFlag("trilat")) { std::cout << "Trilateration" << "\n"; return mainTrilat(argc, argv); } CombinedStats statsAVG; CombinedStats statsMedian; CombinedStats statsSTD; CombinedStats statsQuantil; CombinedStats tmp; std::string evaluationName = "prologic/tmp"; std::vector setupsToRun = { //Settings::data.Path5, //Settings::data.Path7, //Settings::data.Path8, //Settings::data.Path9, //Settings::data.Path10, //Settings::data.Path11 //Settings::data.Path20, //Settings::data.Path21, Settings::data.Path22, }; for (Settings::DataSetup setupToRun : setupsToRun) { Settings::CurrentPath = setupToRun; for (size_t walkIdx = 0; walkIdx < Settings::CurrentPath.training.size(); walkIdx++) { std::cout << "Executing walk " << walkIdx << "\n"; for (int i = 0; i < 1; ++i) { std::cout << "Start of iteration " << i << "\n"; tmp = run(Settings::CurrentPath, walkIdx, evaluationName); statsAVG.ftm.add(tmp.ftm.getAvg()); statsMedian.ftm.add(tmp.ftm.getMedian()); statsSTD.ftm.add(tmp.ftm.getStdDev()); statsQuantil.ftm.add(tmp.ftm.getQuantile(0.75)); statsAVG.rssi.add(tmp.rssi.getAvg()); statsMedian.rssi.add(tmp.rssi.getMedian()); statsSTD.rssi.add(tmp.rssi.getStdDev()); statsQuantil.rssi.add(tmp.rssi.getQuantile(0.75)); std::cout << "Iteration " << i << " completed" << std::endl; } } std::cout << "Results for path " << Settings::CurrentPath.name << std::endl; std::cout << "==========================================================" << std::endl; std::cout << "Average of all statistical data FTM: " << std::endl; std::cout << "Median: " << statsMedian.ftm.getAvg() << std::endl; std::cout << "Average: " << statsAVG.ftm.getAvg() << std::endl; std::cout << "Standard Deviation: " << statsSTD.ftm.getAvg() << std::endl; std::cout << "75 Quantil: " << statsQuantil.ftm.getAvg() << std::endl; std::cout << "==========================================================" << std::endl; std::cout << "==========================================================" << std::endl; std::cout << "Average of all statistical data RSSI: " << std::endl; std::cout << "Median: " << statsMedian.rssi.getAvg() << std::endl; std::cout << "Average: " << statsAVG.rssi.getAvg() << std::endl; std::cout << "Standard Deviation: " << statsSTD.rssi.getAvg() << std::endl; std::cout << "75 Quantil: " << statsQuantil.rssi.getAvg() << std::endl; std::cout << "==========================================================" << std::endl; } //std::vector> error; //std::ofstream error_out; //error_out.open(Settings::errorDir + evaluationName + "_error_path1" + ".csv", std::ios_base::app); //for (kalman_procNoiseDistStdDev = 0.8f; kalman_procNoiseDistStdDev < 1.5f; kalman_procNoiseDistStdDev += 0.1f) //{ // for (kalman_procNoiseVelStdDev = 0.1f; kalman_procNoiseVelStdDev < 0.5f; kalman_procNoiseVelStdDev += 0.1f) // { //for (size_t walkIdx = 0; walkIdx < Settings::data.CurrentPath.training.size(); walkIdx++) //{ // std::cout << "Executing walk " << walkIdx << "\n"; // for (int i = 0; i < 1; ++i) // { // std::cout << "Start of iteration " << i << "\n"; // tmp = run(Settings::data.CurrentPath, walkIdx, evaluationName); // statsMedian.add(tmp.getMedian()); // statsAVG.add(tmp.getAvg()); // statsSTD.add(tmp.getStdDev()); // statsQuantil.add(tmp.getQuantile(0.75)); // std::cout << kalman_procNoiseDistStdDev << " " << kalman_procNoiseVelStdDev << std::endl; // std::cout << "Iteration " << i << " completed" << std::endl; // } //} // error.push_back({{ kalman_procNoiseDistStdDev, kalman_procNoiseVelStdDev, statsAVG.getAvg() }}); // auto minElement = std::min_element(error.begin(), error.end(), [](std::array a, std::array b) { // return a[2] < b[2]; // }); // std::cout << "Current min error " << (*minElement)[2] << "\t Q(0)=\t" << (*minElement)[0] << "\t Q(1)=" << (*minElement)[1] << "\n"; // error_out << kalman_procNoiseDistStdDev << ";" << kalman_procNoiseVelStdDev << ";" << statsAVG.getAvg() << std::endl; // // reset stats // statsAVG.reset(); // statsMedian.reset(); // statsSTD.reset(); // statsQuantil.reset(); // } //} //auto minElement = std::min_element(error.begin(), error.end(), [](std::array a, std::array b) { // return a[2] < b[2]; //}); //std::cout << "Global Min error " << (*minElement)[2] << "\t Q(0)=\t" << (*minElement)[0] << "\t Q(1)=" << (*minElement)[1] << "\n"; //error_out.close(); //std::this_thread::sleep_for(std::chrono::seconds(60)); }