import numpy as np 
from collections import defaultdict 
import matplotlib.pyplot as plt

def simulate(alpha, num_blocks=20000):
    honest_chain = 0  # cumulative sum of weight of honest chain
    attacker_chain = 0  # cumulative sum of weight of attackers chain
    attackers_blocks = 0
    attacker_reward = 0
    honest_reward = 0

    for _ in range(num_blocks):
               
        r = np.random.rand()
        weight = np.log(r)
               
        if r < alpha:  # Attacker mines
            attacker_chain += weight
            attackers_blocks = attackers_blocks + 1
           
            if attacker_chain < honest_chain:
           
                # attacker reveals branch when it's got the least probability of being mined block by block
                attacker_reward += attackers_blocks
                honest_reward -= attackers_blocks
               
                # start again
                attackers_blocks = attacker_chain = honest_chain = 0
               
           
        else:  # Honest miner mines
            honest_chain += weight
            honest_reward += 1
               
    if honest_reward > 0:
        return attacker_reward / honest_reward
    else:
        return float("inf")

# selfish mining as explained in "Majority is not Enough: Bitcoin Mining is Vulnerable" https://arxiv.org/pdf/1311.0243
# note that gamma is optimistic so this is a lower bound
def simulateBitcoin(alpha, num_blocks=20000, gamma=0.5):
    honest_chain = 0  # Length of honest chain
    attacker_chain = 0  # Length of attacker's private fork
    attacker_reward = 0
    honest_reward = 0
    deltaPrev = 0

    for _ in range(num_blocks):
        deltaPrev = attacker_chain - honest_chain
       
        if np.random.rand() < alpha:  # Attacker mines
            attacker_chain += 1
           
            if deltaPrev == 0 and attacker_chain == 2:
                attacker_reward += attacker_chain
                attacker_chain=honest_chain=0
           
        else:  # Honest miner mines
            honest_chain += 1
           
            if deltaPrev<=0:
                honest_reward += honest_chain
                attacker_chain = honest_chain = 0
            elif deltaPrev==1:
                # miners chose at random which fork to mine on based on gamma parameter
                g = np.random.rand()
                if g <= gamma:
                    attacker_reward += 1
                else:
                    honest_reward += 1
               
                attacker_chain=honest_chain=0
            elif deltaPrev>=2:
                attacker_reward += deltaPrev+1
                attacker_chain = honest_chain = 0       

    return attacker_reward / honest_reward
   
alphas = np.linspace(0.15, 0.5, 20)  # 15% to 50% attacker hashing power
results_poem = [simulate(alpha) for alpha in alphas] 
results_bitcoin = [simulateBitcoin(alpha) for alpha in alphas] 

# Plot results 
plt.plot(alphas, results_poem, label="PoEM") 
plt.plot(alphas, results_bitcoin, label="Bitcoin") 
plt.axhline(y=1, color="red", linestyle="--", label="Attack break-even cost") 
plt.xlabel("Attacker Hashrate (α)") 
plt.ylabel("Attacker Reward / Honest Reward") 
plt.ylim([0, 1.2])
plt.legend() 
plt.show()