void netcdf_write_blitz(NcVar *nc_var, blitz::Array<T, rank> const &val)
{
	long counts[rank];
	for (int i=0; i<rank; ++i) counts[i] = val.extent(i);
//printf("netcdf_write_blitz: %p %p\n", nc_var, val.data());
	nc_var->put(val.data(), counts);
}

void GridDomain::global_to_local(
	blitz::Array<double,1> const &global,
	std::vector<blitz::Array<double,1>> &olocal)
{
	if (olocal.size() != this->num_local_indices) {
		fprintf(stderr, "MatrixDomainer::get_rows() had bad dimension 1 = %d (expected %ld)\n", olocal.extent(1), this->num_local_indices);
		throw std::exception();
	}

	for (auto ii = olocal.begin(); ii != olocal.end(); ++ii) {
		// Make sure it has the right dimensions
		if (olocal[i].extent(0) != global.extent(0)) {
			fprintf(stderr, "MatrixDomainer::get_rows() had bad dimension 0 = %d (expected %ld)\n", olocal.extent(0), global.extent(0));
			throw std::exception();
		}
	}


	// Copy out data, translating to local coordinates
	for (int j=0; j < global.extent(0); ++j) {
		int lindex[this->num_local_indices];
		this->global_to_local(global(j), lindex);
		for (int i=0; i<this->num_local_indices; ++i)
			olocal[i](j) = lindex[i];
	}
}

void GmresWrapper<Rank>::ApplyOperator(blitz::Array<cplx, 1> &input, blitz::Array<cplx, 1> &output)
{
	if (Psi == 0)
	{
		throw std::runtime_error("Psi is 0");
	}
	if (TempPsi == 0)
	{
		throw std::runtime_error("TempPsi is 0");
	}

	//Map the 1d vectors to a a blitz array of correct shape
	DataVector shape = Psi->GetData().shape();
	DataVector stride = Psi->GetData().stride();
	DataArray inData(input.data(), shape, stride, blitz::neverDeleteData);
	DataArray inData2(input.data(), shape, stride, blitz::neverDeleteData);
	DataArray outData(output.data(), shape, stride, blitz::neverDeleteData);
	DataArray outData2(output.data(), shape, stride, blitz::neverDeleteData);
	outData = 0;

	//Set psi and tempPsi to point to correct vectorsbuffers
	DataArray oldData = Psi->GetData();
	DataArray oldTempData = TempPsi->GetData();
	Psi->SetData(inData);
	TempPsi->SetData(outData);

	OperatorCallback(Psi, TempPsi);

	//Restore the former buffers
	Psi->SetData(oldData);
	TempPsi->SetData(oldTempData);
}

double bob::learn::em::GMMMachine::logLikelihood(const blitz::Array<double, 1> &x,
  blitz::Array<double,1> &log_weighted_gaussian_likelihoods) const
{
  // Check dimension
  bob::core::array::assertSameDimensionLength(log_weighted_gaussian_likelihoods.extent(0), m_n_gaussians);
  bob::core::array::assertSameDimensionLength(x.extent(0), m_n_inputs);
  return logLikelihood_(x,log_weighted_gaussian_likelihoods);
}

inline double Inverse_Matrix_with_lapack_and_Return_Determinant (blitz::Array<double,2> & M, std::string message="")
{  
  assert (M.rows() ==M.cols());
  int n = M.rows();
  if (n==0) return 1;
  ref<double,2> refM(M,false);//do not check the ordering since inverse and transposition commutes
  double * p = refM;
  return Inverse_Matrix_with_lapack_and_Return_Determinant(p,n,n,message);
}

