static void draw_velocity ( void )
{
	int i, j;
	float x, y, h;

	h = 1.0f/NX;

	glColor3f ( 1.0f, 0.0f, 0.0f );
	glBegin(GL_LINES);
		for ( i=0 ; i<NX ; i++ ) {
			//x = (i-0.5f)*h;
			x = (i+0.5f)*h;

			for ( j=0 ; j<NY ; j++ ) {
				//y = (j-0.5f)*h;
				y = (j+0.5f)*h;
				
				glColor3f(1.0,1.0,1.0);
				glVertex2f(x, y );
				glColor3f(1.0,0.0,0.0);
				glVertex2f(x+g_u[IX(i,j,0)], y+g_v[IX(i,j,0)] );
			}
		}
	glEnd();
}

static A jttayamp(J jt,A w,B nf,A x,A h){A y;B ng=!nf;I j,n;V*v=VAV(h);
 ASSERT(AR(x)<=(nf?v->lr:v->rr),EVRANK);
 switch(v->id){
  case CPLUS:  R tpoly(over(x,one));
  case CMINUS: R tpoly(nf?over(x,num[-1]):over(negate(x),one));
  case CSTAR:  R tpoly(over(zero,x));
  case CDIV:   ASSERT(ng,EVDOMAIN); R tpoly(over(zero,recip(x)));
  case CJDOT:  R tpoly(nf?over(x,a0j1):over(jdot1(x),one));
  case CPOLY:  ASSERT(nf,EVDOMAIN); R tpoly(BOX&AT(x)?poly1(x):x);
  case CHGEOM: ASSERT(nf,EVDOMAIN); RE(j=i0(x)); ASSERT(0<=j,EVDOMAIN);
               y=IX(j);
               R tpoly(divide(hgcoeff(y,h),fact(y)));
  case CBANG:  ASSERT(nf,EVDOMAIN); RE(j=i0(x)); ASSERT(0<=j,EVDOMAIN); 
               R tpoly(divide(poly1(box(iota(x))),fact(x)));
  case CEXP:   if(nf)R eva(x,"(^.x)&^ % !");
               RE(n=i0(x));   
               R 0<=n?tpoly(over(reshape(x,zero),one)):atop(ds(CDIV),amp(h,sc(-n))); 
  case CFIT:   ASSERT(nf&&CPOLY==ID(v->f),EVDOMAIN);
               y=over(x,IX(IC(x)));
               R tpoly(mdiv(df2(x,y,h),atab(CEXP,y,IX(IC(x)))));
  case CCIRCLE:
   switch(i0(x)){
    case 1:    R eval("{&0 1 0 [email protected](4&|) % !");
    case -3:   R eval("{&0 1 0 [email protected](4&|) % ]");
    case 2:    R eval("{&1 0 _1 [email protected](4&|) % !");
    case 5:    R eval("2&|    % !");
    case -7:   R eval("2&|    % ]");
    case 6:    R eval("2&|@>: % !");
    case -1:   R eval("(2&|              % ]) * ([: */ (1&+ % 2&+)@([email protected]<.&.-:))\"0");
    case -5:   R eval("({&0 1 0 [email protected](4&|) % ]) * ([: */ (1&+ % 2&+)@([email protected]<.&.-:))\"0");
 }}
 ASSERT(0,EVDOMAIN);
}

static void get_from_UI(ScalarField *d, VectorField *u_v) {
    int i, j, size = (N + 2) * (N + 2);

    for (i = 0; i < size; i++) {
        (*u_v)[i][0] = (*u_v)[i][1] = (*d)[i] = 0.0f;
    }

    if (!mouse_down[0] && !mouse_down[2]) return;

    i = (int) ((mx / (float) win_x) * N + 1);
    j = (int) (((win_y - my) / (float) win_y) * N + 1);

    if (i < 1 || i > N || j < 1 || j > N) return;

    if (mouse_down[0]) {
        (*u_v)[IX(i, j)][0] = force * (mx - omx);
        (*u_v)[IX(i, j)][1] = force * (omy - my);
    }

    if (mouse_down[2]) {
        (*d)[IX(i, j)] = source;
    }

    omx = mx;
    omy = my;

    return;
}

//-----------------------------------------------------------------------------
inline float StamFluidSolver::getCurl(int i, int j)
{
    float du_dy = (u[IX(i, j + 1)] - u[IX(i, j - 1)]) * 0.5f;
    float dv_dx = (v[IX(i + 1, j)] - v[IX(i - 1, j)]) * 0.5f;

    return du_dy - dv_dx;
}

