匹配子轨迹

<pre name="code" class="cpp">#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <cmath>
#include <dirent.h>
#include <string>
#include <fstream>
#include <cmath>
#include <climits>

using namespace std;


const double pi = acos(-1.0);
const double INF = 1e10;
const int Max = 4010;
int minNum = 500;
int maxNum = 2000;
int num = 0;
string rem[Max];

int numtra = 2000;//the number of trajectories
double coe1 = 0.2;
double coe2 = 0.2;


int len = 80;//the length of the subtrajectories
double down_ratio = 0.25;// the ratio of the down sample
double aggravate = 6.0;// the intensity of the interference
int TYPE = 0; //the type of the noisy
double step = 1.3;//the step of the learning



struct point{
    double x;
    double y;
    point(double x1,double y1){
        x=x1;
        y=y1;
    }
    point(){};
};
const int LEN = 100000;
vector<vector<point> > trajecs(LEN);
vector<vector<point> > subs(LEN);
int be[LEN];
double mean[Max];
int beginA[Max][Max];
int beginB[Max][Max];
double dpA[Max][Max];
double dpB[Max][Max];

double dp[Max][Max];
int begin[Max][Max];

double dis(const point& p1,const point& p2)
{
    return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
double dis(const double& a,const double& b)
{
    return fabs(a - b);
}
double randomDouble()
{
    double den = 2147483647;
    return rand()/den;
}
double randNormal(double std)//the mean value is equal to zero;
{
    double t1 = randomDouble();
    double t2 = randomDouble();
    return std*sqrt(-2.0*log(t1))*cos(2.0*pi*t2);
}
double randUniform(double maxmum)
{
    double section = maxmum*2;
    double int_max = INT_MAX;
    int_max/=section;
    double ans = rand()/int_max - maxmum;
    return ans;
}
int getIntersection(int l1,int r1,int l2,int r2)
{
    if(r2<l1||r1<l2)
      return 0;
    return min(r1,r2)-max(l1,l2)+1;
}
void getRandSubtra(vector<point>& tra,vector<point>& sub,int len,int& be)
{
    len = min(int(tra.size()),len);
    sub.resize(len);
    be = rand()%(tra.size() - len + 1);
    for(int i = 0;i < len;i++)
    {
        sub[i].x = tra[i+be].x ;
        sub[i].y = tra[i+be].y ;
    }


}
void add_noise(vector<point>& vec,double down_ratio,double mean,int type)
{
    int newnum = vec.size() * (1.0 - down_ratio);
    int getout = vec.size() - newnum;
    vector<double> down;
    vector<int> help(vec.size(), 0);
    vector<bool> rem(vec.size(), 1);
    for (int j = 0; j < vec.size(); j++)
    {
        help[j] = j;
    }
    for (int j = 0; j < getout; j++)
    {
        int position = rand() % (vec.size() - j);
        rem[help[position]] = 0;
        swap(help[position], help[vec.size() - j - 1]);
    }
    vector<point> now(newnum);
    int add = 0;
    for (int j = 0; j < vec.size(); j++)
    {
        if (rem[j])
            now[add++] = vec[j];
    }
    for(int i=0;i<now.size();i++)
    {
        if(type == 0)
        {
            now[i].x += randNormal(mean*aggravate);
            now[i].y += randNormal(mean*aggravate);
        }
        else if(type == 1)
        {
            now[i].x += randUniform(mean*aggravate);
            now[i].y += randUniform(mean*aggravate);
        }
        
    }
    swap(vec, now);
}
void getTramean(vector<point>& vec,double& mean)
{
    int n = vec.size();
    mean = 0;
    for(int i = 2;i < n;i++)
    {
        mean += dis(vec[i],vec[i-1]);
    }
    mean /= (n-1);
}

void getNearestSub2(vector<point>& tra,vector<point>& sub,int& beg,int& end)
{
    memset(begin,-1,sizeof begin);
    int n=tra.size();
    int m=sub.size();
    for(int i=0;i<n;i++)
    {
        dp[0][i]=dis(sub[0],tra[i]);
        begin[0][i]=i;
    }
    for(int i=1;i<m;i++)
    {
        dp[i][0]=dp[i-1][0]+dis(sub[i],tra[0]);
        begin[i][0]=0;
    }
    for(int i=1;i<m;i++)
    for(int j=1;j<n;j++)
    {
        if(dp[i-1][j-1]<=dp[i-1][j] && dp[i-1][j-1]<=dp[i][j-1])
        {
            dp[i][j]=dp[i-1][j-1];
            begin[i][j]=begin[i-1][j-1];
        }
        else if(dp[i-1][j]<=dp[i-1][j-1] && dp[i-1][j]<=dp[i][j-1])
        {
            dp[i][j]=dp[i-1][j];
            begin[i][j]=begin[i-1][j];
        }
        else
        {
            dp[i][j]=dp[i][j-1];
            begin[i][j]=begin[i][j-1];
        }
        dp[i][j]+=dis(sub[i],tra[j]);
    }
    double ans=INF;
    for(int i=0;i<n;i++)
    {
        if(ans>dp[m-1][i])
        {
            ans=dp[m-1][i];
            beg=begin[m-1][i];
            end=i;
        }
    }
}
void getNearestSub3(vector<point>& tra,vector<point>& sub,int& beg,int& end)
{
    double ans = INF;
    int n=tra.size();
    int m=sub.size();
    for(int i = 0;i<n-m+1;i++)
    {
        double disnow = 0;
        for(int j = 0;j<sub.size();j++)
        {
            disnow += dis(sub[j],tra[j + i]);
        }
        if(ans > disnow)
        {
            ans = disnow;
            beg = i;
            end = i + m - 1;
        }
    }
}
void getNearestSub(vector<point>& tra,vector<point>& sub,int& beg,int& end)
{
    int n = tra.size();
    int m = sub.size();
    vector<double> prefix(m+1,0);
    for(int i = 1;i < m ;i++)
    {
        prefix[i] = prefix[i - 1] + dis(sub[i],sub[i - 1]);
    }
    dpB[0][0] = dis(sub[0],tra[0])*coe1;
    beginB[0][0] = 0;
    for(int i = 1;i < m;i++)
    {
        double be = dpA[i-1][0]+dis(sub[i],sub[i-1]);
        double af = dpB[i-1][0]+dis(sub[i],tra[0])*coe1;
        if(be <= af)
        {
            dpA[i][0] = be;
            beginA[i][0] = beginA[i-1][0];
        }
        else
        {
            dpA[i][0] = af;
            beginA[i][0] = 0;
        }
        
        dpB[i][0] = prefix[i]+dis(sub[i],tra[0])*coe1;
        beginB[i][0] = 0;
    }
    for(int i = 0;i < n;i++)
    {
        dpA[0][i] = 0;
        beginA[0][i] = -1;
        dpB[0][i] = dis(sub[0],tra[i])*coe1;
        beginB[0][i] = i;
    }
    for(int i = 1;i < m;i++)
    for(int j = 1;j < n;j++)
    {
        double be = dpA[i - 1][j] + dis(sub[i - 1],sub[i]);
        double af = dpB[i - 1][j] + dis(sub[i],tra[j]) * coe1;
        if(be<=af)
        {
            dpA[i][j] = be;
            beginA[i][j] = beginA[i - 1][j];
        }
        else
        {
            dpA[i][j] = af;
            beginA[i][j] = beginB[i - 1][j];
        }


        be = dpA[i][j - 1]+dis(sub[i],tra[j])*coe1;
        af = dpB[i][j - 1]+dis(tra[j],tra[j - 1])*coe1;//ceo2
        if(be<=af)
        {
            dpB[i][j] = be;
            beginB[i][j] = beginA[i][j - 1];
            if(beginA[i][j - 1] == -1)
            {
                beginB[i][j] = j;
            }
        }
        else
        {
            dpB[i][j] = af;
            beginB[i][j] = beginB[i][j - 1];
        }
    }
    double ans = prefix[m - 1];
    for(int i= 0;i < n;i++)
    {
        if(dpA[m - 1][i] < ans)
        {
            ans = dpA[m - 1][i];
            beg = beginA[m - 1][i];
            end = i;
        }
    }
}

double getlen(vector<double>& tra)
{
    double ans = 0;
    for(int i = 2;i<tra.size();i++)
    {
        ans += dis(tra[i],tra[i - 1]);
    }
    return ans;
}

int ma=0;
void readfile(string path)
{
    rem[num] = path;
    ifstream file;
    file.open(path.c_str());
    if(!file.is_open())
    {
        cout<<"not open:  "<<path<<endl;
    }
    char s[5000];
    for(int i = 0;i<6;i++)
        file.getline(s,400);
    double x,y;
    char help[1000];
    vector<point> vec;
    vec.reserve(3000);
    point fir(0,0);
    bool flag = 1;
    while(file.getline(s,1000))  
    {
        sscanf(s,"%lf,%lf%s",&x,&y,help);
        point p(x,y);
        vec.push_back(p);
        if(flag)
        {
            fir.x = x;
            fir.y = y;
            flag = 0;
        }
        vec[vec.size()-1].x = (vec[vec.size()-1].x - fir.x)*1000;
        vec[vec.size()-1].y = (vec[vec.size()-1].y - fir.y)*1000;
        if(vec.size() > maxNum)
         return;
    }
    if(vec.size() >= minNum && vec.size() <= maxNum)
    {
        ma=max(ma,int(vec.size()));
        swap(trajecs[num],vec);
        getTramean(trajecs[num],mean[num]);
        num++;
    }
            
}
void get_file()
{
    DIR * dir;
    struct dirent * ptr;
    int i;
    string path = "/home/xiefubao/myproject/Data";
    char s[10];
    for(int i = 1;i < 20;i++)
    {
        sprintf(s,"/%3d",i);
        if(i < 100)
            s[1] = '0';
        if(i < 10)
            s[2] = '0';
        string nowpath = path +string(s)+"/Trajectory";
        dir = opendir(nowpath.c_str());
        if(dir == NULL)
        {
            cout<< "not exist: " <<nowpath<<endl;
        }
        int pp=0;
        while((ptr = readdir(dir)) != NULL)
        {
            if(ptr->d_type == 4)
                continue;
            char filename[100];
            sprintf(filename,"/%s", ptr->d_name);
            readfile(nowpath + string(filename));
            if(num == numtra)
                 break;
        }
        closedir(dir);
        if(num == numtra)
           break;
    }
    
}
double out=0;
double return_ans()
{
    double md1 = 0;
    double md2 = 0;
    for(int i = 0;i<num;i++)
    {
        int begin = 0,end = 0 ;
        getNearestSub(trajecs[i],subs[i],begin,end);
        int inter = getIntersection(be[i],be[i]+len-1,begin,end);
        md1 += inter*1.0/len;
        md2 += inter*1.0/(end - begin + 1);
    }
    md1/=num;
    md2/=num;
    return  2*md1*md2/(md1+md2);
}
void find_best_coe2(double& c,ofstream& file)
{
    double best_ans = 0;
    double ans = 0;
    int times_of_down = 0;
    double be = 0;
    for(c = 1.0;c>=0.01;c-=0.01)
    {
        double now = return_ans();
        if(best_ans<now)
        {
            best_ans = now;
            out = now;
            ans = c;
        }
        cout<<now<<" ";
        file<<now<<" ";

        if(now>=be)
            times_of_down = 0;
        else 
        {
            times_of_down++;
            if(times_of_down == 300)
            {
                break;
            }
        }
        be = now;
    }
    cout<<endl;
    file<<endl;
    c = ans;
}

void find_best_coe(double& c,ofstream& file)
{
   double left = 0.001;
   double right = 0.8;
   while(right - left>=0.01)
   {
       double sec = (right - left)/3;
       double middle1 = left + sec*1;
       double middle2 = left + sec*2;
       c = middle1;
       double ans_middle1 = return_ans();
       c = middle2;
       double ans_middle2 = return_ans();
       if(ans_middle1<ans_middle2)
       {
           left = middle1;
           c = middle1;
       }
       else
       {
           right = middle2;
           c = middle2;
       }
        cout<<"coe1: "<<coe1<<" now: "<<max(ans_middle1,ans_middle2)<<endl;

        file<<"coe1: "<<coe1<<" now: "<<max(ans_middle1,ans_middle2)<<endl;
        out = max(ans_middle1,ans_middle2);
   }
   c = left;
}
void output(ofstream& file)
{
    double md1 = 0;
    double md2 = 0;

    double dtw1 = 0;
    double dtw2 = 0;

    double ed1 = 0;
    double ed2 = 0;

    for(int t=0;t<1;t++)
    {
        find_best_coe2(coe1,file);
        cout<<"the learning has gone a loop: the following is the result: "<<endl;
        cout<<"mean: "<<out<<endl;
        cout<<"coe1: "<<coe1<<endl<<endl;

        file<<"the learning has gone a loop: the following is the result: "<<endl;
        file<<"mean: "<<out<<endl;
        file<<"coe1: "<<coe1<<endl<<endl;
    }
    for(int i = 0;i<num;i++)
    {
        vector<point> sub;
        int begin = 0,end = 0 ;

        getNearestSub(trajecs[i],subs[i],begin,end);
        int inter = getIntersection(be[i],be[i]+len-1,begin,end);
        md1 += inter*1.0/len;
        md2 += inter*1.0/(end - begin + 1);

        //cout<<be[i]<<" "<<be[i]+len-1<<endl;
        //cout<<begin<<" "<<end<<endl;


        getNearestSub2(trajecs[i],subs[i],begin,end);
        inter = getIntersection(be[i],be[i]+len-1,begin,end);
        //cout<<begin<<" "<<end<<endl;
        dtw1 += inter*1.0/len;
        dtw2 += inter*1.0/(end - begin + 1);

        getNearestSub3(trajecs[i],subs[i],begin,end);
        inter = getIntersection(be[i],be[i]+len-1,begin,end);
        //cout<<begin<<" "<<end<<endl<<endl;
        ed1 += inter*1.0/len;
        ed2 += inter*1.0/(end - begin + 1);
    }
    md1/=num;
    md2/=num;
    cout<<endl<<coe1<<endl<<"md recall rate: "<<md1<<endl;
    cout<<"md accuracy rate: "<<md2<<endl;
    cout<<"md mean rate: "<<2*md1*md2/(md1+md2)<<endl<<endl;
    file<<endl<<coe1<<endl<<"md recall rate: "<<md1<<endl;
    file<<"md accuracy rate: "<<md2<<endl;
    file<<"md mean rate: "<<2*md1*md2/(md1+md2)<<endl<<endl;

    dtw1/=num;
    dtw2/=num;
    cout<<"dtw recall rate: "<<dtw1<<endl;
    cout<<"dtw accuracy rate: "<<dtw2<<endl;
    cout<<"dtw mean rate: "<<2*dtw1*dtw2/(dtw1+dtw2)<<endl<<endl;
    file<<"dtw recall rate: "<<dtw1<<endl;
    file<<"dtw accuracy rate: "<<dtw2<<endl;
    file<<"dtw mean rate: "<<2*dtw1*dtw2/(dtw1+dtw2)<<endl<<endl;

    ed1/=num;
    ed2/=num;
    cout<<"dtw recall rate: "<<ed1<<endl;
    cout<<"dtw accuracy rate: "<<ed2<<endl;
    cout<<"dtw mean rate: "<<2*ed1*ed2/(ed1+ed2)<<endl<<endl;
    file<<"dtw recall rate: "<<ed1<<endl;
    file<<"dtw accuracy rate: "<<ed2<<endl;
    file<<"dtw mean rate: "<<2*ed1*ed2/(ed1+ed2)<<endl<<endl;
}
void read_data()
{
    for(int i=0;i<num;i++)
    {

        getRandSubtra(trajecs[i],subs[i],len,be[i]);
        add_noise(subs[i],down_ratio,mean[i],TYPE);
    }
}
int main()
{
    //srand(time(0));
    get_file();
    ofstream file;
    file.open("rem2.txt");
    
    for(down_ratio = 0.1; down_ratio<=0.6; down_ratio+=0.1)
      for(aggravate = 1;aggravate<=20;aggravate += 2)
        {
            cout<<"****************************************************************************************"<<endl;
            cout<<"****************************************************************************************"<<endl;
            cout<<"down_ratio: "<<down_ratio<<endl;
            cout<<"aggravate: "<<aggravate<<endl;

            file<<"****************************************************************************************"<<endl;
            file<<"****************************************************************************************"<<endl;
            file<<"down_ratio: "<<down_ratio<<endl;
            file<<"aggravate: "<<aggravate<<endl;
            read_data();
            output(file);
        }
    // for(int t=0;t<1;t++)
    // {
    //     find_best_coe(coe1,file);
    //     cout<<"the learning has gone a loop: the following is the result: "<<endl;
    //     cout<<"mean: "<<out<<endl;
    //     cout<<"coe1: "<<coe1<<endl<<endl;
    // }

    
    return 0;
}