void IceModel_Decode::run_timestep(double time_s,
	blitz::Array<int,1> const &indices,
	std::map<IceField, blitz::Array<double,1>> const &vals2)
{
printf("BEGIN IceModel_Decode::run_timestep(%f) size=%ld\n", time_s, indices.size());
	std::map<IceField, blitz::Array<double,1>> vals2d;	/// Decoded fields

	// Loop through the fields we require
	std::set<IceField> fields;
	get_required_fields(fields);
	for (auto field = fields.begin(); field != fields.end(); ++field) {
printf("Looking for required field %s\n", field->str());

		// Look up the field we require
		auto ii = vals2.find(*field);
		if (ii == vals2.end()) {
			fprintf(stderr, "Cannot find required ice field = %s\n", field->str());
			throw std::exception();
		}
		blitz::Array<double,1> vals(ii->second);

		// Decode the field!
		blitz::Array<double,1> valsd(ndata());
		valsd = nan;
		int n = indices.size();
		for (int i=0; i < n; ++i) {
			int ix = indices(i);
			// Do our own bounds checking!
			if (ix < 0 || ix >= ndata()) {
				fprintf(stderr, "IceModel: index %d out of range [0, %d)\n", ix, ndata());
				throw std::exception();
			}

#if 0
			// Sanity check for NaN coming through
			if (std::isnan(vals(i))) {
				fprintf(stderr, "IceModel::decode: vals[%d] (index=%d) is NaN!\n", i, ix);
				throw std::exception();
			}
#endif

			// Add this value to existing field
			double &oval = valsd(ix);
			if (std::isnan(oval)) oval = vals(i);
			else oval += vals(i);
		}

		// Store decoded field in our output
		vals2d.insert(std::make_pair(*field, valsd));
printf("Done decoding required field, %s\n", field->str());
	}

	// Pass decoded fields on to subclass
	run_decoded(time_s, vals2d);
printf("END IceModel_Decode::run_timestep(%ld)\n", time_s);
}

blitz::Array<double,1> bob::example::library::reverse (const blitz::Array<double,1>& array){
  // create new array in the desired shape
  blitz::Array<double,1> retval(array.shape());
  // copy data
  for (int i = 0, j = array.extent(0)-1; i < array.extent(0); ++i, --j){
    retval(j) = array(i);
  }
  // return the copied data
  return retval;
}

  void Machine::setInputDivision (const blitz::Array<double,1>& v) {

    if (m_weight.extent(0) != v.extent(0)) {
      boost::format m("mismatch on the input division shape: expected a vector of size %d, but you input one with size = %d instead");
      m % m_weight.extent(0) % v.extent(0);
      throw std::runtime_error(m.str());
    }
    m_input_div.reference(bob::core::array::ccopy(v));

  }

  void Machine::setBiases (const blitz::Array<double,1>& bias) {

    if (m_weight.extent(1) != bias.extent(0)) {
      boost::format m("mismatch on the bias shape: expected a vector of size %d, but you input one with size = %d instead");
      m % m_weight.extent(1) % bias.extent(0);
      throw std::runtime_error(m.str());
    }
    m_bias.reference(bob::core::array::ccopy(bias));

  }

Mat naiveBlitzToCvMat(blitz::Array<float, 2> a, float multiplier){
	int i, j;
	Mat b = Mat(a.rows(), a.cols(), CV_32FC1);
	for (i=0; i< a.rows(); i++){
		for(j=0; j<a.cols(); j++){
			b.at<float>(i,j) = a(i, j) * multiplier;
		}
	}
	return b;
}

void binary_save(const blitz::Array<T,n> & M,std::string filename)
{
  T z;assert(M.isStorageContiguous());
  std::stringstream f;M.dumpStructureInformation(f);
  std::ofstream out(filename.c_str(),std::ios::binary);
  std::string s = f.str();int i=int(s.length())+1;char c='1';
  out.write( (char*) &i,sizeof(i));
  out.write( (s.c_str()) ,i*sizeof(c));  
  out.write( (char*)M.dataFirst(),M.numElements()*sizeof(z));
}

void convert2rgb( blitz::Array<Imath::Color3f,2>& img)
{
vigra::BasicImageView<vigra::TinyVector<float,3> > img_view( (vigra::TinyVector<float,3> *) img.dataFirst(), img.cols(), img.rows());

	for( int j=0;j<img.rows();++j)
		for( int i=0;i<img.cols();++i)
			img( j, i) *= lab_scale;

	vigra::transformImage( vigra::srcImageRange( img_view), vigra::destImage( img_view), vigra::Lab2RGBFunctor<float>( 1.0f));
}

  inline void for_cons(const blitz::Array<T,n> & M, bool check_order){
    need_copy = (!(M.isStorageContiguous()));
    if (check_order) for (int i=0; i<n;i++) need_copy = (need_copy || (M.ordering(i)!=i));
#ifdef DEBUG_REF_WARNING
    if (need_copy) std::cout<<"WARNING : REF : COPY NEEDED. Performance will be degraded"<<std::endl;
#endif
    Mref = (blitz::Array<T,n> *)&M;
    // The copy has the same shape but is ordered like a fortran array
    if (need_copy) {Mcopy.resize(M.shape());Mcopy=M;}
  }

