percolation without error

//code which creates network for 25 times and does percolation over it 1000times and saves the file

#include<iostream>
#include<vector>
#include<string>
#include<cstdlib>
#include<iomanip>
#include<ctime>
#include<cmath>
#include<random>
#include<fstream>
#include<algorithm>
#include<omp.h>

//#include "openacc_curand.h"
using namespace std;

int EMPTY = 0;
int connectionNumber = 0; // it indexes the connection between different nodes of the network

// this function does the union find work for the percolation

int findroot(int ptr[],int j)
{
    if(ptr[j]<0) return j;

    return ptr[j]=findroot(ptr,ptr[j]);
}

//this function does the same thing but without recursion
/*
int findroot(int ptr[],int i)
{
    //cao = 1;
    int r,s;
    r = s = i;
    while (ptr[r]>=0)
    {
        ptr[s] = ptr[r];
        s = r;
        r = ptr[r];
    }
    return r;
}
*/
int main()
{
    int seed,vertices,m,run,filenum;

    /*

    I am just going to set the initial value for your need

    cout<<"enter seed size: ";
    cin>>seed;
    cout<<endl;
    cout<<"enter vertice number: ";
    cin>>vertices;
    cout<<endl;
    cout<<"order number: ";
    cin>>m;
    cout<<endl;
    cout<<"enter ensemble number: ";
    cin>>run;
    cout<<endl;
    cout<<"enter filenumber: ";
    cin>>filenum;
    cout<<endl;
*/
    seed = 6;
    vertices = 4000;
    m = 5;
    run = 25000;
    filenum = 6;

    //this sets up the connection and initializes the array;
    int con = 0;

    for(int i=1;i<seed;i++)
    {
        con = con + i;
    }

    con = con + (vertices-seed)*m;

    int* netmap1 = new int[con+1]; //node 1 that is connected to a certain connectionNumber
    int* netmap2 = new int[con+1]; //node 2 that is connected to a certain connectionNumber

    for(int i=1;i<=con;i++)
    {
        netmap1[i] = 0;
        netmap2[i] = 0;
    }

    connectionNumber = con;

    srand(time(NULL));
    int A,B,C;
    A = vertices;
    B = run;
    C = filenum;

    //saved filename

    string filename1;
    filename1 = "maxclus_";
    string filename2;
    filename2 = "delmx_";
    string filename3;
    filename3 = "entropy_";

    filename1 = filename1+to_string(A)+"node_"+to_string(m)+"M_"+to_string(B)+"ens"+to_string(C)+".dat";
    filename2 = filename2+to_string(vertices)+"node_"+to_string(m)+"M_"+to_string(run)+"ens"+to_string(filenum)+".dat";
    filename3 = filename3+to_string(vertices)+"node_"+to_string(m)+"M_"+to_string(run)+"ens"+to_string(filenum)+".dat";

    ofstream fout1,fout2,fout3;



    double* maxclus = NULL;
    maxclus = new double[connectionNumber+1];

    double* delmx = NULL;
    delmx = new double[connectionNumber+1];

    double* entropycalc = NULL;
    entropycalc = new double[connectionNumber+1];

    for(int i=0;i<=connectionNumber;++i)
    {
        maxclus[i]=0;
        delmx[i]=0;
        entropycalc[i]=0;
    }

    int* random = new int[connectionNumber+1];

            for(int i=0;i<=connectionNumber;i++)
            {
                random[i] = i;
            }

            random_shuffle(&random[1],&random[connectionNumber]);
            //this is the pointer that needs to be made private for all the parallel loops
            int* ptr = new int[vertices+1];
    for(int i=0;i<vertices+1;i++) ptr[i]=0;

    //the main program starts here



    for(int i=1;i<=25;i++)
    {
        cout<<"run : "<<i<<endl;
        vector<int> network;
        connectionNumber = 0;

        //seeds are connected to the network
        for(int i=1;i<=seed;i++)
        {
            for(int j=1;j<=seed;j++)
            {
                if(j>i)
                {

                    connectionNumber=connectionNumber + 1;
                    netmap1[connectionNumber]=i; //connections are addressed
                    netmap2[connectionNumber]=j;
                    network.push_back(i); // the vector is updated for making connection
                    network.push_back(j);
                }
            }
        }

        int concheck = 0;
        int ab[m]; //this array checks if a node is chosen twice
        int countm = 0;

        for(int i = seed+1;i<=vertices; ++i)
        {
            countm = 0;
            for(int k=0;k<m;k++) ab[k] = 0;
            //if(i%50000==0)cout<<"connection given to : "<<i<<endl;
            for(int j = 1; ;j++)
            {
                concheck = 0;
                int N1=network.size() ;
                int M1=0;
                int u = M1 + rand()/(RAND_MAX/(N1-M1+1) + 1);

                for(int n=0;n<m;n++)
                {
                    if(ab[n] == network[u]) concheck = 1;
                }

                //if condition is fulfilled the connection are given to the nodes
                //the data is saved in the arrays of the connection
                if(concheck == 0 && network[u]!=i)
                {
                    ab[countm] = network[u];
                    countm=countm+1;



                    connectionNumber=connectionNumber+1;
                    netmap1[connectionNumber] = i;
                    netmap2[connectionNumber] = network[u];

                    network.push_back(i);
                    network.push_back(network[u]);
                }

                if(countm==m) break;

            }
        }


        double rand_seed = time(NULL);



        connectionNumber = con;
        //check data
        for(int i=0;i<=connectionNumber;i++)
            {

                cout<<"maxclus: "<<maxclus[i]<<'\t'<<"delmx: "<<delmx[i]<<'\t'<<"entropy: "<<entropycalc[i]<<'\t'<<"net1:"<<netmap1[i]<<'\t'<<"net2:"<<netmap2[i]<<endl;
                cout<<"con: "<<connectionNumber<<'\t'<<"vertices:"<<vertices<<endl;
            }



            //this is where the problem lies
        //basically i want to make all the rx loops parallel in such a way that every parallel loop will have their own copy of ptr[ ] and random[ ] which they can modify themselves
        // this whole part does the 'explosive percolation' and saves the data in maxclus, delmx, entropycalc array of different runs

        //#pragma omp parallel for shared(maxclus,delmx,entropycalc,netmap1,netmap2) schedule(guided)
        for(int rx=1;rx<=10;++rx)
        {
            //EMPTY = 0;
            //connectionNumber = con;
            //int* random = new int[connectionNumber+1];
            //int* ptr = new int[vertices+1];
            //for(int i=0;i<=connectionNumber;++i)
            //{
            //  random[i] = i;
            //}

            random_shuffle(&random[1],&random[connectionNumber]);


            //this is the pointer that needs to be made private for all the parallel loops

            for(int i=0;i<=vertices;++i) ptr[i] = 0;

            //cout<<"run commencing : "<<rx<<endl;


            int index=0,big=0,bigtemp=0,jump=0,en1=0,en2=0;
            int node1=0,node2=0,node3=0,node4=0,nodeA=0,nodeB=0;
            int clus1=0,clus2=0,clus3=0,clus4=0;
            double entropy = log(vertices);


            for(int i=1;i<=connectionNumber;++i)
            {

                if(i!=connectionNumber)
                {

                    node1 = netmap1[random[i]];
                    node2 = netmap2[random[i]];
                    node3 = netmap1[random[i+1]];
                    node4 = netmap2[random[i+1]];

                    if(ptr[node1]==EMPTY) clus1 = 1;
                    else

                    {
                        int x = findroot(ptr,node1);
                        clus1 = -ptr[x];
                    }

                    if(ptr[node2]==EMPTY) clus2 = 1;
                    else

                    {
                        int b = findroot(ptr,node2);
                        clus2 = -ptr[b];
                    }

                    if(ptr[node3]==EMPTY) clus3 = 1;
                    else

                    {
                        int c = findroot(ptr,node3);
                        clus3 = -ptr[c];
                    }

                    if(ptr[node4]==EMPTY) clus4 = 1;
                    else

                    {
                        int d = findroot(ptr,node4);
                        clus4 = -ptr[d];
                    }

                    if(clus1*clus2<clus3*clus4)
                    {
                        nodeA = node1;
                        nodeB = node2;


                        int N1=connectionNumber;
                        int M1=i+1;

                        int u = M1 + (rand()/RAND_MAX)*(N1-M1);
                        //int u = M1 + rand()/(RAND_MAX/(N1-M1+1) + 1);

                        int temp = random[u];
                        random[u] = random[i+1];
                        random[i+1] = temp;
                    }
                    else
                    {
                        nodeA = node3;
                        nodeB = node4;

                        int temp = random[i+1];
                        random[i+1] = random[i];
                        random[i] = temp;
                        int N1=connectionNumber;
                        int M1=i+1;
                        int u = M1 + (rand()/RAND_MAX)*(N1-M1);
                        //int u = M1 + rand()/(RAND_MAX/(N1-M1+1) + 1);


                        temp = random[u];
                        random[u] = random[i+1];
                        random[i+1] = temp;
                    }
                }
                else
                {
                    nodeA = netmap1[random[i]];
                    nodeB = netmap2[random[i]];
                }




                if(ptr[nodeA]==EMPTY && ptr[nodeB]==EMPTY)
                {

                en1=1;
                en2=1;
                ptr[nodeA] = -2;
                ptr[nodeB] = nodeA;
                index = nodeA;
                entropy = (double)(entropy-(-2.0/vertices*log(1.0/vertices))+(-2.0/vertices*log(2.0/vertices)));
                if(entropy<0) entropy = 0;
                }
                else if(ptr[nodeA]==EMPTY && ptr[nodeB]!=EMPTY)
                {

                en1=1;
                int e;

               e = findroot(ptr,nodeB);
                en2 = -(ptr[e]);
                ptr[nodeA] = e;
                ptr[e] += -1;
                index = e;
                entropy = entropy-(-(double)1.0/vertices*log(1.0/(double)vertices))-(-(double)en2/vertices*log((double)en2/vertices))+(-( double)(-ptr[index])/vertices*log((-ptr[index])/(double)vertices));
                if(entropy<0) entropy = 0;


                }
                else if(ptr[nodeA]!=EMPTY && ptr[nodeB]==EMPTY)
                {

                en2 = 1;
                int f;
                f = findroot(ptr,nodeA);
                en1 = -(ptr[f]);
                ptr[nodeB] = f;
                ptr[f] += -1;
                index = f;
                entropy = entropy-(-(double)1.0/(double)vertices*log(1.0/(double)vertices))-(-(double)en1/(double)vertices*log((double)en1/vertices))+(-(double)(-ptr[index])/vertices*log((-ptr[index])/(double)vertices));
                if(entropy<0) entropy = 0;
                }
                else if(ptr[nodeA]!=EMPTY && ptr[nodeB]!=EMPTY)
                {

                    int g,h;

                    g = findroot(ptr,nodeA);
                    en1 = -(ptr[g]);

                    h = findroot(ptr,nodeB);
                    en2 = -(ptr[h]);
                    if(g!=h)
                    {
                        if(ptr[g]<ptr[h])
                        {

                            ptr[g] += ptr[h];
                            ptr[h] = g;
                            index = g;
                        }
                        else
                        {

                            ptr[h] += ptr[g];
                            ptr[g] = h;
                            index = h;
                        }
                        entropy = entropy-(-(double)en1/(double)vertices*log((double)en1/(double)vertices))-(-(double)en2/vertices*log((double)en2/(double)vertices))+(-(double)(-ptr[index])/vertices*log((double)(-ptr[index])/(double)vertices));
                    if(entropy<0) entropy = 0;
                    }
                    else
                    {

                        jump=big-bigtemp;
                        //#pragma omp atomic update
                        maxclus[i] += big;
                        //#pragma omp atomic update
                        delmx[i] += jump;
                        if(entropy<0) entropy = 0;
                        //#pragma omp atomic update
                        entropycalc[i] += entropy;
                        bigtemp = big;

                        continue;
                    }
                }

                if(-ptr[index]>big) big = -ptr[index];

                jump = big - bigtemp;
                //#pragma omp atomic update
                maxclus[i] += big;
                //#pragma omp atomic update
                delmx[i] += jump;
                //#pragma omp atomic update
                entropycalc[i] += entropy;
                bigtemp = big;


            }
            //cout<<"a"<<endl;
            //delete [] ptr;
            //ptr = NULL;
            //delete[] random;
            //random = NULL;
            //#pragma omp barrier
        }

        //delete [] ptr;
            //ptr = NULL;
        //delete[] random;
            //for(int i=0;i<connectionNumber;i++)
            //{

            //  cout<<"maxclus: "<<maxclus[i]<<'\t'<<"delmx: "<<delmx[i]<<'\t'<<"entropy: "<<entropycalc[i]<<'\t'<<endl;

            //}
    }

    fout1.open(filename1.c_str());
    fout2.open(filename2.c_str());
    fout3.open(filename3.c_str());

    connectionNumber = con;

    for(int i=1;i<=connectionNumber;i++)
    {
        fout1<<(double)i/vertices<<'\t'<<(double)maxclus[i]/vertices/run<<endl;
        fout2<<(double)i/vertices<<'\t'<<(double)delmx[i]/run<<endl;
        fout3<<(double)i/vertices<<'\t'<<(double)entropycalc[i]/run<<endl;
    }

    fout1.close();
    fout2.close();
    fout3.close();



    delete[] maxclus;
    maxclus = NULL;
    delete[] delmx;
    delmx = NULL;
    delete[] entropycalc;
    entropycalc = NULL;

    delete [] netmap1;
    netmap1 = NULL;
    delete [] netmap2;
    netmap2 = NULL;
    delete [] ptr;
    ptr = NULL;
    delete[] random;
    random = NULL;


    return 0;
}