==========
新手教程-最大割问题
==========

问题描述
==========

最大割问题是NP完备问题. 给定一张图, 求一种分割方法,
将所有顶点分割成两群, 同时使得被切断的边的数量最大,或边的权重最大.

以无向无权图为例. 在图 :math:`G(V,E)`\ 中, :math:`V`\为图的顶点集合, :math:`E`\为图的边集, :math:`w`\ 为图的邻接矩阵. 
对于 :math:`i,j∈V`\, :math:`w_{ij}`\ 表示顶点 :math:`i`\到顶点 :math:`j`\是否有边, 有连边关系则取 :math:`1`\ ,无连边关系则取 :math:`0`\.
以决策变量 :math:`\sigma_i`\表示顶点 :math:`i`\的分类, 其可能的取值为 :math:`{1,-1}`\ ,分别表示将顶点 :math:`i`\ 分为A类或B类.

则在给定的无向图中，将所有顶点分割成两群的分割方法所对应割取的边的个数为Z，模型表示为:

.. math:: max Z = (\sum_{i<j,i∈V}\sum_{j∈V}w_{ij}-\sum_{i<j,i∈V}\sum_{j∈V}w_{ij}\sigma_i\sigma_j)/2

以一个四顶点实例说明，如下图所示，通过观察可以发现将1、2分为A类，3、4分为B类的“割”法将得到问题的最优解 :math:`Z=4`\.

.. image:: maxcut1.png

通过连边关系可知，邻接矩阵为

.. image:: formula1.png

.. math:: \sum_{i<j,i∈V}\sum_{j∈V}w_{ij} = w_{12}+w_{13}+w_{14}+w_{23}+w_{24}+w_{34}

.. math:: \sum_{i<j,i∈V}\sum_{j∈V}w_{ij}\sigma_i\sigma_j = w_{12}\sigma_1\sigma_2+w_{13}\sigma_1\sigma_3+w_{14}\sigma_1\sigma_4+w_{23}\sigma_2\sigma_3+w_{24}\sigma_2\sigma_4+w_{34}\sigma_3\sigma_4

当顶点1、2为一组，顶点3、4为另一组时, :math:`\sigma_1=\sigma_2=1, \sigma_3=\sigma_4=-1`. 则上式变为

.. math:: \sum_{i<j,i∈V}\sum_{j∈V}w_{ij}\sigma_i\sigma_j = w_{12}-w_{13}-w_{14}-w_{23}-w_{24}+w_{34}

此时目标函数为

.. math:: max Z = (\sum_{i<j,i∈V}\sum_{j∈V}w_{ij}-\sum_{i<j,i∈V}\sum_{j∈V}w_{ij}\sigma_i\sigma_j)/2= w_{13}+w_{14}+w_{23}+w_{24} = 4

最大割数量为4，符合前文通过观察得到的答案。

注意到, :math:`w_{ij}`\ 为输入的常量，并不影响模型的计算，所以上式可以简化为：

.. math:: min H = \sum_{i<j,i∈V}\sum_{j∈V}w_{ij}\sigma_i\sigma_j

其中, :math:`H`\表示哈密尔顿量, :math:`w`\ 为输入的邻接矩阵,决策变量 :math:`\sigma_i`\表示顶点 :math:`i`\的分类,上述式子就是一个最大割问题的Ising模型.


建模代码
==========

输入矩阵
------------

