VECTOR CTempConvs::dotdiv(VECTOR& a, VECTOR& b)
{
	
	///*
	VECTOR div(a.size());
	for (int i = 0; i < (int)a.size(); i++)
	{
		// Armadillo elementwise division gives NaN for 0/0
		if ((a(i) == 0))
			div(i) = 0;
		else if (b(i) == 0)
			div(i) = INF;
		else div(i) = a(i) / b(i);
	}
	//*/

	// Alternatively, the following gives identicle results to Matlab
	/*
	VECTOR div(a.size());
	for (int i = 0; i<(int)a.size(); i++){
		div(i) = 0;
	}
	if (a.size() != b.size())
		return	div;

	for (int i = 0; i<(int)a.size(); i++)
	{
		if (a(i) == NaN)
			div(i) = NaN;
		if (a(i) == INF)
		{
			if (b(i) == INF || b(i) == NaN)
				div(i) = NaN;
			else
				div(i) = INF;
		}
		if (a(i) == 0)
		{
			if (b(i) != 0 || b(i) == INF)
				div(i) = 0;
			else
				div(i) = NaN;
		}
		else
		{
			if (b(i) == 0 || b(i) == NaN)
				div(i) = NaN;
			if (b(i) == INF)
				div(i) = 0;
			else
				div(i) = a(i) / b(i);
		}
	}
	//*/

	//a.print("a:"); b.print("b:"); div.print("div:");

	return div;

}

void System::ROTATE_FRAGMENT(double degree_amount, VECTOR direction,int Gr, VECTOR center){
/**
  \param[in] degree_amount The magnitude of rotation, in degrees
  \param[in] direction The vector definining the axis of rotation. The magnitude of this vector does not matter.
  \param[in] Gr The ID (not index!) of the group/fragment to be rotated
  \param[in] center The vector defining the center of the rotating coordinate system

  Simplest manipulation
  Rotates the fragment with the fragment ID "int Gr" on amount of "double amount"
  around the axis given by "VECTOR direction" around the given center
*/

  int v = get_fragment_index_by_fragment_id(Gr);

  VECTOR dir; dir = center - Fragments[v].Group_RB.rb_cm;
  double amount = dir.length();

  // Translate the fragment's center of mass to the "center" point
  TRANSLATE_FRAGMENT(amount, dir,Gr);

  // Rotate the fragment around new center of mass
  ROTATE_FRAGMENT(degree_amount, direction, Gr);

  // Translate the fragment's center of mass back to the original position
  TRANSLATE_FRAGMENT(-amount, dir,Gr);

}

/*-<==>-----------------------------------------------------------------
/ 
/----------------------------------------------------------------------*/
bool CSphere::hits (const CLine &line, SCALAR &t_hit) {
  // Pendiente de implementar correctamente
	
	// REVISADA

	VECTOR aux = line.loc - loc;
	SCALAR b = aux.dot(line.dir);
	SCALAR c = aux.dot(aux) - radius*radius;
	SCALAR d = b*b - c;
	if (d<0)
			return false; // Noqueremos negativos pq no toca linea con esfera
	else{
		SCALAR tm = -b - sqrt(d); // -b - SQRT(b^2 -4ac)
		SCALAR tM = -b + sqrt(d);// -b + SQRT(b^2 -4ac)
		if(tm > 0)	{
			t_hit = tm;
			return true;
		}
		
		if(tM > 0){
			t_hit = tM;
			return true;
		}
	}
	  return false;
}

void Backpropagation::trainOnlineCV(Mlp& network, 
	MATRIX& trainingInputs, 
	VECTOR& trainingTargets,
	MATRIX& testInputs,
	VECTOR& testTargets)
{
	VECTOR trainingOutputs(trainingTargets.size(), 0.0);
	VECTOR testOutputs(testTargets.size(), 0.0);
	
	while(error > tolerance && testCount < maxTestCount)
	{
		VECTOR::iterator output = trainingOutputs.begin();
		VECTOR::iterator target = trainingTargets.begin();
		for(MATRIX::iterator input = trainingInputs.begin(); 
			input != trainingInputs.end(); 
			++input, ++target, ++output)
		{
			*output = network(*input);
			double err = *output - *target;
			
			getWeightUpdates(network, *input, err);
			
			applyWeightUpdates(network);
			
			++iteration;
			
			if(iteration >= maxIterations)
				break;
		}
		
		++epoch;
		
		error = mse(trainingTargets, trainingOutputs);
		
		// Early-stopping using test (cross-validation) error
		testOutputs = network(testInputs);
		testError = mse(testTargets, testOutputs);
		if(testError < minTestError)
		{
			// Preserve test error and network weights
			minTestError = testError;
			W = network.W;			
			V = network.V;
			biasW = network.biasW;
			biasV = network.biasV;
			testCount = 0;
		}
		else
		{
			++testCount;
		}
	}
	
	network.W = W;
	network.V = V;
	network.biasW = biasW;
	network.biasV = biasV;
	testError = minTestError;
}

