Поправте код программы

[GanJIk]

Здравствуйте.Недавно сел писать программу для решения задач по симплекс методу,вроде что то с чем то но не работает.Если кто может поправте.А если уже кто писал такие программы дайте посмотреть исходный код!

[ACE Valery]

Выложите, пожалуйста, сюда код программы ТЕКСТОМ, а не вложением. И не забудьте оформить его тегом [code].

[GanJIk]

#include <vector>
#include <string> 
#include <math.h>
using namespace std;

#define SIMPLEX_DONE 0		  // оптимизация успешно завершена 
#define SIMPLEX_NO_SOLUTION 1 // задача не имеет решения (не удается найти базис)  
#define SIMPLEX_NO_BOTTOM 2   // решения нет, т.к. линейная форма не ограничена снизу  
#define SIMPLEX_NEXT_STEP 3   // для получения решения нужно сделать еще хотя бы один шаг 
#define MAX_VAL 0.1e-12		  //точность (значение, удовлетворяющее -MAX_VAL < X < MAX_VAL считается нулем)  
#define MaxExtended 1.1e+100  //Максимальное значение
//1.7E+308
//extended -3.6 x 10^4951 до 1.1 x 10^4932
//MaxExtended 1.1 E+ 4932

typedef enum {  Equal,Less,Greater } TOperation;
typedef struct TConstrain 
{ 
	vector<double>    A; 
	double    B; 
	TOperation   Sign; 
	bool      isT; 
} TConstrain; 

class TSimplex
{
public:
	int      M, N;   //  M - число строк, N - число столбцов  
	int     RealN;   // реальное число переменных, изначально вошедших в задачу  
	vector<TConstrain>   Cons; 
	vector<double>	C; 
	double     L; 
	vector<int>      Basis; 
	bool      Max;
	
	TSimplex::TSimplex()
	{
	}

	//Конструктор
	TSimplex::TSimplex(vector<double> _C, bool Maximize)
	{ 
		
		int      j; 

		N = _C.size(); 
		RealN = N; 
		M = 0; 
		C.resize(N); 
		Max = Maximize; 
		if ((!Maximize)) 
			for (j = 0; j < N ; j++) 
				C[j] = -_C[j]; 
		else 
			for (j = 0; j < N ; j++) 
				C[j] = _C[j]; 
		Max = Maximize; 
		L = 0; 
		
	}
	
	//Если num лежит внутри области точности значений, возвращает 0
	double TSimplex::DoPrec(double num)
	{ 
		if ((num < MAX_VAL) && (num > -MAX_VAL)) 
			num = 0; 
		return num; 
	} 

	//Добавление ограничения вида Ai[1]*x1+Ai[2]*x2+....+Ai[n]*xn <=> Bi
	void TSimplex::AddCons(double _B, vector<double> _A, TOperation   Sign)
	{ 
		int	j; 
		if ((_A.size() > N))
			SetAllLengths(_A.size());

		M++; 

		Cons.resize(M); 
		Cons[M - 1].B = _B;
		Cons[M - 1].Sign = Sign;
		Cons[M - 1].A.resize(N);
		for (j = 0; j < _A.size() ; j++) 
			Cons[M - 1].A[j] = _A[j];
		if (_A.size() < N) 
			for (j = _A.size(); j < N ; j++) 
				Cons[M-1].A[j] = 0;
	} 

	//Суммирование строки 1 со строкой 2, домноженной на коэффициент Value   
	void TSimplex::AddString(int Num1,int Num2,double Value)
	{ 
		int      j; 
		for (j = 0; j <= N - 1; j++) 
			Cons[Num1].A[j] = Cons[Num1].A[j] + Cons[Num2].A[j] * Value;
		Cons[Num1].B = Cons[Num1].B + Cons[Num2].B * Value;
	} 

	//
	bool TSimplex::CheckBasis(void)
	{ 
		int      i, j, k; 
		bool      f; 

		Basis.resize(M); 
		for (i = 0; i <M; i++) 
			Basis[i] = -1; 
		for (j = 0; j <= N - 1; j++) 
		{ 
			f = true; 
			k =	-1; 
			i = 0; 
			while ((f && (i < M))) 
			{ 
				if (((Cons[i].A[j] != 0) && (Cons[i].A[j] != 1))) 
					f = false; 
				if ((Cons[i].A[j] == 1)) 
				{ 
					if ((k == -1)) k = i; 
					else f = false; 
				} 
				i+=1; 
			} 
			if ((f && (k != -1)))
				Basis[k] = j; 
		} 
		f = true; 
		for (i = 0; i <= M - 1; i++) 
			f = f && (Basis[i] != -1); 
		return f; 
	}

	

