void testCopy() {
			DataKeyDimension dim(*m_dim1);
			CPPUNIT_ASSERT(dim==*m_dim1);
		}

double & band_matrix::saved_diag(int i)
{
    assert( (i>=0) && (i<dim()) );
    return m_lower[0][i];
}

 const SpdMatrix & IMMGS::Sigma() const {
   sigma_scratch_.resize(dim());
   sigma_scratch_.diag() = unscaled_variance_diagonal();
   sigma_scratch_.diag() *= sigsq();
   return sigma_scratch_;
 }

int main(int argc, char** argv){// 1000G.marker.gz -> rs pos A1 A2
	op = (char*)calloc(100, sizeof(char));
	if(argc==1){fprintf(stderr,"tabix tbxfile CHR | ./ase 1000G.marker.gz [-printFrag] [-overlap] > output\n"); return -1;}
	
	int i;
	char* buf; buf = (char*)calloc(1000,sizeof(char));
	long printFrag=0;
	long overlap=0;
	//long ascount;
	for(i=0; i<argc; i++){if(strcmp(argv[i],"-printFrag")==0){printFrag=1;}}
	for(i=0; i<argc; i++){if(strcmp(argv[i],"-overlap")==0){overlap=1;}}
	//for(i=0; i<argc; i++){if(strcmp(argv[i],"-ascount")==0){ascount=1;}}
	
	
	long N=0;
	gzFile fmarker;
	
	if((fmarker = gzopen(argv[1], "rb6f"))==NULL){fprintf(stderr, "gzread failed. no file!\n"); return -1;}
	long P=0;
	//N = getLine(fmarker);
	dim(fmarker, &buf, 1000*sizeof(char), &N, &P, 0);
fprintf(stderr,"N of markers: %ld %ld\n", N, P);
	
	char** rs;
	long* pos;
	char* cpos;
	char* A1;
	char* A2;
	long* Y;
	double* ds;
	rs = (char**)calloc(N,sizeof(char*));  for(i=0;i<N;i++){ rs[i]=(char*)calloc(100,sizeof(char)); }
	cpos = (char*)calloc(100,sizeof(char));
	pos = (long*)calloc(N,sizeof(long));
	ds = (double*)calloc(N,sizeof(double));
	A1 = (char*)calloc(N,sizeof(char));
	A2 = (char*)calloc(N,sizeof(char));
	Y = (long*)calloc(N*3,sizeof(long));
	
	for(i=0; i<N; i++){ds[i]=0.0;}
	for(i=0; i<N*3; i++){Y[i]=0;}
	int ret;
	//char* buf;
	int l=0,cl=0,ncol=0;
	//buf = (char*)calloc(1000,sizeof(char));
	while((ret = gzread(fmarker, buf, sizeof(buf))) != 0 && ret!=-1){
		for(i=0;i<ret/sizeof(char);i++){
			if(buf[i]=='\n'){
				l++;
				ncol=cl=0;
			}else if(buf[i]==' ' || buf[i]=='\t' || buf[i]=='\n'){
				if(ncol==2){rs[l][cl]='\0';}else if(ncol==1){cpos[cl]='\0'; pos[l]=atoi(cpos);}
				if(buf[i]=='\n'){l++; ncol=cl=0;}else{cl=0; ncol++;}
			}else{
				if(ncol==2){
					rs[l][cl++]=buf[i];
				}else if(ncol==1){
					cpos[cl++]=buf[i];
				}else if(ncol==3){
					if(cl==0){A1[l]=buf[i];}else{A1[l]='0';}
					cl++;
				}else if(ncol==4){
					if(cl==0){
						A2[l]=buf[i];
						if(A2[l]=='-'){A2[l]='1';A1[l]='0';}
					}else{// allele length > 1
						A2[l]='1';A1[l]='0';
					}
					cl++;
				}
			}
		}
	}
	fprintf(stderr, "marker position: [%ld..%ld]\n",pos[0],pos[N-1]);
	//fprintf(stderr, "marker position: [%c %c %c %c..%c %c %c %c]\n",A1[0],A1[1],A1[2],A1[3],A1[N-4],A1[N-3],A1[N-2],A1[N-1]);
	//fprintf(stderr, "marker position: [%c %c %c %c..%c %c %c %c]\n",A2[0],A2[1],A2[2],A2[3],A2[N-4],A2[N-3],A2[N-2],A2[N-1]);
	
	char* chr;
	char* c1;
	char* seq;
	char* num;
	chr = (char*)calloc(100, sizeof(char));
	seq = (char*)calloc(1000, sizeof(char));
	c1 =  (char*)calloc(1000, sizeof(char));
	num = (char*)calloc(100, sizeof(char));
	if(num==NULL || c1==NULL || seq==NULL || chr==NULL){fprintf(stderr, "memory allocation error!\n"); return 0;}
	long left, right;
	long* at;
	at = (long*)calloc(100, sizeof(long));
	long j,j0=0;
	long l0=0, l1=0;
	while(scanf("%s\t%ld\t%ld\t%s\t%s", chr, &left, &right, c1, seq)!=EOF){
		l=0;
		int flag=0;
//fprintf(stderr, "%s\t%ld\t%ld\t%s\t%s\n", chr, left, right, c1, seq);
		for(j=j0; j<N; j++){
			if(pos[j]>right){break;}
			if(pos[j]>=left){l++;}else{j0=j+1;}
		}
		if(l>0){
			// flag: number of snps existing in the fragment
//.........这里部分代码省略.........

/**
 * For that test you would need to luanch camera
 * so it would punch a hole for you.
 */
int main(int, char**)
{
   LOGE("%s start", LOG_TAG);

   overlay2::GenericPipe<ovutils::FB0> pipe;
   overlay2::RotatorBase* r = new overlay2::NullRotator();
   LOGE("Using null rotator");
   bool ret = pipe.open(r);
   OVASSERT(ret, "Failed to open pipe");

   LOGE("start ...");
   ovutils::Whf whf(1024, 600, HAL_PIXEL_FORMAT_YCbCr_420_SP_TILED);
   /* since we are using camera apk for punching a hole for us
    * we need ZOERDER_1 here so it will be on top of camera */
   ovutils::eMdpFlags flags = ovutils::OV_MDP_FLAGS_NONE;
   ovutils::PipeArgs args(flags,
                          ovutils::OVERLAY_TRANSFORM_0,
                          whf,
                          ovutils::WAIT,
                          ovutils::ZORDER_1,
                          ovutils::IS_FG_OFF,
                          ovutils::ROT_FLAG_DISABLED,
                          ovutils::PMEM_SRC_SMI,
                          ovutils::RECONFIG_OFF);
   ret = pipe.start(args);

   OVASSERT(ret, "Failed to start pipe");

   // let's try openning /dev/zero and mmap it.
   LOGE("Open /dev/pmem_adsp");
   overlay2::OvMem m;
   uint32_t bfsz = 1024*600*4;
   ret = m.open(ovutils::PMEM_SRC_ADSP,
                1,     // one buf
                bfsz); //size
   OVASSERT(ret, "Failed to open /dev/pmem_adsp pipe");
   m.dump();
   ::memset(m.addr(), 0, bfsz);
   pipe.setMemoryId(m.getFD());
   ovutils::Dim dim(0, 0, 1024, 600); // x,y,w,h
   ret = pipe.setPosition(dim);
   OVASSERT(ret, "Failed to setPosition");

   LOGE("Before play ...");
   pipe.dump();
   enum { NUM_PLAY = 10 };
   for (uint32_t i=0; i < NUM_PLAY; ++i) {
      ret = pipe.queueBuffer(0); // offset 0
      OVASSERT(ret, "Failed to queueBuffer");
      sleep(1);
   }

   ret = r->close();
   OVASSERT(ret, "Error in rotator close");

   ret = pipe.close();
   OVASSERT(ret, "Error in pipe close");

   ret = m.close();

   delete r;
   r=0;
   LOGE("%s end", LOG_TAG);
   return 0;
}

                }
                pKbdEventRmpCurr->uiVk = pvkMap->uiVkGenerated;
            }
        }
    }

    return cRmpKbdEvents;    
}
    

