import networkx as nx
import random
import math
import numpy as np
#import matplotlib.pyplot as plt                        # The script is usually run using pypy . For this reason, this is commented. Uncomment to be able to plot..
from array import array

'''
This is the code used for the paper

Hipsters on Networks: How a Small Group of Individuals Can Lead to an Anti-Establishment Majority, by Jonas S. Juul & Mason A. Porter.

If you wish to use this script, the paper should be cited. 
Please do not distribute this script without the consent of the authors.

For speed, run this script using pypy.
'''

# Definitions :
def makegraph(N,degree_sequence)  :                     #(number of nodes, [z1, z2], [phi1,phi2], p(z1))
    Graph = nx.configuration_model(degree_sequence)
    Adjacency = np.zeros((N,N))
    for edge in Graph.edges() :
        Adjacency[edge[0]][edge[1]] +=1
        Adjacency[edge[1]][edge[0]] +=1
    return Adjacency


write_to_file = True # Set to "True" if output should be printed to file.

# Now choose which of the following networks the model should be realised on:
#  z1z2 : Network of 2 degrees, z1 and z2. If this code is not changed z1=z2 and this choice results in a z-regular graph.
#  ER   : ER-network with expected mean degree z 
#  Facebook_original : Northwestern25-network from the Facebook100 dataset. To use this network, one should have access to a file which lists the edges of the network (here Northwestern_edges.txt)
network_types = ['z1z2', 'ER', 'Facebook_original']
network = network_types[0] #Choose which of the 3 networks to choose.

N = 10000                                               # number of nodes in synthetic network
p0 = .8                                                 # In networks with 2 kinds of nodes (characterised by degree or threshold) this is the probability of a node being the first kind. 
prob = [p0, 1-p0]                                       # Probabilities of the two kinds of nodes.
pHip = .00                                              # Fraction of hipster-nodes in network. This parameter is also the minimal pHip-value in the pHip-parameter scan beneath
dpHip = .01                                             # Increment of pHip in pHip-parameter scan 
pHip_max = 1.01                                         # pHip-parameter scan stops when this pHip-value (or larger) is reached
z = 5                                                   # Parameter used in z1z2 and ER-network choices. See above
if (network == 'z1z2') :                                # if network choice is 'z1z2', the list of degrees is [z,z]. Can be changed to [z1,z2] in which nodes can have one of two different degrees.
    k_list = [z,z]
    phi_list = [0.1,0.80]  
phi_star = 1/33.                                        # Threshold given to all nodes, if all nodes are assigned the same threshold.
timesteps = 40                                          # number of timesteps in each experiment (here chosen much larger than it takes for steady state to be reached).
Nexp = 200                                              # number of experiments for each pHip-value.
tau = 4                                                 # Delay in the knowledge of total product distribution for hipster nodes.           
pHip_vec = np.zeros(int((pHip_max-pHip)/dpHip+.5))      # Array of pHip-values from pHip-parameter scan.

# Define arrays
Mean_rho_vecs = np.zeros((2,len(pHip_vec)))             # Array for saving mean fraction of nodes which adopted product each product for particular pHip
Var_rho_vecs = np.zeros((2,len(pHip_vec)))              # Array for saving variance of the above mean for particular pHip

Low_rho_vec = array('f',[])                             # In this array the fraction of nodes that has adopted any product is saved, if the realization is discarded. (For Table I in paper)

entry_pHip = 0                                          # Index which changes for each pHip-value. Increases by one for each pHip-value