	void TSimplex::CreateBasis(TSimplex Simplex)
	{ 
		int i, j; 

		M = Simplex.M; 
		N = Simplex.N; 
		RealN = Simplex.RealN; 
		L = 0; 

		Cons.resize(M);
		Basis.resize(M); 
	    C.resize(N); 
		for (i = 0; i < N ; i++) 
			C[i] = 0; 
		for (i = 0; i < M ; i++) 
		{ 
			Cons[i].A.resize(N); 
			for (j = 0; j < N; j++) 
				Cons[i].A[j] = Simplex.Cons[i].A[j]; 
			Cons[i-1].B = Simplex.Cons[i-1].B; 
			Cons[i-1].Sign = Equal; 
			Cons[i-1].isT = false; 
		} 
		for (i = 0; i < M; i++) 
		{ 
			if ((Simplex.Basis[i] != -1)) 
				Basis[i] = Simplex.Basis[i]; 
			else 
			{ 
				SetAllLengths(N + 1); 
				for (j = 0; j <= M - 1; j++) 
					Cons[j].A[N - 1] = 0; 
				Cons[i].A[N - 1] = 1; 
				Cons[i].isT = true; 

				C[N - 1] = 0; 
			    for (j = 0; j <= Simplex.N - 1; j++) 
					C[j] = C[j] + Simplex.Cons[i].A[j]; 
				L = L + Cons[i].B; 
			} 
		} 
	}

	void TSimplex::Free()
	{ 	
		C.resize(0); 
		Basis.resize(0); 
		Cons.resize(0);
		M = 0; 
		N = 0; 
		RealN = 0; 
	}

	double TSimplex::GetMin()
	{ 
		//int      i; 
		if ((Max)) 
			return -L; 
		else 
			return L; 
	} 
vector<double> TSimplex::GetSolution()
	{ 
		vector<double>    Solution; 
		int     i, j; 

		Solution.resize(RealN); 
		for (j = 0; j < RealN; j++) 
		{ 
			Solution[j] = 0; 
			i = 0; 
			while (((i < M) && (Basis[i] != j))) 
				i++; 
			if (((Basis[i] == j) && (i < M))) 
				Solution[j] = Cons[i].B; 
		} 
		return Solution; 
	} 

	void TSimplex::MulString(int Number,double Value)
	{ 
		int      j; 
		for (j = 0; j < N; j++) 
			Cons[Number].A[j] = Cons[Number].A[j] * Value; 
		Cons[Number].B = Cons[Number].B * Value; 
	}

[GanJIk]

void TSimplex::Normalize(void)
	{ 
		int i; 
		for (i = 0; i <= M - 1; i++) 
			if ((Cons[i].Sign != Equal)) 
			{ 
				SetAllLengths(N + 1); 
				if ((Cons[i].Sign == Greater)) 
					Cons[i].A[N - 1] = -1; 
				else 
					Cons[i].A[N - 1] = 1; 
				Cons[i].Sign = Equal; 
			} 
	} 
	
	//Дополнение нулями до нужной длины массива коэффициентов при ограничениях
	void TSimplex::SetAllLengths(int Len)
	{ 
		int      i, j; 
		int      OldN; 

		OldN = N; 
		N = Len; 
		C.resize(N); 
		for (i = 0; i < M; i++) 
			Cons[i].A.resize(N); 
		if ((OldN < N)) 
		{ 
			for (j = OldN; j < N; j++) 
			{ 
				C[j] = 0; 
				for (i = 0; i < M; i++) 
					Cons[i].A[j] = 0; 
			} 
		} 
	}

	//Поиск числа в базисе
	int	TSimplex::FoundInBasis(int num)
	{ 
		int i; 
		bool f; 

		f = false; 
		i = 0; 
		while ((!f && (i < M))) 
		{ 
			f = (Basis[i] == num); 
			i++; 
		} 

		if ((f)) 
			return i - 1; 
		else 
			return -1; 
	}

