Profile_Generator.cpp

/*
 Probabilistic time-dependant speed profile generator     

Vytvořený software implementuje automatické generování pravděpodobnostních rychlostních profilů pro algoritmus Probabilistic Time-Dependent Routing (PTDR).
Účelem PTDR je odhad pravděpodobnostního rozdělení dojezdového času na dané trase. Toto rozdělení je určeno pomocí tzv. pravděpodobnostních rychlostních profilů,
které popisují stochastické vlastnosti hran grafové reprezentace silniční sítě. Cílem tohoto software je generovat tyto profily na základě historických dat o rychlosti
dopravy na hranách silniční sítě. Vzhledem k značnému množství hran v dopravní síti je nutné zpracovat relativně velké množství dat a výsledné profily uložit tak, aby jejich
následné zpracování bylo co nejefektivnější. Přínos výsledku spočívá v inovativním využití datového formátu HDF5 v kombinaci s metodou Aliasů pro vzorkování diskrétních
pravděpodobnostních rozdělení. Tato kombinace umožňuje nejen velmi rychlý výpočet a zápis výsledných profilů ale také zefektivňuje následnou práci s výslednými
profily v rámci algoritmu PTDR díky využitelnosti HDF5 pro High Perforamnce Computing (HPC). Využití metody Aliasů také zároveň eliminuje zaokrouhlovací chyby
často spojené s operacemi využívajícími diskrétní rozdělení pravděpodobnosti. Software je implementován v jazyce C/C++, využívá knihovny pro čtení a zápis HDF5 a byl kompilován kompilátorem Cygwin. 

Vstupy:

Argv1: Cesta ke složce s naměřenými agregovanými rychlostmi pro jednotlivé segmenty. Tato složka musí obsahovat dílčí složky s jednotlivými měřeními rychlosti pro každý interval ve fromátu csv a sloupci SegmentID, Speed (km/h)
Argv2: Počet složek ve složce Argv1 a také počet měření ze kterých se profily počítají ( např v případě týdenních profilů 5 znamemná že se počítají z 5 ti týdnů historických dat)
Argv3: Počet intervalů ve vstupních datech
Argv4: Cesta k výslednému souboru HDF5
Argv5: Délka časového okna profilu v sekundách
Argv6: Cesta k souboru popisujícího silniční síť ve fromátu SegmentID,Length(m),TravelTime(s)

Výstupem je pak HDF5 soubor obsahující profily ze segmentů s naměřenými daty pro počet intervalů ze vstupních dat. Ten obsahuje jak pravděpodobnostní rychlostní profily v Alias formátu tak prosté rychlostní profily.

V případě dotazů kontakutjte na rap027@vsb.cz
 
 
 
 */


#include "hdf5.h"
#include <cstdlib>
#include <dirent.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#define MAX_NAME 1024


#ifdef NAN
/* NAN is supported */
#endif

using namespace std;

struct alias {
    double prob[4];
    double alias[4];
};

struct road_network {
    int gid;
    double length;
    double tt;
};

bool isvalueinarray(int val, int *arr, int size) {
    int i;
    for (i = 0; i < size; i++) {
        if (arr[i] == val)
            return true;
    }
    return false;
}

double *probab(double arr [], double ffs, int len) {
    double *pr = (double *) malloc(4 * sizeof (double));
    for (int i = 0; i < 4; i++)
        pr[i] = 0;
    for (int i = 0; i < len; i++) {
        // printf("ar %f\n",arr[i]);

        if (arr[i] > ffs * 0.95)
            pr[0] += 1;
        else if (arr[i] > ffs * 0.75)
            pr[1] += 1;
        else if (arr[i] > ffs * 0.45)
            pr[2] += 1;
        else
            pr[3] += 1;

    }
    for (int i = 0; i < 4; i++) {
        pr[i] = pr[i] / (double) len;
    }
    return pr;
}