/*-<==>-----------------------------------------------------------------
/ Define the axis of the camera (front, up, left) in world coordinates 
/ based on the current values of the vectors target & loc 
/---------------------------------------------------------------------*/
void CCamera::initAxis() {
	front = target - loc;
	front.normalize();
	VECTOR aux = VECTOR(0,1,0);
	left = aux.cross(front);
	left.normalize();
	up = -left.cross(front);
}

  void
  JacobiPreconditioner<MATRIX,VECTOR>::apply_preconditioner(const VECTOR& Px, VECTOR& x) const
  {
    x = Px;
    unsigned int i(0);
    for (typename VECTOR::iterator it(x.begin()), itend(x.end());
	 it != itend; ++it, ++i)
      *it /= A.get_entry(i, i);
  }

bool insert_keep_sorted_and_unique(typename VECTOR::value_type p, VECTOR& procs)
{
  typename VECTOR::iterator iter = std::lower_bound(procs.begin(), procs.end(), p);
  if (iter == procs.end() || *iter != p) {
    procs.insert(iter, p);
    return true;
  }
  return false;
}

//////////////////////////////////////////////////////////////////////
// subtract two vectors
//////////////////////////////////////////////////////////////////////
VECTOR operator-(const VECTOR& x, const VECTOR& y) 
{
  VECTOR z(x.size());

  for (int i = 0; i < x.size(); i++)
    //z(i) = x(i) - y(i);
    z(i) = x[i] - y[i];

  return z;
}

//#####################################################################
// Function Initialize_Faces
//#####################################################################
void HEXAHEDRON_MESH::
Initialize_Faces()
{
    delete faces;faces=new ARRAY<VECTOR<int,4> >;
    HASHTABLE<VECTOR<int,4> > quad_list(2*elements.m); // list of faces currently found
    for(int h=1;h<=elements.m;h++){
        const VECTOR<int,8>& nodes=elements(h);
        for(int f=0;f<6;f++){
            VECTOR<int,4> face;for(int k=1;k<=4;k++) face[k]=nodes(face_indices[f][k-1]);
            if(quad_list.Set(face.Sorted())) faces->Append(VECTOR<int,4>(nodes(face[1]),nodes(face[2]),nodes(face[3]),nodes(face[4])));}}
}

// Calculate equal categories
void MgFeatureNumericFunctions::GetEqualCategories(VECTOR &values, int numCats, double dataMin, double dataMax, VECTOR &distValues)
{
    // Expected categories should be more than zero
    if (numCats <= 0)
    {
        STRING message = MgServerFeatureUtil::GetMessage(L"MgInvalidComputedProperty");

        MgStringCollection arguments;
        arguments.Add(message);
        throw new MgFeatureServiceException(L"MgServerSelectFeatures.GetEqualCategories", __LINE__, __WFILE__, &arguments, L"", NULL);
    }

    // find the range of the data values
    double min = DoubleMaxValue;
    double max = -DoubleMaxValue;

    int cnt = (int)values.size();
    if (cnt <= 0) { return; } // Nothing to do, we just send back Property Definition to clients from reader

    for (int i=0; i < cnt; i++)
    {
        double val = values[i];

        if (val > max)
            max = val;
        if (val < min)
            min = val;
    }

    // expand the range a little to account for numerical instability
    double delta = 0.0001 * (max - min);
    min -= delta;
    max += delta;
    // but don't let the values extend beyond the data min/max
    if (min < dataMin)
        min = dataMin;
    if (max > dataMax)
        max = dataMax;

    // This method ignores dataMin and dataMax.  A different "Equal" distribution
    // might ignore the actual data values and create categories based on dataMin
    // and dataMax when those values are not +/- infinity.

    // fill in the categories
    distValues.push_back(min);
    delta = (max - min) / (double)numCats;
    for (int i=1; i<numCats; i++)
    {
        double nextval = distValues[i-1] + delta;
        distValues.push_back(nextval);
    }
    distValues.push_back(max);
}