.. code:: python

   import numpy as np
   import kaiwu as kw
 
   # Import the plotting library
   import matplotlib.pyplot as plt
 
   # invert input graph matrix 
   matrix = -np.array([
                   [0, 1, 0, 1, 1, 0, 0, 1, 1, 0],
                   [1, 0, 1, 0, 0, 1, 1, 1, 0, 0],
                   [0, 1, 0, 1, 1, 0, 0, 0, 1, 0],
                   [1, 0, 1, 0, 0, 1, 1, 0 ,1, 0],
                   [1, 0, 1, 0, 0, 1, 0, 1, 0, 1],
                   [0, 1, 0, 1, 1, 0, 0, 0, 1, 1],
                   [0, 1, 0, 1, 0, 0, 0, 0, 0, 1],
                   [1, 1, 0, 0, 1, 0, 0, 0, 1, 0],
                   [1, 0, 1, 1, 0, 1, 0, 1, 0, 1],
                   [0, 0, 0, 0, 1, 1, 1, 0, 1, 0]])

使用CIM模拟器查看量子比特演化
--------------------------
查看单次CIM模拟的量子比特演化过程用到了\ ``kaiwu.cim.simulator_core``\ 模拟器,
该模拟器输入CIM-Ising矩阵、pump功率、噪声强度、圈数、圈步长, 输出量子比特演化的模拟量, 需要注意的是, 使用这个模拟器前一般会使用\ ``kaiwu.cim.normalizer``\ 对输入矩阵的特征值进行归一化:

.. code:: python

   matrix_n = kw.cim.normalizer(matrix, normalization=0.5)
   output = kw.cim.simulator_core(
               matrix_n,
               c = 0,
               pump = 0.7,
               noise = 0.01,
               laps = 1000,
               dt = 0.1)

计算哈密顿量, 一般使用未进行归一化的矩阵进行计算:

.. code:: python

   h = kw.sampler.hamiltonian(matrix, output)

绘制量子比特演化图与哈密顿量图：

.. code:: python

   plt.figure(figsize=(10, 10))

   # pulsing diagram
   plt.subplot(211)
   plt.plot(output, linewidth=1)
   plt.title("Pulse Phase")
   plt.ylabel("Phase")
   plt.xlabel("T")

   # Energy diagram
   plt.subplot(212)
   plt.plot(h, linewidth=1)
   plt.title("Hamiltonian")
   plt.ylabel("H")
   plt.xlabel("T")

   plt.show()

.. image:: e1p1.png

查看最优解
查看最优解前需要将\ ``kaiwu.cim.simulator_core``\ 模拟器输出的数据使用\ ``kaiwu.sampler.binarizer``\ 进行二值化.
关于\ ``kaiwu.sampler.binarizer``\ 我们在本示例中只输入第一个参数表示输入解集, 
其余参数保持默认即可.

.. code:: python

   c_set = kw.sampler.binarizer(output)

最优解采样,
``kaiwu.sampler.optimal_sampler``\ 将对二值的解集进行按能量排序,
越靠前的解能量越低, 即解越优.

.. code:: python

   opt = kw.sampler.optimal_sampler(matrix, c_set, 0)

取最优解, 其中 opt = (解集, 能量), 下面第一个索引取解集,
第二个索引取解集第一个解.

.. code:: python

   best = opt[0][0]
   print(best)

模拟计算
--------

使用CIM模拟器进行计算
^^^^^^^^^^^^^^^^^^^^^

若像大规模计算, 我们用到了\ ``kaiwu.cim.simulator``\ 模拟器,
该模拟器输入CIM-Ising矩阵、pump功率、噪声强度、圈数、圈步长、归一化因子、独立运行次数, 输出量子比特经过去重的二值化值,
使用这个模拟器前就不需要使用\ ``kaiwu.cim.normalizer``\ 对输入矩阵的特征值进行归一化了:

.. code:: python

   output = kw.cim.simulator(
               matrix,
               pump = 0.7,
               noise = 0.01,
               laps = 50,
               dt = 0.1,
               normalization = 0.3,
               iterations = 10)

输出结果
--------

.. code:: python

   opt = kw.sampler.optimal_sampler(matrix, output, 0)
   best = opt[0][0]
   max_cut = (np.sum(-matrix)-np.dot(-matrix,best).dot(best))/4
   print("The obtained max cut is "+str(max_cut)+".")