example_sat_problem = \
  [[{"sign": True, "variable": "A"}, 
    {"sign": False, "variable": "B"}, 
    {"sign": False, "variable": "D"}],
   [{"sign": False, "variable": "B"}, 
    {"sign": True, "variable": "C"}, 
    {"sign": False, "variable": "D"},
    {"sign": True, "variable": "E"}],
   [{"sign": True, "variable": "C"}, 
    {"sign": False, "variable": "C"}, 
    {"sign": True, "variable": "A"}]]

example_assignment = {"B": False, "E": True, "C": False}

def satisfied(SAT_problem, variable_assignment):
    from collections import defaultdict
    variable_assignment = defaultdict(lambda: False, **variable_assignment)
    return any(variable_assignment[literal["variable"]] == literal["sign"]
               for clause in SAT_problem
               for literal in clause)

satisfied(example_sat_problem, example_assignment)

def next_random(n):
    return (abs(n) * 3**160 + 7**80) % (2**256 - 4294968273)

def generate_literals(n, generators):
    return [{"variable": n % g, "sign": ((n >> idx) % 2) == 1, "generator": g} 
            for idx,g in enumerate(generators)]

generate_literals(3681122926843, [9,13,14,21,33,35])

next_random(3681122926843+9+21+33)

def sat_problem_generated(header, generators, sat_problem):
    n = header
    for clause in sat_problem:
        possible_literals = generate_literals(n, generators)
        if not all(literal in possible_literals for literal in clause):
            return False
        n = next_random(n + sum(literal["generator"] for literal in clause))
    return True

header = 5751887518028651109

example_generators = [5,7,8,9,11,13,17,19]

example_generated_sat_problem = \
 [({'generator': 5, 'sign': True, 'variable': 4},
   {'generator': 7, 'sign': False, 'variable': 5},
   {'generator': 9, 'sign': False, 'variable': 6}),
  ({'generator': 5, 'sign': False, 'variable': 3L},
   {'generator': 7, 'sign': True, 'variable': 4L},
   {'generator': 8, 'sign': False, 'variable': 2L}),
  ({'generator': 5, 'sign': False, 'variable': 0L},
   {'generator': 7, 'sign': True, 'variable': 6L},
   {'generator': 8, 'sign': True, 'variable': 6L})]

sat_problem_generated(5751887518028651109, example_generators, example_generated_sat_problem)

def number_of_clauses(n, clauses):
    return n == len(clauses)

def clause_variables(clause):
    "Returns a frozenset of variables in a clause"
    return frozenset(literal["variable"] for literal in clause)

def clause_width(n, clause):
    return n == len(clause_variables(clause))

def all_clauses_width(n, clauses):
    return all(clause_width(n, clause) for clause in clauses)

def no_repeated(clauses):
    visited = set([])
    for clause in clauses:
        variables = clause_variables(clause)
        if variables in visited:
            return False
        visited.add(variables)
    return True

example_unique = \
 [({'sign': True, 'variable': 4L},
   {'sign': False, 'variable': 6L},
   {'sign': True, 'variable': 5L}),
  ({'sign': True, 'variable': 3L},
   {'sign': False, 'variable': 3L},
   {'sign': True, 'variable': 5L}),
  ({'sign': True, 'variable': 3L},
   {'sign': True, 'variable': 2L},
   {'sign': False, 'variable': 3L}),
  ({'sign': False, 'variable': 4L},
   {'sign': False, 'variable': 6L},
   {'sign': True, 'variable': 4L})]

no_repeated(example_unique)

example_not_unique = \
 [({'sign': True, 'variable': 0L},
   {'sign': False, 'variable': 0L},
   {'sign': False, 'variable': 7L}),
  ({'sign': False, 'variable': 0L},
   {'sign': False, 'variable': 1L},
   {'sign': False, 'variable': 0L}),
  ({'sign': True, 'variable': 3L},
   {'sign': True, 'variable': 2L},
   {'sign': True, 'variable': 7L}),
  ({'sign': True, 'variable': 4L},
   {'sign': False, 'variable': 3L},
   {'sign': True, 'variable': 5L}),
  ({'sign': True, 'variable': 4L},
   {'sign': True, 'variable': 1L},
   {'sign': False, 'variable': 3L}),
  ({'sign': True, 'variable': 0L},
   {'sign': True, 'variable': 0L},
   {'sign': True, 'variable': 7L})]

no_repeated(example_not_unique)