VECTOR<TYPE>::VECTOR(const VECTOR& copy)
{
	VECTOR();
	if(this != &copy)
	{
		this->reserve(copy.capacity());
		for(unsigned int i = 0; i < copy.size(); i++)
		{
			this->mpData[i] = copy[i];
		}
	}
}

void Cell::brute_force(VECTOR& rij, int degree, vector<triple>& res,triple& central_translation){
/**
  \brief Brute force generation of the neighbor list for a given pair of atoms

  \param[in] rij  Vector connecting the two atoms for which we want to construct neighbor list
  \param[in] degree The maximal number of unit cells to consider in all directions (usually, if the
             cell is large and the cutoff distance is not, it may be sufficient to have degree=1). Degree = 0 
             implies only the original unit cell with no periodic re[licas. Set degree to a larger value, if
             the simulation cell is small.
  \param[out] res The list of triples with each triple describing the integer translations (for given specific pair of atoms)
             of the original cell needed to account for all neighbors
  \param[out] central_translation The triple that makes the two atoms to be minimally separated (be in the central cell). For instance, if the two
              atoms are near the opposite sides of the box, the 1 box translation will put them together (nearby). That translation is 
              then the central translation.
    
*/


  double Roff2 = Roff * Roff;
  VECTOR r;

  if(res.size()>0){  res.clear(); }    
  central_translation.n1 = -degree;
  central_translation.n2 = -degree;
  central_translation.n3 = -degree;

  r = rij;
  central_translation.n1 = 0;
  central_translation.n2 = 0;
  central_translation.n3 = 0;
  double min_dist = r.length2();
  
  for(int n1=-degree;n1<=degree;n1++){
    for(int n2=-degree;n2<=degree;n2++){
      for(int n3=-degree;n3<=degree;n3++){

        r = (rij - n1*t1 - n2*t2 - n3*t3);
        double d = r.length2(); 
        if(d<=Roff2){
          triple t;
          t.n1 = n1;
          t.n2 = n2;
          t.n3 = n3;          
          if(d<=min_dist){ central_translation = t; min_dist = d;}
          res.push_back(t);
        }
      }
    }
  }

}

// Calculate Standard Deviation for the values
void MgFeatureNumericFunctions::GetStandardDeviation(VECTOR &values, VECTOR &distValues)
{
    double mean = 0;

    int cnt = (int)values.size();
    if (cnt <= 0) { return; } // Nothing to do, we just send back Property Definition to clients from reader

    double min = DoubleMaxValue;
    double max = -DoubleMaxValue;

    for (int i=0; i < cnt; i++)
    {
        double val = values[i];

        if (val > max)
            max = val;

        if (val < min)
            min = val;

        mean += val;
    }

    // expand min and max a little to account for numerical instability
    double delta = 0.0001 * (max - min);
    min -= delta;
    max += delta;

    // compute the mean, variance and standard deviation
    double count = (double)cnt;  // (guaranteed to be > 0)
    mean /= count;
    double variance = 0;

    for (int i=0; i < cnt; i++)
    {
        variance += (values[i] - mean) * (values[i] - mean);
    }

    double deviation = sqrt((double)(variance / count));

    // Set the base date as min date
    if (m_type == MgPropertyType::DateTime)
    {
        deviation += min;
    }

    distValues.push_back(deviation);

    return;
}