while (pHip < pHip_max) :                               # This is pHip-parameter scan. Minimum pHip-value is above defined pHip
    r_sing = np.zeros((2,Nexp))                         # Array with final fraction of nodes adopting each producted for all realisations and a single pHip-value. From this, mean adopted fraction and variance is calculated.
    r_var = 0                                           # Entry in r_sing which a value is safed into. Increases by one for each succesful (non-discarded) realisation

    rho_fin = np.zeros((2,timesteps))                   # Fraction of nodes which has adopted each product in a single realisation (timeseries array plotted in FIG 2.)

    exp = 0                                             # Counter of relisations. Increases by one if realisation is not discarded. 
    while (exp < Nexp) :                                # pHip-value increases by one when exp = Nexp
        #print pHip, exp                                # Uncomment to get output for each realisation. Nice for long simulations.

        Total_adopted = 1                                # Number of nodes that adopted a product (are "active")

        if (network == 'z1z2') :                        # If chosen network is 'z1z2', define an array with degrees and another with thresholds
            deg = array('i', [0]*N)                     # array with degrees of nodes
            phi = np.ones(N)                            # array with thresholds of nodes

        if (network == 'Facebook_original') :           # If chosen network is 'Facebook_original', define the following
            N = 10567                                   # Number of nodes in specific network
            k_max = 2105                                # Maximal degree in this network
            phi = np.ones(N)*phi_star                   # Array of thresholds. Every node has the same threshold 'phi_star'.
            deg = array('i', [0]*N)                     # Array with degrees of nodes
            Ndegree = np.zeros(k_max)                   # Array which on entry k with contain an integer: the number of nodes with degree k+1

            A = np.zeros((10567,10567))                 # Make adjacency matrix (brute force an ugly method)
            c0 = array('f',[])                          # Temporary array. If Northwestern25 edge E connects (i,j), 'i' will be appended to this array
            c1 = array('f',[])                          # Temporary array. If Northwestern25 edge E connects (i,j), 'j' will be appended to this array

            with open('Facebook100/Northwestern_edges.txt') as source_file: # File with edges in Northwestern25. Format: node1\tnode2
                for line in source_file:                # Save edges in the two temporary 'column-arrays' c0 and c1
                    cols = [float(x) for x in line.split()]
                    c0.append(cols[0])
                    c1.append(cols[1])

            for i in range (len(c0)) :                  # Make adjacency matrix for imported edges. In file, smallest node index is 1. Here it is 0.
                A[int(c0[i]-1)][int(c1[i]-1)] = 1
                A[int(c1[i]-1)][int(c0[i]-1)] = 1
            for node_num in range (N) :                 # Run the following for all nodes
                node_degree = sum(A[:][node_num]) +0    # Find degree of node
                deg[node_num] += node_degree            # Save the degree in degree array
                Ndegree[int(node_degree) - 1] +=1       # Add onde to entry k-1 in array Ndegree. This keeps track of the number of nodes having degree k.

        if (network == 'ER') :                          # If chosen network type is ER
            G = nx.fast_gnp_random_graph(N,z/float(N))  # Make an ER-graph with N nodes and mean degree z
            A = np.zeros((N,N))                         # Make an adjacency matrix too.
            edges = G.edges()                           # Edges of created ER-graph
            for edge in range (np.shape(G.edges())[0]) :# Now make corresponding adjacency matrix
                A[edges[edge][0]][edges[edge][1]] +=1   # Save 1 on entry [i,j] if edge (i,j) exists
                A[edges[edge][1]][edges[edge][0]] +=1   # Save 1 on entry [i,j] if edge (i,j) exists (symmetric adjacency)
            deg = sum(A)                                # Make array with degrees of nodes
            phi = np.ones(N)*phi_star                   # Make array with thresholds of nodes (all have same threshold: 'phi_star')
            Ndegree = np.zeros(int(max(deg)))           # Array which on entry k contains an integer: the number of nodes with degree k+1 
            for index in range (len(Ndegree)) :
                Ndegree[index] = sum(deg==index+1) # Insert number of nodes with degree k on entry k of array Ndegree

        if (network == 'z1z2') :                        # If network choice is 'z1z2'
            for index in range (len(deg)) :             # Above deg was defined.

                if (random.uniform(0,1)<prob[0]) :      # Draw number uniformly at random. If less than prob[0]: Make first kind of node
                    deg[index] = int(k_list[0])         # Give node first kind of degree
                    phi[index] = phi_list[0]            # Give node first kind of threshold
                else :                                  # if random number was not less than prob[0], make second kind of node.
                    deg[index] = int(k_list[1])         # Give node second kind of degree
                    phi[index] = phi_list[1]            # Give node second kind of treshold
            A = makegraph(N,deg)                        # Make configuration-model network with given degrees

        S = np.zeros(N)                                 # Array with "states". S[i] = 0 if node i is inactive (have not adopted no producted), S[i] = 1 if adopted product 1, S[i] if adopted product 2.
        H = np.zeros(N)                                 # Array with "hipster values". H[i] = 0 if nodes i is conformist, H[i] = 1 if Hipster
        for i in range (N) :                            # For all nodes, now assign "hipster-values"
            if (random.uniform(0,1) < pHip) :           # Draw random number uniformly at random. If less than pHip, node is hipster
                H[i] = 1.                               
        change_status = np.zeros(N)                     # Will be used in network dynamics - will node change status from S[i]=0 to S[i] \ge 1 or not?

        seed1 = int(random.uniform(0,N))                # Choose seed uniformly at random. "Seed" refers to the node which starts out with product 1 in this realisation
        S[seed1] = 1.                                   # activate the above chosen seed at t=0
        rho = np.zeros((2,timesteps))                   # array containing Fraction of active nodes (one row for each product)

        rho[0][0] += 1/float(N*Nexp)                    # One node already adopted product 1

        for t in range (1,timesteps) :                  # Run the model for timesteps 1 to 'timesteps'

            for index in range (np.shape(rho)[0]) :     # At this point, the fraction of nodes that have adopted product m is equal to the number of nodes that had adopted it at last time step.
                rho[index][t] += rho[index][t-1]
            tr = 0.                                     # node index in the following loop
            while (tr<N) :
                node = int(tr)                          # In this loop, examine whether node = int(tr) will change adopt a product in this time step
                tr +=1                                      
                     
                if (S[node] == 0.) :                    # If node in question has not adopted a product, it is now examined if it will adopt one in this time step
                    active_neighbour = 0                # number of active neighbours of node (will be counted beneath)
                    nb_1 = 0.                           # NUmber of active neighbours that have adopted "product 1" (will be counted)
                    nb_2 = 0.                           # NUmber of active neighbours that have adopted "product 2" (will be counted)
                    for neighbour in range (N) :        # Loop through all nodes in the network
                        if (S[neighbour] == 1. and A[node][neighbour] >0) : # 'neighbour' has adopted product 1 and is a neighbour of 'node'
                            active_neighbour +=A[node][neighbour]   # add one to "active_neighbour" counting
                            nb_1 +=1                    # Add one to counting of number of neighbours that have adopted product 1
                        if (S[neighbour] == 2. and A[node][neighbour] >0) :# 'neighbour' has adopted product 2 and is a neighbour of 'node'
                            active_neighbour +=A[node][neighbour]   # add one to "active_neighbour" counting
                            nb_2 +=1                    # Add one to counting of number of neighbours that have adopted product 2

                    if (active_neighbour >= phi[node]*deg[node]) : #Check if the number of fraction of active neighbours is at least the threshold of the nodes
                        Total_adopted +=1               # Add 1 to the number of nodes that are active
                        if (H[node] == 0.) :            # Check if node is conformist or hipster. If conformist:
                            if (nb_1> nb_2) :           
                                change_status[node] = 1.# Adopt product 1 if more neighbours adopted product 1 than product 2
                            elif (nb_1 == nb_2) :       # Adopt either product if equally many neighbours adopted product 1 and product 2
                                change_status[node] = int(1+int(random.uniform(0,1)+0.5))
                            else :
                                change_status[node] =2. # Adopt product 2 if more neighbours adopted product 2 than product 1
                        if (H[node] == 1.) :            # If the node is hipster...
                                if (t>=tau) :           # Check if time is larger than delay for hipsters
                                    t_hip = int(t-tau)  # If time is larger than delay, define time at which hipsters look check population product distribution as t-tau
                                else :
                                    t_hip = 0           # If time is NOT larger than delay, define time at which hipsters look check population product distribution as 0
                                if (rho[0][t_hip] >rho[1][t_hip]) : # If more larger fraction in population adopted product 1 at time at which hipsters know product distribution. Hipster adopts product 2
                                    change_status[node] = 2.
                                elif (rho[0][t_hip] <rho[1][t_hip]) :# If more larger fraction in population adopted product 2 at time at which hipsters know product distribution. Hipster adopts product 1
                                    change_status[node] = 1.
                                else :                  # If equally many adopted each product at time at which hipsters know product distribution. Hipster adopts each product with equal probability
                                    change_status[node] = int(1+int(random.uniform(0,1)+0.5))
                
            for node in range (N) :                     # Now: If node was set to adopt a product. Adopt the product
                if (S[node] == 0.) :                   
                    S[node] = change_status[node]       # If node had adopted no product before, change S[node] to equal product which it now adopts.

                    if (change_status[node] == 1.) :    # If node adopted product 1 this time step, add 1/(N*Nexp) to the fraction of nodes that had adopted product 1 this time step.
                        rho[0][t] += 1/float(N*Nexp)

                    if (change_status[node] == 2.) :    # If node adopted product 2 this time step, add 1/(N*Nexp) to the fraction of nodes that had adopted product 2 this time step.
                        rho[1][t] += 1/float(N*Nexp)



        if (Total_adopted/float(N) > 0.1) :             # Only go on to next experiment if larger fraction than 0.1 chose to adopt a product. 
            for t in range (timesteps) :                # Save fraction of nodes that had adopted a product at each time step. This is not used when a pHip-parameter scan is made.
                for index in range (2) :
                    rho_fin[index][t] += rho[index][t]
            exp +=1                                     # Increase experiment-counter by 1
            r_sing[0][r_var] += rho[0][-1]*Nexp         # Save final fraction of nodes that adopted product 1 in this realisation. Mean and variance of all entries in the array will be calculated when parameters are changed.
            r_sing[1][r_var] += rho[1][-1]*Nexp         # Save final fraction of nodes that adopted product 2 in this realisation. Mean and variance of all entries in the array will be calculated when parameters are changed.
            r_var +=1                                   # Increase index used when saving in the array 'r_sing' by one
        else :
            Low_rho_vec.append(Total_adopted/float(N))  # If less than 0.1 adopted a product in this realisation, save how large a fraction actually adopted a product (for Table I)



    pHip_vec[entry_pHip] = pHip +0                      # Save current pHip-value to array

    Mean_rho_vecs[0][entry_pHip] = np.mean(r_sing[0][:]) +0 # Save mean of fraction of nodes that adopted product 1 in realisations with these parameter.
    Mean_rho_vecs[1][entry_pHip] = np.mean(r_sing[1][:]) +0 # Save mean of fraction of nodes that adopted product 2 in realisations with these parameter.
    Var_rho_vecs[0][entry_pHip] = np.var(r_sing[0][:])+0# Save variance of fraction of nodes that adopted product 1 in realisations with these parameter.
    Var_rho_vecs[1][entry_pHip] = np.var(r_sing[1][:])+0# Save variance of fraction of nodes that adopted product 2 in realisations with these parameter.
    pHip += dpHip                                       # Increase fraction of hipster nodes ('dpHip') by amount 'dpHip'
    entry_pHip += 1                                     # Increase the entry used for saving in array pHip_vec by one.