def connected(c1, c2):
    ls1 = clause_variables(c1)
    ls2 = clause_variables(c2)
    return any(l in ls2 for l in ls1)

def clauses_connected_by_traversal(clauses, traversal):
    if len(traversal) != len(clauses) or \
       not (0 <= traversal[0]['child clause index'] < len(clauses)):
        return False
    visited = set([traversal[0]['child clause index']])
    for i in range(1,len(traversal)):
        child_clause_index = traversal[i]['child clause index']
        parent_clause_index = traversal[i]['parent clause index']
        if child_clause_index in visited or \
           not (parent_clause_index in visited) or \
           not (0 <= child_clause_index < len(clauses)) or \
           not connected(clauses[child_clause_index],
                         clauses[parent_clause_index]):
               return False
        visited.add(child_clause_index)
    return True

connected_example_sat_problem = \
 [({'sign': True, 'variable': 4L},
   {'sign': False, 'variable': 4L},
   {'sign': True, 'variable': 5L}),
  ({'sign': True, 'variable': 3L},
   {'sign': False, 'variable': 4L},
   {'sign': False, 'variable': 1L}),
  ({'sign': True, 'variable': 2L},
   {'sign': True, 'variable': 6L},
   {'sign': False, 'variable': 3L}),
  ({'sign': True, 'variable': 3L},
   {'sign': False, 'variable': 0L},
   {'sign': True, 'variable': 5L}),
  ({'sign': True, 'variable': 4L},
   {'sign': False, 'variable': 0L},
   {'sign': False, 'variable': 1L}),
  ({'sign': False, 'variable': 2L},
   {'sign': False, 'variable': 1L},
   {'sign': True, 'variable': 4L})]
    
connected_example_spanning_tree_traversal = \
 [{'child clause index': 0},
  {'child clause index': 1, 'parent clause index': 0},
  {'child clause index': 2, 'parent clause index': 1},
  {'child clause index': 3, 'parent clause index': 2},
  {'child clause index': 4, 'parent clause index': 1},
  {'child clause index': 5, 'parent clause index': 2}]

clauses_connected_by_traversal(connected_example_sat_problem, connected_example_spanning_tree_traversal)

%load_ext hierarchymagic

%%dot -f svg
graph G {
    0--1; 1--2; 2--3; 1--4; 2--5;
}

def number_of_generators(n, generators):
    return n == len(generators)

def min_generator(m, generators):
    return m <= min(*generators)

def sha256_less_than_target(block_header, variable_assignment, target):
    import hashlib
    import json
    s = hashlib.sha256()
    s.update(str(block_header))
    s.update(json.dumps(variable_assignment, sort_keys=True))
    return int('0x' + s.hexdigest(), 16) < target

example_block_header = 28850183455121774801176625632015427996404687904360779760892949139292104258030L
example_variable_assignment = {
  1L: True,
  3L: True,
  6L: False,
  8L: False,
  9L: True,
  15L: True,
  17L: True,
  18L: False,
  19L: True,
  23L: True,
  28L: False,
  30L: False,
  35L: True,
  36L: False,
  37L: True }

sha256_less_than_target(example_block_header, example_variable_assignment, 2**256 / 3)

sha256_less_than_target(example_block_header, example_variable_assignment, 2**256 / 4)

def proof_of_work_validate(solution, parameters):
    return (satisfied(solution["clauses"], solution["variable assignment"]) and
            sat_problem_generated(parameters["block header"], solution["generators"], solution["clauses"]) and
            number_of_clauses(parameters["number of clauses"], solution["clauses"]) and
            all_clauses_width(parameters["clause width"], solution["clauses"]) and
            no_repeated(solution["clauses"]) and
            clauses_connected_by_traversal(solution["clauses"], solution["spanning tree traversal"]) and
            number_of_generators(parameters["number of generators"], solution["generators"]) and
            min_generator(parameters["minimal generator"], solution["generators"]) and
            sha256_less_than_target(parameters["block header"], solution["variable assignment"], parameters["target"]))

def mine(params):
    i = 0
    while True:
        start = params["minimal generator"] + i
        generators = range(start, start + parameters["number of generators"])
        for state in initial_states(generators, params):
            for solution in search(state,params):
                if sha256_less_than_target(params["block header"], solution["variable assignment"], params["target"]):
                    yield solution
                else:
                    break
        i += 1
    
