//分支限界法求解TSP問題的程序
#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
using namespace std;

//定義相關結構體
typedef struct node_data      //結點數據結構
{
	float c[20][20];          /*費用矩陣*/
	int row_init[20];         /*費用矩陣的當前行映射為原始行*/
	int col_init[20];         /*費用矩陣的當前列映射為原始列*/
	int row_cur[20];          /*費用矩陣的原始行映射為當前行*/
	int col_cur[20];          /*費用矩陣的原始列映射為當前列*/
	int ad[20];               /*回路頂點鄰接表*/
	int k;                    /*當前費用矩陣的階*/
	float w;                  /*結點的下界*/
}NODE;

typedef struct                /*堆結構數據*/
{
	NODE *p;                  /*指向結點元素的指針*/
	float w;                  /*所指向結點的下界，堆元素的關鍵字*/ 
}HEAP;

//主要函數申明部分
float row_min(NODE *node,int row,float &second);
float col_min(NODE *node,int col,float &second);
float array_red(NODE *node);
float edge_sel (NODE *node,int &vk,int &vl);
void del_rowcol(NODE *node,int vk,int vl);
void edge_byp(NODE *node,int vk, int vl);
NODE *initial(float c[20][20],int n);
void insert(HEAP *&heap,int &n_heap,HEAP z);
HEAP delete_min(HEAP *&heap,int &n_heap);

//主函數開始
int main()
{
    int n=20;
	int ad[20],path[20];
	float c[20][20];
	float max_distance, min_distance;
	string filename="distance.txt";
	ifstream sensor;

	sensor.open(filename.c_str());             //打開文件讀取距離數據信息
	if(sensor.fail())
	{
		cout<<"Error opening input file\n";
	}
	else 
	{
		sensor>>max_distance>>min_distance;    //沒有到達文件末尾時循環讀取數據

		while(!sensor.eof())
		{
			for(int i=0;i<20;i++)
			{
				for(int j=0;j<20;j++)
				{
					sensor>>c[i][j];
				}
				c[i][i]=max_distance+1;
			}
		}
	}                                         //讀取距離信息結束

	int i,j,vk,vl;
    float d,w=0;
	NODE *xnode;
	NODE *ynode;
	NODE *znode;
	int n_heap=0;
    HEAP *heap=new HEAP[50*30];
    HEAP  x,y,z;

    xnode=initial(c,n);
    xnode->w=array_red(xnode);
    while(xnode->k!=0)
	{
		d=edge_sel(xnode,vk,vl);
        znode=new NODE;
        *znode =*xnode;
        znode->c[vk][vl]=max_distance+1;
        array_red(znode);
        znode->w=xnode->w+d;
        z.w=znode->w;
        z.p=znode;
        insert(heap,n_heap,z);  
        ynode=new NODE;
        *ynode=*xnode;
        edge_byp(ynode,vk,vl);
        del_rowcol(ynode,vk,vl);
        ynode->w=array_red(ynode);
        ynode->w+=xnode->w;
        y.w=ynode->w;
        y.p=ynode;

        if(ynode->k==2)
		{
			if(ynode->c[0][0]==0)
			{
				ynode->ad[ynode->row_init[0]]=ynode->col_init[1];
                ynode->ad[ynode->row_init[1]]=ynode->col_init[0];
			}
            else
			{
			    ynode->ad[ynode->row_init[0]]=ynode->col_init[0];
                ynode->ad[ynode->row_init[1]]=ynode->col_init[1];
			}
            ynode->k=0;
		}
        insert (heap,n_heap,y);
        delete xnode;
        x=delete_min(heap,n_heap);  
        xnode=x.p;
	}
    
    w=xnode->w;

    for(i=0;i<n;i++)               //保存鄰接表
		ad[i]=xnode->ad[i];
	
    delete xnode;

	for (i=0;i<n_heap;i++)         //釋放結點數據空間
	{
		delete heap->p;
		heap--;
	}
	heap++;

    delete []heap;                //釋放堆空間

	int t=0,k=1;                  //將鄰接表轉化為遍歷城市的順序列表，保存在path[]中
	path[0]=0;
	while(k<n)
	{
		j=t;
		t=ad[j];
		path[k]=t;
		k++;
	}                            

	cout<<"the result of the path:"<<endl;        //輸出最優遍歷城市順序的結果，轉化為A、B的形式

	for(j=0;j<20;j++)
	{
		cout<<char(path[j]+65)<<' ';
	}

	cout<<endl<<"the total length is: "<<w<<endl; //輸出最短路程的值


	return 0;
}                   //主函數結束

/**********************************函數定義部分***********************************/