static void get_from_UI ( float * d, float * u, float * v )
{
	int i, j, size = (N+2)*(N+2);

	for ( i=0 ; i<size ; i++ ) {
		u[i] = v[i] = d[i] = 0.0f;
	}

	if ( !mouse_down[0] && !mouse_down[2] ) return;

	i = (int)((       mx /(float)win_x)*N+1);
	j = (int)(((win_y-my)/(float)win_y)*N+1);

	if ( i<1 || i>N || j<1 || j>N ) return;

	if ( mouse_down[0] ) {
		u[IX(i,j)] = force * (mx-omx);
		v[IX(i,j)] = force * (omy-my);
	}

	if ( mouse_down[2] ) {
		d[IX(i,j)] = source;
	}

	omx = mx;
	omy = my;

	return;
}

static void draw_density ( void )
{
	int i, j;
	float x, y, h, d00, d01, d10, d11;

	h = 1.0f/N;

	glBegin ( GL_QUADS );

		for ( i=0 ; i<=N ; i++ ) {
			x = (i-0.5f)*h;
			for ( j=0 ; j<=N ; j++ ) {
				y = (j-0.5f)*h;

				d00 = dens[IX(i,j)];
				d01 = dens[IX(i,j+1)];
				d10 = dens[IX(i+1,j)];
				d11 = dens[IX(i+1,j+1)];

				glColor3f ( d00, d00, d00 ); glVertex2f ( x, y );
				glColor3f ( d10, d10, d10 ); glVertex2f ( x+h, y );
				glColor3f ( d11, d11, d11 ); glVertex2f ( x+h, y+h );
				glColor3f ( d01, d01, d01 ); glVertex2f ( x, y+h );
			}
		}

	glEnd ();
}

///////////////////////////////////////////////////////////////////////////////
/// Calcuate source term in the difussion equation
///
///\param para Pointer to FFD parameters
///\param var Pointer to FFD simulation variables
///\param var_type Type of variable 
///\param index Index of trace substances or species
///
///\return 0 if no error occurred
///////////////////////////////////////////////////////////////////////////////
int source_diff(PARA_DATA *para, REAL **var, int var_type, int index) {
  int i, j, k;  
  int imax = para->geom->imax, jmax = para->geom->jmax; 
  int kmax = para->geom->kmax;
  int IMAX = imax+2, IJMAX = (imax+2)*(jmax+2);
  REAL *b = var[B];
 
  FOR_EACH_CELL
    switch(var_type) {
      case VX:
        b[IX(i,j,k)] += var[VXS][IX(i,j,k)];   
        break;
      case VY: 
        b[IX(i,j,k)] += var[VYS][IX(i,j,k)];  
        break;
      case VZ: 
        b[IX(i,j,k)] += var[VZS][IX(i,j,k)]; 
        break;
      case TEMP: 
        b[IX(i,j,k)] += var[TEMPS][IX(i,j,k)];
        break;
      case TRACE:  
        b[IX(i,j,k)] += var[TRACE+para->bc->nb_Xi+index][IX(i,j,k)];
        break;
    }
  END_FOR

  return 0;
} // End of source_diff()

static void draw_velocity ( void )
{
	int i, j;
	float x, y, h, linescale;

	h = 1.0f/N;
  linescale = 0.2f;;

	glColor3f ( 1.0f, 1.0f, 1.0f );
	glLineWidth ( 1.0f );

	glBegin ( GL_LINES );

		for ( i=1 ; i<=N ; i++ ) {
			x = (i-0.5f)*h;
			for ( j=1 ; j<=N ; j++ ) {
				y = (j-0.5f)*h;

        

				glVertex2f ( x, y );
				glVertex2f ( x+(u[IX(i,j)]*linescale), y+(v[IX(i,j)]*linescale) );
			}
		}

	glEnd ();
}

void
dump_array (int step_count, char* name, int N, float scale, float* d)
{
    int i, j;
    char str[256];

    // find maximum value in image.
    float max = 0;
    for (j = 1; j < N+1; j++) {
        for (i = 1; i < N+1; i++) {
            float d00       = d[IX(i, j)] * 255 * scale;
            if (d00 >= max) max = d00;
        }
    }

    // Write out the file
    snprintf(str, 256, "out/%04d-%s.ppm", step_count, name);
    FILE* file = fopen (str, "w+");
    fprintf (file, "P2\n");
    fprintf (file, "%d %d\n", N, N);
    fprintf (file, "%d\n", (int)max);


    // Write out image file.
    for (j = N; j >= 1; j--) {
        for (i = 1; i < N+1; i++) {
            float d00       = d[IX(i, j)] * 255 * scale;
            fprintf (file, "% 4d ", (int)d00);
        }

        fprintf(file, "\n");
    }

    fclose(file);
}

void print_vel(float *u, float *v, int Nx, int Ny)
{
    int i, j;
    for (j=1;j<=Ny;j++)
        for (i=1;i<=Nx;i++)
            fprintf(file, "%f %f ", u[IX(i, j)], v[IX(i, j)]);
    fprintf(file, "\n");
}