def search(state, params):
    if (len(state["clauses"]) == params["number of clauses"]):
        yield {"clauses": state["clauses"], 
               "spanning tree traversal": state["spanning tree traversal"],
               "variable assignment": state["variable assignment"],
               "generators": state["generators"]}
    else:
        for new_state in possible_new_states(state, params):
            for solution in search(new_state, params):
                yield solution

def possible_new_states(state, params):
    from itertools import combinations
    for clause in combinations(generate_literals(state["current random"], 
                                                 state["generators"]), 
                               params["clause width"]):
        parent_clause_index = None
        for literal in clause:
            if literal["variable"] in state["variable spanning tree indexes"]:
                parent_clause_index = state["variable spanning tree indexes"][literal["variable"]]
                break
        if (parent_clause_index == None or
            not clause_width(params["clause width"], clause) or
            clause_variables(clause) in state["visited clause variables"]):
            continue
        for pivot in clause:
            if not (pivot["variable"] in state["variable assignment"] and
                    state["variable assignment"][pivot["variable"]] != pivot["sign"]):
                yield update_state(state, clause, parent_clause_index, pivot)

def initial_states(generators, params):
    from itertools import combinations
    for clause in combinations(generate_literals(params["block header"], generators),
                               params["clause width"]):
        if not clause_width(params["clause width"], clause):
            continue
        for pivot in clause:
            yield {"current random": next_random(params["block header"] + 
                                                 sum(literal["generator"] 
                                                     for literal in clause)),
                   "clauses": [clause],
                   "variable assignment": {pivot["variable"]: pivot["sign"]},
                   "spanning tree traversal": [{ "child clause index": 0 }],
                   "visited clause variables": set([clause_variables(clause)]),
                   "variable spanning tree indexes": { literal["variable"]: 0 
                                                       for literal in clause },
                   "generators": generators }

def update_state(state, clause, parent_clause_index, pivot):
    from copy import deepcopy
    index = len(state["clauses"])
    new_state = deepcopy(state)
    new_state["clauses"].append(clause)
    new_state["visited clause variables"].add(clause_variables(clause))
    new_state["variable assignment"][pivot["variable"]] = pivot["sign"]
    for l in clause:
        new_state["variable spanning tree indexes"][l["variable"]] = index
    new_state["spanning tree traversal"].append({"parent clause index": parent_clause_index,
                                                 "child clause index": index})
    new_state["current random"] = \
        next_random(state["current random"] + 
                    sum(literal["generator"] for literal in clause))
    return new_state

block_header = 59807446017497444054206621583131621044750062436043876690313232144928752617878L
parameters = {
    "block header": 59807446017497444054206621583131621044750062436043876690313232144928752617878L,
    "number of clauses": 50,
    "clause width": 3,
    "target": 2**256 / 2,
    "minimal generator": 40,
    "number of generators": 8
}
solution = mine(parameters).next()
proof_of_work_validate(solution, parameters)

def bitcoin_mine(block_header, target):
    import hashlib
    nonce = 0
    while int('0x' + hashlib.sha256(hashlib.sha256(str(block_header + nonce)).digest()).hexdigest(), 16) >= target:
        nonce += 1
    return nonce

from timeit import timeit

max_val = 2**256
def time_bitcoin_difficulty(d):
    return timeit('bitcoin_mine(randint(0,max_val), %d)' % (max_val / d), 
                  'from random import randint ; from __main__ import bitcoin_mine, max_val', 
                  number=1)

bitcoin_difficulties = range(100,1000,10)
bitcoin_sample_times = [[time_bitcoin_difficulty(d) for _ in range(500)] for d in bitcoin_difficulties]

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import numpy as np
from pylab import polyfit
import matplotlib.pyplot as plt
average_times = map(np.mean, bitcoin_sample_times)
C,k = polyfit(bitcoin_difficulties, average_times, 1)
plt.plot(bitcoin_difficulties, average_times)
plt.plot(bitcoin_difficulties, [C*d + k for d in bitcoin_difficulties], '-r')
plt.title("Average Mining Time (Linear)")
plt.show()
print "Estimated Ĉ:", C, "secs"
print "Overhead:", k, "secs"

variances = map(np.var, bitcoin_sample_times)
C2,C1,k = polyfit(bitcoin_difficulties, variances, 2)
plt.plot(bitcoin_difficulties, variances)
plt.plot(bitcoin_difficulties, [C2*d*d + C1*d + k for d in bitcoin_difficulties], '-r')
plt.title("Mining Time Variance (Quadratic)")
plt.show()
print "C₁:", C1, "C₂:", C2
print "Estimated Ĉ²:", C**2