	int TSimplex::SimplexStep()
	{ 
		int i,j; 
		bool f,opt; 
		int x,y;    //координаты опорного элемента  
		double CurMax; 
		vector<TConstrain>   temp; 
		vector<double>    tempC; 
 
		opt = true; 
		CurMax = -1; 
		for (i = 0; i <= N - 1; i++) 
		{ 
		//проверка на разрешимость  
			if ((C[i] > 0)) 
			{ 
				opt = false;   //попутная проверка на оптимальность 
				if ((C[i] > CurMax))    //поиск ведущего столбца (максимальный элемент в C[i])  
				{ 
					CurMax = C[i]; 
					x = i; 
				} 
				f = true; 
				for (j = 0; j <= M - 1; j++) 
					f = f && (Cons[j].A[i] < 0); 
				if ((f)) 
				{ 
					return SIMPLEX_NO_BOTTOM; 
				} 
			} 

		} 
		if (opt) 
			return SIMPLEX_DONE; 
		else 
		{ 
		 //	зная номер ведущего столбца, ищем номер ведущей строки  
			CurMax = MaxExtended;    //на самом деле тут будем искать минимум, а не Max  
			for (j = 0; j <= M - 1; j++) 
				if ((Cons[j].A[x] > 0)) //идем только по положительным элементам  
					if ((Cons[j].B / Cons[j].A[x] < CurMax)) 
					{ 
						CurMax = Cons[j].B / Cons[j].A[x]; 
						y = j; 
					} 
					else 
						if ((DoPrec(Cons[j].B / Cons[j].A[x] - CurMax) == 0)) 
							if ((Cons[j].isT)) 
								y = j; 

			//сохраняем текущие значения 
			temp.resize(M);
			for (j = 0; j <= M - 1; j++) 
			{ 
				temp[j].A.resize(N); 
				for (i = 0; i <= N - 1; i++) 
					temp[j].A[i] = Cons[j].A[i]; 
				temp[j].B = Cons[j].B; 
			} 
			tempC.resize(N); 
			for (i = 0; i <= N - 1; i++) 
				tempC[i] = C[i]; 
			//делаем пересчет таблицы
			//строка делиться на ведущий элемент 
			MulString(y,1 / Cons[y].A[x]); 
			//преобразование остальных элементов
			for (j = 0; j <= M - 1; j++) 
			{ 
				if ((j != y)) 
				{ 
					for (i = 0; i <= N - 1; i++) 
					{ 
						Cons[j].A[i] = DoPrec(temp[j].A[i] - temp[j].A[x] * temp[y].A[i] / temp[y].A[x]); 
					} 
					Cons[j].B = DoPrec(temp[j].B - temp[j].A[x] * temp[y].B / temp[y].A[x]); 
				} 
				else 
				{ 
					for (i = 0; i <= N - 1; i++) 
						Cons[j].A[i] = DoPrec(Cons[j].A[i]); 
				} 
			} 

			//и строка с коэффициентами функции 
			for (i = 0; i <= N - 1; i++) 
			{ 
				C[i] = DoPrec(tempC[i] - tempC[x] * temp[y].A[i] / temp[y].A[x]); 
			} 

			Basis[y] = x; 

			//и сама функция: 
			L = DoPrec(L - tempC[x] * temp[y].B / temp[y].A[x]); 

			for (i = 0; i <= M - 1; i++) 
				temp[i].A.resize(0); 
			
			temp.resize(0); 
			tempC.resize(0); 

			return SIMPLEX_NEXT_STEP; 

		} 
	} 

	//Поиск решения для заданной функции с учетом ограничений
	//Возвращаемые значения (тип integer):
	//SIMPLEX_DONE - решение успешно найдено
	//SIMPLEX_NO_SOLUTION - решения нет, или невозможно найти опорный базис
	//SIMPLEX_NO_BOTTOM - решения нет в силу неограниченного убывания линейной формы
	int	TSimplex::Solve()
	{ 
		int      i, j; 
		TSimplex     Simplex; 
		bool      f; 
		int      Step; 
		double    cc; 
 
		Normalize(); 
		f = false; 
		if ((!CheckBasis())) 
		{ 
			CreateBasis(*this); 
			Simplex.Solve(); 
			f = Simplex.GetMin() != 0; 
			if ((!f)) for (i = 0; i <= M - 1; i++) 
			{ 
				for (j = 0; j <= N - 1; j++) 
					Cons[i].A[j] = Simplex.Cons[i].A[j]; 
				Cons[i].B = Simplex.Cons[i].B; 
				Cons[i].isT = false; 
				Basis[i] = Simplex.Basis[i]; 
				cc = C[Basis[i]]; 
				for (j = 0; j <= N - 1; j++) 
					C[j] = DoPrec(C[j] - cc * Cons[i].A[j]); 

				L = DoPrec(L - cc * Cons[i].B); 

			} 
			Simplex.Free(); 
		} 
		if ((f)) Step = SIMPLEX_NO_SOLUTION; 
		else 
			do 
			{ 
				Step = SimplexStep(); 
			} 
			while (Step == SIMPLEX_NEXT_STEP);
			
		return Step; 
	}
};