C++ DblVector类代码示例

// x given 
int FitParabola(const DblVector& PointsY, const DblVector& PointsX, double& a, double& b, double& c)
{
	// fit parabola
	assert(PointsX.size()==PointsY.size());
	DblMatrix S;
	S.Assign(4, 3, 0.0);
	for(unsigned int i=0; i<PointsX.size(); i++)
	{
		S[0][0] += PointsX[i] * PointsX[i] * PointsX[i] * PointsX[i];
		S[0][1] += PointsX[i] * PointsX[i] * PointsX[i];
		S[0][2] += PointsX[i] * PointsX[i];
		S[0][3] += PointsX[i] * PointsX[i] * PointsY[i];
		S[1][2] += PointsX[i];
		S[1][3] += PointsX[i] * PointsY[i];
		S[2][3] += PointsY[i];
		
	}
	S[1][1] = S[0][2];
	S[2][2] = (double)PointsX.size();
	S.Diagonalize();
	DblVector Solutions;
	int Ret = S.SolveLinearCramer(Solutions);
	if(Ret==-1) return -1;
	else
	{
		a = Solutions[0];
		b = Solutions[1];
		c = Solutions[2];
	}
	return Ret;
}

//------------------------------------------------------------------------------
//  void build(const NodePairDeque &nodepairs, Topology* topology)
//------------------------------------------------------------------------------
void TopologyFactory::build(const NodePairDeque &nodepairs, Topology* topology)
{
    TRACE("TopologyFactory::build -->");
    for(NodePairDeque::const_iterator iter = nodepairs.begin();
        iter != nodepairs.end(); ++iter)
    {
		if (iter->source != iter->destination)
		{
	        double link_cap = link_capacity == (RandomVar*) 0 ?
	            default_capacity : link_capacity->generate();
	        DblVector d;
	        d.push_back(iter->source);             // source
	        d.push_back(iter->destination);        // dest
	        d.push_back(0.0);                      // core
	        d.push_back(link_cap); // cap
	        topology->link_list->insert(d);
	        if (!topology->is_directed)
	        {
	            DblVector d2;
	            d2.push_back(iter->destination);  // source
	            d2.push_back(iter->source);       // dest
	            d2.push_back(0.0);                // core
	            d2.push_back(link_cap); // cap
	            topology->link_list->insert(d2);
	        } //end: if (!topology->is_directed)
		} //end: if (iter->source != iter->destination)
    }
    TRACE("TopologyFactory::build <--");
}

void DblMatrix::Assign(int w, int h, double Val)
{
	DblVector v;
	v.assign(w, Val);
	for(int j=0; j<h; j++)
		this->push_back(v);
}

void DblVector::MulElem(const DblVector& B, DblVector& Res) const
{
	assert(this->size()==B.size());
	Res.clear();
	for(unsigned int i=0; i<this->size(); i++)
		Res.push_back((*this)[i]*B[i]);
}

// x given
int FitLine(const DblVector& PointsY, const DblVector& PointsX, double& a, double& b)
{
	assert(PointsX.size()==PointsY.size());
	DblMatrix S;
	S.Assign(3, 2, 0.0);
	for(unsigned int i=0; i<PointsX.size(); i++)
	{
		S[0][0] += PointsX[i] * PointsX[i];
		S[0][1] += PointsX[i];
		S[0][2] += PointsX[i] * PointsY[i];
		S[1][2] += PointsY[i]; 
	}
	S[1][1] = (double)PointsX.size();
	S.Diagonalize();
	DblVector Solutions;
	int Ret = S.SolveLinearCramer(Solutions);
	if(Ret) 
	{
		a = Solutions[0];
		b = Solutions[1];
	}
	else std::cout << "Error in line fitting" << std::endl;
	
	return Ret;
}

int FitLineIt(const DblVector& PointsY, const DblVector& PointsX,
			  double& a, double& b, int It, double ErrWidth)
{
	// base case
	if(It==0) return 0;
	std::cout << "Line fitting no samples: " << PointsY.size() << std::endl;

	// fitting
	FitLine(PointsY, PointsX, a, b);

	double Error = 0.0;
	for(int i=0; i<(int)PointsY.size(); i++)
	{
		double Val0 = a*PointsX[i]+b;
		double Val1 = PointsY[i]; 
		Error += dblsqr(Val1-Val0);
		//std::cout << Val1 << " " << Val0 << std::endl;
	}
	Error /= (double)PointsY.size();
	Error = sqrt(Error);

	DblVector px, py; int Cnt=0;
	for(int k=0; k<(int)PointsX.size(); k++)
	{
		double Val = a*PointsX[k]+b;
		double Dist = sqrt(dblsqr(PointsY[k] - Val));
		if(Dist<Error*ErrWidth)
		{
			px.push_back(PointsX[k]);
			py.push_back(PointsY[k]);
			Cnt++;
		}
	}
	return FitLineIt(py, px, a, b, It-1, ErrWidth);
}