if (write_to_file == True) :                            # If 'write_to_file' is set to True, Print the following to a file.
    file = open("Hipsters_network_%s_phi_%s_Nexp_%s_tau=_%s_pHip=%s.txt"%(network,phi_star,Nexp,tau,pHip-dpHip),"w")
    file.write("Array of fraction that adopted product 1 is \n frac_1 = %s \n"%(rho_fin[0][:]))
    file.write("Array of fraction that adopted product 2 \n frac_2 = %s \n"%(rho_fin[1][:]))
    file.write("Prob[0] was \n p0 = %s \n"%(p0))
    file.write("The total number of experiments was Nexp = %3.0i \n"%(Nexp))
    file.write("The pHip-values were = %s \n"%(pHip_vec))
    file.write("Runs that were not counted: full array %s \n \n"%(Low_rho_vec)) 
    file.write("Mean of fraction of product 1 adopters array and variance of product 1 adopters array are \n Mean_0 = %s \n Var_0 =  %s \n \n"%(Mean_rho_vecs[0][:], Var_rho_vecs[0][:])) 
    file.write("Mean of fraction of product 2 adopters array and variance of product 2 adopters array are \n Mean_1 = %s \n Var_1 = %s \n \n"%(Mean_rho_vecs[1][:], Var_rho_vecs[1][:])) 
    file.close()