/*
--------------------------------------------------------------------------------------------------
- check collision with stairs
--------------------------------------------------------------------------------------------------
*/
bool PlanesPhysicHandler::ColisionWithStair(const AABB & actorBB, const VECTOR &Speed, VECTOR &ModifiedSpeed)
{
    float moveX = Speed.x;
    float moveZ = Speed.z;

    // calculate norm of speed
    VECTOR speedNorm = Speed.unit();

    float startX = (actorBB.P.x+actorBB.E.x)/2.0f;
    float startZ = (actorBB.P.z+actorBB.E.z)/2.0f;

    std::vector<StairPlane>::const_iterator it = _stairs.begin();
    std::vector<StairPlane>::const_iterator end = _stairs.end();

    // for each stairs
    for(int i=0; it != end; ++it, ++i)
    {
        // project point to plane and check if we cross it
        float DotProduct=speedNorm.dot(it->Normal);

        // Determine If Ray Parallel To Plane
        if (abs(DotProduct) > 0.000001f)
        {
            // Find Distance To Collision Point
            float l2=(it->Normal.dot(it->C1-VECTOR(startX, actorBB.P.y, startZ)))/DotProduct;

            // Test If Collision Behind Start or after end
            if (l2 > 0 && l2 < Speed.length())
            {
                float collionsX = startX + (speedNorm.x * l2);
                float collionsZ = startZ + (speedNorm.z * l2);

                if((collionsX >= it->minX) && (collionsX <= it->maxX))
                {
                    if((collionsZ >= it->minZ) && (collionsZ <= it->maxZ))
                    {
                        VECTOR spmY(Speed.x, 0, Speed.z);
                        VECTOR Vt(it->Normal.dot(spmY)*it->Normal);
                        VECTOR Vn(spmY - Vt);
                        ModifiedSpeed = Vn;
                        return true;
                    }
                }
            }
        }
    }

    return false;
}

void MgFeatureNumericFunctions::GetMaximum(VECTOR &values, VECTOR &distValues)
{
    // TODO: Change this algorithm to take reader directly instead of vector

    // find the range of the data values
    distValues.push_back(MgServerFeatureUtil::Maximum(values));
}

bool CheckEqualVector(
    const VECTOR& rActualVector
  , const VECTOR& rExpectedVector)
{
  // Check that the two array are of equal length
  auto ret = rActualVector.size() == rExpectedVector.size();
  if (ret) {
    for (auto i = 0u; i < rExpectedVector.size(); ++i) {
      ret = ret && rActualVector[i] == rExpectedVector[i];
      if (!ret) break;
    }
  }
  else {
    std::cout << "Vectors have different size!" << std::endl;
  }
  return ret;
}

		/**
		* Re-order a vector of labels according to the permutation.
		* perm[i] = k means that label L(i) initially at position i must be put at pos k.
		*
		* see getPermute() for the opposite transformation.
		**/
		template<typename VECTOR> VECTOR getAntiPermute(const VECTOR & labels) const
			{
			const size_t l = labels.size();
			MTOOLS_INSURE(_perm.size() == l);
			VECTOR res(l);
			for (size_t i = 0; i < l; i++) { res[_perm[i]] = labels[i]; }
			return res;
			}

Rectangle::VECTOR Rectangle::End() const
{
	VECTOR RET;
	if(Array[0]==Array[1])
		RET=Array[0];
	else
	{
		if(Array[0].x() > Array[1].x())
			RET.x()=Array[0].x();
		else
			RET.x()=Array[1].x();		
		if(Array[0].y() > Array[1].y())
			RET.y()=Array[0].y();
		else
			RET.y()=Array[1].y();
	}
	return RET;
}

//#####################################################################
// Function Find_And_Append_Adjacent_Elements
//#####################################################################
// find adjacent simplices that contains face and append them to the adjacency list
template<int d> void SIMPLEX_MESH<d>::
Find_And_Append_Adjacent_Elements(const int simplex,const VECTOR<int,d>& face)
{
    int first_node=face[1];VECTOR<int,d-1> other_nodes=face.Remove_Index(1);
    for(int t=1;t<=(*incident_elements)(first_node).m;t++){
        int simplex2=(*incident_elements)(first_node)(t);
        if(simplex2!=simplex && Nodes_In_Simplex(other_nodes,simplex2))
            (*adjacent_elements)(simplex).Append_Unique(simplex2);}
}

bool is_sorted_and_unique(const VECTOR& vec, COMPARE compare)
{
    bool sorted_and_unique = true;
    for(size_t i=1; i<vec.size(); ++i) {
        if (!compare(vec[i-1],vec[i])) {
            sorted_and_unique = false;
        }
    }
    return sorted_and_unique;
}