bool getentboundingbox(extentity &e, ivec &o, ivec &r)
{
    switch(e.type)
    {
        case ET_EMPTY:
            return false;
        case ET_MAPMODEL:
        {
            model *m = loadmodel(NULL, e.attr[1]);
            if(m)
            {
                vec center, radius;
                m->boundbox(center, radius);
                rotatebb(center, radius, e.attr[0]);

                o = e.o;
                o.add(center);
                r = radius;
                r.add(1);
                o.sub(r);
                r.mul(2);
                break;
            }
        }
        // invisible mapmodels use entselradius
        default:
            o = e.o;
            o.sub(entselradius);
            r.x = r.y = r.z = entselradius*2;
            break;
    }
    return true;
}

std::complex<double> operator*(const ivec &a, const cvec &b)
{
  it_assert_debug(a.size() == b.size(), "operator*(): sizes does not match");
  std::complex<double> temp = 0;
  for (int i = 0;i < a.size();i++) {temp += (double)a(i) * b(i);}
  return temp;
}

bool getentboundingbox(extentity &e, ivec &o, ivec &r)
{
    switch(e.type)
    {
        case ET_EMPTY:
            return false;
        case ET_MAPMODEL:
        {
            LogicEntityPtr entity = LogicSystem::getLogicEntity(e); // INTENSITY
            model *m = entity.get() ? entity->getModel() : NULL; // INTENSITY
            if(m)
            {
                vec center, radius;
                m->boundbox(0, center, radius, entity.get()); // INTENSITY: entity
                rotatebb(center, radius, e.attr1);
                o = e.o;
                o.add(center);
                r = radius;
                r.add(1);
                o.sub(r);
                r.mul(2);
                break;
            }
        }
        // invisible mapmodels use entselradius
        default:
            o = e.o;
            o.sub(GETIV(entselradius));
            r.x = r.y = r.z = GETIV(entselradius)*2;
            break;
    }
    return true;
}

static
void assert_ivec(const ivec &expected, const ivec &actual)
{
  ASSERT_EQ(expected.length(), actual.length());
  for (int n = 0; n < expected.length(); ++n)
  {
    ASSERT_EQ(expected[n], actual[n]);
  }
}

cvec operator+(const std::complex<double> &s, const ivec &v)
{
  it_assert_debug(v.size() > 0, "operator+(): Vector of zero length");

  cvec temp(v.size());
  for (int i = 0;i < v.size();i++) {
    temp(i) = s + std::complex<double>(v(i));
  }
  return temp;
}

	void r_search(Node *root, const interval<T>& i, ivec& vec, int mode) {
		if (root == NULL) return;

		if (mode == OVERLAP && root->m_interval->overlaps(i)) vec.push_back(*(root->m_interval));
		if (mode == CONTAIN && root->m_interval->contains(i)) vec.push_back(*(root->m_interval));
		if (mode == CONTAINED && i.contains(*(root->m_interval))) vec.push_back(*(root->m_interval));

		if (root->left != NULL && i.overlaps(root->left->min, root->left->max))
			r_search(root->left, i, vec, mode);
		if (root->right != NULL && i.overlaps(root->right->min, root->right->max))
			r_search(root->right, i, vec, mode);
	}

