int counts(struct node *p)
{
				int l,r;
				if(p) {
								 l = 1 + counts(p -> left);
								 r = 1 + counts(p -> right);
								 return l>r?l:r;
				} else
					 return 0;
}

size_t getImportList(
   const PartitioningSolution<SolutionAdapter> &solution,
   const DataAdapter * const data,
   ArrayRCP<typename DataAdapter::zgid_t> &imports // output
)
{
  typedef typename PartitioningSolution<SolutionAdapter>::part_t part_t;
  typedef typename PartitioningSolution<SolutionAdapter>::gno_t gno_t;
  typedef typename DataAdapter::zgid_t zgid_t;

  size_t numParts = solution.getActualGlobalNumberOfParts();
  int numProcs = solution.getCommunicator()->getSize();

  if (numParts > size_t(numProcs)) {
    char msg[256];
    sprintf(msg, "Number of parts %lu exceeds number of ranks %d; "
                 "%s requires a MappingSolution for this case\n",
                  numParts, numProcs, __func__);
    throw std::logic_error(msg);
  }

  size_t localNumIds = data->getLocalNumIDs();
  const zgid_t *gids = NULL;
  data->getIDsView(gids);

  const part_t *parts = solution.getPartListView();

  // How many ids to each process?
  Array<int> counts(numProcs, 0);
  for (size_t i=0; i < localNumIds; i++)
    counts[parts[i]]++;

  Array<gno_t> offsets(numProcs+1, 0);
  for (int i=1; i <= numProcs; i++){
    offsets[i] = offsets[i-1] + counts[i-1];
  }

  Array<typename DataAdapter::zgid_t> gidList(localNumIds);
  for (size_t i=0; i < localNumIds; i++) {
    gno_t idx = offsets[parts[i]];
    gidList[idx] = gids[i];
    offsets[parts[i]] = idx + 1;
  }

  Array<int> recvCounts(numProcs, 0);
  try {
    AlltoAllv<zgid_t>(*(solution.getCommunicator()),
                      *(solution.getEnvironment()),
                      gidList(), counts(), imports, recvCounts());
  }
  Z2_FORWARD_EXCEPTIONS;

  return imports.size();
}

int counts(int tot, int* coins, int highest) {
	calls++;

	if( tot <  0) return 0; 
	if( tot == 0) return 1;
	
	int sum = counts(tot - coins[highest], coins, highest );
	
	if(highest > 0)
			sum += counts(tot, coins, highest-1 );

	return sum;
}

        void score_value(const Model & model, AlignedFloats scores) const {
            if (DIST_DEBUG_LEVEL >= 1) {
                DIST_ASSERT_EQ(scores.size(), counts().size());
            }

            const size_t size = counts().size();
            const float shift = -fast_log(sample_size() + model.alpha);
            const float * __restrict__ in = VectorFloat_data(shifted_scores_);
            float * __restrict__ out = VectorFloat_data(scores);

            for (size_t i = 0; i < size; ++i) {
                out[i] = in[i] + shift;
            }
        }

    void ArithmeticUtilDecoder::decode(std::vector<byte>& data) {
        long p = in->pos;
        uint64 length = in->read48bits();
        bits_left = in->read48bits()*8;
        //       std::cout<<"Current position: "<<in->pos<<" - "<<" bits to read: "<<bits_left<<"\n";
        data.resize(length);
        std::vector<uint64> counts(256);
        std::vector<SYMBOL> symbols(256);
        utils::gammaDecode(counts,in);
        in->flushBuffer();
        uint64 cumul=0;
        for(int i=0;i<256;i++) {
            symbols[i].low_count=cumul;
            symbols[i].high_count = cumul = cumul + counts[i];
        }
        byte aa = in->readByte();
        byte bb = in->readByte();
        byte cc = in->readByte();
        byte dd = in->readByte();
        bits_left-=32;
        code = (aa<<24)|(bb<<16)|(cc<<8)|dd;
        for(int i=0;i<length;i++) {
            uint64 v = get_value();
            byte b = get_byte(symbols,v);
            if(i+1<length)    del_symbol(symbols[b]);
            data[i]=b;
        }
        in->flushBuffer();

    }