import scipy.stats
import pandas

data = []
for d, samples in zip(bitcoin_difficulties[::5], bitcoin_sample_times[::5]):
    cdf = lambda t: 1 - (1 - (1. / d))**(t / C)
    ks_statistic, p_value = scipy.stats.kstest(samples, cdf)
    data.append({"Difficulty D": d,
                 "KS Statistic": ks_statistic,
                 "P-Value": p_value})


pandas.DataFrame.from_records(data)

def random_parameters_with_difficulty(d, number_of_clauses, min_generator):
    from random import randint
    return {
        "block header": randint(0,2**256),
        "number of clauses": number_of_clauses,
        "clause width": 3,
        "target": 2**256 / d,
        "minimal generator": min_generator,
        "number of generators": 8
    }

def time_sat_difficulty(d,n,m):
    return timeit('next(mine(random_parameters_with_difficulty(%d,%d,%d)))' % (d,n,m), 
                  'from __main__ import random_parameters_with_difficulty, mine',number=1)

sat_difficulties = range(1,50)
sat_sample_times = [[time_sat_difficulty(d,10,5) for _ in range(100)] for d in sat_difficulties]

plt.plot(sat_difficulties, map(np.mean, sat_sample_times))
C,k = polyfit(sat_difficulties, map(np.mean, sat_sample_times), 1)
plt.plot(sat_difficulties, [C*d + k for d in sat_difficulties], '-r')
plt.title("Average Mining Time (Linear)")
plt.show()
print "Estimated Ĉ:", C, "secs"
print "Overhead:", k, "secs"

variances = map(np.var, sat_sample_times)
C2,C1,k = polyfit(sat_difficulties, variances, 2)
plt.plot(sat_difficulties, variances)
plt.plot(sat_difficulties, [C2*d*d + C1*d + k for d in sat_difficulties], '-r')
plt.title("Mining Time Variance (Quadratic)")
plt.show()
print "C₁:", C1, "C₂:", C2
print "Estimated Ĉ²:", C**2

data = []
for d, samples in zip(sat_difficulties[4::5], sat_sample_times[4::5]):
    cdf = lambda t: 1 - (1 - (1. / d))**(t / C)
    ks_statistic, p_value = scipy.stats.kstest(samples, cdf)
    data.append({"Difficulty D": d,
                 "KS Statistic": ks_statistic,
                 "P-Value": p_value})


pandas.DataFrame.from_records(data)

clause_lengths = range(10,100,10)
sat_clause_sample_times = [[time_sat_difficulty(1,l,5) for _ in range(100)] for l in clause_lengths]

plt.plot(clause_lengths, map(np.mean, sat_clause_sample_times))
C,k = polyfit(clause_lengths, map(lambda x: np.log(np.mean(x)), sat_clause_sample_times), 1)
plt.plot(clause_lengths, [np.exp(C*d + k) for d in clause_lengths], '-r')
plt.yscale('log')
plt.title("Average Mining Time (Exponential)")
plt.show()

plt.plot(clause_lengths, map(np.var, sat_clause_sample_times))
C,k = polyfit(clause_lengths, map(lambda x: np.log(np.var(x)), sat_clause_sample_times), 1)
plt.plot(clause_lengths, [np.exp(C*d + k) for d in clause_lengths], '-r')
plt.yscale('log')
plt.title("Mining Time Variance (Exponential)")
plt.show()

minimal_generators = range(100,1000,100)
sat_min_gen_sample_times = [[time_sat_difficulty(1,5,g) for _ in range(100)] for g in minimal_generators]

plt.plot(minimal_generators, map(np.mean, sat_min_gen_sample_times))
C,k = polyfit(minimal_generators, map(lambda x: np.log(np.mean(x)), sat_min_gen_sample_times), 1)
plt.plot(minimal_generators, [np.exp(C*d + k) for d in minimal_generators], '-r')
plt.yscale('log')
plt.title("Average Mining Time (Exponential)")
plt.show()

plt.plot(minimal_generators, map(np.var, sat_min_gen_sample_times))
C,k = polyfit(minimal_generators, map(lambda x: np.log(np.var(x)), sat_min_gen_sample_times), 1)
plt.plot(minimal_generators, [np.exp(C*d + k) for d in minimal_generators], '-r')
plt.yscale('log')
plt.title("Mining Time Variance (Exponential)")
plt.show()