void ClassifierMLP::Train(std::vector< fvec > samples, ivec labels)
{
	u32 sampleCnt = samples.size();
	if(!sampleCnt) return;
	DEL(mlp);
	dim = samples[0].size();

	CvMat *layers;
//	if(neuronCount == 3) neuronCount = 2; // don't ask me why but 3 neurons mess up everything...

	if(!layerCount || neuronCount < 2)
	{
		layers = cvCreateMat(2,1,CV_32SC1);
		cvSet1D(layers, 0, cvScalar(dim));
		cvSet1D(layers, 1, cvScalar(1));
	}
	else
	{
		layers = cvCreateMat(2+layerCount,1,CV_32SC1);
		cvSet1D(layers, 0, cvScalar(dim));
		cvSet1D(layers, layerCount+1, cvScalar(1));
		FOR(i, layerCount) cvSet1D(layers, i+1, cvScalar(neuronCount));
	}

	u32 *perm = randPerm(sampleCnt);

	CvMat *trainSamples = cvCreateMat(sampleCnt, dim, CV_32FC1);
	CvMat *trainLabels = cvCreateMat(labels.size(), 1, CV_32FC1);
	CvMat *sampleWeights = cvCreateMat(samples.size(), 1, CV_32FC1);
	FOR(i, sampleCnt)
	{
		FOR(d, dim) cvSetReal2D(trainSamples, i, d, samples[perm[i]][d]);
		cvSet1D(trainLabels, i, cvScalar(labels[perm[i]]));
		cvSet1D(sampleWeights, i, cvScalar(1));
	}

void cofdm_map::set_data(ivec x)
{
cvec qv;
bvec ce;
int K = x.length();
int i;

	ce.set_length(K); ce.ones();
	#if (DEBUG_LEVEL == 3)
	cout << "***** cofdm_map::set_data *****" << endl;	
	cout << "K=" << K << endl;
	cout << "ce=" << ce << endl;
	cout << "x=" << x << endl;
	cout << "data_carriers=" << data_carriers << endl;
	#endif

	if( K == data_carriers.length() ) {
		qv = qammod.process(ce,x);
		#if (DEBUG_LEVEL == 3)
		cout << "qv=" << qv << endl;
		#endif
		for (i=0; i<K; i++) {
			y0(data_carriers(i))=qv(i);
		}
	}
	else {
		throw sci_exception("cofdm_map::set_data - x.size() <> data_carriers.size()=", data_carriers.length());
	}

	#if (DEBUG_LEVEL == 3)
	cout << "y0(piltos)  y0(zeros) migh have rubbish" << endl;
	cout << "y0=" << y0 << endl;
	cout << "+++++ cofdm_map::set_data +++++" << endl;	
	#endif
}

ivec TSP(const ivec& purchases)
{
	if (purchases.size()<2) return ivec(purchases);
	ivec cur = purchases;
	ivec best = cur;
	double bdist = pathCost(cur);
	double cdist = bdist;
	for(double t=10; t>.1; t*=.999) {
		size_t a = rand()%(cur.size()-1);
		size_t b = rand()%(cur.size()-1);
		if (b>=a) ++b;
		else swap(a,b);

		int ca=cur[a], cb=cur[b];
		double da = a ? itemDist[cur[a-1]][cb] -  itemDist[cur[a-1]][ca] : startDist[cb] - startDist[ca];
		double db = b<cur.size()-1 ? itemDist[ca][cur[b+1]] - itemDist[cb][cur[b+1]] : endDist[ca] - endDist[cb];
		double d = da+db;
		if (d<0 || randf()<exp(-d/t)) {
			reverse(cur.begin()+a,cur.begin()+b+1);
			cdist = pathCost(cur);
			if (cdist < bdist) {
				best = cur;
				bdist = cdist;
			}
		}
	}
	return best;
}

int bpsk_sym_err(bvec bv_src, ivec iv_dem) {
	int err = 0;
	for(int i=0; i<iv_dem.size(); i++) {
		if(bv_src(i)!=iv_dem(i))
			err++;
	}
	return err;
}

double expectedTotalCost(const ivec& path)
{
    double r=0;
    for(size_t i=1; i<path.size(); ++i) {
//		r += dist[path[i-1]][path[i]];
        int a = path[i-1], b = path[i];
        int n = lower_bound(conn[a].begin(),conn[a].end(),b)-conn[a].begin();
        assert(conn[a][n]==b);
        r += edgeDist[a][n];
//		if (final) cout<<"lol "<<a<<' '<<b<<' '<<n<<' '<<edgeDist[a][n]<<'\n';
    }
    if (final) cout<<path<<'\n'<<"length cost "<<r<<'\n';
    r *= LENGTH_FACTOR;

    double rr=0;

    vector<double> pdist;
    for(size_t k=0; k<samples.size(); ++k) {
        vector<int>& v = samples[k];
        int n = v.size();
        pdist.resize(n);
        for(int i=0; i<n; ++i) {
            int a = v[i];
            double d=1e100;
            for(size_t j=0; j<path.size(); ++j)
                d = min(d, dist[a][path[j]]);
            pdist[i] = d;
//			if (final) cout<<"pdist "<<k<<' '<<a<<": "<<d<<'\n';
        }

        for(int i=0; i<n; ++i)
            for(int j=0; j<i; ++j) {
                int x=v[i], y=v[j];
//				cout<<"asdasd "<<tsps[a].dists<<'\n';
                tsps[k].dists[i][j] = tsps[k].dists[j][i] = min(dist[x][y], pdist[i]+pdist[j]);
            }

        double cr = tsps[k].calc();
        if (robustOpt) rr = max(rr, cr);
        else rr += probs[k] * cr;

        if (final) cout<<"tsp cost "<<k<<": "<<cr<<'\n';
    }
    return r + rr;
}

void gfx::set(int qvalue, const ivec &invalues)
{
  //    it_assert0(qvalue>0, "gfx::set, out of range");
  degree=invalues.length()-1;
  coeffs.set_size(degree+1);
  for (int i=0;i<degree+1;i++)
	coeffs[i].set(qvalue,invalues[i]);
  q=qvalue;
}

int qam_sym_err(bvec bv_src, ivec iv_dem) {
	int err = 0;
	for(int i=0, k=0; i<iv_dem.size(); i++, k+=2) {
		int src = (bv_src(k).value()<<1) + bv_src(k+1).value();
		if(src!=iv_dem(i))
			err++;
	}
	return err;
}

CFix operator*(const cfixvec &a, const ivec &b)
{
  it_assert_debug(a.size() == b.size(), "operator+(): sizes do not match");
  CFix temp(0);
  for (int i = 0; i < a.size(); i++) {
    temp += a(i) * b(i);
  }
  return temp;
}

ivec max_to_all(const ivec &pt) {
  int in[5], out[5];
  for (int i=0; i<5; ++i) in[i] = out[i] = pt.in_direction(direction(i));
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,5,MPI_INT,MPI_MAX,mycomm);
#endif
  ivec ptout(pt.dim);
  for (int i=0; i<5; ++i) ptout.set_direction(direction(i), out[i]);
  return ptout;
}