static DEVICE_LAYOUT dlMatrixEngUs =
{
    sizeof(DEVICE_LAYOUT),
    MATRIX_PDD,
    rgscvkMatrixEngUSTables,
    dim(rgscvkMatrixEngUSTables),
    MatrixUsRemapVKey,
};

extern "C"
BOOL
Matrix(
    PDEVICE_LAYOUT pDeviceLayout
    )
{
    SETFNAME(_T("Matrix"));

    BOOL fRet = FALSE;

    if (pDeviceLayout == NULL)
    {

 // second diag (used in LU decomposition), saved in m_lower
 double CBandMatrix::saved_diag(int i) const
 {
     assert( (i>=0) && (i<dim()) );
     return m_lower[0][i];
 }

void Action_AtomicCorr::Print() {
  int idx, idx3, vec1size;
  Vec3 V1, V2;
  mprintf("    ATOMICCORR: Calculating correlations between %s vectors:\n",
          ModeString[acorr_mode_]);
  if (atom_vectors_.empty()) {
    mprinterr("Error: No vectors calculated.\n");
    return;
  }
  DataSet_MatrixFlt* tmatrix = static_cast<DataSet_MatrixFlt*>( dset_ );
  tmatrix->AllocateTriangle( atom_vectors_.size() );
  // Calculate correlation coefficient of each atomic vector to each other
  ACvector::iterator av_end = atom_vectors_.end();
  ACvector::iterator av_end1 = av_end - 1;
  ProgressBar progress( tmatrix->Size() );
  int iprogress = 0;
  for (ACvector::const_iterator vec1 = atom_vectors_.begin(); vec1 != av_end1; ++vec1)
  {
    for (ACvector::const_iterator vec2 = vec1 + 1; vec2 != av_end; ++vec2)
    {
      double corr_coeff = 0.0;
      progress.Update( iprogress++ );
      // If vectors are too close, skip. vec2 always > vec1
      if ( (*vec2) - (*vec1) > min_ ) {
        // Make sure same # of frames in each and not empty
        if ( vec1->empty() || vec2->empty() ) {
          mprintf("Warning: A vector is empty: Vec%zu=%zu, Vec%zu=%zu\n",
                  vec1 - atom_vectors_.begin(), vec1->size(),
                  vec2 - atom_vectors_.begin(), vec2->size());
        } else if ( vec1->size() != vec2->size() ) {
          mprintf("Warning: Vec %zu and Vec %zu do not have same # of frames.\n",
                  vec1 - atom_vectors_.begin(), vec2 - atom_vectors_.begin());
        } else {
          vec1size = (int)(vec1->size() / 3);
#ifdef _OPENMP
#pragma omp parallel private(idx, idx3, V1, V2) reduction(+: corr_coeff)
{
#pragma omp for
#endif
          for (idx = 0; idx < vec1size; ++idx) {
            idx3 = idx * 3;
            V1 = vec1->VXYZ(idx3);
            V2 = vec2->VXYZ(idx3);
            V1.Normalize();
            V2.Normalize();
            corr_coeff += (V1 * V2);
          }
#ifdef _OPENMP
} // END pragma omp parallel
#endif
          corr_coeff /= (double)vec1size;
          if (fabs(corr_coeff) <= cut_)
            corr_coeff = 0.0;
        }
      }
      tmatrix->AddElement( corr_coeff );
    } // END inner loop
  } // END outer loop

  if (acorr_mode_ == ATOM) {
    Dimension dim(1.0, 1.0, "Atom");
    dset_->SetDim(Dimension::X, dim);
    dset_->SetDim(Dimension::Y, dim);
  } else {
    Dimension dim(1.0, 1.0, "Residue");
    dset_->SetDim(Dimension::X, dim);
    dset_->SetDim(Dimension::Y, dim);
  }
  std::string labels;
  for (ACvector::const_iterator atom = atom_vectors_.begin();
                                atom != atom_vectors_.end(); ++atom)
    labels += ( atom->Label() + "," );
  if (outfile_ != 0)
    outfile_->ProcessArgs( "xlabels " + labels + " ylabels " + labels );
}

separation find_wheel_minor(matroid_permuted <MatroidType>& permuted_matroid, matrix_permuted <MatrixType>& permuted_matrix,
                            matroid_element_set& extra_elements)
{
    assert (permuted_matrix.size1() >= 3 && permuted_matroid.size2() >= 3);

    /// Permute 1's in first row to the left.
    matroid_reorder_columns(permuted_matroid, permuted_matrix, 0, 1, 0, permuted_matroid.size2(), std::greater <int>());
    size_t count_first_row_ones = matrix_count_property_column_series(permuted_matrix, 0, 1, 0, permuted_matrix.size2(), is_non_zero());

    /// Trivial 1-separation
    if (count_first_row_ones == 0)
    {
        return separation(std::make_pair(1, 0));
    }

    matroid_reorder_rows(permuted_matroid, permuted_matrix, 1, permuted_matroid.size1(), 0, 1, std::greater <int>());
    size_t count_first_column_ones = matrix_count_property_row_series(permuted_matrix, 0, permuted_matroid.size1(), 0, 1, is_non_zero());

    /// 1- or 2-separation
    if (count_first_row_ones == 1)
    {
        if (count_first_column_ones == 0)
        {
            return separation(std::make_pair(1, 1));
        }
        else
        {
            return separation(std::make_pair(1, 1), std::make_pair(1, 0));
        }
    }
    else if (count_first_column_ones == 1)
    {
        return separation(std::make_pair(1, 1), std::make_pair(0, 1));
    }

    assert ((permuted_matrix(0,0) == 1) && (permuted_matrix(1,0) == 1) && (permuted_matrix(0,1) == 1));

    /// Ensure we have a 2x2 block of ones
    if (permuted_matrix(1, 1) != 1)
    {
        matroid_binary_pivot(permuted_matroid, permuted_matrix, 0, 0);
        extra_elements.insert(permuted_matroid.name1(0));
        extra_elements.insert(permuted_matroid.name2(0));
    }

    assert ((permuted_matrix(0,0) == 1) && (permuted_matrix(1,0) == 1) && (permuted_matrix(0,1) == 1) && (permuted_matrix(1,1) == 1));

    /// Grow the block to a set-maximal one.
    matroid_reorder_columns(permuted_matroid, permuted_matrix, 0, 2, 2, permuted_matroid.size2(), std::greater <int>());
    size_t block_width = 2 + matrix_count_property_column_series(permuted_matrix, 0, 2, 2, permuted_matrix.size2(), is_all_ones());

    matroid_reorder_rows(permuted_matroid, permuted_matrix, 2, permuted_matroid.size1(), 0, block_width, std::greater <int>());

    size_t block_height = 2 + matrix_count_property_row_series(permuted_matrix, 2, permuted_matrix.size1(), 0, block_width, is_all_ones());

    /// Search for a path in BFS graph
    zero_block_matrix_modifier modifier(block_height, block_width);
    matrix_modified <matrix_permuted <MatrixType> , zero_block_matrix_modifier> modified_matrix(permuted_matrix, modifier);
    bipartite_graph_dimensions dim(permuted_matrix.size1(), permuted_matrix.size2());
    std::vector <bipartite_graph_dimensions::index_type> start_nodes(block_height);
    for (size_t i = 0; i < block_height; i++)
        start_nodes[i] = dim.row_to_index(i);
    std::vector <bipartite_graph_dimensions::index_type> end_nodes(block_width);
    for (size_t i = 0; i < block_width; i++)
        end_nodes[i] = dim.column_to_index(i);

    std::vector <bipartite_graph_bfs_node> bfs_result;
    bool found_path = bipartite_graph_bfs(modified_matrix, dim, start_nodes, end_nodes, false, bfs_result);

    /// There must be a 2-separation.
    if (!found_path)
    {
        std::pair <size_t, size_t> split(0, 0);

        /// Swap unreachable rows to top
        std::vector <int> reachable(permuted_matrix.size1());
        for (size_t i = 0; i < permuted_matrix.size1(); ++i)
        {
            const bipartite_graph_bfs_node& node = bfs_result[dim.row_to_index(i)];
            int value = (node.is_reachable() ? (node.distance > 0 ? 2 : 1) : 0);
            reachable[permuted_matrix.perm1()(i)] = value;
            if (value == 0)
                split.first++;
        }
        vector_less <int, std::less <int> > less(reachable, std::less <int>());

        sort(permuted_matrix.perm1(), less);
        permuted_matroid.perm1() = permuted_matrix.perm1();

        /// Swap unreachable columns to left
        reachable.resize(permuted_matrix.size2());
        for (size_t i = 0; i < permuted_matrix.size2(); ++i)
        {
            const bipartite_graph_bfs_node& node = bfs_result[dim.column_to_index(i)];
            int value = (node.is_reachable() ? 2 : (node.distance == -2 ? 1 : 0));
            reachable[permuted_matrix.perm2()(i)] = value;
            if (value < 2)
                split.second++;
        }

//.........这里部分代码省略.........

    /*--------------------------------------------------------------------------------
        Function    : InventoryScreen::refresh
        Description : draws the InventoryScreen to the game window
        Inputs      : None
        Outputs     : updated InventoryScreen
        Return      : void
    --------------------------------------------------------------------------------*/
    void InventoryScreen::refresh()
    {
        Point dim(screen.getWidth(), screen.getHeight());

        // set the colorscheme of the selected line
        TCODConsole::setColorControl(TCOD_COLCTRL_1, TCODColor::white, TCODColor::blue);

        // create a vector of the different item symbols
        vector< pair<int, TCODColor> > tabChars;
        tabChars.push_back(pair<int, TCODColor>('\\', UI_FORE_COLOR));
        tabChars.push_back(pair<int, TCODColor>('[', UI_FORE_COLOR));
        tabChars.push_back(pair<int, TCODColor>('!', UI_FORE_COLOR));
        tabChars.push_back(pair<int, TCODColor>('*', UI_FORE_COLOR));

        // draw interface
        int x=0, y=2;
        //screen.printEx(dim.X()/2, 0, TCOD_BKGND_SET, TCOD_CENTER, 
        //               (client->personalStats.name() + string("'s Inventory")).c_str());
        screen.printEx(dim.X()/2, 0, TCOD_BKGND_SET, TCOD_CENTER, "Inventory");

        // 'start' is the point at which the folder contents are drawn
        Point start = drawFolders(x, y, dim, tabChars);

        // draw the information line
        x = start.X();
        y = start.Y();
        screen.print(x+1, y++, "Qt.  Item");
        screen.hline(x, y++, dim.X()-3);

        // update start position to where items will be drawn
        start.setY(y);
        
        linesPerPage = dim.Y()-y;
        currentPage = (selectedLine==0) ? 0 : selectedLine / (linesPerPage-1);

        int numItems = client->inventory.getIndexSize(selectedTab);

        if(numItems <= 0) return;

        drawScrollBar(Point(dim.X()-2, y), linesPerPage, currentPage, 
                      (numItems/linesPerPage)+1);

        // determine at what line should be highlighted
        int hlOffset = selectedLine % ((currentPage==0) ? linesPerPage : linesPerPage-1); 

        // set an iterator to the lowest spot that will show the selected line
        vector<ItemPtr>::const_iterator it;
        it = client->inventory.getIndexBegin(selectedTab) + (currentPage*linesPerPage) - currentPage;

        for(; it!=client->inventory.getIndexEnd(selectedTab); ++it, ++y)
        {
            // ignore lines that are off the display area
            if(y >= start.Y()+linesPerPage-1) break;

            // construct the display string
            // three spaces are allocated to the item count number.  The buffer string
            // ensures that this number is right aligned in this space.
            string buffer = "";
            unsigned int c = (*it)->getCount();
            while(c<100) { buffer += " "; c *= 10;}
            string disp = buffer + boost::lexical_cast<string>((*it)->getCount()) + "  " + (*it)->getName();

            // if this is the selected line on the current page, highlight it
            if(y == (hlOffset) + start.Y()) 
            {
                disp = "%c" + disp + "%c";
                screen.printEx(start.X(), y, TCOD_BKGND_SET, TCOD_LEFT, 
                               disp.c_str(), TCOD_COLCTRL_1, TCOD_COLCTRL_STOP);
            }
            else // display the line normally
            {
                screen.printEx(start.X(), y, TCOD_BKGND_SET, TCOD_LEFT, disp.c_str()); 
            }
        }
    }

/// \brief Sets the matrix to the \p nrows x \p ncols zero matrix.
void matrix::zero(int nrows, int ncols) {

    dim(nrows, ncols);
    for(int i=0;i<nf*nc;i++) p[i]=0;

}

int main(int, char**)
{
   LOGE("OverlayCtrlTest start");

   overlay2::Ctrl ctrl;
   overlay2::RotatorBase* rot = new overlay2::NullRotator();

   // open rot here
   bool ret = rot->open();
   OVASSERT(ret, "Rotator failed to open");
   ovutils::PipeArgs args(ovutils::OV_MDP_FLAGS_NONE,//flags PIPE SHARED
                          ovutils::OVERLAY_TRANSFORM_0,
                          ovutils::Whf(),
                          ovutils::WAIT,
                          ovutils::ZORDER_0,
                          ovutils::IS_FG_OFF,
                          ovutils::ROT_FLAG_ENABLED,
                          ovutils::PMEM_SRC_ADSP,
                          ovutils::RECONFIG_OFF);
   ret = rot->remap(ROT_NUM_BUF, args);
   OVASSERT(ret, "Rotator failed to remap");

   ret = ctrl.open(ovutils::FB0, rot);
   OVASSERT(ret, "Ctrl failed to open");
   ovutils::Whf whf(100, 200, HAL_PIXEL_FORMAT_RGBA_8888);
   ovutils::eMdpFlags flags = ovutils::OV_MDP_FLAGS_NONE;
   ret = ctrl.start(args);
   OVASSERT(ret, "Ctrl failed to start");
   ctrl.dump();

   LOGE("setPosition ...");
   ovutils::Dim dim(0, 0, 50, 100);
   ret = ctrl.setPosition(dim);
   OVASSERT(ret, "Ctrl failed to setPosition");
   ctrl.dump();

   LOGE("setParameter ...");
   ovutils::Params p ( ovutils::OVERLAY_DITHER, 0 );
   ret = ctrl.setParameter(p);
   OVASSERT(ret, "Ctrl failed to setParameter");
   ctrl.dump();

   LOGE("setSource ...");
   ovutils::PipeArgs args1(flags,
                           ovutils::OVERLAY_TRANSFORM_0,
                           ovutils::Whf(75, 150, HAL_PIXEL_FORMAT_RGBA_8888),
                           ovutils::WAIT,
                           ovutils::ZORDER_0,
                           ovutils::IS_FG_OFF,
                           ovutils::ROT_FLAG_DISABLED,
                           ovutils::PMEM_SRC_SMI,
                           ovutils::RECONFIG_OFF);
   ret = ctrl.setSource(args1);
   OVASSERT(ret, "Ctrl failed to setSource");
   ctrl.commit();
   ctrl.dump();

   ctrl.close();
   rot->close();
   delete rot;

   LOGE("OverlayCtrlTest end");
   return 0;
}

void SPxSolver::qualRedCostViolation(Real& maxviol, Real& sumviol) const
{   
   maxviol = 0.0;
   sumviol = 0.0;

   int i;
   // TODO:   y = c_B * B^-1  => coSolve(y, c_B)
   //         redcost = c_N - yA_N 
   // solve system "x = e_i^T * B^-1" to get i'th row of B^-1
   // DVector y( nRows() );
   // basis().coSolve( x, spx->unitVector( i ) );
   // DVector rdcost( nCols() );
#if 0 // un-const
   if (lastUpdate() > 0)
      factorize();

   computePvec();

   if (type() == ENTER)
      computeTest();
#endif
   if (type() == ENTER)
   {
      for(i = 0; i < dim(); ++i)
      {
         Real x = coTest()[i];
         
         if (x < 0.0)
         {
            sumviol -= x;
            
            if (x < maxviol)
               maxviol = x;
         }
      }
      for(i = 0; i < coDim(); ++i)
      {
         Real x = test()[i];
         
         if (x < 0.0)
         {
            sumviol -= x;
            
            if (x < maxviol)
               maxviol = x;
         }
      } 
   }
   else
   {
      assert(type() == LEAVE);

      for(i = 0; i < dim(); ++i)
      {
         Real x = fTest()[i];
         
         if (x < 0.0)
         {
            sumviol -= x;
            
            if (x < maxviol)
               maxviol = x;
         }
      }
   }
   maxviol *= -1;
}

void RandomPCA::pca(MatrixXd &X, int method, bool transpose,
   unsigned int ndim, unsigned int nextra, unsigned int maxiter, double tol,
   long seed, int kernel, double sigma, bool rbf_center,
   unsigned int rbf_sample, bool save_kernel, bool do_orth, bool do_loadings)
{
   unsigned int N;

   if(kernel != KERNEL_LINEAR)
   {
      transpose = false;
      verbose && std::cout << timestamp()
	 << " Kernel not linear, can't transpose" << std::endl;
   }

   verbose && std::cout << timestamp() << " Transpose: " 
      << (transpose ? "yes" : "no") << std::endl;

   if(transpose)
   {
      if(stand_method != STANDARDIZE_NONE)
	  X_meansd = standardize_transpose(X, stand_method, verbose);
      N = X.cols();
   }
   else
   {
      if(stand_method != STANDARDIZE_NONE)
	 X_meansd = standardize(X, stand_method, verbose);
      N = X.rows();
   }

   unsigned int total_dim = ndim + nextra;
   MatrixXd R = make_gaussian(X.cols(), total_dim, seed);
   MatrixXd Y = X * R;
   verbose && std::cout << timestamp() << " dim(Y): " << dim(Y) << std::endl;
   normalize(Y);
   MatrixXd Yn;

   verbose && std::cout << timestamp() << " dim(X): " << dim(X) << std::endl;
   MatrixXd K; 
   if(kernel == KERNEL_RBF)
   {
      if(sigma == 0)
      {
	 unsigned int med_samples = fminl(rbf_sample, N);
      	 double med = median_dist(X, med_samples, seed, verbose);
      	 sigma = sqrt(med);
      }
      verbose && std::cout << timestamp() << " Using RBF kernel with sigma="
	 << sigma << std::endl;
      K.noalias() = rbf_kernel(X, sigma, rbf_center, verbose);
   }
   else
   {
      verbose && std::cout << timestamp() << " Using linear kernel" << std::endl;
      K.noalias() = X * X.transpose() / (N - 1);
   }

   //trace = K.diagonal().array().sum() / (N - 1);
   trace = K.diagonal().array().sum();
   verbose && std::cout << timestamp() << " Trace(K): " << trace 
      << " (N: " << N << ")" << std::endl;

   verbose && std::cout << timestamp() << " dim(K): " << dim(K) << std::endl;
   if(save_kernel)
   {
      verbose && std::cout << timestamp() << " saving K" << std::endl;
      save_text("kernel.txt", K);
   }

   for(unsigned int iter = 0 ; iter < maxiter ; iter++)
   {
      verbose && std::cout << timestamp() << " iter " << iter;
      Yn.noalias() = K * Y;
      if(do_orth)
      {
	 verbose && std::cout << " (orthogonalising)";
	 ColPivHouseholderQR<MatrixXd> qr(Yn);
	 MatrixXd I = MatrixXd::Identity(Yn.rows(), Yn.cols());
	 Yn = qr.householderQ() * I;
	 Yn.conservativeResize(NoChange, Yn.cols());
      }
      else
	 normalize(Yn);

      double diff =  (Y -  Yn).array().square().sum() / Y.size(); 
      verbose && std::cout << " " << diff << std::endl;
      Y.noalias() = Yn;
      if(diff < tol)
	 break;
   }

   verbose && std::cout << timestamp() << " QR begin" << std::endl;
   ColPivHouseholderQR<MatrixXd> qr(Y);
   MatrixXd Q = MatrixXd::Identity(Y.rows(), Y.cols());
   Q = qr.householderQ() * Q;
   Q.conservativeResize(NoChange, Y.cols());
   verbose && std::cout << timestamp() << " dim(Q): " << dim(Q) << std::endl;
   verbose && std::cout << timestamp() << " QR done" << std::endl;

   MatrixXd B = Q.transpose() * X;
//.........这里部分代码省略.........

//! constructor
CGUIMessageBox::CGUIMessageBox(IGUIEnvironment* environment, const wchar_t* caption,
                               const wchar_t* text, s32 flags,
                               IGUIElement* parent, s32 id, core::rect<s32> rectangle)
	: CGUIWindow(environment, parent, id, rectangle), StaticText(0),
	OkButton(0), YesButton(0), NoButton(0), CancelButton(0)
{
	#ifdef _DEBUG
	setDebugName("CGUIMessageBox");
	#endif

	// remove focus
	Environment->setFocus(0);

	// remove buttons

	getMaximizeButton()->remove();
	getMinimizeButton()->remove();

	if (caption)
		setText(caption);

	IGUISkin* skin = Environment->getSkin();

	s32 buttonHeight = skin->getSize(EGDS_BUTTON_HEIGHT);
	s32 buttonWidth = skin->getSize(EGDS_BUTTON_WIDTH);
	s32 titleHeight = skin->getSize(EGDS_WINDOW_BUTTON_WIDTH)+2;
	s32 buttonDistance = skin->getSize(EGDS_WINDOW_BUTTON_WIDTH);

	// add static multiline text

	core::dimension2d<s32> dim(AbsoluteClippingRect.getWidth() - buttonWidth,
	                           AbsoluteClippingRect.getHeight() - (buttonHeight * 3));
	core::position2d<s32> pos((AbsoluteClippingRect.getWidth() - dim.Width) / 2,
	                          buttonHeight / 2 + titleHeight);

	StaticText = Environment->addStaticText(text,
	                                        core::rect<s32>(pos, dim), false, false, this);
	StaticText->setWordWrap(true);
	StaticText->grab();

	// adjust static text height

	s32 textHeight = StaticText->getTextHeight();
	core::rect<s32> tmp = StaticText->getRelativePosition();
	tmp.LowerRightCorner.Y = tmp.UpperLeftCorner.Y + textHeight;
	StaticText->setRelativePosition(tmp);
	dim.Height = textHeight;

	// adjust message box height

	tmp = getRelativePosition();
	s32 msgBoxHeight = textHeight + (s32)(2.5f * buttonHeight) + titleHeight;

	// adjust message box position

	tmp.UpperLeftCorner.Y = (parent->getAbsolutePosition().getHeight() - msgBoxHeight) / 2;
	tmp.LowerRightCorner.Y = tmp.UpperLeftCorner.Y + msgBoxHeight;
	setRelativePosition(tmp);


	// add buttons

	s32 countButtons = 0;
	if (flags & EMBF_OK) ++countButtons;
	if (flags & EMBF_CANCEL) ++countButtons;
	if (flags & EMBF_YES) ++countButtons;
	if (flags & EMBF_NO) ++countButtons;

	core::rect<s32> btnRect;
	btnRect.UpperLeftCorner.Y = pos.Y + dim.Height + buttonHeight / 2;
	btnRect.LowerRightCorner.Y = btnRect.UpperLeftCorner.Y + buttonHeight;
	btnRect.UpperLeftCorner.X = (AbsoluteClippingRect.getWidth() -
	                             ((buttonWidth + buttonDistance)*countButtons)) / 2;
	btnRect.LowerRightCorner.X = btnRect.UpperLeftCorner.X + buttonWidth;

	// add ok button
	if (flags & EMBF_OK)
	{
		OkButton = Environment->addButton(btnRect, this);
		OkButton->setText(skin->getDefaultText(EGDT_MSG_BOX_OK));
		OkButton->grab();

		btnRect.LowerRightCorner.X += buttonWidth + buttonDistance;
		btnRect.UpperLeftCorner.X += buttonWidth + buttonDistance;

		Environment->setFocus(OkButton);
	}

	// add yes button
	if (flags & EMBF_YES)
	{
		YesButton = Environment->addButton(btnRect, this);
		YesButton->setText(skin->getDefaultText(EGDT_MSG_BOX_YES));
		YesButton->grab();

		btnRect.LowerRightCorner.X += buttonWidth + buttonDistance;
		btnRect.UpperLeftCorner.X += buttonWidth + buttonDistance;
	}

//.........这里部分代码省略.........

// ----------------------------------------------------------------------
int Dense3d::setup(map<string,string>&)
{
  iA(_srcPos!=NULL && _srcNor!=NULL && _trgPos!=NULL);
  iA((*_srcPos).m()==dim() && (*_trgPos).m()==dim());  //nothing to do
  return 0;
}

bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
{
    scalarField dphi(phiq-phi());
    bool isMechRedActive = chemistry_.mechRed()->active();
    label dim(0);
    if (isMechRedActive)
    {
        dim = nActiveSpecies_;
    }
    else
    {
        dim = completeSpaceSize() - nAdditionalEqns_;
    }

    scalar epsTemp = 0;
    List<scalar> propEps(completeSpaceSize(), scalar(0));

    for (label i=0; i<completeSpaceSize()-nAdditionalEqns_; i++)
    {
        scalar temp = 0;

        // When mechanism reduction is inactive OR on active species multiply L
        // by dphi to get the distance in the active species direction else (for
        // inactive species), just multiply the diagonal element and dphi
        if
        (
            !(isMechRedActive)
          ||(isMechRedActive && completeToSimplifiedIndex_[i] != -1)
        )
        {
            label si=(isMechRedActive) ? completeToSimplifiedIndex_[i] : i;

            for (label j=si; j<dim; j++)// LT is upper triangular
            {
                label sj=(isMechRedActive) ? simplifiedToCompleteIndex_[j] : j;
                temp += LT_(si, j)*dphi[sj];
            }

            temp += LT_(si, dim)*dphi[idT_];
            temp += LT_(si, dim+1)*dphi[idp_];
            if (variableTimeStep())
            {
                temp += LT_(si, dim+2)*dphi[iddeltaT_];
            }
        }
        else
        {
            temp = dphi[i]/(tolerance_*scaleFactor_[i]);
        }

        epsTemp += sqr(temp);

        if (printProportion_)
        {
            propEps[i] = temp;
        }
    }

    // Temperature
    if (variableTimeStep())
    {
        epsTemp +=
            sqr
            (
                LT_(dim, dim)*dphi[idT_]
               +LT_(dim, dim+1)*dphi[idp_]
               +LT_(dim, dim+2)*dphi[iddeltaT_]
            );
    }
    else
    {
        epsTemp +=
            sqr
            (
                LT_(dim, dim)*dphi[idT_]
               +LT_(dim, dim+1)*dphi[idp_]
            );
    }

    // Pressure
    if (variableTimeStep())
    {
        epsTemp +=
            sqr
            (
                LT_(dim+1, dim+1)*dphi[idp_]
               +LT_(dim+1, dim+2)*dphi[iddeltaT_]
            );
    }
    else
    {
        epsTemp += sqr(LT_(dim+1, dim+1)*dphi[idp_]);
    }

    if (variableTimeStep())
    {
        epsTemp += sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
    }

    if (printProportion_)
//.........这里部分代码省略.........

HashGrid3<Vector3> Docking::createHashGrid(){
	
	
	float min_x = 1000;
	float min_y = 1000;
	float min_z = 1000;

	float max_x = -1000;
	float max_y = -1000;
	float max_z = -1000;
	
	//iterate over all atoms of protein A to calc dimesion (s. calcSASSurface)
	for (vector<Vector3>::iterator a_it = positionsA.begin(); a_it != positionsA.end(); a_it++){
			
		Vector3 vec = *a_it;
					
		if (vec.x > max_x){
			max_x = vec.x;
		}
		if (vec.x < min_x){
			min_x = vec.x;
		}

		if (vec.y > max_y){
			max_y = vec.y;
		}

		if (vec.y < min_y){
			min_y = vec.y;
		}

		if (vec.z > max_z){
				max_z = vec.z;
		}

		if (vec.z < min_z){
			min_z = vec.z;
		}


	}
	
	float x = max_x - min_x;
	float y = max_y - min_y;
	float z = max_z - min_z;


	x = ceil(x);
	y = ceil(y);
	z = ceil(z);
	

	Vector3 dim(x, y, z);
	Vector3 origin(min_x, min_y, min_z);

	HashGrid3<Vector3> pointGrid(origin,dim, 4.0);

	Vector3 vec2;


	for (vector<Vector3>::iterator A_it= positionsA.begin(); A_it != positionsA.end(); ++A_it){
		
		vec2 = *A_it;
		pointGrid.insert(vec2,vec2);

	}
	
	return pointGrid;

	
}

void CameraTestTool::leftButtonDrag(const TPointD &pos, const TMouseEvent &e) {
  m_dragged = true;
  CleanupParameters *cp =
      CleanupSettingsModel::instance()->getCurrentParameters();

  TPointD dp(pos - m_lastPos);
  if (m_scaling != eNoScale) {
    TDimensionD dim(cp->m_camera.getSize());
    TPointD delta;

    // Mid-edge cases
    if (m_scaling == e1M)
      delta = TPointD(dp.x / Stage::inch, 0);
    else if (m_scaling == e0M)
      delta = TPointD(-dp.x / Stage::inch, 0);
    else if (m_scaling == eM1)
      delta = TPointD(0, dp.y / Stage::inch);
    else if (m_scaling == eM0)
      delta = TPointD(0, -dp.y / Stage::inch);
    else {
      // Corner cases
      if (e.isShiftPressed()) {
        // Free adjust
        if (m_scaling == e11)
          delta = TPointD(dp.x / Stage::inch, dp.y / Stage::inch);
        else if (m_scaling == e00)
          delta = TPointD(-dp.x / Stage::inch, -dp.y / Stage::inch);
        else if (m_scaling == e10)
          delta = TPointD(dp.x / Stage::inch, -dp.y / Stage::inch);
        else if (m_scaling == e01)
          delta = TPointD(-dp.x / Stage::inch, dp.y / Stage::inch);
      } else {
        // A/R conservative
        bool xMaximalDp = (fabs(dp.x) > fabs(dp.y));

        if (m_scaling == e11)
          delta.x = (xMaximalDp ? dp.x : dp.y) / Stage::inch;
        else if (m_scaling == e00)
          delta.x = (xMaximalDp ? -dp.x : -dp.y) / Stage::inch;
        else if (m_scaling == e10)
          delta.x = (xMaximalDp ? dp.x : -dp.y) / Stage::inch;
        else if (m_scaling == e01)
          delta.x = (xMaximalDp ? -dp.x : dp.y) / Stage::inch;

        // Keep A/R
        delta.y = delta.x * dim.ly / dim.lx;
      }
    }

    TDimensionD newDim(dim.lx + 2.0 * delta.x, dim.ly + 2.0 * delta.y);
    if (newDim.lx < 2.0 || newDim.ly < 2.0) return;

    cp->m_camera.setSize(newDim,
                         true,    // Preserve DPI
                         false);  // A/R imposed above in corner cases
  } else {
    if (e.isShiftPressed()) {
      TPointD delta = pos - m_firstPos;
      cp->m_offx    = m_firstCameraOffset.x;
      cp->m_offy    = m_firstCameraOffset.y;
      if (fabs(delta.x) > fabs(delta.y)) {
        if (!cp->m_offx_lock) cp->m_offx += -delta.x * (2.0 / Stage::inch);
      } else {
        if (!cp->m_offy_lock) cp->m_offy += -delta.y * (2.0 / Stage::inch);
      }
    } else {
      if (!cp->m_offx_lock) cp->m_offx += -dp.x * (2.0 / Stage::inch);
      if (!cp->m_offy_lock) cp->m_offy += -dp.y * (2.0 / Stage::inch);
    }
  }

  m_lastPos = pos;

  CleanupSettingsModel::instance()->commitChanges();
  invalidate();
}

// Remapping function for the matrix keyboard
static
UINT
WINAPI
MatrixUsRemapVKey(
    const KEYBD_EVENT *pKbdEvents,
    UINT               cKbdEvents,
    KEYBD_EVENT       *pRmpKbdEvents,
    UINT               cMaxRmpKbdEvents
    )
{
    SETFNAME(_T("MatrixUsRemapVKey"));

    static BOOL fFnDown = FALSE;

    UINT cRmpKbdEvents = 0;
    UINT ui;

    if (pRmpKbdEvents == NULL)
    {
        // 1 to 1 mapping
        if (cMaxRmpKbdEvents != 0)
        {
            DEBUGMSG(ZONE_ERROR, (_T("%s: cMaxRmpKbdEvents error!\r\n"), pszFname));
        }
        return cKbdEvents;
    }

    if (pKbdEvents == NULL)
    {
        DEBUGMSG(ZONE_ERROR, (_T("%s: pKbdEvents error!\r\n"), pszFname));
    }

    if (cMaxRmpKbdEvents < cKbdEvents)
    {
        DEBUGMSG(ZONE_ERROR, (_T("%s: Buffer is not large enough!\r\n"), pszFname));
        return 0;
    }

    for (ui = 0; ui < cKbdEvents; ++ui)
    {
        const KEYBD_EVENT *pKbdEventCurr = &pKbdEvents[ui];
        KEYBD_EVENT *pKbdEventRmpCu