struct alias compalias(double *probb) {
    struct alias a;
    int n = 5;
    double prob[5];
    double probbm[5];
    double alias[5];
    int small[5];
    int large [5];
    for (int i = 0; i < 4; i++) {
        prob[i + 1] = NAN;
        alias[i + 1] = NAN;
        probbm[i + 1] = probb[i] * ((n - 1) / 1);
        small[i + 1] = NAN;
        large[i + 1] = NAN;
    }
    int ns = 0;
    int nl = 0;
    for (int j = 1; j < n; j++) {
        if (probbm[j] > 1) {
            nl = nl + 1;
            large[nl] = j;
        } else {
            ns = ns + 1;
            small[ns] = j;
        }
    }
    while (ns != 0 && nl != 0) {
        int j = small[ns];
        ns = ns - 1;
        int k = large[nl];
        nl = nl - 1;
        prob[j] = probbm[j];
        alias[j] = k;
        probb[k] = probbm[k] + probbm[j] - 1;
        if (probbm[k] > 1) {
            nl = nl + 1;
            large[nl] = k;
        } else {
            ns = ns + 1;
            small[ns] = k;
        }
    }
    while (ns != 0) {
        prob[small[ns]] = 1;
        ns = ns - 1;
    }
    while (nl != 0) {
        prob[large[nl]] = 1;
        nl = nl - 1;
    }
    for (int i = 1; i < 5; i++) {
        alias[i] = alias[i] - 1;
    }
    for (int i = 1; i < 5; i++) {
        if (isnan(alias[i]) == 1) {
            //printf("fdsf");
            alias[i] = -1;
        }
    }
    double aliasm[4];
    double probm[4];
    for (int i = 1; i < 5; i++) {
        aliasm[i - 1] = alias[i];
        probm[i - 1] = prob[i];
    }
    memcpy(&a.alias, &aliasm, sizeof a.alias);
    memcpy(&a.prob, &probm, sizeof a.prob);
    return a;
}