int maxneb_disc () {
	Rcpp::IntegerVector targs(targets);
	Rcpp::IntegerVector nebs(neighbors);
	Rcpp::IntegerVector plot(plots);
	Rcpp::CharacterVector status(stat);
	Rcpp::IntegerVector xx(bqudx);
	Rcpp::IntegerVector yy(bqudy);
	Rcpp::IntegerVector tt(time);
	int rad = Rcpp::as<int>(sr);
	Rcpp::IntegerVector counts(targs.size());
	int max_neb = 0;

	for (int i = 0, n = targs.size(); i < n; i++) {
		int targ = index(targs[i], nebs);
		int num_nebs= 0;
		for (int neb = 0, n2 = nebs.size(); neb < n2; neb++) {
			if (targ != neb && plot[targ] == plot[neb] &&
				(std::strcmp("ALIVE", status[neb]) == 0) &&
				(xx[targ] + sr > xx[neb]) && (xx[targ] - sr < xx[neb]) &&
				(yy[targ] + sr > yy[neb]) && (yy[targ] - sr < yy[neb]) &&
				tt[targ] == tt[neb])
				num_nebs++;
		}
		counts[i] = num_nebs;
	}

	return std::max_element(counts.begin(), counts.end(), myfn);

}

// bootstrapping
std::vector<float> Classifier::getBootstrap(const KmerSequence& querySeq, RandomGen &generator, const searchHit& hit)
{
    KmerSequence bootstrap;
    std::vector<kmerSize_t>kmers = querySeq.kmers;
    std::vector<float> counts(referenceData_.numLevels(), 0.f);
    int bootstrap_size = kmers.size() / subsampleSize_;
    
    for (unsigned int bs=0; bs<numBootstrap_; bs++)
    {
        std::random_shuffle(kmers.begin(), kmers.end(), generator);
        bootstrap.kmers = std::vector<kmerSize_t>(kmers.begin(), kmers.begin() + bootstrap_size);
        searchHit bsHit = referenceData_.search(bootstrap);
        
        for (unsigned int i=0; i<referenceData_.numLevels(); i++)
        {
            if (hit.annotationIds[i].size() == 1)
            {
                if(std::find(bsHit.annotationIds[i].begin(), bsHit.annotationIds[i].end(), *hit.annotationIds[i].begin()) != bsHit.annotationIds[i].end())
                {
                    counts[i] += static_cast<float>(std::count(bsHit.annotationIds[i].begin(), bsHit.annotationIds[i].end(), *hit.annotationIds[i].begin())) / static_cast<float>(bsHit.annotationIds[i].size());
                }
            }
        }
    }

    std::vector<float> retVec;
    for (unsigned int i=0; i<referenceData_.numLevels(); i++)
    {
        retVec.push_back(counts[i] / static_cast<float>(numBootstrap_));
    }

    return retVec;
}

TGraph* 
makeGraph(const TArrayI& adcs, Int_t rate) 
{
  Int_t    last = adcs.fArray[0];
  TArrayI  counts(4);
  TGraph*  graph = new TGraph(rate * adcs.fN);
  graph->SetLineColor(rate);
  graph->SetMarkerColor(rate);
  graph->SetMarkerStyle(20+rate);
  graph->SetLineStyle(rate);
  graph->SetName(Form("rate%d", rate));
  graph->SetTitle(Form("Rate %d", rate));
  for (Int_t i = 0; i < adcs.fN; i++) { 
    counts.Reset(-1);
    convert(rate, adcs.fArray[i], last, counts);
    
    for (Int_t j = 0; j < rate; j++) { 
      Int_t    idx = (i * rate + j);
      Double_t x   = (i + (rate > 1 ? Float_t(j+1) / rate-1 : 0));
      graph->SetPoint(idx, x, counts[j]);
    }
    last = counts[rate - 1];
  }
  return graph;
}

char MultiSequence::GetAnnotationChar(const std::vector<char> &column) const
{
    std::vector<int> counts(256, 0);
    char annot[] = ":::::::::...........";
    const char *groups[20] = { "STA", "NEQK", "NHQK", "NDEQ", "QHRK",
                               "MILV", "MILF", "HY", "FYW", "CSA", 
                               "ATV", "SAG", "STNK", "STPA", "SGND", 
                               "SNDEQK", "NDEQHK", "NEHQRK", "FVLIM", "HFY" };
    
    for (size_t i = 0; i < column.size(); i++)
        counts[BYTE(toupper(column[i]))]++;
    
    for (int i = 0; i < 256; i++)
        if (char(i) != '-' && counts[i] == int(column.size()))
            return '*';
    
    for (int i = 0; i < 20; i++)
    {
        int total_counts = 0;
        for (size_t j = 0; j < strlen(groups[i]); j++) 
            total_counts += counts[BYTE(groups[i][j])];
        if (total_counts == int(column.size()))
            return annot[i];
    }

    return ' ';
}

