Tricky C++ Code

[ASQ19]

I've the following C++ code that calculate a variable t, the code is dealing with a 2D problem but it was modified to have the output printed in 1D array.

I want to add a line in the code to calculate the #of grid at a specific time or let's the # of grids with respect to time profile. For example: at t=1 or t=<1 we've three grid blocks, at t=<2 we've 5 grid block ,...etc

There are three source codes associated with it, and I'm only attaching part of the main code, the part that needs to be modified, and this part is divided int two parts


//--------Part-1----------

//modification starts from her

int main()
{
clock_t startClock,finishClock,currentClock;
vector<int>::iterator pos;
double timeCount,temp_time;
startClock = clock();

//-----your sort function goes here here--------
// Declare size of two dimensional array and initialize.
int nx,ny,nz,nt,i,j,k,initial_zero,index,inf;
cout << endl << "enter nx == " ;
cin>> nx;

cout << endl << "enter ny == " ;
cin>> ny;

cout << endl << "enter nz == " ;
cin>> nz;

nt=nx*ny*nz;

startClock = clock();
//ofstream timefile;
cout<<"t1= ";
cout<<startClock/1000<<endl;

int s,x,y,z,XYZ_index,w,IJK_index,itt,npixels,q;
double Tt;
double *F, *SourcePoints;
int ne[18]={-1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1};
/* Neighbour list */
int neg_free;
int neg_pos;
double *neg_listv;
double *neg_listx;
double *neg_listy;
double *neg_listz;
double *neg_listo;
/* Matrix containing the Frozen Pixels" */
bool *Frozen;
/* Size of input image */
int *dims_c;
int dims[3];
/* Size of SourcePoints array */
int num_source;

cout << endl << "enter num_source == " ;
cin>> num_source;


// assign the dimensions
dims[0]=nx;
dims[1]=ny;
dims[2]=nz;

npixels=dims[0]*dims[1]*dims[2];


//double* dx = NULL;
dx = new double[nt];

//double* dy = NULL;
dy = new double[nt];

//double* dz = NULL;
dz = new double[nt];

for (int ii=0;ii<npixels;ii++){
dx[ii]=1;
dy[ii]=1;
dz[ii]=1;
}


/* The output distance image */
double *T;
int *listprop;
double **listval;
double *Final_T;
double *perm;
double *permx;
double *permx_orig;
double *time_orig;
double *dtc;
double *porosity;
double *porosity_orig;
double *source;

permx_orig=new double[npixels];
time_orig=new double[npixels];
dtc=new double[npixels];
porosity_orig=new double[npixels];
Final_T=new double[npixels];
perm=new double[npixels];
permx=new double[npixels];
porosity=new double[npixels];

source=new double[num_source*num_source*3];

int kk;
kk=-1;

for(int i=0;i<num_source;i++){
for(int j=0;j<num_source;j++){
k=0;
kk=kk+1;

source[3*kk]=nx/2+1;
source[3*kk+1]=ny/2+1;
source[3*kk+2]=nz/2+1;
cout<<source[3*(kk)]<<"-"<<source[3*kk+1]<<"-"<<source[3*kk+2]<<endl;

}
}

for (int ii=0;ii<npixels;ii++){
permx[ii]=1;
permx_orig[ii]=1;
porosity[ii]=1;
porosity_orig[ii]=1;
}
int ii;





T= Final_T;
F=permx;
SourcePoints=source;

/* Pixels which are processed and have a final distance are frozen */
Frozen = new bool[npixels];
for(q=0;q<npixels;q++){Frozen[q]=0; T[q]=-1;}


/*Free memory to store neighbours of the (segmented) region */
neg_free = 1000000;
neg_pos=0;

neg_listx = new double [neg_free];
neg_listy = new double [neg_free];
neg_listz = new double [neg_free];

/* List parameters array */
listprop=new int [3];
/* Make jagged list to store a maximum of 2^64 values */
//listval= (double **)malloc( 64* sizeof(double *) );/// num of rows=64?
listval= new double* [64];
/* Initialize parameter list */
initialize_list(listval, listprop);
neg_listv=listval[listprop[1]-1];


for (s=0; s<num_source*num_source; s++) {
/*starting point */
x= (int)SourcePoints[0+s*3]-1;
y= (int)SourcePoints[1+s*3]-1;
z= (int)SourcePoints[2+s*3]-1;
XYZ_index=mindex3(x, y, z, dims[0], dims[1]);

Frozen[XYZ_index]=1;
T[XYZ_index]=0;
}

for (s=0; s<num_source*num_source; s++) {
/*starting point */
x= (int)SourcePoints[0+s*3]-1;
y= (int)SourcePoints[1+s*3]-1;

z= (int)SourcePoints[2+s*3]-1;


XYZ_index=mindex3(x, y, z, dims[0], dims[1]);
for (w=0; w<6; w++) {
/*Location of neighbour */
i=x+ne[w];
j=y+ne[w+6];
k=z+ne[w+12];

IJK_index=mindex3(i, j, k, dims[0], dims[1]);

/*Check if current neighbour is not yet frozen and inside the */

/*picture */
if(isntfrozen3d(i, j, k, dims, Frozen)) {
Tt=(1/(max(F[IJK_index],eps)));
/*Update distance in neigbour list or add to neigbour list */
if(T[IJK_index]>0) {
if(neg_listv[(int)T[IJK_index]]>Tt) {
listupdate(listval, listprop, (int)T[IJK_index], Tt);
}
}
else {
/*If running out of memory at a new block */
//if(neg_pos>=neg_free) {
// neg_free+=100000;
// neg_listx = (double *)realloc(neg_listx, neg_free*sizeof(double) );
// neg_listy = (double *)realloc(neg_listy, neg_free*sizeof(double) );
// neg_listz = (double *)realloc(neg_listz, neg_free*sizeof(double) );
//}
list_add(listval, listprop, Tt);
neg_listv=listval[listprop[1]-1];
neg_listx[neg_pos]=i;
neg_listy[neg_pos]=j;
neg_listz[neg_pos]=k;
T[IJK_index]=neg_pos;
neg_pos++;
}
}
}
}

/*Loop through all pixels of the image */
for (itt=0; itt<(npixels); itt++) /* */ {
/*Get the pixel from narrow list (boundary list) with smallest */
/*distance value and set it to current pixel location */
index=list_minimum(listval, listprop);
neg_listv=listval[listprop[1]-1];
/* Stop if pixel distance is infinite (all pixels are processed) */
if(IsInf(neg_listv[index])) { break; }

[ASQ19]

//----------------Part-2------------------

/*index=minarray(neg_listv, neg_pos); */
x=(int)neg_listx[index]; y=(int)neg_listy[index]; z=(int)neg_listz[index];
XYZ_index=mindex3(x, y, z, dims[0], dims[1]);
Frozen[XYZ_index]=1;
T[XYZ_index]=neg_listv[index];
//cout<< T[XYZ_index]<<endl;
//show_list(listval, listprop);
/*Remove min value by replacing it with the last value in the array */
//1- remove the min value and go up in tree by comparing it with its adjacent bothers and parent
//2- remove the last value in the list and again go up in tree
//3- replace the minimum by the last value and again go up in the tree
//- if we have extra tree delete them
list_remove_replace(listval, listprop, index) ;
//show_list(listval, listprop);
neg_listv=listval[listprop[1]-1];
if(index<(neg_pos-1)) {
neg_listx[index]=neg_listx[neg_pos-1];
neg_listy[index]=neg_listy[neg_pos-1];
neg_listz[index]=neg_listz[neg_pos-1];
T[(int)mindex3((int)neg_listx[index], (int)neg_listy[index], (int)neg_listz[index], dims[0], dims[1])]=index;
}
neg_pos =neg_pos-1;

/*Loop through all 6 neighbours of current pixel */
for (w=0;w<6;w++) {
/*Location of neighbour */
i=x+ne[w]; j=y+ne[w+6]; k=z+ne[w+12];
IJK_index=mindex3(i, j, k, dims[0], dims[1]);

/*Check if current neighbour is not yet frozen and inside the */
/*picture */
if(isntfrozen3d(i, j, k, dims, Frozen)) {
//Tt=CalculateDistance(T, F[IJK_index], dims, i, j, k, usesecond, usecross, Frozen);
Tt=CalculateDistance(T, F[IJK_index], dims, i, j, k, Frozen);

/*Update distance in neigbour list or add to neigbour list */
IJK_index=mindex3(i, j, k, dims[0], dims[1]);
if((T[IJK_index]>-1)&&T[IJK_index]<=listprop[0]) {
if(neg_listv[(int)T[IJK_index]]>Tt) {
listupdate(listval, listprop, (int)T[IJK_index], Tt);
}
}
else {
list_add(listval, listprop, Tt);
//show_list(listval, listprop);
neg_listv=listval[listprop[1]-1];
neg_listx[neg_pos]=i; neg_listy[neg_pos]=j; neg_listz[neg_pos]=k;

T[IJK_index]=neg_pos;
neg_pos++;
}
}
}
}
/* Free memory */
/* Destroy parameter list */
destroy_list(listval, listprop);
free(neg_listx);
free(neg_listy);
free(neg_listz);
free(Frozen);
finishClock=clock();
timeCount = finishClock - startClock ;
cout<<"run time = ";
cout<<timeCount/1000/60<<endl;

//==================================================================================
ofstream myfile1;
myfile1.open ("fine-time.txt");
//myfile << "Writing this to a file.\n";

double temp=0;
myfile1 << setw(16) ;
for (int k=1; k<nz+1;k++){
for (int j=1; j<ny+1;j++){
for (int i=1; i<nx+1;i++){
index= (k-1) * nx * ny +(j-1) * (nx) + (i-1) ;
/*if(Final_T[index]==inf){
Final_T[index]=1.;
}*/
myfile1 <<Final_T[index]<< setw(16) ;
// save fine scale time
time_orig[index]= Final_T[index];
//myfile1 <<F[index]<< setw(16) ;
temp=temp+Final_T[index];

}

myfile1 << endl ;
}
}
myfile1.close();
cout<<"Summation-time = ";
cout<<temp/1000;





//=====================================================================================
}

