C++ help

Форум » C/C++ для начинающих (C/C++ for beginners) » C++ help » Ответить

C++ help

asq19: I've attached a link to my problem that I need to code http://shot.qip.ru/00bNfp-4mIrA3itZ/

Ответов - 9

Сыроежка: It would be better if you will describe your problem here without any attachment. If the problem is big enough then try to split in into smaller problems.

asq19: The code is shown below: It's working perfect, What I'm trying to do is code what I attached in the excel sheet into this one below: // Upscaling-FMM.cpp : Defines the entry point for the console application. // 28jan-2012.cpp : Defines the entry point for the console application. // 27jan2012.cpp : Defines the entry point for the console application. // // 11-dec-noon.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include <iostream> #include <vector> #include <math.h> #include <time.h> #include <ctime> #include <fstream> #include <iomanip> //#include <vector> //#include <algorithm> #include <stdio.h> #include <stdlib.h> using namespace std; //========================================================================================================= double *dxc; double *dyc; double *dzc; double *dx; double *dy; double *dz; //========================================================================================================= //#define eps 2.2204460492503131e-16 #define eps 3.e-16 #define doublemax 1e50 #define INF 2e50 //#define INF 111 #define listINF 2.345e50 //#define listINF 111 #ifndef min #define min(a,b) ((a) < (b) ? (a): (b)) #endif #ifndef max #define max(a,b) ((a) > (b) ? (a): (b)) #endif /* Find minimum value of an array and return its index */ int minarray(double *A, int l) { int i; int minind=0; for (i=0; i<l; i++) { if(A<A[minind]){ minind=i; } } return minind; } double pow2(double val) { return val*val; } int iszero(double a){ return a*a<eps; } int isnotzero(double a){ return a*a>eps; } void roots(double* Coeff, double* ans) { double a=Coeff[0]; double b=Coeff[1]; double c=Coeff[2]; double r1, r2; double d; d=max(pow2(b)-4.0*a*c,0); if(isnotzero(a)) { ans[0]= (-b - sqrt(d)) / (2.0*a); ans[1]= (-b + sqrt(d)) / (2.0*a); } else { r1=(-b - sqrt(d)); r2=(-b + sqrt(d)); if(isnotzero(r1)) { if(isnotzero(r2)) { ans[0]= (2.0*c)/r1; ans[1]= (2.0*c)/r2; } else { ans[0]= (2.0*c)/r1; ans[1]= (2.0*c)/r1; } } else if(isnotzero(r2)) { ans[0]= (2.0*c)/r2; ans[1]= (2.0*c)/r2; } else { ans[0]=0; ans[1]=0; } } } int maxarray(double *A, int l) { int i; int maxind=0; for (i=0; i<l; i++) { if(A>A[maxind]){ maxind=i; } } return maxind; } __inline int mindex3(int x, int y, int z, int sizx, int sizy) { return x+y*sizx+z*sizx*sizy; } __inline bool IsFinite(double x) { return (x <= doublemax && x >= -doublemax ); } __inline bool IsInf(double x) { return (x >= doublemax ); } __inline bool IsListInf(double x){ return (x == listINF ); } __inline bool isntfrozen3d(int i, int j, int k, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(k>=0)&&(i<dims[0])&&(j<dims[1])&&(k<dims[2])&&(Frozen[mindex3(i, j, k, dims[0], dims[1])]==0); } __inline bool isfrozen3d(int i, int j, int k, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(k>=0)&&(i<dims[0])&&(j<dims[1])&&(k<dims[2])&&(Frozen[mindex3(i, j, k, dims[0], dims[1])]==1); } int p2x(int x) /* 2^x */ { /* return pow(2,x); */ int y=1; int p2x[16] ={1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}; while(x>15) { x=x-15; y=y*32768; } return y*p2x[x]; } void show_list(double **listval, int *listprop) { int z, k; double mm; for(z=0;z<listprop[1]; z++) { for(k=0;k<p2x(z+1); k++) { mm=k; if((z>0)&&(listval[z-1][(int)floor(mm/2)]>listval[z][k])) { printf("*%6.2f", listval[z][k]); } else { printf(" %6.2f", listval[z][k]); } } printf("\n"); } } void initialize_list(double ** listval, int *listprop) { /* Loop variables */ int i; /* Current Length, Orde and Current Max Length */ listprop[0]=0; listprop[1]=1; listprop[2]=2; /* Make first orde storage of 2 values */ listval[0]=new double[2]; /* Initialize on infinite */ for(i=0;i<2;i++) { listval[0]=listINF; } } void destroy_list(double ** listval, int *listprop) { /* Loop variables */ int i, list_orde; /* Get list orde */ list_orde=listprop[1]; /* Free memory */ for(i=0;i<list_orde;i++) { free(listval); } free(listval); free(listprop); } void list_add(double ** listval, int *listprop, double val) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i, j; /* Temporary list location */ int listp; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* If list is full expand list */ if(list_length==list_lengthmax) { list_lengthmax=list_lengthmax*2; for (i=list_orde; i>0; i--) { listval=listval[i-1]; listval = (double *)realloc(listval, p2x(i+1)*sizeof(double)); for(j=p2x(i); j<p2x(i+1); j++) { listval[j]=listINF; } } listval[0]=new double[2]; listval[0][0]=min(listval[1][0], listval[1][1]); listval[0][1]=listINF; list_orde++; } /* Add a value to the list */ listp=list_length; list_length++; listval[list_orde-1][listp]=val; /* Update the links minimum */ for (i=list_orde-1; i>0; i--) { listp=(int)floor(((double)listp)/2); if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; } } /* Set list parameters */ listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax; } int list_minimum(double ** listval, int *listprop) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i; /* Temporary list location */ int listp; /* Index of Minimum */ int minindex; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Follow the minimum through the binary tree */ listp=0; for(i=0;i<(list_orde-1);i++) { /* <= ????? */ if(listval[listp]<=listval[listp+1]) { listp=listp*2; } else { listp=(listp+1)*2; } } i=list_orde-1; if(listval[listp]<=listval[listp+1]){minindex=listp; } else { minindex=listp+1; } return minindex; } void list_remove(double ** listval, int *listprop, int index) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i; /* Temp index */ int index2; double val; double valmin; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Temporary store current value */ val=listval[list_orde-1][index]; valmin=listINF; /* Replace value by infinite */ listval[list_orde-1][index]=listINF; /* Follow the binary tree to replace value by minimum values from */ /** the other values in the binary tree. */ i=list_orde-1; while(true) { if((index%2)==0) { index2=index+1; } else { index2=index-1; } if(val<listval[index2]) { index=(int)floor(((double)index2)/2.0); if(listval[index2]<valmin) { valmin=listval[index2]; } if(i==0) { break; } listval[i-1][index]=valmin; i--; if(i==0) { break; } } else { break; } } } void list_remove_replace(double ** listval, int *listprop, int index) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i, listp; /* Temporary store value */ double val; int templ; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Remove the value the minimum number and goup in tree by comparision its adjacent number to the parents */ list_remove(listval, listprop, index); /* Replace the removed value by the last in the list. (if it was */ /* not already the last value) */ if(index<(list_length-1)) { /* Temporary store last value in the list */ val=listval[list_orde-1][list_length-1]; /* Remove last value in the list */ list_remove(listval, listprop, list_length-1); /* Add a value to the list */ listp=index; listval[list_orde-1][index]=val; /* Update the links minimum */ for (i=(list_orde-1); i>0; i--) { listp=(int)floor(((double)listp)/2); if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; } } } /* List is now shorter */ list_length--; /* Remove trunk of binary tree / Free memory if list becomes shorter */ if(list_orde>2&&IsListInf(listval[0][1])) { //cout<<"truncate-start"<<endl; //show_list(listval, listprop); list_orde--; list_lengthmax=list_lengthmax/2; /* Remove trunk array */ free(listval[0]); /* Move the other arrays one up */ templ=2; for (i=0; i<list_orde; i++) { listval=listval[i+1]; /* Resize arrays to their new shorter size */ listval = (double *)realloc(listval, templ*sizeof(double)); templ*=2; } } /* Set list parameters */ listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax; //show_list(listval, listprop); //cout<<"truncate-end"<<endl; } void listupdate(double **listval, int *listprop, int index, double val) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* loop variable */ int i, listp; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Add a value to the list */ listp=index; listval[list_orde-1][index]=val; /* Update the links minimum */ for (i=(list_orde-1); i>0; i--) { listp=(int)floor(((double)listp)/2); if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; } } /* Set list parameters */ listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax; } __inline int mindex2(int x, int y, int sizx) { return x+y*sizx; } __inline bool isntfrozen2d(int i, int j, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(i<dims[0])&&(j<dims[1])&&(Frozen[i+j*dims[0]]==0); } __inline bool isfrozen2d(int i, int j, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(i<dims[0])&&(j<dims[1])&&(Frozen[i+j*dims[0]]==1); } //========================================================================================================= double second_derivative(double Txm1, double Txm2, double Txp1, double Txp2) { bool ch1, ch2; double Tm; Tm=INF; ch1=(Txm2<Txm1)&&IsFinite(Txm1); ch2=(Txp2<Txp1)&&IsFinite(Txp1); if(ch1&&ch2) { Tm =min( (4.0*Txm1-Txm2)/3.0 , (4.0*Txp1-Txp2)/3.0);} else if(ch1) { Tm =(4.0*Txm1-Txm2)/3.0; } else if(ch2) { Tm =(4.0*Txp1-Txp2)/3.0; } return Tm; } double CalculateDistance(double *T, double Fijk, int *dims, int i, int j, int k, bool *Frozen) { /* Loop variables */ int q, t; /* Current location */ int in, jn, kn; /* Derivatives */ double Tm[3]={0, 0, 0}; double Coeff[3]; /* local derivatives in distance image */ double Txm1, Txm2, Txp1, Txp2; double Tym1, Tym2, Typ1, Typ2; double Tzm1, Tzm2, Tzp1, Tzp2; double Tt, Tt2,xm1,xp1,ym1,yp1,zm1,zp1,x1,x2,y1,y2,z1,z2; /* Return values root of polynomial */ double ansroot[2]={0, 0}; /* Order derivatives in a certain direction */ int Order[3]={0, 0, 0}; /* Neighbours 4x2 */ int ne[18]={-1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1}; /* Stencil constants */ double G1[3]={1, 1, 1}; x1 = dx[mindex3(i, j, k, dims[0], dims[1])]; y1 = dy[mindex3(i, j, k, dims[0], dims[1])]; z1 = dz[mindex3(i, j, k, dims[0], dims[1])]; /*Get First order derivatives (only use frozen pixel) */ in=i-1; jn=j+0; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Txm1=T[mindex3(in, jn, kn, dims[0], dims[1])];xm1=dx[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Txm1=INF; xm1=0; } in=i+1; jn=j+0; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Txp1=T[mindex3(in, jn, kn, dims[0], dims[1])];xp1=dx[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Txp1=INF; xp1=0; } in=i+0; jn=j-1; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tym1=T[mindex3(in, jn, kn, dims[0], dims[1])];ym1=dy[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tym1=INF; ym1=0;} in=i+0; jn=j+1; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Typ1=T[mindex3(in, jn, kn, dims[0], dims[1])];yp1=dy[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Typ1=INF; yp1=0;} in=i+0; jn=j+0; kn=k-1; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tzm1=T[mindex3(in, jn, kn, dims[0], dims[1])];zm1=dz[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tzm1=INF; zm1=0;} in=i+0; jn=j+0; kn=k+1; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tzp1=T[mindex3(in, jn, kn, dims[0], dims[1])];zp1=dz[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tzp1=INF; zp1=0;} /*Make 1e order derivatives in x and y direction */ Tm[0] = min( Txm1 , Txp1); if(IsFinite(Tm[0])){ Order[0]=1; } else { Order[0]=0; } Tm[1] = min( Tym1 , Typ1); if(IsFinite(Tm[1])){ Order[1]=1; } else { Order[1]=0; } Tm[2] = min( Tzm1 , Tzp1); if(IsFinite(Tm[2])){ Order[2]=1; } else { Order[2]=0; } // X2=XM1 IF CONDITIONAL EXPRESSION IS TRUE AND XP1 IF IT IS FALSE x2 = Tm[0] == Txm1 ? xm1 : xp1; y2 = Tm[1] == Tym1 ? ym1 : yp1; z2 = Tm[2] == Tzm1 ? zm1 : zp1; G1[0] = abs(x2 + x1)/2; G1[1] = abs(y2 + y1)/2; G1[2] = abs(z2 + z1)/2; /*Calculate the distance using x and y direction */ Coeff[0]=0; Coeff[1]=0; Coeff[2]=-1/(max(pow2(Fijk),eps)); for (t=0; t<3; t++) { switch(Order[t]) { case 1: Coeff[0]+=G1[t]; Coeff[1]+=-2.0*Tm[t]*G1[t]; Coeff[2]+=pow2(Tm[t])*G1[t]; break; } } roots(Coeff, ansroot); Tt=max(ansroot[0], ansroot[1]); return Tt; } //starting your modificatiob //============================================================================== 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; } /*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; //===================================================================================== }