void threads::output() {
    set_cpu(1);

    /*fill the indexes*/
    std::vector<size_t> indexes(blocks.size());
    for (size_t i = 0; i < blocks.size(); i++) { indexes[i] = i; }

    /*initialize stored last-output counts*/
    std::vector<size_t> counts(blocks.size(), 0);

#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wmissing-noreturn"
    while (true) {

        /*mix around the indexes*/
        std::shuffle(indexes.begin(), indexes.end(), primary_rng);

        for (size_t i = 0; i < blocks.size(); i++) {

            /*if the block is new*/
            if (blocks[indexes[i]].count > counts[indexes[i]]) {

                /*update the block count*/
                counts[indexes[i]] = blocks[indexes[i]].count;
                output_stream.write((const char *) blocks[indexes[i]].data,
                                    blocksize);
                output_stream.flush();

            }
        }
    }
#pragma clang diagnostic pop
}

    static void gather(const communicator& comm, const std::vector<T>& in, std::vector< std::vector<T> >& out, int root)
    {
      std::vector<int>  counts(comm.size());
      Collectives<int,void*>::gather(comm, (int) in.size(), counts, root);

      std::vector<int>  offsets(comm.size(), 0);
      for (unsigned i = 1; i < offsets.size(); ++i)
        offsets[i] = offsets[i-1] + counts[i-1];

      std::vector<T> buffer(offsets.back() + counts.back());
      MPI_Gatherv(Datatype::address(const_cast<T&>(in[0])),
                  in.size(),
                  Datatype::datatype(),
                  Datatype::address(buffer[0]),
                  &counts[0],
                  &offsets[0],
                  Datatype::datatype(),
                  root, comm);

      out.resize(comm.size());
      size_t cur = 0;
      for (unsigned i = 0; i < (unsigned)comm.size(); ++i)
      {
          out[i].reserve(counts[i]);
          for (unsigned j = 0; j < (unsigned)counts[i]; ++j)
              out[i].push_back(buffer[cur++]);
      }
    }

 Vector DPMM::allocation_counts() const {
   Vector counts(mixture_components_.size());
   for (int i = 0; i < mixture_components_.size(); ++i) {
     counts[i] = mixture_components_[i]->suf()->n();
   }
   return counts;
 }

BOOST_AUTO_TEST_CASE_TEMPLATE(GoodnessOfFit, RandomGenerator, RandomGenerators)
{
  // Kolmogorov-Smirnov Goodness-of-Fit Test
  // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35g.htm

  const size_t MAX_BINS = 32;
  const uint32_t MAX_ITERATIONS = 35;

  std::vector<uint32_t> counts(MAX_BINS, 0);

  for (uint32_t i = 0; i < MAX_ITERATIONS; i++) {
    counts[RandomGenerator::generate() % MAX_BINS]++;
  }

  std::vector<double> edf(MAX_BINS, 0.0);
  double probability = 0.0;
  for (size_t i = 0; i < MAX_BINS; i++) {
    probability += 1.0 * counts[i] / MAX_ITERATIONS;
    edf[i] = probability;
  }

  double t = 0.0;
  for (size_t i = 0; i < MAX_BINS; i++) {
    t = std::max(t, std::abs(edf[i] - (i * 1.0 / MAX_BINS)));
  }

  // Check if it is uniform distribution with confidence 0.95
  // http://dlc.erieri.com/onlinetextbook/index.cfm?fuseaction=textbook.appendix&FileName=Table7
  BOOST_WARN_LE(t, 0.230);
}

	std::vector<std::pair<double, double>> calculateProbabilityHist(const std::vector<int>& data, int bin_size) {
		// Check for valid input data
		if ((int)data.size() == 0) {
			cout << "Error! Cannot calculate probability histogram because the input data vector is empty." << endl;
			throw invalid_argument("Error! Cannot calculate probability histogram because the input data vector is empty.");
		}
		// Determine the starting bin position
		int min_val = *min_element(data.begin(), data.end());
		int max_val = *max_element(data.begin(), data.end());
		// Determine number of bins
		int num_bins = (int)((double)(max_val - min_val + 1) / (double)bin_size);
		// Calculate bins
		vector<pair<double, double>> hist(num_bins, make_pair(0.0, 0.0));
		for (int i = 0; i < num_bins; i++) {
			hist[i].first = min_val + 0.5*(bin_size - 1) + (bin_size - 1) * i;
		}
		// Calculate histogram
		vector<int> counts(num_bins, 0);
		int index;
		for (int i = 0; i < (int)data.size(); i++) {
			index = (data[i] - min_val) / bin_size;
			counts[index]++;
		}
		// total counts
		int total_counts = accumulate(counts.begin(), counts.end(), 0);
		// Normalized histogram to get probability
		for (int i = 0; i < num_bins; i++) {
			hist[i].second = (double)counts[i] / (double)(total_counts);
		}
		return hist;
	}

  int maximalRectangle(vector<vector<char> > &matrix) {
    if (0 == matrix.size()) {
      return 0;
    }

    int result = INT_MIN;
    int height = matrix.size();
    int width = matrix[0].size();
    int row;
    int col;
    vector<vector<int> > counts(height, vector<int>(width, 0));

    for (row = 0; row < height; ++row) {
      for (col = 0; col < width; ++col) {
        if (matrix[row][col] != '0') {
          if (0 == row) {
            counts[row][col] = 1;
          }
          else {
            // 每一行往上数，包括自己有几个1。
            counts[row][col] = counts[row - 1][col] + 1;
          }
        }
      }
    }

    for (row = 0; row < height; ++row) {
      result = max(result, maximal_rectangle(counts[row]));
    }

    return result;
  }

      void readTokenMatrix(Dataset& data, const char* filename) {
	std::vector<string> lines;
	loadFile(lines, filename);

	assert(lines.size() > 1);
	std::vector<int> counts(2);
	tokenizeLineIntoIntVector(lines[1], counts);
	assert(counts.size() == 2);
	//cout << "Found: " << counts[0] << " and " << counts[1] << endl;
	data.numDocs = counts[0];
	data.numTokens = counts[1];

	vec trainClasses(data.numDocs);
	mat features(data.numDocs, data.numTokens);
	rowvec row(data.numTokens);

	for(uint i = 3; i < data.numDocs + 3; i++) {
	  std::vector<double> tokens(getNumTokens(lines[i]));
	  tokenizeLineIntoDoubleVector(lines[i], tokens);
	  trainClasses(i-3) = tokens[0];
	  //cout << "class[" << tokens[0] << "]";
	  row.zeros();
	  int cumsum = 0;
	  for(uint j = 1; j < tokens.size() - 1; j+=2) {
	    cumsum += tokens[j];
	    //cout << "cumsum[" << cumsum << "]token[" << tokens[j+1] << "]" << flush ;
	    row[cumsum] = tokens[j+1];
	  }
	  Matrix::setMatrixRowToVector(features, i-3 ,row);
	}
	
	data.classifications = trainClasses;
	data.features = features;

      }

    virtual double             error(void) const{
        // Compute uncertainty
        double error = sqrt(counts()+1.0e-50);

        // Return error
        return error;
    }