[JLBorges]

> I want to add a line in the code to calculate the #of grid at a specific time
> or let's the # of grids with respect to time profile.
> For example: at t=1 or t=<1 we've three grid blocks, at t=<2 we've 5 grid block ,...etc

You have an array; and the contents of the array are modified at certain points in time.
Now, what does "# of grid" mean? What is a "grid block"?

[ASQ19]

Let's say we've some results printed in a 5x5 or 200x200 matrix,...etc . I need to sort the numeric results printed in the matrix by dividing the results, which is the variable t, into intervals.

Check the following example:

I run the same code ,that I posted part of it earlier, and I choose a 5x5 system

Here is the results of "t", after imported to matlab:

5x5 1 2 3 4 5
1 3.2524 2.5453 2 2.5453 3.2524
2 2.5453 1.7071 1 1.7071 2.5453
3 2.0000 1.0000 0 1.0000 2.0000
4 2.5453 1.7071 1 1.7071 2.5453
5 3.2524 2.5453 2 2.5453 3.2524

So, since it's a 5x5 I did the sorting of the results manually as follow:

I looked at each number of t and count how many times ( or how many grids) are repeated and record the data as shown below:

t matrix element (i,j, i-1,j,...ect) = # of grid
0 0
1 4
1.7 8
2 12
2.54 20
3.25 25