void linear_combine_sub(int N, float * f_out, float * f, float * f0)
{
	for (int i = 1; i <= N; i++){
		for (int j = 1; j <= N; j++){
			f_out[IX(i, j)] = f[IX(i, j)] - f0[IX(i, j)];
		}
	}
}

//static void get_from_UI ( float * d, float * u, float * v, float * heat, float * compressibility, CombustionParticleSystem* cps )
static void get_from_UI ( float * d, float * u, float * v, float * heat, float * compressibility)
{
	int i, j, size = NX*NY*NZ;

	/*
	for ( i=0 ; i<size ; i++ ) {
		u[i] = v[i] = w[i] = d[i] = 0.0f;
	}
	*/

	if ( !mouse_down[0] && !mouse_down[2] ) return;

	i = (int)((       mx /(float)win_x)*NX);
	j = (int)(((win_y-my)/(float)win_y)*NY);

	if ( i<2 || i>=NX-1 || j<2 || j>=NY-1 ) return;

	if ( mouse_down[0] ) {
		/*
		compressibility[IX(i,j)] = 1.0f;
		compressibility[IX(i+1,j)] = 1.0f;
		compressibility[IX(i,j+1)] = 1.0f;
		compressibility[IX(i+1,j+1)] = 1.0f;
		heat[IX(i,j)] = 0.05f;
     */
		d[IX(i,j,0)] = source;

		float diffx = mx - omx;
		float diffy = win_y - my - omy;
		
		if(abs(diffx) > 0.000001)
			diffx = diffx/diffx * diffx<0.0f?-1.0f:1.0f;
		if(abs(diffy) > 0.000001)
			diffy = diffy/diffy * diffy<0.0f?-1.0f:1.0f;
		
		u[IX(i,j,0)] = 0.1f * force * diffx;
		v[IX(i,j,0)] = 0.1f * force * diffy;

//		u[IX(i-1,j,0)] = 0.075f *  force;
//		u[IX(i,j,0)] = 0.075f *  force;
//		u[IX(i+1,j,0)] = 0.075f *  force;
//		v[IX(i,j+1,0)] = 0.5f * force * diffy;
	}

	/*
	if ( mouse_down[2] ) {
		cps->trigger(mx, my, 100);
		//d[IX(i,j)] = source;
	}
	*/

	omx = mx;
	omy = win_y-my;

	return;
}

/*
	Add buoyancy to the density
	- Bouyancy adds a directional velocity to the current velocity of the cell
	- Buoyance is equal to timestep * (strength of the buoyance) * (direction of buoyance) * (density amount at cell)
	- Velocity gets higher the more density there is in a cell
*/
void addDensityBuoyancy (int N, float *velx, float *vely, float *dens, float bdirx, float bdiry, float bstrength, float dt) {
{	
	for (int j=1; j<=N; j++) {
		for (int i=1; i<=N; i++) {
			velx[IX(i,j)] += dt * bstrength * bdirx * dens[IX(i,j)];
			vely[IX(i,j)] += dt * bstrength * bdiry * dens[IX(i,j)];
		}
	}
}
}

void get_from_UI ( float * d, float * u, float * v, float elapsed, TCOD_key_t k,TCOD_mouse_t mouse ) {
	int i,j;
	float vx=0.0f,vy=0.0f;

	stepDelay -= elapsed;
	if ( stepDelay < 0.0f ) {
		if ( TCODConsole::isKeyPressed(TCODK_UP) && playery > 0 ) {
			playery--;
			vy -= force;
		} 
		if ( TCODConsole::isKeyPressed(TCODK_DOWN) && playery < N/2-1 ) {
			playery++;
			vx += force;
		} 
		if ( TCODConsole::isKeyPressed(TCODK_LEFT) && playerx > 0 ) {
			playerx--;
			vx -= force;
		}
		if ( TCODConsole::isKeyPressed(TCODK_RIGHT) && playerx < N/2-1 ) { 
			playerx++;
			vx += force;
		}
		stepDelay=0.2f; // move 5 cells per second
		// try to move smoke when you walk inside it. doesn't seem to work...
		u[IX(playerx*2,playery*2)] = 5*vx;
		v[IX(playerx*2,playery*2)] = 5*vy;
	}

	for ( i=0 ; i<SIZE; i++ ) {
		u[i] = v[i] = d[i] = 0.0f;
	}

	if ( !mouse.lbutton && !mouse.rbutton ) return;

	i = mouse.cx*2;
	j = mouse.cy*2;
	if ( i<1 || i>N || j<1 || j>N ) return;

	if ( mouse.lbutton ) {
		float dx,dy,l;
		dx=mouse.cx-playerx;
		dy=mouse.cy-playery;
		l=sqrt(dx*dx+dy*dy);
		if ( l > 0 ) {
			l = 1.0f/l;
			dx*=l;
			dy*=l;
			u[IX(playerx*2,playery*2)] = force * dx;
			v[IX(playerx*2,playery*2)] = force * dy;
			d[IX(playerx*2,playery*2)] = source;
		}
	}

}