//返回由指針node所指向節點費用矩陣中第row行最小值，并把次小值回送與引用變量second中
float row_min(NODE *node,int row,float &second)
{
	float temp;
    int i;
    if (node->c[row][0],node->c[row][1])
	{
		temp =node->c[row][0];   second=node->c[row][1];
    }
    else 
	{
		temp=node->c[row][1]; second=node->c[row][0];
    }
    for (i=2;i<node->k;i++)
	{
		if(node->c[row][i]<temp)
		{
			second=temp;temp=node->c[row][i];
		}
        else if(node->c[row][i]<second)
			second=node->c[row][i];
	}
    return temp;
}
//返回由指針node所指向節點費用矩陣中第col列最小值，并把次小值回送與引用變量second中
float col_min(NODE *node,int col, float &second)
{
	float temp;
    int i;
    if (node->c[0][col],node->c[1][col])
	{
		temp =node->c[0][col]; second=node->c[1][col];
    }
    else
	{
		temp=node->c[1][col]; second=node->c[0][col];
    }
    for (i=2;i<node->k;i++)
	{
		if (node->c[i][col]<temp )
		{
			second =temp; temp =node->c[i][col];
		}
        else if (node->c[i][col]<second)
			second=node->c[i][col];
     }
  return temp;
}
//歸約指針node指向的節點的費用矩陣 返回值為歸約常數
float array_red(NODE *node)
{
  int i,j;
  float temp,temp1,sum=0;
  for(i=0;i<node->k;i++)
  {
	  temp=row_min(node,i,temp1);
      for (j=0;j<node->k;j++)
		  node->c[i][j] -=temp;
      sum +=temp;
  }
  for(j=0;j<node->k;j++)
  {
	  temp=col_min(node,j,temp1);
      for (i=0;i<node->k;i++)
		  node->c[i][j] -=temp;
      sum +=temp;
  }
  return sum;
}
//計算dkl 并且搜索分支的邊 返回dkl的值
float edge_sel(NODE *node,int &vk,int &vl)
{
  int i,j;
  float temp,d=0;
  float *row_value=new float[node->k];
  float *col_value=new float[node->k];
  for (i=0;i<node->k;i++)
      row_min(node,i,row_value[i]);
  for (i=0;i<node->k;i++)
	  col_min(node,i,col_value[i]);

  for (i=0;i<node->k;i++)
  {
	  for(j=0;j<node->k;j++)
	  {
		  if (node->c[i][j]==0)
		  {
			  temp=row_value[i]+col_value[j];
		      if(temp>d)
			  {
				  d=temp; vk=i;vl=j;
			  }
		  }
	  }
  }
  delete row_value;
  delete col_value;
  return d;
}
//del_rowcol(NODE *node,int vk,int vl)刪除費用矩陣當前vk行 第vl列的所有元素
void del_rowcol(NODE *node,int vk,int vl)
{
	int i,j,vk1,vl1;
    for (i=vk;i<node->k-1;i++)
		for (j=0;j<vl;j++)
			node->c[i][j]=node->c[i+1][j];
	for (j=vl;j<node->k-1;j++)
		for (i=0;i<vk;i++)
			node->c[i][j]=node->c[i][j+1];
    for (i=vk;i<node->k-1;i++)
		for (j=vl;j<node->k-1;j++)
			node->c[i][j]=node->c[i+1][j+1];
    vk1=node->row_init[vk];
    node->row_cur[vk1]=-1;
    for (i=vk1+1;i<20;i++)
		node->row_cur[i]--;//
    vl1=node->col_init[vl];
    node->col_cur[vl1]=-1;
	for (i=vl1+1;i<20;i++)
		node->col_cur[i]--;
    for (i=vk;i<node->k-1;i++)
		node->row_init[i]=node->row_init[i+1];
    for(i=vl;i<node->k-1;i++)
		node->col_init[i]=node->col_init[i+1];
    node->k--;
}
//把vk vl表示的邊登記到頂點鄰接表 并旁路矩陣中有關的邊
void edge_byp(NODE *node,int vk, int vl)
{
	int i,j,k,l;
    vk=node->row_init[vk];
    vl=node->col_init[vl];
    node->ad[vk]=vl;
    for (i=0;i<20;i++)
	{
		j=i;
        while (node->ad[j]!=-1)
        j=node->ad[j];
        if(i!=j)
		{
			l=node->row_cur[j];
            k=node->col_cur[i];
            if((k>0)&&(l>0))
            node->c[l][k]=100;   //以后要改進的地方 100是個特殊值
		}
	}
}
//初始化函數
NODE *initial(float c[20][20],int n)
{
	int i,j;
    NODE *node=new NODE;
    for (i=0;i<n;i++)
		for (j=0;j<n;j++)
			node->c[i][j]=c[i][j];
    for (i=0;i<n;i++)
	{
		node->row_init[i]=i;
        node->col_init[i]=i;
        node->row_cur[i]=i;
        node->col_cur[i]=i;
    }
    for (i=0;i<n;i++)
		node->ad[i]=-1;
    node->k=n;
    return node;
}
//向最小堆棧插入結點函數
void insert(HEAP *&heap,int &n_heap,HEAP z)
{
	HEAP *temp;
	int k=n_heap;
	if(k==0) 
	{
		*heap=z;
	}
	else
	{
		temp=heap;
	    while(heap->w<z.w&&k>0)
		{
			*(heap+1)=*heap;
		    k--;
		    if(k>0)  heap--;
		}
	    if(k>0)  heap++;
	    *heap=z;
        heap=temp+1;
	}
    n_heap++;
}

//刪除堆棧頭結點函數
HEAP delete_min(HEAP *&heap,int &n_heap)
{
	HEAP temp;
	temp=(*heap);
	heap--;
	n_heap--;
	return temp;
}