[JLBorges]

> So, since it's a 5x5 I did the sorting of the results manually
> I looked at each number of t and count how many times are repeated and record the data

So do the same thing in code now; sort the sequence and make a single pass through it, counting the number of repeats as you move through it.

#include <vector>
#include <algorithm>
#include <utility>
#include <iostream>

inline bool almost_equals( double a, double b )
{
    static constexpr double EPSILON = 1.0e-6 ;
    return std::fabs(a-b) < EPSILON ;
}

std::vector< std::pair<double,std::size_t> > counts( std::vector<double> t )
{
    std::vector< std::pair<double,std::size_t> > result ;

    if( !t.empty() )
    {
        std::sort( t.begin(), t.end() ) ;

        double previous = t.front() ;

        for( std::size_t i = 1 ; i < t.size() ; ++i )
        {
            if( !almost_equals( t[i], previous ) ) result.emplace_back( previous, i ) ;
            previous = t[i] ;
        }

        result.emplace_back( t.back(), t.size() ) ;
    }

    return result ;
}

int main()
{
    std::vector<double> t = { 1.1, 2.2, 3.3, 2.2, 4.4,
                                1.1, 5.5, 3.3, 4.4, 3.3,
                                8.8, 9.9, 7.7, 3.3, 5.5,
                                4.4, 2.2, 7.7, 5.5, 1.1,
                                6.6, 6.6, 7.7, 8.8, 2.2 } ;
    auto c = counts(t) ;
    for( const auto& p : c ) std::cout << p.first << ' ' << p.second << '\n' ;
}