int main(int argc, char *argv[]) {
    struct dirent *de; // Pointer for directory entry 
    struct dirent *de2; // Pointer for directory entry 
    int num_per = atoi(argv[2]);
    int sim_length = atoi(argv[3]);
    int dim1, dim2, dim3;
    int *prof = (int *) malloc(100000 * sizeof (int));
    int *len = (int *) malloc(100000 * sizeof (int));
    int numprof = 4;
    int timestep = atoi(argv[5]);
    int size = 0;
    double *ffs = (double *) malloc(100000 * sizeof (double));
    dim1 = num_per;
    dim2 = sim_length + 1;
    dim3 = 100000;
    int siz, flag;
    int i, j, k;
    char ver[] = "1.0";
    char mystring1 [100];
    char *str1;
    double f;
    FILE* my_file1 = fopen(argv[6], "r");
    if (my_file1 == NULL) {
        perror(argv[1]);
        exit(1);
    }
    struct road_network *cz = (road_network*) malloc(sizeof (struct road_network)*5000000);
    size_t count1 = 0;
    while (fgets(mystring1, 100, my_file1) != NULL) {
        str1 = strtok(mystring1, ",");
        cz[count1].gid = atoi(str1);
        str1 = strtok(NULL, ",");
        cz[count1].length = atof(str1);
        str1 = strtok(NULL, ",");
        cz[count1].tt = atof(str1);
        str1 = strtok(NULL, ",");
        ++count1;
    }
    hid_t file_id, dataset_id, dataspace_id, group_id, plist, aid2, attr2, ret, atype, aid3, attr3; /* identifiers */
    hsize_t dims[2];
    hsize_t cdims[2];
    herr_t status;
    file_id = H5Fcreate(argv[4], H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    aid2 = H5Screate(H5S_SCALAR);
    attr2 = H5Acreate1(file_id, "IntervalsPerSegment", H5T_NATIVE_INT, aid2, H5P_DEFAULT);
    ret = H5Awrite(attr2, H5T_NATIVE_INT, &sim_length);
    ret = H5Sclose(aid2);
    aid2 = H5Screate(H5S_SCALAR);
    attr2 = H5Acreate1(file_id, "TimeStep", H5T_NATIVE_INT, aid2, H5P_DEFAULT);
    ret = H5Awrite(attr2, H5T_NATIVE_INT, &timestep);
    ret = H5Sclose(aid2);
    aid2 = H5Screate(H5S_SCALAR);
    attr2 = H5Acreate1(file_id, "NumberOfProfiles", H5T_NATIVE_INT, aid2, H5P_DEFAULT);
    ret = H5Awrite(attr2, H5T_NATIVE_INT, &numprof);
    ret = H5Sclose(aid2);
    aid3 = H5Screate(H5S_SCALAR);
    atype = H5Tcopy(H5T_C_S1);
    H5Tset_size(atype, 4);
    attr3 = H5Acreate1(file_id, "Version", atype, aid3, H5P_DEFAULT);
    ret = H5Awrite(attr3, atype, ver);
    ret = H5Sclose(aid3);
    double *** array = (double ***) malloc(dim1 * sizeof (double**));
    for (i = 0; i < dim1; i++) {
        array[i] = (double **) malloc(dim2 * sizeof (double *));
        for (j = 0; j < dim2; j++) {
            array[i][j] = (double *) malloc(dim3 * sizeof (double));
        }
    }
    char *record = (char *) malloc(1000 * sizeof (char));
    char *line = (char *) malloc(1000 * sizeof (char));
    char *file2 = (char *) malloc(1000 * sizeof (char));
    char *file = (char *) malloc(1000 * sizeof (char));
    DIR *dr = opendir(argv[1]);
    if (dr == NULL) 
    {
        printf("Could not open current directory");
        return 0;
    }
    i = 0;
    while ((de = readdir(dr)) != NULL) {
        if (i > 1) {
            strcpy(file, argv[1]);
            printf("%s\n", de->d_name);
            DIR *dr2 = opendir(strcat(strcat(file, "\\"), de->d_name));
            if (dr2 == NULL) // opendir returns NULL if couldn't open directory 
            {
                printf("Could not open current directory");
                return 0;
            }
            j = 0;
            siz = 0;
            while ((de2 = readdir(dr2)) != NULL) {
                if (j > 1) {
                    printf("%s\n", de2->d_name);
                    strcpy(file2, file);
                    strcat(file2, "\\");
                    strcat(file2, de2->d_name);
                    FILE* my_file1 = fopen(file2, "r");
                    if (my_file1 == NULL) {
                        perror(argv[1]);
                        exit(1);
                    }
                    while (fgets(record, 100000, my_file1) != NULL) {
                        line = strtok(record, ";");
                        if (siz == 0) {
                            array[i - 2][0][siz] = atof(line);
                            if (isvalueinarray(atoi(line), prof, size) == false) {
                                prof[size] = atoi(line);
                                for (int y = 0; y < count1; y++) {
                                    if (prof[size] == cz[y].gid) {
                                        ffs[size] = (cz[y].length / (cz[y].tt / 1000))*3.6;
                                        break;
                                    }
                                }
                                size++;
                            }
                            for (int z = 0; z < size + 1; z++) {
                                if (array[i - 2][0][siz] == prof[z]) {
                                    f = ffs[z];
                                }
                            }
                            for (int l = 1; l < dim2; l++) {

                                array[i - 2][l][siz] = f;
                            }
                            line = strtok(NULL, ";");
                            array[i - 2][j - 1][siz] = atof(line);
                            siz++;
                        } else {
                            flag = 0;
                            for (int l = 0; l < siz; l++) {
                                if (array[i - 2][0][l] == atof(line)) {
                                    line = strtok(NULL, ";");
                                    array[i - 2][j - 1][l] = atof(line);
                                    flag = 1;
                                    break;
                                }
                            }
                            if (flag == 0) {
                                array[i - 2][0][siz] = atof(line);
                                if (isvalueinarray(atoi(line), prof, size) == false) {
                                    prof[size] = atoi(line);
                                    for (int y = 0; y < count1; y++) {
                                        if (prof[size] == cz[y].gid) {

                                            ffs[size] = (cz[y].length / (cz[y].tt / 1000))*3.6;
                                            break;
                                        }
                                    }
                                    size++;
                                }
                                for (int z = 0; z < size + 1; z++) {
                                    if (array[i - 2][0][siz] == prof[z]) {
                                        f = ffs[z];                                     
                                    }
                                }
                                for (int l = 1; l < dim2; l++) {
                                    array[i - 2][l][siz] = f;
                                }
                                line = strtok(NULL, ";");
                                printf("siz:%d\n",siz);
                                array[i - 2][j - 1][siz] = atof(line);
                                siz++;
                            }
                        }
                    }
                    fclose(my_file1);
                }
                j++;
            }
            len[i - 2] = siz;
            closedir(dr2);
        }
        i++;
    }
    closedir(dr);
    free(cz);
    free(file);
    free(record);
    free(file2);
    for (i = 0; i < size; i++) {
        printf("iii %d\n",i);
        double *sp_p[dim2 - 1];
        for (int i = 0; i < dim2 - 1; i++)
            sp_p[i] = (double *) malloc(1000 * sizeof (double));
        int pr = 0;
        for (j = 0; j < dim1; j++) {
            for (int k = 0; k < len[j]; k++) {
                if (array[j][0][k] == prof[i]) {
                    for (int l = 1; l < dim2; l++) {
                        sp_p[pr][l-1] = array[j][l][k];
                        //printf("ar %f\n",sp_p[pr][l-1]);
                    }
                    pr++;
                    break;
                }
            }
        }
        double *sppr = (double *) malloc((dim2 - 1) * sizeof (double));
        struct alias a;
        double avg = 0;
        double probb[pr];
        double *probb2 = (double *) malloc(4 * sizeof (double));
        double *sp_pro[8];
        for (int ii = 0; ii < 8; ii++)
            sp_pro[ii] = (double *) malloc((dim2 - 1) * sizeof (double));
        double *al_alias[4];
        for (int ii = 0; ii < 4; ii++)
            al_alias[ii] = (double *) malloc((dim2 - 1) * sizeof (double));
        double *al_pro[4];
        for (int ii = 0; ii < 4; ii++)
            al_pro[ii] = (double *) malloc((dim2 - 1) * sizeof (double));
        for (int w = 0; w < dim2 - 1; w++) {
            avg = 0;
            for (int j = 0; j < pr; j++) {
                avg += sp_p[j][w];               
                probb[j] = sp_p[j][w];
            }
            probb2 = probab(probb, ffs[i], pr);
            for (int q = 0; q < 4; q++) {
                sp_pro[(2 * q) + 1][w] = probb2[q];
                switch (q) {
                    case 0:
                        sp_pro[2 * q][w] = ffs[i];
                        break;
                    case 1:
                        sp_pro[q * 2][w] = 0.75 * ffs[i];
                        break;

                    case 2:
                        sp_pro[2 * q][w] = 0.5 * ffs[i];
                        break;

                    case 3:
                        sp_pro[2 * q][w] = 0.1 * ffs[i];
                        break;
                }
            }
            a = compalias(probb2);
            for (int r = 0; r < 4; r++) {
                al_alias[r][w] = a.alias[r];
                al_pro[r][w] = a.prob[r];
            }
            sppr[w] = avg / pr;
        }
        double *sp_prol = (double *) malloc((8*sim_length) * sizeof (double));
        double *al_aliasl = (double *) malloc((4*sim_length) * sizeof (double));
        double *al_prol = (double *) malloc((4*sim_length) * sizeof (double));
        int counter=0;
        for (int q=0;q<sim_length;q++){
            for (int s=0;s<8;s++){
                sp_prol[counter++]=sp_pro[s][q];
            }
        }
        counter=0;
        for (int q=0;q<sim_length;q++){
            for (int s=0;s<4;s++){
                al_aliasl[counter]=al_alias[s][q];
                al_prol[counter++]=al_pro[s][q];
            }
        }       
        char str[20];
        sprintf(str, "/%d", prof[i]);
        group_id = H5Gcreate1(file_id, str, 0);
        dims[0] = sim_length;
        dims[1] = 4;
        cdims[0] = 4;
        cdims[1] = 4;
        dataspace_id = H5Screate_simple(2, dims, NULL);
        plist = H5Pcreate(H5P_DATASET_CREATE);
        H5Pset_chunk(plist, 2, cdims);
        H5Pset_deflate(plist, 6);
        char str1[80];
        strcpy(str1, str);
        strcat(str1, "/alias");
        dataset_id = H5Dcreate1(file_id, str1, H5T_NATIVE_DOUBLE, dataspace_id, plist);
        dataset_id = H5Dopen(file_id, str1, H5P_DEFAULT);
        status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, al_aliasl);
        status = H5Dclose(dataset_id);
        strcpy(str1, str);
        strcat(str1, "/probability");
        dataset_id = H5Dcreate1(file_id, str1, H5T_NATIVE_DOUBLE, dataspace_id, plist);
        dataset_id = H5Dopen(file_id, str1, H5P_DEFAULT);
        status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, al_prol);
        status = H5Dclose(dataset_id);
        dims[0] = sim_length;
        dims[1] = 8;
        cdims[0] = 4;
        cdims[1] = 4;
        dataspace_id = H5Screate_simple(2, dims, NULL);
        plist = H5Pcreate(H5P_DATASET_CREATE);
        H5Pset_chunk(plist, 2, cdims);
        H5Pset_deflate(plist, 0);
        strcpy(str1, str);
        strcat(str1, "/prob_profile");
        dataset_id = H5Dcreate1(file_id, str1, H5T_NATIVE_DOUBLE, dataspace_id, plist);
        dataset_id = H5Dopen(file_id, str1, H5P_DEFAULT);
        status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, sp_prol);      
        status = H5Dclose(dataset_id);
        dims[0] = 1;
        dims[1] = sim_length;
        cdims[0] = 1;
        cdims[1] = 4;
        dataspace_id = H5Screate_simple(2, dims, NULL);
        plist = H5Pcreate(H5P_DATASET_CREATE);
        H5Pset_chunk(plist, 2, cdims);
        H5Pset_deflate(plist, 6);
        strcpy(str1, str);
        strcat(str1, "/speed_profile");
        dataset_id = H5Dcreate1(file_id, str1, H5T_NATIVE_DOUBLE, dataspace_id, plist);
        dataset_id = H5Dopen(file_id, str1, H5P_DEFAULT);
        status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, sppr);
        status = H5Dclose(dataset_id);
    }
    status = H5Fclose(file_id);
    return 0;
}