# Graph Theory
###### tags: `GraphTheory`

- [Google 表單](https://docs.google.com/forms/d/e/1FAIpQLSdI5q_QOfxGXrC-DXSiZ1w4fIsGtTNGkk-5cwfIlY6OndRvtw/viewform)
- [其他組的 Paper](https://docs.google.com/spreadsheets/d/1rqvpjmUlqdxjt12ZBXTQGgKwWhIeyPx2nAyxiY8gNqk/edit#gid=96068209)
- 30~40 mins for each presentation, 5 mins for Q&A
- The paper you select must published in 2015~2021
- [CAD contest](http://iccad-contest.org/2021/tw/)
- [Paper 1 簡報共筆](https://docs.google.com/presentation/d/1rBfCTFF-9d4RHfFD0clj4T2PwCIhLKff8Cxp6WnMgCQ/edit?usp=sharing)
- [Cad 雲端](https://drive.google.com/drive/folders/1MTdTuuoYQxMV0McV_0kTKrgt5HBJVdW2?usp=sharing)
- [Paper II Report](https://hackmd.io/V_RwJzcpSTGUzKHJdqEM8A)
- [Cad Contest PPT](https://docs.google.com/presentation/d/18fZN017ETy_ACC6brOKtXVhd1nHXG7_oJya56ACcy7I/edit?usp=sharing)
- [Cad Contest report](https://hackmd.io/A7xfwCNoQbGhp8iK3iv54w)

## 重要日程
4/2: Online form filling: [Google 表單](https://docs.google.com/forms/d/e/1FAIpQLSdI5q_QOfxGXrC-DXSiZ1w4fIsGtTNGkk-5cwfIlY6OndRvtw/viewform)
4/9: Paper rejection announcement
4/16: Rejected group resubmit new paper

Online presentation
4/30、5/7: GT on Security Design and Verification
5/28、6/4: GT on Quantum Computing

## Security Design and Verification

- [x] [An Approach for Internal Network Security Metric Based on Attack Probability](https://www.researchgate.net/publication/324752954_An_Approach_for_Internal_Network_Security_Metric_Based_on_Attack_Probability)
- ~~[Efficient Dynamic Updates of Distributed Components Through Version Consistency](https://ieeexplore.ieee.org/document/7516718)~~
--在不損害安全性情形下提高動態更新服務效率，並通過圖形象化此方法，也有驗證此方法正確性
- [Efficient and Scalable Integrity Verification of Data and Query Results for Graph Databases](https://ieeexplore.ieee.org/document/8118163)
-此文章所建的圖形數據庫具有高安全性，且圖的edge note線性複雜度為最佳之外，計算速度還比普通數字簽呈快。

### An Approach for Internal Network Security Metric Based on Attack Probability
**mail address: `[e-mail]`**

Questions to ask:
Why internal network?

---
### Replenish

- About psudocode 1 IsAviailableMonitor:
- Where are EO sets came from? If the system's normal data flow are known,why should we have to create the attack path rather than just attack these normal path?
Attack action node set are comparing with many set of EO separatly. By the way, EO set are known just like the training data. This algorithm is to make sure the attck data flows don't against some known path that normal system pass through.
Also some of the normal data path are not known for hacker. As a "white attacker",we must find out all of posible path that can invade the system. Hence using known EO set for attacking are not enough at all.

- Assume AO={ba7,ba1,ba2,ba3,ba5,ba4,ba6}、EO={ba1,ba2,ba3,ba4}.If ba5、ba6、ba7 are somehow been offseted by at least one of ba1、ba2、ba3、ba4. Is that posible for this condition? Also, coverage rate would not be "1"?
The meaning of EO set is the "observation record" of the normal data flow. That is, if you had been through the offseted note, this action would be record whatever it been offseted or not. Thus, the circumstance described above are invalid.
- About calculation of posibility:
- Base on the difference of CVSS parameter and the calculation, why can we compared with reference[5]? Is that fair?
We think it's ambiguous too......
- How to calculate Key value pairs $P_{AS}$ K M N ？ Why use number "12" as $P_A$'s denominator ？
Maybe author hope to constrain posibility must less than 1. Because no one can garentee the node will "100%" been triggered.
- In figure 10, why delete attack node7 rather than node8?
There are same in that case because once node7 or node8 been triggered, the attack "success". Hence, we guess author tried to show this concept in this case.

1. 是否為好幾組 attack action node set 跟好幾組 EO 做比較，以確認這些 AO set 是否真正 match 實際情況？
2. EO是怎麼來的？是否已知？
3. 若EO為已知，為何不直接使用EO實現attack graph?以及如此一來attacker直接利用已知EO做攻擊即可？
4. 在procedure IsAvailableMonitor中，假設AO={ba7,ba1,ba2,ba3,ba5,ba4,ba6}、EO={ba1,ba2,ba3,ba4}，若ba5、ba6、ba7中至少有一個會跟ba1、ba2、ba3、ba4抵銷，則這樣coverage rate就不會是1，請問這種情況會發生嗎?

第二是與機率的計算方法有關
1. 與 reference [5] 使用之 cvss 參數以及 or 算法不同，為何可以進行比較，以及為何使用 max 概念做 or 運算?
2. Key value pairs 之 P_{AS} 的 K M N 參數從何算得？ P_A 之分母為何要除以 12？
3. 在 figure 10 中，我們發現 attack node a7 被刪掉了，但是我們以您的思路邏輯認為應該是 a8 應被刪除?
4. 另外我們也想詢問，為何 User6 及 Root 6 可以被簡化?

關於結論
1.

---

Two types of analysis method(in network security)
- unknown vulnerabilities: prevention measures.
- known vulnerabilities: repair and improve the security.

#### Attack graph
Building the actual network into a theoretical graph so as to simplify or idealize the problem.
- Cons: It cannot quantitatively describe the network security.

Proposed solution: Introducing the **cumulatice reachable probability**

#### An approach for internal network security metric based on attack probability

#### Temporal difference relationshop
(of atomic attack event / or other)

### Summary
- **yo**: Introduction: 該技術的應用領域介紹 (現有方法跟方式)
> As the rapid development of network and information technology, the role of information system in enterprise becomes important.
> The number of attacks from internal network has also increased, so it's necessary to build a internal network security system.
> Network security analysis method can be divided into two types:
>> 1. **unknown vulnerabilities**, considering the prevention measures.
>> 2. **known vulnerabilities**, repairing the weak parts and improve the security of whole network.

> Fixing unknown vulnerabilities is an effective but abstract way, and it's hard to implement.Since we will start from known network security vulnerabilities.
> For example, attack graph is a kind of graph theory method to judge the network security of known vulnerabilities.
> In this paper, we construct an attack graph model from known aspects to **simplify or idealize the actual elements**. Nevertheless, attack graph model cannot quantitatively describe the network security. In order to conduct quantitative analysis of the possible attacks, we introduce the cumulative reachable probability for each node.
> Above all, we combine attack graph model and attack probability to solve the network security problem.

- **yo、pei**: Related works:
- [3] 貝氏網路
> [3] is centered around near real-time security analysis such as intrusion response. By identifying the important types of uncertainty and using Bayesian networks, [3] capture what is really happening in the scope and severity level, possible consquences, and potential countermeasures.
> Three sources of uncertainty including:
>>1. ***the uncertainty on attack success***
>>2. ***the uncertainty of attacker choice***
>>3. ***the uncertainty from imperfect IDS sensors***

> Although [3] provides an effective way to combine attack graph and Bayesian network, it did not calculate the probability of uncertainty testing data in the final derivation process.
<!--
- 說什麼：
1. 考慮的點（deterministic $\rightarrow$ real-time）
2. IDS（這個跟monitoring event node很像）
3. CPT（這個跟CRP類似，都是計算攻擊成功的機率）
4. summary
- reference 3考慮attack success/attacker choice/IDS sensor的不確定性
- reference 3是首次想將BN和attack graph整合起來的paper
- paper作者認為reference 3沒有將不確定性的數值做量化
- 修正部分：
- Introduction
- deterministic(所有預先條件滿足攻擊才會發生), limited their use(consider real-time intrusion response)
- real-time security analysis is imprecise(compare to deterministic reason)
- intrusion detection system(IDS)
- false positives：答案是假的測出來是真的
- false negative ：答案是真的測出來是假的
- Bayesuan network(BN)
- cause and effect relationships
- Directed Acyclic Graph(DAG)
- strength of influence represented by conditional probability tables(CPT)
- BN很有用但是要從attack graph建構BN不容易（考慮現實因素）
- likelihood，CPT中的attack-choice很難被分離出來
- BN的建構方法
1. 圖形結構要基準化
2. CPT要可以從現實中自動計算
3. BN model不能太容易受CPT parameter影響
- Modularized CPT computation
- Noisy-And/Noisy-Or
- Evaluation Methodology
- comparing outputs with a Referee's(標準答案)
- BN-based tool所看到的會合ground truth有落差，因為IDS sensor's uncertainty
- 除了IDS uncertainty之外，errors in firewall rules/vulnerability reports也會造成ground truth's distortion
- pre-deployment
- 專注在pre-deployment vulnerabilities
- post-deployment
- 專注在post-deployment attack events
- 包含attack events還有good events
- How Sensitive is the BN-based tool to its CPT tables?
- BN-based tool受到CPT的影響應該要很小（穩定）
- 這篇paper只會同時考慮一個參數的影響（理由是多參數的穩定性太難做）
- attack graph model (three sources of uncertainty)
1. ***the uncertainty on attack success***
2. ***the uncertainty of attacker choice***
3. ***the uncertainty from imperfect IDS sensors***
- uncertainty testing data is not calculated
- Related work
- 他自吹說：
1. 利用了intrusion detectors
2. incorporate a holistic security analysis framework
- Conclusions
1. attack graph and Bayesian networks will not work nicely together
2. to use BN correctly $\rightarrow$ identify and represent relevant uncertainties
3. improve enterprise security analysis through BN
4. first effort to combine attack graph and BN-->
- [7] reference7中有關attack graph的定義與本篇類似，較為不同的是沒有考慮多個resource state指向同一個attack node的Or Relationship，以及多個attack node指向同一個resource state的And Relationship

- [8]
- [5] 實作方法
- **wendy**: Attack Graph
- Atomic attack
An atomic attack is the minimum set of basic actions that an attacker needs from a resource state node to another resource state node.
For example, There are two methods to open a word file.Method 1,open the file by double clicking it. Method 2, click on the file, and then click "Enter" to open it.
- Monitoring Event Node
This paper place monitoring event node to each source node for observating whether the node has been triggered or not. Whenever the monitoring event nodes capture the actions, these actions will be record into alarm log(attack evidence).
- Attack Graph
Attack graph is an directed acyclic graph which represent the relation of attacked source state nodes and the attack actions.Containing components are as follow:
- S : the set of resource state nodes that attacker has occupied.
- A : the set of attack action nodes
- O : the set of monitoring event nodes
- E : directed edges between nodes.
- P : the probability of attack
$P = P_A \cup P_{AS}$
$P_A$ the probability that attack action node being attacked(triggered)
$P_{AS}$ the success probability of attack action
- Directed Edges Relationship
- **AND** : Attacker needs to have all of the current sources to carry out the next attack
- **OR** : Attacker needs to have one of the current sources to carry out the next attack
- **pei**: Strategies to simplfy attack graph
figure 3: the schematic diagram
- introducing temperal difference relationship
- if the monitoring event node a~i~ is detected by the security monitoring system earlier than the monitoring event node a~j~ all the time, then we could say the monitoring events a~i~ and a~j~ have a temporal diference
relationship
- figure 2
- **wendy**: is available monitor (動畫): pseudocode I
- This procedure IsAvailableMonitor is to check wheather the attack evidences of a monitoring event nodes in the alarm clock trigger's order are correct or not.Assume AO is the set of attack evidences in alarm log.EO is the set of corresponding atomatic attack action.
- While the current attack evidence does match the current atomic action, means that the current attack evidence confirms the actual system dataflow. In the mean time, counter which record the matched pairs of current attack evidence and atomic attack action will added one. Otherwise, keep comparing the next atomatic action until all of element in the set of EO or AO has been checked. The calculation of coverage rate, which representing the completeness of attack evidence set matches actral system dataflow.
- **pei**: is temporal difference relationship: pseudocode II
根據monitor system所提供的alarm log，可觀察攻擊者的行為存在時間的先後關係，在attack graph中因此存在一些不可能發生的path，不利於計算和分析系統的安全性，因此需要將這些path刪掉，pseudocode2的作法為任選兩個monitor event node O1和O2，若O1和O2存在temporal difference relationship，會驅動O1和O2的action nodes之集合分別定義為Wo1和Wo2，在Wo1和Wo2中各挑一個action node A1和A2，若在攻擊圖中A2到A1存在一條path，則將違反temporal difference relation ship，該條path就要刪掉，
- Calculation of probability
- **wendy**: CVSS (Common Vulnerability Scoring System)
- Proposed by National Infrastructure Advisory Council(NIAC)
- A self-test scale measurement that evaluate the system's security and risk level
- AC(Access Complexity):
High(Indicator value:1) : Exist highly specific access condition.For example, if an attacker has to attacked DNS and have more higher privilige to cheat other host, and these complexible condition we called that a high access complexity.
Medium(Indicator value:2) : Needs a little bit specific access condition.That is, attacker must to some extent's privilige to attack system.
Low(Indicator value:3) : Don't need specific access condition to attack system.
- AU(Authentication):
Multiple(Indicator value:1) : Exploring this vulnerability required attacker authenticate two or more times.Even if the same credentials are used each time.For example, Attacker authenticating to an operating system in addition to providing credentials to access an application host on that system.
Single(Indicator value:2) : Exploring this vulnerability required attacker authenticate just one time.
None(Indicator value:3) : Exploring this vulnerability does not required authentication.
- AV(Access Vector):
Local(Indicator value:1) : Only physical access(for example: sudo,COM port)can activate attacked nodes
Adjacent(Indicator value:2) : Only local IP subnet, Bluetooth, local Ethernet access can activate attacked nodes
Network(Indicator value:3) : Remote access(which do not require local access)can activate attacked nodes
- Key-value pairs : Divide simplified attack graph into the form of "resource state node, attack action, and next resource state node."
- probability of Key-value pairs: $𝑃 = (𝑃_𝐴, 𝑃_{AS})$
- $P_A$ : probability of attack action
$P_A = (V_s + H_i)/12$.
$H_i$ : the weight where the host is located
Office network(Indicator value:2) : Host is located on the office network, easy to be used.
Core network(Indicator value:1) : Host is located on the core network, difficult to be used

$V_s$ : pair's own probability(the posability that attackers triggered such a pairs)
- $P_{AS}$ : probability of successful attack
$P_{AS} = K + M + N$
K : indicating whether the vulnerability information is published
K = 0.1,The vulnerability has been issued
K = 0 ,The vulnerability has not been issued
M : indicating whether the atomic attack method or step is currently available.
M = 0.4,the vulnerability has a detailed attack step scheme
M = 0.2,simple attack scheme
M = 0 ,extremly not available
N : indicating whether the attack tools are required in the atomic attack
N = 0.4,vulnerability is not required to use the attack tools
N = 0.2,vulnerability exploited needs to use the attack tools and the corresponding attack tools are available
N = 0,need to use the attack tools, but there are no available attack tools
[CVSS_webside](https://nvd.nist.gov/vuln-metrics/cvss/v2-calculator?calculator&.0#score)
- Cumulative Reachable Probability
- **AND** : need both resource/attack action nodes been activated
<a href="https://imgbb.com/"><img src="https://i.ibb.co/dWg0Ddw/and.png" alt="and" border="0"></a>
$P_{AC} = P_{S1} × P_{S2}× P_A$
- **OR** : need one of resource/attack action nodes been activated
<a href="https://imgbb.com/"><img src="https://i.ibb.co/nD0k5vp/or.png" alt="or" border="0"></a>
$P_{AC} = (P_{S1} + P_{S2}) × P_A$
:::info
關於那個特殊的 operator 這邊好像有說明 [An Approach for Security Assessment of Network Configurations using Attack Graph](https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=5383993) 但是不知道是不是同一個
:::
- **yen、yo**: Experiments, comparing with [5]
- the difference compared with [5]
- The cumulative model.
In [5], the author use $p(t)$ to represent the probabilities, and use $P(t)$ to represent the cumulative probabilities.
==The individual probabilities==: The attacking event probability of any atomic attack:
$$
p(t)=AV*AU*AC
$$
(為什麼是乘法?????, P.S. 在原本的 paper 裡面是加法)

:::info
- $t \in T$, $T$ is the set of atomic attacking nodes ($\approx a$ in this paper)
- $c \in C$, $C$ is the set of states ($\approx S$ in this paper)
- $C_0$: initial state set
- $C_d$: reachable state set
:::
==The cumulative probabilities==
1. Initial nodes: $\forall c \in C_0, P(c) = p(c)$
2. Other nodes: $\forall c \in C_d, P(c) = p(c)*(P(t_1)+P(t_2)-P(t_1)*P(t_2))$, where $t_1, t_2 \in pre(c)$
- The $p(c)$ means the probability of the state $c$ to be obtained.
- The $P(t)$ means the probability of the attacking event $t$ to be performed.
4. The attack event happening probabilities: $\forall t \in T, P(t) = p(t) * \prod_{c \in pre(t)}P(c)$ (為什麼不是 $P(c)$ 的聯集????)
:::info
Assuming that the attacking event $t$ requires all the condition $c$, that is,
$$
P\left(\bigcap_{c \in pre(t)} c_i\right) = \prod_{c \in pre(t)} P(c_i)
$$
:::
- Atomic attacking risk $R(t)$ (原子攻擊影響???)
1. The influence of attacking event $t$ : $M(t)$
$$
M(t)=10.41* \Big(1-(1-CI)(1-II)(1-AI)\Big)
$$
2. The atomic attacking risk $R(t)$
$$
R(t) = P(t)*M(t)
$$
3. The attacking risk of the entire graph
$$
R(all) = \sum_{t \in T}P(t)*M(t)
$$
- environments
- topology (fig. 8)
- 用圖去表示 topology
- initial
- simplified
- The procedure to simplify the attack graph

        P1={Root0, a1, User1, a3, User3, a7, User4}
        P2={Root0, a1, User1, a4, User4, a8, User3}
        P3={Root0, a2, User2, a5, User3, a7, User4}
        P4={Root0, a5, User2, a6, User4, a8, User3}

We can find a bijection function to map the different path into a unique form by introducing the temporal difference relationship. The author did this by deleting attack nodes and path. So the resulting attack path can be expressed as:

        P1={Root0, a1, User1, a3, User3, a7, User4}

- Inspect the experiment
- Solve the probability of the root node
$P_{S0}$ can be obtained by solving the following equations:
$$
\left\{\begin{matrix}
P_{S1}=p_{AC1}\times P_{AS1}=(P_{S0}\times P_{A1}\times P_{AS1}) \\
P_{S2}=p_{AC2}\times P_{AS2}=(P_{S0}\times P_{A2}\times P_{AS2})
\end{matrix}\right.
$$
We found that $P_{S0} = 1$.
We calculated the probability manually to check the consistency of our result and the author's result. However, the auther doesn't explain the (+) operator that they adopted clearly, so we can only guess the operation by reverse engineering the author's calculation result. The following equation can be inferred:
$$
P_1 \oplus P_2 = \max(P_1, P_2)
$$
Nonetheless, the operator is quite confussing. We lookuped some papers such as the author's reference in [4], but the $\oplus$ operator defined in [4] is $P_1 \oplus P_2= P_1 + P_2 - P_1 \cap P_2$
- Method by the author
- User 3, User 4
$$
\begin{aligned}
P_{S4}&=(P_{S1}*0.67*0.3) \oplus (P_{S2}*0.67*0.3) \\
&=0.116781 \oplus 0.09246 \\
&=\max(0.1167, 0.09246) \\
&=0.1167
\end{aligned}
$$
$$
\begin{aligned}
P_{S3}&=(P_{S1}*0.67*0.3) \oplus (P_{S2}*0.67*0.3) \oplus (P_{S4}*0.58*0.5) \\
&=0.116781 \oplus 0.09246 \oplus 0.03386\\
&=\max(0.1167, 0.09246, 0.03386) \\
&=0.1167
\end{aligned}
$$
- User 5, Root 5
$$
\begin{aligned}
P_{S5}&=(P_{S3}*0.83*0.7) \oplus (P_{S4}*0.83*0.7) \\
&=0.06784 \oplus 0.06784 \\
&= 0.06784
\end{aligned}
$$
$$
\begin{aligned}
P_{SR5} &= (P_{S5}*P_{A5})*P_{AS5} \\
&= 0.0395
\end{aligned}
$$
- Matching the author's results in table 5:
- Method by [4]
- User 3, User 4
$$
\begin{aligned}
P_{S4}&=(P_{S1}*0.67*0.3) \oplus (P_{S2}*0.67*0.3) \\
&=0.116781 \oplus 0.09246 \\
&=0.1167 + 0.09246 - 0.1167 * 0.09246 \\
&=0.19844
\end{aligned}
$$
$$
\begin{aligned}
P_{S3}&=(P_{S1}*0.67*0.3) \oplus (P_{S2}*0.67*0.3) \oplus (P_{S4}*0.58*0.5) \\
&=0.116781 \oplus 0.09246 \oplus 0.03386\\
&=0.1167 + 0.09246 + 0.03386 \\
& \quad - 0.1167 * 0.09246 - 0.09246 * 0.03386 - 0.1167 * 0.03386 \\
& \quad + 0.09246 * 0.03386 * 0.1167\\
&=0.24457
\end{aligned}
$$
- User 5, Root 5
$$
\begin{aligned}
P_{S5}&=(P_S3*0.83*0.7) \oplus (P_S4*0.83*0.7) \\
&=0.142095 + 0.11529 - 0.142095 * 0.11529 \\
&=0.21729
\end{aligned}
$$
$$
\begin{aligned}
P_{SR5} &= (P_{S5}*P_{A5})*P_{AS5} \\
&= 0.12624
\end{aligned}
$$

- **yen、yo**: Conclusion
- Assumption (Domain of this research): internal network
- Compare with [5]
- improvement from [5] to this paper
- fig. 11
- use monitoring event node and temporal difference relationship to simplify the attack graph.
- has a significant improvement in efficiency due to the smaller graph size
- propose the calculation method to cumulative reachable probability.
- having more smooth fluctuation of cumulative probability compared with [5]

- **yen、yo**: Issues
> In this paper, we discover several issues which author didn't mention and can't be explain clearly. The following are those issues:
>> 1. In fig.9, we can't explain why Root 6 and User 6 were deleted. Since author had mentioned that the attack would terminate when it reaches User 5 or User 6, considering either User 5 or User 6 is acceptable and more effective.
>> 2. From fig.9 to fig.10, directed edge a7 was deleted and we can't give this action an appropriate reason. Because of edge a7 conforms temporal difference relationship we have mentioned, there's no reason to delete this path. We believe that deleting edge a8 and keep edge a7 is more properly.
>> 3. The CVSS parameters used in this paper is different from the parameters on the CVSS official website. Futhermore, this paper did not mention the precise CVSS parameters.
>> 4. Conclusions in this paper noted that it had a better performance when comparing to reference [5]. However, both reference [5] and this paper didn't mentioned the enviroment of experiment. We couldn't confirm whether the results are convicible.
<!--
//fig. 9, Root 6、User 6 為什麼刪掉了
//- CVSS 的定義跟別人不一樣
//- fig. 9, 為什麼刪掉 $a_7$
//- fig. 11, node 的生成方法不明確
//- table 4, 沒有說明他的 CVSS 參數選擇(跟環境, 漏洞有關)-->
### Reference
- [reference 3](https://www.researchgate.net/publication/220957917_Using_Bayesian_Networks_for_Cyber_Security_Analysis)
- [reference 5](http://caod.oriprobe.com/articles/32076521/Network_security_assessment_based_on_probabilities_of_attack_graph_nod.htm)
- [reference 8](https://ieeexplore.ieee.org/document/5936075)

### QA
1. 怎麼知道 monitor event node 對應的 node (如何配置)
2. 為什麼可以用 Cumulative probability 去分析 attack graph
3. 本篇的 attack graph 為什麼長這樣

# Quantum Computing

### ~~[Graph-theoretic Simplification of Quantum Circuits with the ZX-calculus](https://quantum-journal.org/papers/q-2020-06-04-279/pdf/)~~
- Authors:
Hayato Ushijima-Mwesigwa, Christian F. A. Negre, Susan M. Mniszewski
- Publication:PMES'17: Proceedings of the Second International Workshop on Post Moores Era Supercomputing, November 2017, Pages 22–29, https://doi.org/10.1145/3149526.3149531
- quantum circuit optimisation based on the ZX-calculus (文中有介紹啥是 ZX-calculus)
- ZX-diagrams (describing quantum computations graphically)
- show that the resulting reduced diagram can be transformed back into a quantum circuit
- [ ] 看起來很多像是證明跟推導的部分

### ~~[Graph Partitioning using Quantum Annealing on the D-Wave System](https://arxiv.org/ftp/arxiv/papers/1705/1705.03082.pdf)~~
2019, cite 8
- Quantum Graph Neural Networks (QGNN), representing quantum processes which have a graph structure, and are particularly suitable to be executed on distributed quantum systems over a quantum network
- From classical neural network to Quantum Graph Neural Network

## [Graph Partitioning using Quantum Annealing on the D-Wave System](https://arxiv.org/ftp/arxiv/papers/1705/1705.03082.pdf)
- Authors: Hayato Ushijima-Mwesigwa, Christian F. A. Negre, Susan M. Mniszewski
- Publication:PMES'17: Proceedings of the Second International Workshop on Post Moores Era Supercomputing, November 2017, Pages 22–29, https://doi.org/10.1145/3149526.3149531

### ==[雲端](https://drive.google.com/drive/folders/1-oAwOAZKSJuXx7f4TpeeCGkot3RoS22b?usp=sharing)== ==[PPT](https://docs.google.com/presentation/d/1ZZRAAKep4B5pwMrUzI8DvlMjLEjrIPbzFHQTWLcjbJs/edit?usp=sharing)==

### 名詞解釋參考：
#### Quantum
- Hamiltonian problem:
- quantum annealing:use quantum-mechanical property such as tunneling and entanglement to minimize energy-based models(such as Ising model、QUBO model)
- markov chain:下一狀態的機率分布只能由當前狀態決定，在markov chain 的每一步，可以從一個狀態轉變到另一個狀態，也可以保持當前狀態
- Ising model、QUBO model:provide the objective function for quantum annealing
- Ising objective function:
$O(h,J,s)=\sum_{i} h_is_i + \sum_{i<j}J_{ij}s_is_j$
- $h_i$: certain external forces applied to individual particles
- $J_{ij}$: interaction forces between grid neighbors，if particular i and j are not adjacent we assume $J_{ij}$ = 0
- $s_i$: magnetic dipole moments of atomic "spins"that can be in one of two states (+1 or −1)
- [Ising problem](https://en.wikipedia.org/wiki/Ising_model): A mathematical model of ferromagnetism in statistical mechanics. The "spins" $s=\{-1, +1\}$
- the system tends to the lowest energy but heat disturbs this tendency, thus creating the possibility of different structural phases.
- [Chimera graph](https://support.dwavesys.com/hc/en-us/articles/360003695354-What-Is-the-Chimera-Topology-)
- The connectivity graph of the D-wave system (a kind of sparse graph))

#### GP (Graph Partition)
- [modularity](https://en.wikipedia.org/wiki/Modularity_(networks)): measures the strength of division of a network into modules (communities)
- hight modularity: dense within same module and sparse between different modeules
- [Laplacian matrix (of a graph)](https://en.wikipedia.org/wiki/Laplacian_matrix)
- Graph Laplacian, Harmonic matrix(調和矩陣), Admittance matrix, Kichhoff matrix, Discrete Laplacian,
- 是一種用矩陣來表示圖的方法,
- the Laplacian matrix is indeed the ‘discrete’ version of the Laplacian operator over graphs, which means the divergence of the gradient of a graph. [A good article](https://www.quora.com/Whats-the-intuition-behind-a-Laplacian-matrix-Im-not-so-much-interested-in-mathematical-details-or-technical-applications-Im-trying-to-grasp-what-a-laplacian-matrix-actually-represents-and-what-aspects-of-a-graph-it-makes-accessible/answer/Muni-Sreenivas-Pydi)
- That is, the Laplacian of a graph function determines how “smooth” the graph function is.
- 
$L$: Laplacian matrix
$D$: Degree matrix
$A$: Adjacent matrix
- Note, the Dirichlet energy form
$$f^TLf=\frac{1}{2}\sum_{u, v \in G}w(uv)(f(u)-f(v))^2$$
$w(uv)$ is the weight of edge $uv$, if there is no connection, the weight is 0.
if $f(u)$ is defined as which part node $u$ belongs to, that is
$$f(u)=\left\{
\begin{array}
+1, & part A \\
-1, & part B
\end{array}
\right. =s_u$$
- cut-edge
- A cut edge is defined as an edge whose end points are in different parts.

### Q
- ==yen==why $min(s^TLs)$ can be equivalent to $max(s^TAs)$?
- ==yen, yo==how does the $\alpha$ and $\beta$ be chosen?
- ==yo==explain the equation (11)
- ==yin==how to use the complement of the graph to handle a dense graph
- ==pan==What is the purpose of modularity with thresholding (2.1)
- Hint: Preprocessing of the modularity matrix to produce a **sparse version** that is still **representative** of the original allows for larger problems to be solved efficiently by community clustering using quantum annealing. ???????????????????????
- ~~table 3 graph??~~
- what is super edge、super node?
---

### 分工
- Introduction
- ==Wendy==Application
- ==Wendy==Purpose
- Background knowledge
- ==pan==Quantum annealing
- ==pan==Ising model and QUBO
- Purposed method
- ==pan==community detection use modularity
- ==yen==binary partition
- ==yen==k-concurrent partition
- ==yo==Result
- ==yo==Conclusion

---

### Representation structure

#### Introduction

- Some aplications about graph partitioning(GP):

- Social Network : PageRank
A PageRank vector is a weighted sum of the pobability distributions(click into this website) obtained by taking a sequence of random walk steps(click the seperlink inside original websites) starting from a specified initial distribution(original webside).
The weight placed on the distribution obtained after t walk steps decreases exponentially in t, with the rate of decay determined by a parameter.
GP mathed could reduce the complexity of PageRank vector community and finding a rank and importance of a web page.
<a href="https://ibb.co/7rDWp8Y"><img src="https://i.ibb.co/TPG80pv/Websites-interlinking-to-illustrate-Page-Rank-percents.png" alt="Websites-interlinking-to-illustrate-Page-Rank-percents" border="0"></a>
- VLSI design
Using GP method to partite circuit into smaller component in order to reduce the complexity of VLSI design and decrease the number of connection among those components(that is, decrease cut edges).
- Biological Networks: protien-protien interaction
Using GP to recognize the active or effector site in comarison to original protien in order to exchage or combine different component to others(That is the way that emerge allosteric mechanism of protein)
<a href="https://ibb.co/tLy62gf"><img src="https://i.ibb.co/5Kq3jmD/2021-06-02-181605.png" alt="2021-06-02-181605" border="0"></a>)
- Purpose:
Graph partitioning applications are ubiquitous through out many applications. The purpose of the paper is to explore graph partitioning problem using quantum annealing on the D-wave machine. In that case, we can reducing the calculation of the density matrix(original large graph) into smaller subsystem(partitioned smaller graph) rendering the calculation more computationally efficient.
#### Background knowledge
- Quantum computer
A quantum computer is a type of computer that uses quantum mechanics so that it can perform certain kinds of computation more efficiently than a regular computer can.
The method it works is to minimize the objective function(Quantum annealing,explained at next below) to reach the minimum energy summation of the system(Global best solution).

Why use quantum computer?
GOOGLE announce that thier quantum computer(Sycamore 無花果量子計算機)completing a complexible problem with only 200 seconds. Compared with the fastest computation speed in the world("Summit "from IBM)which needs more than 10000 years to complete the same problem.

Because of the high speed calculation of quantum computer, author gives an method to fit the GP problem into quantum machine to inhance computing efficient.
- Qubit
Each elementary particles would have a spin(kinds of direction of angular momentum), which depends on the observation direction we decides. If the spin direction equals to the observation's, than the spin of particle would be "+1", on the contrary it would be "-1". In D-wave machine, we called these directions be "state".
<center><img src="https://i.ibb.co/bPh7k5m/spin1.png" width="500px"></center>
<center><img src="https://i.ibb.co/4FnsPFW/spin2.png" width="500px"></center>

Because qubit is quantum object, it would be in a superposition of state "-1" and "+1" could be at the same time. To begin, there is just one valley, with a single minimum.

<a href="https://ibb.co/Fzqn7Pq"><img src="https://i.ibb.co/M759NX5/superposition1.png" alt="superposition1" border="0"></a>

After quantum annealing runs, the barrier is raised and energy diagram turned to 2 valley.The qubit end up one of these valleys at the end of the anneal.
We can control the probability the valley which qubit would drop in by applying external magnetic field called "bias".

<a href="https://ibb.co/FnBJJkm"><img src="https://i.ibb.co/LvC66DY/superposition3.png" alt="superposition3" border="0"></a>

Each qubit would influence every others. Apply to graph, we simply giving an interaction force between neighbers, these pairs called "couplers". Otherwise set thier interaction force to 0.(That is, assume two qubits are far away).
For a graph system, there is a set of qubit(node). Quantum machine would track all permutation of the qubit state to minimize the system energy.
<a href="https://ibb.co/kXth2VX"><img src="https://i.ibb.co/RcqCTGc/superposition4.png" alt="superposition4" border="0"></a>

- Quantum annealing:
QA is a combinatorial optimization technique meant to exploit quantum-mechanical effects to minimize and sample from energy-based models(such as Ising or QUBO model),output of QA is a low energy ground state

- Ising model:
$$O(h,J,s) = \displaystyle\sum_{i}h_is_i + \displaystyle\sum_{i<j}J_{ij}s_is_j$$
above is the objective function of Ising model which
$h_i$ :certain external forces applied to individual particles
$J_{ij}$:interaction forces between grid neighbors, if particular i and j are not adjacent we assume $J_{ij}$=0
$s_i$ :magnetic dipole moments of atomic "spins"that can be in one of two states (+1 or −1)

- QUBO model:
$$O(Q,x) = \displaystyle\sum_{i}Q_{ii}x_i + \displaystyle\sum_{i<j}Q_{ij}x_ix_j$$
above is the objective function of QUBO model which
$Q_{ii}$ :is an analog to the Ising $h_i$
$Q_{ij}$ :is an analog to the Ising $J_{ij}$
$x_i$ : $s_i$ = $2x_i - 1$
$x_i$ is related to $s_i$ in Ising model, so $x_i$ can be one of two states (1 or 0)
- D-wave system:
computers implementing specialized quantum annealing

#### Flow

#### Methods
##### Community detection use modularity
Our goal is high intraconnectivity and low interconnectivity .This division of a graph into communities differs from the usual GP problem in that the size of the communities are not normally known in advance.
- community
When GP is unconstrained (with no limitation on partition size) and communication volume is minimized, the resulting natural parts are called communities.
- modularity matrix
denoted as $B$ ,is a NxN(N is number of nodes) matrix, and $B_{ij} = A_{ij} - k_ik_j/2m$, which $A_{ij}$ is an element in adjacency matrix； $k_i$ 、 $k_j$ are degrees of vertices i and j； $m = \frac{1}{2}\sum_i{k_i}$ is the total number of edges. $k_ik_j/2m$ is the expected number of edges between vertices i and j.If $B_{ij} >0$,then the number of edges between the node pair is more than we expected by chance
- modularity
One metric that can quantify the quality of a community structure.Performs a comparison of the connectivity of edges within communities with the connectivity of a random network where edges are placed randomly (with same number of nodes and edges in original one).If the number of edges between groups is significantly less than we expect by chance, or equivalent if the number within groups is significantly more, then we can say it's high-qualitied community structure.

- Graph clustering
The modularity Q is given by the sum of $A_{ij} - k_ik_j/2m$ over all pairs of vertices i,j that fall in the same group.First of all, we cluster a graph into two groups, labeled as -1 and +1,and the modularity is
$$Q = \sum_{ij}(A_{ij} - k_ik_j/2m)(s_is_j+1) = \sum_{ij}(A_{ij} - k_ik_j/2m)(s_is_j)$$
in matrix form modularity is:
$$Q(s)=S^TBS$$
the problem of finding the optimal modularity requires solving $max_s(Q(s))$, or $min_s(−Q(s))$.In order to achieve the goal, we need to formulate equation above as an Ising objective function, and $J_{ij}$ and $h_i$ in Ising will be set to $-2B_{ij}$ and $-B_{ii}$ to match the equation；In doing so, we can use D-wave system to find a high-qualitied community structure.

:::info
==prove $k_ik_j/2m$:==
The term $k_ik_j/2m$ is the expected number of edges between vertices i and j if edges are placed at random.The following is the derivation process of it.It will use the concept "half-edge",$k_i$(the degree of i) stands for the number of "half-edges" starting at i. Then, to make a full edge, you need to bind two half-edges in the graph. If you did it by picking half-edges uniformly, you would find that the probability of a randomly-picked edge being between i and j is $k_ik_j/4{m^2}$ + $k_jk_i/4{m^2}$ = $k_ik_j/2{m^2}$ ,because the probability of picking an half-edge of node i is $k_i/2m$. Define $X_k$ is a Bernoulli random variable describing whether the kth edge we picked is an edge between vertices i and j,then $S = \sum_{1}^{m}{X_k}$ is the number of edges between vertices i and j,so $E(S) = m * k_ik_j/2{m^2} = k_ik_j/2m$ is proved.
==prove $\sum_{ij}(A_{ij} - k_ik_j/2m)(s_is_j)$:==

:::

##### graph clustering
###### Background knowledge
- Laplacian matrix
The laplacian matrix has several name such as Harmonic matrix(調和矩陣), Admittance matrix, Kichhoff matrix, Discrete Laplacian. The physical meaning of laplacian of a graph is similart to the laplacian operator in continuous domain. It is indeed the ‘discrete’ version of the Laplacian operator over graphs, which means the divergence of the gradient of a graph. That is, the Laplacian of a graph function determines how “smooth” the graph function is.
Laplacian matrix is defined as follow:
$$
L_{i,j}=\left\{
\begin{array}
deg(v_i) & if\ i=j \\
-1 & if \neq j \ and \ v_i \ is \ adjacent \ to \ v_j \\
0 & otherwise
\end{array}\right.
$$
Or, in the matrix form:
$$
L=D-A
$$
where $L$ is the Laplacian matrix; $D$ is the degree matrix; $A$ is the adjacent matrix.
- Dirichlet energy
We can get more insight from the Dirichlet energy form:
$$f^TLf=\frac{1}{2}\sum_{u, v \in G}w(uv)(f(u)-f(v))^2$$
where $w(uv)$ is the weight of edge $uv$, if there is no connection, the weight is 0. if $f(u)$ is defined as which part node $u$ belongs to, that is
$$f(u)=\left\{
\begin{array}
+1, & part A \\
-1, & part B
\end{array}
\right. =s_u$$
then we can obtain the number of the edges that go across different group by calculating $f^TLf$, therefore, we can find the minimal cut-edge by $\min_f(f^TLf)$.

###### Binary partition

The D-wave system will test different composition of $s$ with quantum annealing algorithm and converge to the optimal answer in the end. Finally, we can get the best solution of the partition by the output $s$.

Our goal is to describe the graph partition problem as a QUBO form, so we can perform the calculation by the D-wave system.

The objective is to find the minimal number of cut-edge, so we group the input graph vertexs into two parts, +1 and -1. After previous introduction of Laplacian matrix, we can write down the objective function, $\min_s(s^TLs)$, where $s$ is the vector with size $N$ indicating the group of tht vetex belonging.

By $L=D-A$, we can simplify the original function by
$$
\begin{align}
\min_s(s^TLs) &= \min_s(s^T(D-A)s) = \min_s(s^TDs - s^TAs)
\end{align}
$$
Note that $D$ is the degree matrix, which is a diagonal matrix, and the physical meaning of $s^TDs$ is twice the number of edges in the given graph. We can inspect it by each term: $s_i D_i s_i = D_i s_i^2 = D_i$, and after the production being performed, the summation is $\sum_i D_i=2m$

In the same graph, no matter how the partition it is, the degree won't be affected, therefore, the term $s^TDs$ is constant. And we can rewrite the original objective function as
$$
\begin{align}
\max_s(s^TAs)
\end{align}
$$
which objective has been targeted to twice the difference of the number of edges($s^TDs$) and fourth the number of cut-edges($s^TLs$).

Next, we need to transform to the function into the QUBO form, which has no constraint and has all variables in the set of $\{0, 1\}$. we introduce the penalty $\alpha (\sum s_i)^2$ to the objective function:
$$
\max_s(\beta s^T A s-\alpha (\sum s_i)^2), s\in \{-1, 1\}^n
$$
By $(\sum s_i)^2=\sum si^2 +2\sum s_is_j=s^T1_{n\times n}s$ and $s=2x-1_n, x\in \{0, 1\}^n$, we can further obtain the following final result:
$$
\min_x\left({x^T (\alpha 1_{n \times n} - \beta A)x - x^T (\alpha 1_{n \times n} - \beta A) 1_n}\right), x\in \{0, 1\}^n
$$

About the choice of the $\alpha$ and $\beta$, the author referred from *Ising formulations of many NP problems* by Andrew Luca. It stated the choice would be determined according to the importance of the two term, that is the objective and the constraint, however, there exists some restrictions on the two terms. Let us look into the original objective: $\max(\beta (H_B + \Delta H_B) - \alpha (H_A + \Delta H_A))$, where $\alpha, \beta >0$, $H_B$ is the term of objective, $H_A$ is the term of constraint and the $Delta$ terms indicate the difference caused by any vertex moving. We consider the condition that the minimal constraint might have the maximal objective contribution, and the constraint term should not be sacrificed. The overall difference $\beta H_B - \alpha H_A$ should not be larger than 0, therefore, $\alpha \Delta H_A > \beta \Delta H_B$, and
$$
\frac{\alpha}{\beta}>\frac{\Delta H_B}{\Delta H_A}
=\frac{\min(\Delta, N/2)}{4}
=\frac{\min(2\Delta, N)}{8}
$$

The minimum is intuitive, considering the perfect balance condition such as $\{1, 1, -1, -1\}$, and we toggle any bit: $\{1, 1, 1, -1\}$, then we get $(\sum s)^2=4$.

To obtain the maximum $\Delta H_B$, we consider two conditions:
1. moving vertex from larger group to smaller group
The worst case is that all vertexes in the larger group have an edge with $u$ and no edge with the smaller one, so the $\Delta H_B$ is the maximal degree, $\Delta$.
2. moving vertex from smaller group to larger group
The worst case is that all vertexes in smaller group have an edge with $u$ and no edge with the larger one, so the $\Delta H_B$ is the maximal size of the smaller group, $N/2$.

Considering that the quantum annealing algorithm would evolve to the optimal case, so the maximal of $\Delta H_B$ is $\min(\Delta, N/2)$.

###### k-concurrent partition

Once we understand how the binary partition works, we can extend the concept to k-concurrent partition by considering that each vertex can be or not be group $k_j$, for $k_j \in k_j,..,k$.

Firstly, we need to have k copies of the original graph and group the same vertexes on single super-node, $v_i$. There are N black lines, connecting the posibilities of the vertexes to be group $k_j$, which indicating the relationship of the existence of single vertex in $V(G)$. There must be only one vertex $x_{i,j}$ representing the vertex in the original graph.

<center><img src='https://i.imgur.com/b8K4KCr.png =300x' style='width: 300px'></center>

Therefore, one of the constraints is that the summation of the probability of the same vertexes should be 1:
$$
\sum_{j=1}^k x_{i,j} = 1
$$

After performing the above procedures, we can obtain $n$ super-nodes, where $n=|V(G)|$.

Secondly, we create super-edges, connecting those vertexes $x_{i,j}$ that belonging to $k_j$ groups, and if we consider the perfect balancing case, the sum of those vertexes should be $n/k$.

$$
\sum_{i=1}^n x_{i,j} = \frac{n}{k}
$$
Without this constraint, the partition result might approach to the extreme solutions, such as single vertex in a group, and only one large partition exist.

Finally, the objective function is the minimization of the number of cut-edges. By the conclusion of the former discussion, we have known that the number of cut-edge between group $k_j$ and other group is $\min_{x_j} (x_j^T L x_j)$, which can also be visualized in the above figure, by summing up those vertexes on super-edge $j$. And the overall objective function is the summation on each group:

$$
\sum_{j=1}^k\min_{x_j} (x_j^T L x_j) = \min_x\sum_{j=1}^k (x_j^T L x_j)
$$

subjects to
$$
\sum_{j=1}^k x_{i,j} = 1 \\
\sum_{i=1}^n x_{i,j} = \frac{n}{k}
$$

with domain
$$
x_{i, j} \in \{0, 1\},\quad j=1,...,k,\quad i=1,...,n
$$

In order to fit the requirements of D-wave system, we have to transform the above optimization problem into QUBO form. so we replace the constraints with penalties:
$$
\min_{x_{i,j}}\left(\beta(\sum_{j=1}^k (x_j^T L x_j)) + \sum_{j=1}^k \alpha_j (\sum_{i=1}^nx_{i,j}-\frac{n}{k})^2 + \sum_{i=1}^n \gamma_i (\sum_{j=1}^k x_{i,j}-1)^2 \right)
$$

#### Result and Discussion
Our experiments on D-Wave machine used software tools such as $sapi$ Python and the hybrid classical-quantum $qbsolv$. We used $sapi$ for small graphs(up to 70 vertices), and for larger graphs we used hybrid classical-quantum $qbsolv$(up to 9000 vertices).

NetworkX was used for generating and processing graphs. Random graphs model, large graphs from the Walshaw Archive, and QMD molecule electronic structure graphs were used to evaluate our methods.

We compared the result of our methods to existing multi-level GP frameworks, METIS and KaHIP family (including KaFFPa and KaFFPaE). The quality of GP was evaluated by a comparison metric as the number of $cut \space edges$ between partition (smaller is better).

CD(Community Detection) with Thresholding
- In this case we applied community detection with thresholding to the karate club graph(34 vertices) and used $qbsolv$ to solve the problem. We set threshold value to 0.12 which means $B_{i,j}=0$ if $B_{i,j}<0.12$. Using this method, modularity reduced but the number of edges reduced significantly.
>>
>> The number of edges has direct consequence on the number of qubits/couplers required on D-Wave machine. Preprocessing of the modularity matrix produced a sparse version that still representative, also enhanced efficiency.
>> 

> Result of 2-Partitioning using $sapi$
>> Table2 shows the result of partitioning at most 70 vertices graphs(random graphs) using $sapi$ Python. The largest complete graph that is fully embeddable on the D-Wave 2X machine has approximately 45 vertices. However, partitioning up to 70 vertices is possible due to our use of the graph complement in the QUBO formulation.
>>
>>Our results gave solutions with a comparable quality to METIS and KaHIP(KaFFPa). Since we used $sapi$ as interface, 100 percent of the partitioning was carried out by the quantum annealer.
>>

> Results of 2-Partitioning using $qbsolv$
>> We partitioned the graphs using the hybrid classical-quantum $qbsolv$, which and partition graphs with over 100 vertices (up to 9000 vertices). Our solutions gave high quality partitions, in all cases having a cut size smaller than METIS and comparable to KaHIP. For the last two graphs, our solutions matched the best known results.
>>

>Results of $k$-Concurrent Graph Partitioning
>> $k$-Concurrent GP requires that a graph of size $N$ (number of nodes) $\times$ $k$ (number of partitions). This quickly uses up the available qubits.
>>
>> Table4 shows the results of running conccurent GP on small random graphs using $sapi$. The limitation of running directly on the D-Wave machine using $sapi$ is ~45 nodes. A 15 nodes graph partitioned into 4 parts used almost all available qubits. Similarly for 20 nodes split into 3 parts.
>> 
>>
>>
>> $k$-Concurrent GP on large graphs requires the use of the hybrid classicla-quantum $qbsolv$. Table5 shows the results of splitting 250, 500, and 1000 nodes graphs into 2, 4, 8, and 16 parts. The quality is comparable to METIS, and the numbers of $cut \space edges$ is reduced by tens to hundreds.
>> 
>>
>>
>> For partitioning a graph into $k=2^i$ for $i>1$, we used molecule electronics structure graphs as density matrices generated from QMD simulation. Figure below shows the phenyl dendrimer molecule. The electronics structure of this molecule consists of 730 orbitals, resulting in a graph of 730 vertices.
>>
>> Electronic structure graph $G_\tau(\rho)$ constructed from the QMD density matrix $\rho$ using a threshold $\tau$ that leads to 730 vertices. In the end we used D-Wave machine to achieve 4-Concurrent GP.
>> 
>>
>>
>> In table6 the results for $k$-Concurrent GP of the molecule electronics structure graphs are shown. This is a difficult case for GP due to the fractal structure. The $qbsolv$ results are comparable or better than METIS for the Phenyl dendrimer and the Peptide laft structure graphs. When 4-partitioning the Phenyl Dendrimer we see a very large number of $cut \space edges$ difference.
>> 

#### Conclusion

> We have shown that GP framed as a QUBO problem is a natural fit for the D-Wave quantum annealer. Results using quantum($sapi$) and hybrid classical-quantum($qbsolv$) approaches are comparable or better than existing traditional GP tools(METIS, KaHIP)
We also showned that solving CD(community detecting) using a threshold modularity matrix did not change the community results and could reduce qubit/coupler.
In 2-partitioning GP problem, our results equaled or outperformed existing GP tools and in some cases the best known results.
>
>In k-Concurrent GP problem, using super-node concept to partition random graphs and molecule electron structure graphs into 2, 4, 8, and 16 parts. Our results improved the quality and reduced the number of $cut edges$. $K$-concurrent GP could run on QPU directly for small graphs and using hybrid classical-quantum $qbsolv$ for large graphs.
>
>

:::info
==未解之謎==
<font color='#f00'>ask: 為甚麼KaHIP只出現在某幾個table?是有甚麼限制嗎?</font>
<font color='#f00'>ask: hybrid classical-quantum approach</font>
<font color='#f00'>ask: QPU?</font>
<font color='#f00'>ask: $sapi$是python下的? $qbsolv$是D-wave的api?</font>
<font color='#f00'>ask: 使用$sapi$ 的好處是甚麼?為甚麼不要都用$qbsolv就好$反正都能解的話</font>
**是考慮全quantum跟hybrid approach的比較?**
<font color='#f00'>ask: table2所指的simulator是?"complement in the QUBO formulation"演算法意義?</font>
<font color='#f00'>ask: table5好像沒有提到是用甚麼graph??</font>
<font color='#f00'>ask: Phenyl dendrimer molecular structure/ QMD density matrix+threshold ==(reference[13])==</font>
<font color='#f00'>ask: the largest complete graph that is fully embeddable on the D-Wave 2X ~ 45 vertices.
是因為45 vertices是考慮complete graph的情況所以如果考慮fraction $p$ 可以做超過45 vertices?
這樣的話造成運算量過大的其實是edges數量不是vertices有幾個?</font>

ask: modularity threshold limit, modularity下降是不好的嗎?
ans: modularity threshold評估graph connect的程度，對於後面GP的關聯是甚麼?edge直接不看?

==已解決==
- ask: QMD?
ans: quantum molecule dynamic。
- ask: graph size N$\times$k 所表示的是甚麼?
ans: 將 N nodes partition成 k parts，所使用的qubits數量為N$\times$k。

:::

- 名詞解釋
- [karate club graph](https://en.wikipedia.org/wiki/Zachary%27s_karate_club)
- visualize:
- [Erdos-Renyi random graph](https://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93R%C3%A9nyi_model)
- [Dendrimer](https://en.wikipedia.org/wiki/Dendrimer)
- dendrimer molecular structure: 樹枝狀分子結構
- KaHIP
- KaHIP is a family of graph partitioning programs. In this work, we used KaFFPaE with its preconfiguration variant set to strong and KaFFPaE with the setting designed for $perfectly balanced$ partitioning.

- softwares tools
1. $sapi$: up to 70 vertices
2. $qbsolv$: up to 9000 vertices
3. NetworkX: gernerating and processing graphs

- quantity
- *cut edges number*(smaller is better)

- compared
- multi-level GP frameworks
- METIS
- KHIP

- Table1
we can set the thresholds for modularity matrix .For instance, with a threshold value of 0.12 meaning that we set Bij = 0 if Bij < 0.12, so we can ignore the edge between the pair of nodes;if all edges of a node are ignored,we can ignore the node,

above the table, we use community dection to the karate club graph(34 vertices). We see that the larger a threshold is, the smaller a modularity is, while number of edges we need decreases. In this case we can see that with a threshold value of 0.12
modularity is reduced by ∼ 30 % (we think it can be accepted )but the required number of edges is significantly less (∼ 65 % less)In this way, we can produce a sparse version that is still representative of the original allows for larger problems to be solved efficiently by community clustering using quantum annealing.

- Table2
result for 4-clustering using the recursive method, using random graphs with powerlaw degree distribution(100 vertices).In this case, quality of the clustering is significantly reduced.


==**$\rightarrow$ both cases reduce the number of edges**==

<font color='#f00'>ask: change of vertice number</font>

- Table3
study two sets of graph chosen in order to determine the limits of the different partitioning methods used.First set had at most 70 vertices, and the second set consisted of graph with at most 9000 vertices

For the first set, we partitioned the graph using $sapi$. Thus 100 percent of the graph was carried out by the quantum annealer. The largest complete graph that is fully embeddable on the D-Wave 2X has approximately 45 vertices.
Our results demonstrate that for large values of $p$($p$ means the fraction of total possible edges), we can successfully partition such graphs even if they have more than 45 vertices.

- Table4
For the second set, we partitioned the graph using $qbsolv$. Since $qbsolv$ is a hybrid classical-quantum approach, we could patition graphs with over 100 vertices.

- Table5
For partitioning a graph into $k=2^i$ parts for $i>1$, we performed our experiments on electronic structure graphs as density matrices generated from QMD simulation.
Table5 shows the result of partitioning a 730 vertices graph(Phenyl dendrimer) into 2, 4 and 8 parts.
Partitioning each molecule electronic structure graph into 2 parts using $qbsolv$ results in equal or reduced number of $edge cuts$ compared to METIS and KaHIP. However the number of $edge cuts$ increased when partitioning into 4 and 8 parts through recursive bisection for the Phenyl dendrimer. Comparing the result of Peptide Laft protein, we have better results for all cases.(2, 4, or 8 parts)

- Table6
k-Concurrent GP requires that a graph of size $N$ (numbers of nodes) $\times$ $k$ (number of partitioning) be embedded in the $Chimera graph$.

Table6 shows the result of partitioning $N \times k$ graph, using $sapi$, METIS, $qbsolv$ three methods. Results using $sapi$ are comparable to METIS and $qbsolv$. A 15 node graph patitioned into 4 parts used almost all available qubits. Similarly for 20 nodes split into 3 parts.

- Table7
Results for random dense graphs of 250, 500, and 1000 nodes split into 2, 4, 8, and 16 parts are shown in Table7.

- Table8
In Table8 the results for k-Concurrent GP of the molecule electronics structure graphs are shown.
The $qbsolv$ results are comparable or better than METIS for the Phenyl dendrimer and Peptide laft. We see a very large number of $cut edges$ for METIS when 4-partitioning the Phenyl Dendrimer.

- reference
- [reference 28](https://link.springer.com/chapter/10.1007/978-3-642-38527-8_16) : p.164-p.175
- [modularity 參考講義](https://www.cs.cmu.edu/~ckingsf/bioinfo-lectures/modularity.pdf)

#### Conclusion
We have shown that GP framed as a QUBO problem is a natural fit for the D-Wave quantum annealer. Results using quantum($sapi$) and hybrid classical-quantum($qbsolv$) approaches are comparable or better than existing traditional GP tools(METIS, KaHIP)
We also showned that solving CD(community detecting) using a threshold modularity matrix did not change the community results and could reduce qubit/coupler.
In 2-partitioning GP problem, our results equaled or outperformed existing GP tools and in some cases the best known results.
In k-Concurrent GP problem, using super-node concept to partition random graphs and molecule electron structure graphs into 2, 4, 8, and 16 parts. Our results improved the quality and reduced the number of $cut edges$. $K$-concurrent GP could run on QPU directly for small graphs and using hybrid classical-quantum $qbsolv$ for large graphs.

---

## 參考資料

---

# Cad Contest

1. 讀檔案
3. 資料結構
- 建成 graph
- 辨識 drill 跟 pad 的關係
- 畫圖
4. 偵測違法
- 切割成違法 / 沒違法
6. 嘗試修復
- 找最佳修復方法
- 移動
- 變形
- 切割
7. 如果還有違法，回到 4

### how
- 要 solve 的問題要講清楚
- Graph 的結構
### 演算法(暫定)
找出每個物件(在graph上代表node:pth、smd)的bounding box(以下簡稱bbox)，判斷其是否重疊，若有重疊則在graph上代表node之間有edge，edge 的權重代表重疊面積，
- step1:
將 node 分成 connected (有和其他物件重疊) 和 disconnected (沒有和其他物件重疊) 兩堆
- step2:
在有重疊那堆點中找出edge數最少的pth，
- step3:
將該pth的bbox依北、東北、東、東南、南、西南、西、西北、北移動，若在某方向上移動能使其與相鄰物件重疊面積減少，則繼續移動直到與drill的bbox的間距為min width，計算此時的重疊面積，若此時重疊面積比原本小，則記住此時的移動方向和偏移距離，八個方向均跑完，可知最佳移動方向和偏移距離
- step4:
將pth依照該最佳移動方向和偏移距離移動
- step5:
回到step1

### Problem Description

#### Goal

設計一自動調整PCB銅箔分布之作業，且為避免銅箔過近可能造成上錫時電路短路、銅箔面積不足/分布不均導致零件鬆脫，因此需同時考慮多種銅箔限制。
![](https://i.ibb.co/jRRjhQx/image.png")

#### Constraint
- mingap : 銅箔與鑽孔最小之距離(如下圖左)
- minwidth : 銅箔之間最小之寬度(如下圖左)
- geometry : 外圍輪廓變化量，計算如下:
若為一圓 : Arc(弧)為1，Line(直線)為0，若調整後的銅箔圖形為長方形，Arc為0，Line為4，則外圍輪廓變化量(geometry)為Arc變化量|1-0|與Line變化量|0-4|之加總，即為1+4=5。

#### Concept
常見的銅箔調整方法有以下兩種 :
Offset case : 直接偏移銅箔，以滿足規定
Cutting case : 將不符合規定的銅箔部分進行裁切
#### Evaluation
- 符合問題基本要求 : (若無法符合以 0 分計算)
- 調整後銅箔數量必須與調整前相等
- 描述圖形之資料內容須滿足簡單多邊形(Simple polygon)之定義
- 所有銅箔間距均大於等於銅箔最小間距限制值
- 所有銅箔寬度均大於等於銅箔最小寬度限制值
- 需在 300 秒內完成執行完畢
- 計分方式 :
- 銅箔外觀相似度評分 (30%)

- 配置位置評分 (30%)
- 銅箔面積評分 (40%)
- 若評分相同者者，則以程式執行效能分數高者為佳
### Data Structure
#### Input file
在ProblemA.txt、ProblemB.txt、ProblemC.txt三組公開測資中，分別標示著銅箔和鑽孔的位置和幾何形狀以及問題的限制，其資料標籤和幾何資訊如下圖所示
其中mingap、minwidth、geometry為問題的constraint，smd、pth、drill為擺放在電路板上的銅箔和鑽孔標籤，其中需要注意的是smd和drill是不可移動的，本題僅對pth做位置上的調整。arc和line描述了smd、pth、drill的幾何形狀，下圖為幾何資訊和其形狀的對應
#### Graph structure
為了方便pth的移動以及判斷pth之間的重疊情形，我們將Input file的資訊轉為Graph的方式儲存
- Node
node 的類型有三種，分別為smd、pth、drill，每個node的屬性包含了bounding box以及geometries，後者描述了此node的幾何結構(由幾條arc、line或circle構成)，前者是根據此幾何結構框出的bounding box，目的是為了檢察node和node之間有沒有碰撞(區域重疊)
- Edge
兩個node之間存在著edge代表兩個node有區域重疊，edge的權重代表兩個node之間重疊的面積

#### 演算法
- Related Work :
- Stephen G. Kobourov, "Chapter 12, Force-Directed Drawing Algorithms", Handbook of Graph Drawing and Visualization, 2013
- Chun-Cheng Lin and Hsu-Chun Yen, "A new force-directed graph drawing method based on edge-edge repulsion," Ninth International Conference on Information Visualisation (IV'05), 2005, pp. 329-334, doi: 10.1109/IV.2005.10.method based on edge-edge repulsion

將圖形中 node 與 node 之間具 egde 連接者，施予一排斥力以分離兩 node，經由迭代不斷使兩兩 node 越來越遠。過程中 edge 越來越長，而排斥力越來越小，最終整個圖形會呈現穩定狀態(或是收斂成幾乎穩定)，觸發終止條件(排斥力為 0)。

- Proposed Design Flow:
- step 1:
進行檔案預處理，將題目之 arc 與 line 等資訊建立成自定義graph。
- step 2:
檢查 graph 中是否有碰撞到的點，若有，進入 step 3；否則，表示整張圖已達終止條件進入 step 6。
- step 3:
將有碰撞/無碰撞之 node 進行 partition。
- step 4:
計算 connected component 內的 directed forces 及其合力，並計算出碰撞各點之位移向量。
- step 5:
根據step 4 之位移向量，更新每個 node 的新位置，更新完後回 step 2 重新檢查是否有碰撞。
- step 6:
將調整後之結果轉回題目規定之格式(arc、line等)，輸出結果。


- Collision Check
- Detection Zone :
首先設計 detection zone 以計算是否有碰撞。detection zone 規劃如下:
- drill detection zone :
黑色部分為鑽孔面積(drill)。
灰色部分為銅箔與drill之最小寬度(min width)。
- pth detection zone :
紅色部分為銅箔面積(pth)，可調整之區域。
淡紅色部分為$\frac{min gap}{2}$，因為mingap為兩銅箔之間最短距離，因此分給兩邊銅箔之距離需除2。

註: 未有重疊的 node 並不會影響其他點，因此只有 connected component 才會受計算的力作用。
- Bounding Box :
CGAL函式庫計算不同形狀圖形之交點連線(intersection)，將圖形和圖形之間的每個區段(在此為geometry)做intersection判斷，因此intersection計算複雜度為 $O((N$ x $N$~G~$)^2)$ ，其中$N$為點數量，$N_G$ 為各點geometry值的平均。
因intersection計算量較為龐雜，我們設計了Bounding Box手法。
Bounding Box(BBox)為不規則多邊形之外接長方形。將各點轉換成的BBox，再使用intersection初步判斷是否兩BBox之間detection zone有無重疊(overlap)，若有則代表collision，可以得知多邊形之間初步的collision數$\gamma$，再針對這些有collision的多邊形做較仔細的intersection判斷

此法為多邊形轉為四邊形後才進行overlap判斷，因此可大幅降低計算複雜度。
- Graph Bipartitioning :
根據overlap與否，將點分為overlapping為一堆，並且將有重疊之兩點用edge做連線(如下圖點2,3,4)。其他沒有overlapping的點為另一堆，沒有edge連接(如下圖點1,5)。
此步驟目的為忽略上階段BBox就沒有overlap的點，由BBox有overlap之點再進行更精準的intersection計算。

- Force Calculation
- Repulsive forces(排斥力):
當兩個nodes發生交疊時會違反constraint，因此施予一penalty distance，意即讓此兩nodes產生排斥力。又此排斥力分為以下兩種 :
- Repulsive between pth :
兩個pth不可靠太近，意即mingap constraint之penalty。
下圖所示，有兩小於mingap之點$u$與點$v$，intersection中點$cI$位置計算，為所有$u，v$交集點座標加總。而對點$u$之repulsive force為點$cI$至點$u$之向量，反之點$v$即為反作用力。
- Repulsive between pth and drill
pth在移動時可能會使drill與pth之間距離太小，意即違反minwidth constraint。
下圖所示，repulsive force為pth圓心至drill圓心之向量。
- Attractive force (吸引力):
因題目其一評分標準為移動後的位置與原本的位置離越近越好($Placement$)，因此設計一attractive force使得移動後的node可以往原本node的位置形成一吸引力，方向指向原本位置。

- Position Updating :
計算上述各力之加權總和，即可更新此次graph所有點的位置。
$A_u$為node $u$面積，類比為下式之質量，意即面積越大，移動位置需越小，$Placement$評分項分數才會越高。
$\Delta t$ 為每次疊代之相隔時間，以套用於衝量概念
- Terminating Condition
疊代終止條件即為整張圖都沒有碰撞產生的穩定狀態，意即 $\vec{F}_{Rep} = 0$
### Result
#### Data Visualization
- Generat visualized graph from .txt file.Graph
- Geometry labels (line, arc)
- Data labels (smd, pth, drill)

  

上方三張圖為我們利用CAD contest所提供之A/B/C測試資料所得到的視覺化後資料，後續我們會去判斷當前的資料關係進行最佳化。

#### Original Graph Expansion
 

左圖為視覺化後所得之圖片，在graph expansion這部份我們會依照演算法所述將他向外擴展，並利用擴展後的資料去實行我們的演算法。右圖為向外擴展後所得之相對位置。
右圖向外拓展的部分為Detection zone，我們會藉由偵測Detection zone的交疊與否判斷是否需要對此物件進行移動。而圖中藍線表示的是pth和smd的輪廓，綠線表示的是drill的輪廓。
#### Check Collision
- Bounding box
- Curve intersection


從Graph Extension後的圖我們可以去尋找各contour的交界，偵測出那些是我們需要移動的部分，而哪些是我們不需要移動的部分。若是將交錯部分放大(下圖)，可以觀察到紅色點的部分為我們偵測到的Detection zone有交疊區域。


#### Configuration of copper foil location
在偵測到交疊區域後我們會移動有交疊的部分，透過演算法所提及的iteration算法，重複移動並偵測不符合我們要求的部分，直到所有Detection zone均沒有重疊為止。完成後的視覺化資料如下圖所示，到此步驟後可以確認所有條件均符合題目要求。


#### 銅箔位置評分
進行銅箔擺放位置最佳化時就可以得到center移動路徑

#### 銅箔面積評分
- step1:
產生出兩張圖，分別將其針對不同element著色，區分出想要判斷重疊面積的部分。
- step2:
裡用opencv估算出重疊面積(兩張圖片做and運算)，精讀由取的pixel大小決定。
- And operation
- Find contour
- Calculate contour area

#### 銅箔相似度

### 參考資料
- [List of graph layout algorithm](https://stats.stackexchange.com/questions/51519/list-of-graph-layout-algorithms/383565)
- [A New Force-Directed Graph Drawing Method Based on Edge-Edge Repulsion](https://acms.nctu.edu.tw/papers/JVLC2011-EERepulsion.pdf)
- Steps:
1. compute the repulsive forces
2. compute the spring force of each edge
- [Force-Directed Drawing Algorithms](https://cs.brown.edu/people/rtamassi/gdhandbook/chapters/force-directed.pdf)
- [Boost forced directed model](https://www.boost.org/doc/libs/1_33_1/libs/graph/doc/fruchterman_reingold.html)