http://liveworkspace.org/code/1I6udJ$0

[ASQ19]

Thank you JLBorges,

I did not understand the last part of the code:

std::vector<double> t = { 1.1, 2.2, 3.3, 2.2, 4.4,
1.1, 5.5, 3.3, 4.4, 3.3,
8.8, 9.9, 7.7, 3.3, 5.5,
4.4, 2.2, 7.7, 5.5, 1.1,
6.6, 6.6, 7.7, 8.8, 2.2 } ;
Is it for another example ??

Is this good for 2000x2000 system?

If you inspect my code, there is a line where the 2D array is converted to 1D and at the very end it's converted to 2D when the results are printed

[JLBorges]

> Is it for another example ??

Yes, it is just an example.


> Is this good for 2000x2000 system?

The code is independent of the size of the sequence.
Subject to the proviso that a larger sequence would require more memory.

For machines these days, 2000x2000 is not a big size.

[ASQ19]

I did figure out how to do it, it's as follow:

ofstream myfile1;
ofstream grid;
myfile1.open ("fine-time.txt");
grid.open ("grid.txt");
//myfile << "Writing this to a file.\n";

double temp=0;
myfile1 << setw(16) ;
for (int k=1; k<nz+1;k++){
for (int j=1; j<ny+1;j++){
for (int i=1; i<nx+1;i++){
index= (k-1) * nx * ny +(j-1) * (nx) + (i-1) ;
/*if(Final_T[index]==inf){
Final_T[index]=1.;
}*/
myfile1 <<Final_T[index]<< setw(16) ;
// save fine scale time
time_orig[index]= Final_T[index];
//myfile1 <<F[index]<< setw(16) ;
temp=temp+Final_T[index];

}

myfile1 << endl ;
}
}


//variables for tempery values and to count # of grids
double temp1,val=0.0;
int count=0;

//sorting Final_T array using buble sort algorithm
for(i=0;i<npixels;i++){
for(j=0;j<npixels-1;j++){

if(Final_T[j]>Final_T[j+1]){

temp1=Final_T[j+1];
Final_T[j+1]=Final_T[j];
Final_T[j]=temp1;
}


}
}

//write in to gird.txt
for(i=0;i<npixels;i++){

if(val<Final_T[i]){
grid<<Final_T[i-1]<<" --> "<<count-1<<endl;
val=Final_T[i];
}
count++;

}
grid<<Final_T[i-1]<<" --> "<<count<<endl;





myfile1.close();
grid.close();
cout<<"Summation-time = ";
cout<<temp/1000;

----------------


However, the grid.txt file will not show any result for systems beyond 200 by 200

In order for somebody to help he needs to run the file first, let me know if you need the source files