int FitQuadraticSurfaceIt(const DblVector& PX, const DblVector& PY,
							const DblVector& PZ, DblVector& Parameters,
							int It, double ErrWidth)
{

	// base case
	if(It==0) return 0;
	std::cout << "Quad surface fitting no samples: " << PX.size() << std::endl;

	// fitting
	FitQuadraticSurface(PX, PY, PZ, Parameters);
	double Error = GetQuadraticSurfaceError(PX, PY, PZ, Parameters);
	
	DblVector px, py, pz; int Cnt=0;
	for(int k=0; k<(int)PX.size(); k++)
	{
		double Val = QuadraticSurfaceVal(PX[k], PY[k], Parameters);
		double Dist = sqrt((PZ[k] - Val) * (PZ[k] - Val));
		if(Dist<Error*ErrWidth)
		{
			px.push_back(PX[k]);
			py.push_back(PY[k]);
			pz.push_back(PZ[k]);
			Cnt++;
		}
		
	}
	return FitQuadraticSurfaceIt(px, py, pz, Parameters, It-1, ErrWidth);
}

void DblMatrix::MakeVec(DblVector& v) const
{
	v.clear();
	for(int j=0; j<GetHeight(); j++)
		for(int i=0; i<GetWidth(); i++)
			v.push_back((*this)[j][i]);
}

// fit a two dimensional quadratic surface
int FitQuadraticSurface(const DblVector& PX, const DblVector& PY,
							const DblVector& PZ, DblVector& Parameters)
{
	// fit parabola
	assert(PX.size()==PY.size());
	DblMatrix S;
	S.Assign(7, 6, 0.0);
	for(unsigned int i=0; i<PX.size(); i++)
	{
		double X2 = PX[i] * PX[i];
		double Y2 = PY[i] * PY[i];
		double X3 = X2 * PX[i];
		double Y3 = Y2 * PY[i];
		double X4 = X3 * PX[i];
		double Y4 = Y3 * PY[i];
		
		S[0][1] += PY[i];
		S[1][1] += Y2;
		
		S[0][2] += PX[i];
		S[1][2] += PX[i]*PY[i];
		S[2][2] += X2;

		S[0][3] += PX[i]*PY[i];
		S[1][3] += PX[i]*Y2;
		S[2][3] += X2*PY[i];
		S[3][3] += X2*Y2;

		S[0][4] += Y2;
		S[1][4] += Y3;
		S[2][4] += PX[i]*Y2;
		S[3][4] += PX[i]*Y3;
		S[4][4] += Y4;

		S[0][5] += X2;
		S[1][5] += X2*PY[i];
		S[2][5] += X3;
		S[3][5] += X3*PY[i];
		S[4][5] += X2*Y2;
		S[5][5] += X4;

		// function values
		S[0][6]+=PZ[i];
		S[1][6]+=PY[i]*PZ[i];
		S[2][6]+=PX[i]*PZ[i];
		S[3][6]+=PX[i]*PY[i]*PZ[i];
		S[4][6]+=Y2*PZ[i];
		S[5][6]+=X2*PZ[i];

	}
	S[0][0] = (double)PX.size();
	S.Diagonalize();
	return S.SolveLinearCramer(Parameters);
}

double GetQuadraticSurfaceError(const DblVector& PX, const DblVector& PY,
							const DblVector& PZ, DblVector& Parameters)
{
	double Sum=0.0;
	for(int i=0; i<(int)PX.size(); i++)
	{
		double Val = QuadraticSurfaceVal(PX[i], PY[i], Parameters);
		Sum += (PZ[i] - Val) * (PZ[i] - Val);
	}
	return sqrt(Sum/(double)PX.size());
}