/* A linear solver*/
void lin_solve(int N, int b, float *x, float *x0, float a, float c)
{
    int i, j, k;
    
    for(k=0;k<20;k++)
    {
        FOR_EACH_CELL
        x[ IX(i,j) ] = ( x0[ IX(i,j) ] + a*(x[ IX(i-1,j) ] + x[ IX(i+1, j) ] + x[IX(i, j-1)] + x[ IX(i, j+1) ])) /c;
        END_FOR
        set_bnd(N, b, x);
    }
}

Matrix & Matrix::homogenize(void) {

  int m = this->_m, n = this->_n;

  float factor = this->matrix[ IX(m-1,n-1) ];

  for (int i=0; i<m; i++) {
    for (int j = 0; j < n; j++) {
      this->matrix[ IX(i,j)] /= factor;
    }
  }

  return *this;

}

///////////////////////////////////////////////////////////////////////////////
/// Check the maximum value of the scalar psi at (ci,cj,ck) and its surrounding
/// cells
///
///\param para Pointer to FFD parameters
///\param psi Pointer to the variable
///\param ci Index in x direction
///\param cj Index in y direction
///\param ck Index in z direction
///
///\return 0 if no error occurred
///////////////////////////////////////////////////////////////////////////////
REAL check_max(PARA_DATA *para, REAL *psi, int ci, int cj, int ck) {
  int imax = para->geom->imax, jmax = para->geom->jmax; 
  int i, j, k;
  int IMAX = imax+2, IJMAX = (imax+2)*(jmax+2);
  REAL tmp = psi[IX(ci,cj,ck)];

  for(i=0;i<=1;i++)
    for(j=0;j<=1;j++)
      for(k=0;k<=1;k++) {
        if(tmp<psi[IX(ci+i,cj+j,ck+k)]) tmp=psi[IX(ci+i,cj+j,ck+k)];
      }
    
return tmp;

}// End of check_max( )

//-----------------------------------------------------------------------------
void StamFluidSolver::linearSolve ( int b, float * x, float * x0, float a, float c )
{
    // NOTE that I don't use linear solve in the new code

    for (int k=0 ; k<6 ; k++ )
    {
        for (int j=1 ; j < _NY-1 ; j++ ) {
            for (int i=1 ; i < _NX-1; i++ )
            {
                x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]))/c;
            }
        }
        setBound ( b, x );
    }
}

    void eval(MonteCarloNode& n) {
        auto& b = n.board;
        auto sym = n.orientation;

        feats.fill(Features());
        oriented.fill(Features());
        extractFeatures(b,feats);
        dihedralTranspose(feats,oriented,sym);

        inVec.clear();
        convertToTCNNInput(oriented,inVec);

        result = &nn.fprop(inVec);
        std::vector<size_t> inds;
        inds.assign(n.moves.size(),0);

        double total = 0;
        for(size_t i = 0; i < n.moves.size(); ++i) {
            inds[i] = /*IX(n.moves[i]);*/IX(dihedral(n.moves[i],sym));
            total += (*result)[inds[i]];
        }

        for(size_t i = 0; i < n.moves.size(); ++i) {
            n.probabilities[i].store((*result)[inds[i]]/total);
        }

    }

/******************************************************************************
  Read the data send from the other program
******************************************************************************/
int read_cosimulation_data(PARA_DATA *para, REAL **var)
{
  int j, k;
  int imax = para->geom->imax, jmax = para->geom->jmax;
  int kmax = para->geom->kmax;
  int IMAX = imax+2, IJMAX = (imax+2)*(jmax+2);
  otherSharedData data;
  REAL feak[1];

  if(read_from_shared_memory(&data))
    exit(1);
  else
    ffd_log("cosimulation.c: read data from shared memory.", FFD_NORMAL);

  feak[0] = data.arr[0];

  /*--------------------------------------------------------------------------
  | The following code is to be modified by the users
  --------------------------------------------------------------------------*/
  for(j=0; j<=jmax+1; j++)
    for(k=0; k<=kmax+1; k++)
      var[TEMPBC][IX(imax+1,j,k)]= feak[0];

   
  printf("\ntime=%f, status=%d\n", data.t, data.status);
  printf("arr[0]=%f, arr[1]=%f, arr[2]=%f\n", data.arr[0], data.arr[1], data.arr[2]);

   //printf("message=%s\n",data.message); 

  return 0;
} // End of read_cosimulation_data()