void SparseMsg::setPattern(UInt num, const UInt *proc) {

  UInt csize = Par::Size();

  // Set dest proc
  nsend = num;

  std::vector<int> sendto(nproc, 0);
  std::vector<int> counts(nproc, 1);
  for (UInt i = 0; i < num; i++) {
    ThrowRequire(proc[i] < csize);
    sendto[proc[i]] = 1;
    if (proc[i] == (UInt) rank) {
      sendself = true;
      self_idx = i;
    }
  }

  !Par::Serial() ? MPI_Reduce_scatter(&sendto[0], &num_incoming, &counts[0], MPI_INT, MPI_SUM, comm)
    : num_incoming = sendto[0];
//std::cout << "Proc:" << rank << "to receive " << num_incoming << " messages" << std::endl;

  // Set up send buffers (so we don't have save the proc ids
  if (nsend > 0) outBuffers.resize(nsend); else outBuffers.clear();
  for (UInt i = 0; i < nsend; i++) {
    outBuffers[i].proc = proc[i];
  }
}

template<class T> static bool load_nc_array(const NcFile& ncf, const string& name, vector<T>& dest, bool required = true, int offset = 0, int count = -1)
{
	NcVar *v = load_nc_variable(ncf, name.c_str(), required);
	if (v)
	{
		vector<long> offsets = list_of(offset).repeat(v->num_dims()-1, 0);
		v->set_cur(&offsets.front());
		vector<long> counts (v->num_dims());
		long* shape = v->edges();
		transform(shape, shape + v->num_dims(), offsets.begin(), counts.begin(), minus<long>());
		delete shape;
		if (count > 0)
		{
			counts[0] = count;
		}
		dest.resize(product(counts));
		bool success = v->get(&dest.front(), &counts.front());
		if (!success)
		{
			dest.resize(0);
			check(!required, string("NetcdfDataset::load_nc_array<") + typeid(T).name() + "> " + name + '\n' + "failed with offset " + str(offsets) + ", counts " + str(counts));
		}
		return success;
	}
	return false;
}

radixSort(std::vector<T> v)
{

    std::vector<T> tmp;

    size_t radix_count = sizeof(T); // we process elements in 256 counting system. Each radix = byte.


    for (size_t radix = 0; radix < radix_count; ++radix)
    {
        std::vector<size_t> counts(radix_count);

        for (auto x : v)
        {
            // get radix
            auto radix = getRadix(x);
            
            // increase radix count at its position


        }

        // calculate cumulative summs
        // write elements from V to TMP according to cum_sums positions. Modify cum_sums after writing an element

        v.swap(tmp);
    }


    

}