Topology* TopologyFactory::createTopologyAdjacency(const TString &description)
{
	TRACE("TopologyFactory::createTopologyAdjacency -->");

	Topology* result = new Topology();
	result->input_type = 0;
	result->is_directed = false;

	{	// First find number of nodes, which requires to read entire file
		FileReader fr(description.at(1));
		std::set<int> nodes;
		TString line;

		// read #nodes
		for(bool cont = fr.firstLine(line); cont; cont = fr.nextLine(line))
		{
			nodes.insert(atoi(line.at(0).c_str()));
			nodes.insert(atoi(line.at(1).c_str()));
		}
		result->number_of_nodes = nodes.size();
		result->edge_nodes.assign(nodes.begin(), nodes.end());
	} //end: block
    result->number_of_qos     = 0;
    result->link_list         = new LinkList(result->number_of_nodes);

	{	// Build links
		FileReader fr(description.at(1));
		TString line;

		// read #nodes
		for(bool cont = fr.firstLine(line); cont; cont = fr.nextLine(line))
		{
			DblVector link;
			link.push_back(atof(line.at(0).c_str()));	// src
			link.push_back(atof(line.at(1).c_str()));	// dst
			link.push_back(0.0);						// core
			link.push_back(atof(line.at(2).c_str()));	// cap
			result->link_list->insert(link);
			if  (!result->is_directed)
			{
				DblVector rlink(link);
				rlink[0] = link[1];
				rlink[1] = link[0];
				result->link_list->insert(rlink);
			} // end: if
		} //end:  for
	} //end: block
	TRACE("TopologyFactory::createTopologyAdjacency <--");
	return result;
}

SortedVector::SortedVector(const DblVector& v)
{
	clear();
	for(unsigned int i=0; i<v.size(); i++)
		push_back(SortStruct(v[i], i));
	std::sort(begin(), end(), SortStructGreater);
}

unsigned long FitNormalizedPolynomial(const DblVector& y, const DblVector& x, int degree, double* coefficients, double* max, double* min)
{

	if(y.size() != x.size())
	{
		std::cerr << "ipa_Utils::FitPolynomial: Error" << std::endl;
		std::cerr << "\t ... vector 'y' and vector 'x' must have same size." << std::endl;
		return ipa_Utils::RET_FAILED;
	}

	/// Extract min and max from 32bit distance values
	/// to scale x-values
	(*min)=DBL_MAX;
	(*max)=-DBL_MAX;
	for(unsigned int k=0; k<x.size(); k++)
	{
		double d = x[k];
		if(d < (*min)) (*min)=d;
		if(d > (*max)) (*max)=d;
	}
	double diffMaxMin = (*max) - (*min);

	/// Normalize data points, to assert numeric stability during Matrix inversion
	std::vector<double> y_Normalized;
	std::vector<double> x_Normalized;
	for(unsigned int k=0; k<x.size(); k++)
	{
		x_Normalized.push_back((x[k]-(*min))/diffMaxMin); 
		y_Normalized.push_back(y[k]/diffMaxMin);
	}

	/// Compute normalized polynomial coefficients through least squares fitting
	double* c;
	c = Wm4::PolyFit2<double>((int)x.size(), &x_Normalized[0], &y_Normalized[0], degree);

	for (int i = 0; i <= degree; i++)
	{
		coefficients[i] = c[i];
	}

	delete[] c;

	return ipa_Utils::RET_OK;
}

unsigned long FitPolynomial(const DblVector& y, const DblVector& x, int degree, double* coefficients)
{

	if(y.size() != x.size())
	{
		std::cerr << "ipa_Utils::FitPolynomial: Error" << std::endl;
		std::cerr << "\t ... vector 'y' and vector 'x' must have same size." << std::endl;
		return ipa_Utils::RET_FAILED;
	}

	/// Compute normalized polynomial coefficients through least squares fitting
	double* c;
	c = Wm4::PolyFit2<double>((int)x.size(), &x[0], &y[0], degree);

	for (int i = 0; i <= degree; i++)
	{
		coefficients[i] = c[i];
	}
	delete[] c;

	return ipa_Utils::RET_OK;
}

// solve linear NxN equation system (cramer)
int DblMatrix::SolveLinearCramer(DblVector& Solution)
{
	// get main determinant
	DblMatrix DetMat;
	double Det=0.0;
	int w = GetWidth();
	int h = GetHeight();
	for(int i=0; i<h; i++)
	{
		DetMat.push_back(DblVector());
		for(int j=0; j<w-1; j++)
			DetMat[i].push_back((*this)[i][j]);
	}
	Det = DetMat.GetDeterminantNxN();
	if(fabs(Det)<THE_EPS_DEF) return 0;
	// get solutions
	Solution.clear();
	for(int k=0; k<h; k++)
		Solution.push_back(this->GetHelpDeterminant(k)/Det);

	return 1;
}