void binary_load(blitz::Array<T,n> & M,std::string filename)
{
  assert(M.isStorageContiguous());T z;
  std::ifstream out(filename.c_str(),std::ios::binary);
  int i;char c='1';
  out.read( (char*) &i,sizeof(i));
  char *st = new char[i+1];
  out.read( st ,i*sizeof(c)); std::string s(st); 
  std::stringstream f;  M.dumpStructureInformation(f);
  if (f.str() != s) FATAL("Can not load binary : array do not conform. Structure (file vs array)"<<s<<"----"<<f.str());
  out.read( (char*)M.dataFirst(),M.numElements()*sizeof(z));
}

 Machine::Machine(const blitz::Array<double,2>& weight)
   : m_input_sub(weight.extent(0)),
   m_input_div(weight.extent(0)),
   m_bias(weight.extent(1)),
   m_activation(boost::make_shared<bob::learn::activation::IdentityActivation>()),
   m_buffer(weight.extent(0))
 {
   m_input_sub = 0.0;
   m_input_div = 1.0;
   m_bias = 0.0;
   m_weight.reference(bob::core::array::ccopy(weight));
 }

double bob::learn::em::GMMMachine::logLikelihood(const blitz::Array<double, 2> &x) const {
  // Check dimension
  bob::core::array::assertSameDimensionLength(x.extent(1), m_n_inputs);
  // Call the other logLikelihood_ (overloaded) function


  double sum_ll = 0;
  for (int i=0; i<x.extent(0); i++)
    sum_ll+= logLikelihood_(x(i,blitz::Range::all()));

  return sum_ll/x.extent(0);  
}

void check_dimensions(
    std::string const &vname,
    blitz::Array<T, rank> const &arr,
    std::vector<int> const &dims)
{
    for (int i=0; i<rank; ++i) {
        if (dims[i] >= 0 && arr.extent(i) != dims[i]) {
            fprintf(stderr,
                    "Error in %s: expected dimension #%d = %d (is %d instead)\n",
                    vname.c_str(), i, dims[i], arr.extent(i));
            throw std::exception();
        }
    }
}

  void Machine::forward (const blitz::Array<double,1>& input, blitz::Array<double,1>& output) const {

    if (m_weight.extent(0) != input.extent(0)) { //checks input dimension
      boost::format m("mismatch on the input dimension: expected a vector of size %d, but you input one with size = %d instead");
      m % m_weight.extent(0) % input.extent(0);
      throw std::runtime_error(m.str());
    }
    if (m_weight.extent(1) != output.extent(0)) { //checks output dimension
      boost::format m("mismatch on the output dimension: expected a vector of size %d, but you input one with size = %d instead");
      m % m_weight.extent(1) % output.extent(0);
      throw std::runtime_error(m.str());
    }
    forward_(input, output);

  }

  void Machine::setWeights (const blitz::Array<double,2>& weight) {

    if (weight.extent(0) != m_input_sub.extent(0)) { //checks 1st dimension
      boost::format m("mismatch on the weight shape (number of rows): expected a weight matrix with %d row(s), but you input one with %d row(s) instead");
      m % m_input_sub.extent(0) % weight.extent(0);
      throw std::runtime_error(m.str());
    }
    if (weight.extent(1) != m_bias.extent(0)) { //checks 2nd dimension
      boost::format m("mismatch on the weight shape (number of columns): expected a weight matrix with %d column(s), but you input one with %d column(s) instead");
      m % m_bias.extent(0) % weight.extent(1);
      throw std::runtime_error(m.str());
    }
    m_weight.reference(bob::core::array::ccopy(weight));

  }

void bob::learn::em::GMMMachine::setMeans(const blitz::Array<double,2> &means) {
  bob::core::array::assertSameDimensionLength(means.extent(0), m_n_gaussians);
  bob::core::array::assertSameDimensionLength(means.extent(1), m_n_inputs);
  for(size_t i=0; i<m_n_gaussians; ++i)
    m_gaussians[i]->updateMean() = means(i,blitz::Range::all());
  m_cache_supervector = false;
}