运输问题可行解 差额法（Vogl近似法）

在前一篇文章已经介绍了如何使用最小元素法求解，这里介绍一下差额法的求解。

差额法的思想是在行和列里找到费用单价最小的两个值，得到这两个值的差值。

在所有行和列的结果做找到差值最大的那行或者那一列，然后先处理这一行或列里面费用最小的那个地方。

其他和最小元素法一样。算法过程也类似.

下面给出算法（可以优化的地方：对于差值相等的情况的处理，算法增加变量记录一行或者一列是否Ok减少比较次数。如果你写成了更好的，希望也告诉我一下。互相学习）

	static void Cal(VecDoubleType vecS, VecDoubleType vecR
		, const VecDoubleType &vecD, VecDoubleType &vecX)
	{
		const VecDoubleType::size_type iSSize = vecS.size(), iRSize = vecR.size(), iDSize = vecD.size();
		assert(iDSize == iSSize * iRSize);

		// 数据初始化
		vecX.clear();
		VecBoolType vecOk;
		VecDoubleType::size_type i = 0, j = 0;
		for (i = 0; i != iSSize; i++)
		{
			for (j = 0; j != iRSize; j++)
			{
				vecX.push_back(0.0);
				vecOk.push_back(false);
			}
		}
		// 如果很慢后面可以优化的地方: 添加行和列的布尔变量，指示这一行或者这一列是否已经Ok，遍历的时候
		// 就不用全部遍历了
		while(true)
		{
			// 处理，
			// 如果某行只有一个或者某列只有一个了，那么他的值就能直接确定了
			CheakOnlyOne(vecS, vecR, vecOk, vecX);

			// 如果某行的发货量为0 那么这行就不用考虑，都OK了
			// 如果某列的收货量为0 那么这列也不用考虑了。
			CheakOk(vecS, vecR, vecOk);
			if (AllOk(vecOk))// 全部都Ok
			{
				break;
			}

			// 找行和列的剩下的 最小两个距离的差值， 记录差值最大的那个行或者列
			// 先重行里找到最大的那行
			VecDoubleType::size_type iRowIndex = iSSize;
			double dRowMaxBetweenDis = -1;
			for (i = 0; i != iSSize; i++)
			{
				double dFirstMin = -1, dSecondMin = -2;// 因为距离一定大于0
				for (j = 0; j != iRSize; j++)
				{
					if (vecOk[i * iRSize + j])
					{
						continue;
					}
					double dThis = vecD[i * iRSize + j];
					if (dFirstMin < 0)
					{
						dFirstMin = dThis;
					}
					else if (dSecondMin < 0)
					{
						if (dThis - dFirstMin < g_dTol)
						{
							dSecondMin = dFirstMin;
							dFirstMin = dThis;
						}
						else
						{
							dSecondMin = dThis;
						}
					}
					else
					{
						if (dThis - dFirstMin < g_dTol)
						{
							dSecondMin = dFirstMin;
							dFirstMin = dThis;
						}
						else if (dThis - dSecondMin < g_dTol)
						{
							dSecondMin = dThis;
						}
					}
				}
				if (dFirstMin > -1 && dSecondMin > -1)
				{
					double dBetwnntDis = dSecondMin - dFirstMin;
					if (dBetwnntDis - dRowMaxBetweenDis > g_dTol)
					{
						dRowMaxBetweenDis = dBetwnntDis;
						iRowIndex = i;
					}
				}
			}

			// 同样来找列的
			VecDoubleType::size_type iColIndex = iRSize;
			double dColMaxBetweenDis = -1;
			for (j = 0; j != iRSize; j++)
			{
				double dFirstMin = -1, dSecondMin = -2;// 因为距离一定大于0
				for (i = 0; i != iSSize; i++)
				{
					if (vecOk[i * iRSize + j])
					{
						continue;
					}
					double dThis = vecD[i * iRSize + j];
					if (dFirstMin < 0)
					{
						dFirstMin = dThis;
					}
					else if (dSecondMin < 0)
					{
						if (dThis - dFirstMin < g_dTol)
						{
							dSecondMin = dFirstMin;
							dFirstMin = dThis;
						}
						else
						{
							dSecondMin = dThis;
						}
					}
					else
					{
						if (dThis - dFirstMin < g_dTol)
						{
							dSecondMin = dFirstMin;
							dFirstMin = dThis;
						}
						else if (dThis - dSecondMin < g_dTol)
						{
							dSecondMin = dThis;
						}
					}
				}
				if (dFirstMin > -1 && dSecondMin > -1)
				{
					double dBetwnntDis = dSecondMin - dFirstMin;
					if (dBetwnntDis - dColMaxBetweenDis > g_dTol)
					{
						dColMaxBetweenDis = dBetwnntDis;
						iColIndex = j;
					}
				}
			}

			assert(iRowIndex != iSSize && iColIndex != iRSize);
			bool bInRow = true;// 在那行里面找
			if (dColMaxBetweenDis - dRowMaxBetweenDis > g_dTol)
			{
				bInRow = false;// 在列里面找
			}

			VecDoubleType::size_type iMin = iSSize, jMin = iRSize;// 位置
			if (bInRow)
			{
				iMin = iRowIndex;
				double dMinDis = -1;
				for (j = 0; j != iRSize; j++)
				{
					if (vecOk[iRowIndex * iRSize + j])
					{
						continue;
					}
					double dThis = vecD[iRowIndex * iRSize + j];
					if (dMinDis < 0)
					{
						dMinDis = dThis;
						jMin = j; 
					}
					else if (dMinDis - dThis > g_dTol)
					{
						dMinDis = dThis;
						jMin = j;
					}
				}
			} 
			else
			{
				jMin = iColIndex;
				double dMinDis = -1;
				for (i = 0; i != iSSize; i++)
				{
					if (vecOk[i * iRSize + iColIndex])
					{
						continue;
					}
					double dThis = vecD[i * iRSize + iColIndex];
					if (dMinDis < 0)
					{
						dMinDis = dThis;
						iMin = i; 
					}
					else if (dMinDis - dThis > g_dTol)
					{
						dMinDis = dThis;
						iMin = i;
					}
				}
			}
			
			// 这个位置先处理
			assert(iMin != iSSize && jMin != iRSize);
			double dMinItem  = min(vecS[iMin], vecR[jMin]);
			vecX[iMin * iRSize + jMin] = dMinItem;
			vecS[iMin] -= dMinItem;
			vecR[jMin] -= dMinItem;
			vecOk[iMin * iRSize + jMin] = true;
		}


	}


	static void CheakOnlyOne(VecDoubleType &vecS, VecDoubleType &vecR,
		VecBoolType &vecOk, VecDoubleType &vecX)
	{
		const VecDoubleType::size_type iSSize = vecS.size(), iRSize = vecR.size();
		VecDoubleType::size_type i = 0, iOnlyOne = 0, j = 0, jOnlyOne = 0;
		
		// 判断行
		for ( i = 0; i != iSSize; i++)
		{
			int iNotOkNum = 0;
			for (j = 0; j != iRSize; j++)
			{
				if (!vecOk[i * iRSize + j])
				{
					iNotOkNum++;
					if (iNotOkNum > 1)
					{
						break;
					}
					iOnlyOne = i;
					jOnlyOne = j;
				}
			}
			if (1 == iNotOkNum)
			{
				// 这一行只有一个了
				double dMin = min(vecS[iOnlyOne], vecR[jOnlyOne]);
				vecX[iOnlyOne * iRSize + jOnlyOne] = dMin;
				vecS[iOnlyOne] -= dMin;
				vecR[jOnlyOne] -= dMin;
				vecOk[iOnlyOne * iRSize + jOnlyOne] = true;
			}
		}
		// 判断列
		for (j = 0; j != iRSize; j++)
		{
			int iNotOkNum = 0;
			for (i = 0; i != iSSize; i++)
			{
				if (!vecOk[i * iRSize + j])
				{
					iNotOkNum++;
					if (iNotOkNum > 1)
					{
						break;
					}
					iOnlyOne = i;
					jOnlyOne = j;
				}
			}

			if (1 == iNotOkNum)
			{
				// 这一列只有一个了
				double dMin = min(vecS[iOnlyOne], vecR[jOnlyOne]);
				vecX[iOnlyOne * iRSize + jOnlyOne] = dMin;
				vecS[iOnlyOne] -= dMin;
				vecR[jOnlyOne] -= dMin;
				vecOk[iOnlyOne * iRSize + jOnlyOne] = true;
			}
		}
	}