asq19: I've the 3 source codes as well if needed

asq19: This code is 2D but it was modified to have the output printed in 1D array

asq19: #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, previous ) ) result.emplace_back( previous, i ) ; previous = t ; } 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' ; }

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 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 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

Сыроежка: Well, I would like to make some comments. I did not look through all your code but even the very beginning of the code arises questions. For example let consider inclusing of headers #include "stdafx.h" #include <iostream> #include <vector> #include <math.h> #include <time.h> #include <ctime> #include <fstream> #include <iomanip> //#include <vector> //#include <algorithm> #include <stdio.h> #include <stdlib.h> using namespace std; Standard C headers in C++ are named with prefix 'c' and without any extension. So It would be correctly to write #include "stdafx.h" #include <iostream> #include <vector> #include <cmath> #include <ctime> //#include <ctime> This header duplicates the header above so it may be removed. #include <fstream> #include <iomanip> //#include <vector> //#include <algorithm> #include <cstdio> #include <cstdlib> Further it is a bad idea to use the usong directive. Try to write code without this directive. So I would remove the directive using namespace std; As I said I did not look through all your code but it vseems that defining global variables double *dxc; double *dyc; double *dzc; double *dx; double *dy; double *dz; is not a good idea. I think that the following code #ifndef min #define min(a,b) ((a) < (b) ? (a): (b)) #endif #ifndef max #define max(a,b) ((a) > (b) ? (a): (b)) #endif /* Find minimum value of an array and return its index */ int minarray(double *A, int l) { int i; int minind=0; for (i=0; i<l; i++) { if(A<A[minind]){ minind=i; } } return minind; } should be removed. As you write the program on C++ then you should use C++ features. There std::min, std::max, std::min_element, and std::max_element algorrithms in C++. You should use them instead of macros and user defined functions. I do not like for example the following function inline bool almost_equals( double a, double b ) { static constexpr double EPSILON = 1.0e-6 ; return std::fabs(a-b) < EPSILON ; } First of all C defines already the appropriate constant as DBL_EPSILON. And its value usually less than or equal to 1E-9. The same value can be obbtained in C++ by using member function epsilon of standard class std::numeric_limits. In my computer it is equal to 2.22045e-016 The standard defiens this value as Machine epsilon: the difference between 1 and the least value greater than 1 that is representable. Meaningful for all floating point types But as it occured it does not mean that x == x + some value less than 2.22045e-016 I run the following program in my computer tthat you also can test [pre2]#include "stdafx.h" #include <iostream> #include <climits> #include <limits> main() { std::cout << std::numeric_limits<double>::epsilon() << std::endl; std::cout << DBL_EPSILON << std::endl; double x = 1.0; double y = x + std::numeric_limits<double>::epsilon(); if ( x != y ) std::cout << "x != y\n"; else std::cout << "x == y \n"; double eps = 2.22045e-016; do { y = x + eps; std::cout << eps << std::endl; eps -= 0.00001e-016; } while ( x != y ); std::cout << "x = " << x << ", y = " << y << std::endl; return 0; }[/pre2] The result value of eps when x == y + eps is equal to 1.11022e-016.

asq19: I'll take your comments and apply them to the code and get back to you shortly Thanks alot

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){ grid<<Final_T[i-1]<<" --> "<<count-1<<endl; val=Final_T; } 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

полная версия страницы