//gene_idx is for a single gene (otherwise it should be imat rather than ivec)
void extract_score(fvec &S, const fvec &score0, 
const ivec &gene_idx){
  
	int ns = gene_idx.size();//number of SNPs in a gene
		
	S = fvec(ns);
	for(int j = 0; j < ns; ++j){
		int J = gene_idx[j];
		S[j] = score0[J];
	}
	
}

void Triangles::setNodes(ivec n)
{
    int size = n.size();
    nodes = ivec(size);

    for (int i = 0; i < size; i++)
    {
        nodes[i] = n[i];
    }

    return;
}

cvec qam_mod::process(bvec ce, ivec x)
{
 cvec y;	
 int N;
 ivec iv;
 cvec cv;

	#if (DEBUG_LEVEL==3)
	cout << "***** qam_mod::process *****" << endl;	
	cout << "ce=" << ce << endl;
	cout << "x=" << x << endl;
	sleep(1000);
	#endif

	iv.set_length(1);
	cv.set_length(1);
	N=ce.length();
	y.set_length(N);
	if (x.length()!=N) {
		throw sci_exception("qam_mod::process - ce.size <> x.size", x.length() );
	}	
	for (int i=0; i<N; i++) {
		if ( bool(ce[i])) {
			iv[0]  = x[i];
			cv = modulate(iv);
			y0 = scale * cv[0];
		}
		y[i]=y0;
	}

	#if (DEBUG_LEVEL==3)
	cout << "y=" << y << endl;
	cout << "+++++ qam_mod::process +++++" << endl;	
	sleep(1000);
	#endif

	return (y);
}

void match_one (svec a, svec b, ivec& ret) {

  si_map A;
  svec inter;

  for(int i=0; i < a.size(); i++) A[a[i]] = i;

  std::sort(a.begin(), a.end());
  std::sort(b.begin(), b.end());

  std::set_intersection(a.begin(), a.end(), b.begin(), b.end(),
        std::back_inserter(inter));
  ret.resize(inter.size());
  for(int i=0; i < inter.size(); i++) ret[i] = A[inter[i]];
}

void match_both (svec a, svec b, ivec& ai, ivec& bi) {

  si_map A;
  si_map B;
  svec inter;

  int i = 0;

  for(i=0; i < a.size(); i++) A[a[i]] = i;
  for(i=0; i < b.size(); i++) B[b[i]] = i;

  std::sort(a.begin(), a.end());
  std::sort(b.begin(), b.end());

  std::set_intersection(a.begin(), a.end(), b.begin(), b.end(),
        std::back_inserter(inter));

  ai.resize(inter.size());
  for(i=0; i < inter.size(); i++) ai[i] = A[inter[i]];

  bi.resize(inter.size());
  for(i=0; i < inter.size(); i++) bi[i] = B[inter[i]];

}