import scipy.sparse as sparse
import numpy as np
import matrix_functions as mf
import scipy.misc as misc
import networkx as nx
import numpy.random as rand
import dwave_networkx as dnx
import minorminer as mm
import time as time

#written by Nicholas Chancellor
#Based on the work in Nicholas Chancellor 2019 Quantum Sci. Technol. 4 045004 https://doi.org/10.1088/2058-9565/ab33c2
#feel free to resuse/modify but please attribute the source and cite the paper in any published work

def make_var_encode_core(qubits,Ntot,encType):
	if encType=='one_hot':
		# makes a one hot Ising encoding for a single variable over the qubits 'qubits' given Nqubits total qubits
		Nencode=len(qubits) # length of the list of qubits determines the size of the variable
		h=np.zeros(Ntot) # empty list for putting fields in
		h[qubits]=-(Nencode-1) # local fields cancel precisely when one qubit is in the one configuration
		J=np.zeros([Ntot,Ntot]) # make matrix of zeros
		J[np.kron(qubits,np.ones(len(qubits),dtype=int)),np.kron(np.ones(len(qubits),dtype=int),qubits)]=1 # connect all elements
		J=np.triu(J) # remove lower triangular part
		J=J-np.diag(np.diag(J)) # remove diagonal part
	elif encType=='domain_wall':
		# makes a domain wall encoding for a single variable over the qubits 'qubits'
		h=np.zeros(Ntot) # empty list for putting fields in
		J=np.zeros([Ntot,Ntot]) # make matrix of zeros
		if len(qubits)>1: # only if there is more than one qubit in the variable		
			h[qubits[0]]=1 # encourage first variable to be one
			h[qubits[len(qubits)-1]]=-1 # encourages last variable to be zero
			J[qubits[0:(len(qubits)-1)],qubits[1:len(qubits)]]=-1 # ferromagnetic interactions
			J[qubits[1:len(qubits)],qubits[0:(len(qubits)-1)]]=-1 # ferromagnetic interactions
			J=np.triu(J) # upper triangular part
	else:
		print('encoding type not recognized, returinging NaN values, this will cause errors downstream')
		J=np.nan
		h=np.nan
	return (J,h)


def two_var_interaction(qubits1,qubits2,val1,val2,Ntot,encType):
	if encType=='one_hot': # one hot case
		# encodes an interaction penalizing the simultaneous case where variable 1 takes val1 and variable 2 takes val 2
		h=np.zeros(Ntot) # no fields needed for encoding
		J=np.zeros([Ntot,Ntot]) # initialize empty J
		J[qubits1[val1],qubits2[val2]]=1 # coupling in one direction
		J[qubits2[val2],qubits1[val1]]=1 # coupling in the other
		J=np.triu(J) # only keep upper triangular part
	elif encType=='domain_wall': # domain wall case
		# encodes an interaction penalizing the simultaneous case where variable 1 takes val1 and variable 2 takes val 2
		h=np.zeros(Ntot) # empty list for putting fields in
		J=np.zeros([Ntot,Ntot]) # make matrix of zeros
		if val1==0: # special case where some couplers are replaced by fields
			if val2==0:	# special case where some coupers are replaced by fields
				J[qubits1[0],qubits2[0]]=1 # anti-ferro coupling
				h[qubits1[0]]=1 # encourages qubit to be 1
				h[qubits2[0]]=1 # encourages qubit to be 1
			elif val2==len(qubits2):	# special case where some coupers are replaced by fields
				J[qubits1[0],qubits2[len(qubits2)-1]]=-1 # ferro coupling
				h[qubits1[0]]=1 # encourages qubit to be 1
				h[qubits2[len(qubits2)-1]]=-1 # encourages qubit to be 0
			else: # case where no interactions with val2 are replaced by fields
				J[qubits1[0],qubits2[val2-1]]=-1 # ferro coupling
				J[qubits1[0],qubits2[val2]]=1 # anti-ferro coupling
				h[qubits2[val2-1]]=-1 # encourages qubit to be 0
				h[qubits2[val2]]=1 # encourages qubit to be 1
		elif  val1==len(qubits1): # special case where some couplers are replaced by fields
			if val2==0:	# special case where some coupers are replaced by fields
				J[qubits1[len(qubits1)-1],qubits2[0]]=-1 # ferro coupling
				h[qubits1[len(qubits1)-1]]=-1 # encourages qubit to be 0
				h[qubits2[0]]=1 # encourages qubit to be 1
			elif val2==len(qubits2):	# special case where some coupers are replaced by fields
				J[qubits1[len(qubits1)-1],qubits2[len(qubits2)-1]]=1 # anti-ferro coupling
				h[qubits1[len(qubits1)-1]]=-1 # encourages qubit to be 0
				h[qubits2[len(qubits2)-1]]=-1 # encourages qubit to be 0
			else: # case where no interactions with val2 are replaced by fields
				J[qubits1[len(qubits1)-1],qubits2[val2-1]]=1 # anti-ferro coupling
				J[qubits1[len(qubits1)-1],qubits2[val2]]=-1 # ferromagnetic coupling
				h[qubits2[val2-1]]=-1 # encourages qubit to be 0
				h[qubits2[val2]]=1 # encourages qubit to be 1
		else: # val1 is not at either end point of the chain
			if val2==0:	# special case where some coupers are replaced by fields
				J[qubits1[val1],qubits2[0]]=1 # anti-ferro coupling
				J[qubits1[val1-1],qubits2[0]]=-1 # ferro coupling
				h[qubits1[val1]]=1 # encourages qubit to be 1
				h[qubits1[val1-1]]=-1 # encourages qubit to be 0
			elif val2==len(qubits2):	# special case where some coupers are replaced by fields
				J[qubits1[val1],qubits2[len(qubits2)-1]]=-1 # ferro coupling
				J[qubits1[val1-1],qubits2[len(qubits2)-1]]=1 # anti-ferro coupling
				h[qubits1[val1]]=1 # encourages qubit to be 1
				h[qubits1[val1-1]]=-1 # encourages qubit to be 0
			else: # case where no values are at the end of chains
				J[qubits1[val1],qubits2[val2]]=1 # anti-ferro coupling
				J[qubits1[val1-1],qubits2[val2-1]]=1 # anti-ferro coupling
				J[qubits1[val1],qubits2[val2-1]]=-1 # ferro coupling
				J[qubits1[val1-1],qubits2[val2]]=-1 # ferro coupling	
		J=J+J.T # make symetric
		J=np.triu(J) # make upper triangular
		J=J/4 # factor of two between Ising and Qubo
		h=h/4 # factor of two between Ising and Qubo	
	else:
		print('encoding type not recognized, returinging NaN values, this will cause errors downstream')
		J=np.nan
		h=np.nan
	return (J,h)

def generate_color_encoding(graph,nColor):
	Ntot_oh=len(graph)*nColor # total number of qubits for one hot
	J_oh_core=np.zeros([Ntot_oh,Ntot_oh]) # J matrix for one hot cores
	h_oh_core=np.zeros(Ntot_oh) # h vector for one hot cores
	J_oh_int=np.zeros([Ntot_oh,Ntot_oh]) # J matrix for one hot interactions
	h_oh_int=np.zeros(Ntot_oh) # h vector for one hot interactions
	Ntot_dw=len(graph)*(nColor-1) # total number of qubits for domain wall
	J_dw_core=np.zeros([Ntot_dw,Ntot_dw]) # J matrix for domain wall cores
	h_dw_core=np.zeros(Ntot_dw) # h vector for domain wall cores
	J_dw_int=np.zeros([Ntot_dw,Ntot_dw]) # J matrix for domain wall interactions
	h_dw_int=np.zeros(Ntot_dw) # h vector for domain wall interactons
	for iCore in range(len(graph)): # build varaible cores for each encoding type
		(J_oh1,h_oh1)=make_var_encode_core(list(range((iCore*nColor),(iCore+1)*nColor)),Ntot_oh,'one_hot') # core for one hot encoding
		J_oh_core=J_oh_core+J_oh1 # add to couplings
		h_oh_core=h_oh_core+h_oh1 # add to fields
		(J_dw1,h_dw1)=make_var_encode_core(list(range((iCore*(nColor-1)),(iCore+1)*(nColor-1))),Ntot_dw,'domain_wall') # core for domain wall encoding
		J_dw_core=J_dw_core+J_dw1 # add to couplings
		h_dw_core=h_dw_core+h_dw1 # add to fields
	for iCoup in range(len(graph)): # loop over first graph index
		for jCoup in range(iCoup+1,len(graph)): # loop over second graph index
			if graph[iCoup,jCoup]==1: # if the edge is included in the graph
				qubits1_oh=list(range((iCoup*nColor),((iCoup+1)*nColor))) # qubits involved in first one hot variable
				qubits2_oh=list(range((jCoup*nColor),((jCoup+1)*nColor))) # qubits involved in second one hot variable
				qubits1_dw=list(range((iCoup*(nColor-1)),((iCoup+1)*(nColor-1)))) # qubits involved in first one hot variable
				qubits2_dw=list(range((jCoup*(nColor-1)),((jCoup+1)*(nColor-1)))) # qubits involved in second one hot variable
				for iColor in range(nColor): # loop over colors
					(J_oh1,h_oh1)=two_var_interaction(qubits1_oh,qubits2_oh,iColor,iColor,Ntot_oh,'one_hot') # generate J and h values for one hot
					J_oh_int=J_oh_int+J_oh1 # add to couplings
					h_oh_int=h_oh_int+h_oh1 # add to fields
					(J_dw1,h_dw1)=two_var_interaction(qubits1_dw,qubits2_dw,iColor,iColor,Ntot_dw,'domain_wall') # generate J and h values for domain wall
					J_dw_int=J_dw_int+J_dw1 # add to couplings
					h_dw_int=h_dw_int+h_dw1 # add to fields
	one_hot_enc={} # dictionary for storing one hot encoding
	one_hot_enc['J_core']=J_oh_core # core coupling
	one_hot_enc['h_core']=h_oh_core # core fields
	one_hot_enc['J_int']=J_oh_int # interaction coupling
	one_hot_enc['h_int']=h_oh_int # interaction fields
	domain_wall_enc={} # dictionary for storing domain wall encoding
	domain_wall_enc['J_core']=J_dw_core # core coupling
	domain_wall_enc['h_core']=h_dw_core # core fields
	domain_wall_enc['J_int']=J_dw_int # interaction coupling
	domain_wall_enc['h_int']=h_dw_int # interaction fields
	return (one_hot_enc,domain_wall_enc)

def J_to_edges(J): # converts interaction matrix J to a list of edges
	edges=[] # empt list of edges
	for iJ in range(len(J)): # loop over fist index
		for jJ in range(iJ+1,len(J)): # loop over second index
			if not J[iJ,jJ]==0: # if an interaction is present
				edges=edges+[(iJ,jJ)] # add edge to list
	return edges # return the list of edges

def make_graph_from_edges(edges): # make a networkx graph from a list of edges (useful for plotting, etc...)
	G=nx.Graph()
	G.clear()
	G.add_edges_from(edges)
	return(G)

def find_smallest_square_chimera_emb(J,start_size,*args,**kwargs): # finds the smallest square chimera which the coupling matrix J is embedable into
	edges=J_to_edges(J) # create list of edges 
	lin_size=start_size # start with given linear size
	chi=dnx.chimera_graph(lin_size,lin_size,4) # create chimera graph
	if 'chainlength_patience' in kwargs: # externally set chain length patience
		emb=mm.find_embedding(edges,chi,chainlength_patience=kwargs['chainlength_patience']) # try to create embedding
	else:
		emb=mm.find_embedding(edges,chi) # try to create embedding
	if emb=={}: # if the embedding has failed increment size until succes
		while emb=={}: # while the embedding has failed
			lin_size=lin_size+1 # increment linear size
			chi=dnx.chimera_graph(lin_size,lin_size,4) # create chimera graph
			if 'chainlength_patience' in kwargs: # externally set chain length patience
				emb=mm.find_embedding(edges,chi,chainlength_patience=kwargs['chainlength_patience']) # try to create embedding
			else:
				emb=mm.find_embedding(edges,chi) # try to create embedding
	else: # otherwise decrement linear size until failure occurs or single unit cell
		emb1=emb # dummy variable to prevent empty embedding from being output	
		if lin_size==1: # let us know if we started with a linear size of 1
			start1=1
		else:
			start1=0	
		while not ((emb1=={}) or (lin_size==1)): # stop if embedding fails or single unit cell
			emb=emb1 # update emb with non-failed value
			lin_size=lin_size-1 # decrement linear size
			chi=dnx.chimera_graph(lin_size,lin_size,4) # create chimera graph
			if 'chainlength_patience' in kwargs: # externally set chainlength patience
				emb1=mm.find_embedding(edges,chi,chainlength_patience=kwargs['chainlength_patience']) # try to create embedding
			else:
				emb1=mm.find_embedding(edges,chi) # try to create embedding
		if start1==0: # only undo decrement if did not start with size 1
			lin_size=lin_size+1 # undo last decrement
	return (lin_size,emb) # return linear size and smallest embedding


def count_qubits(emb): # function for counting total number of qubits involved in an embedding
	nEmb=0 # start with no auxilliaries for embedding
	for ichain in range(len(emb)):
		if ichain in emb: # make sure qubit has been emebedded
			nEmb=nEmb+len(emb[ichain]) # add to number of auxilliaries
	return(nEmb)


def smallest_chimera_embedding_independent_random_color(edgeProb,nColor,sizeList):
	one_hot_min_sizes=np.zeros(len(sizeList))
	domain_wall_min_sizes=np.zeros(len(sizeList))
	one_hot_counts=np.zeros(len(sizeList))
	domain_wall_counts=np.zeros(len(sizeList))
	one_hot_counts_noEmb=np.zeros(len(sizeList))
	domain_wall_counts_noEmb=np.zeros(len(sizeList))
	lin_size_oh=1 # start with chimera unit cell
	lin_size_dw=1 # start with chimera unit cell
	for iSize in range(len(sizeList)): # iterate over different graph sizes
		colorGraph_adj=np.array(rand.random([sizeList[iSize],sizeList[iSize]])<edgeProb,dtype=int) # create random graph to try to color
		(one_hot_enc,domain_wall_enc)=generate_color_encoding(colorGraph_adj,nColor) # create domain wall and one hot encodings
		startEmb_time=time.time() # time when embedding started
		(lin_size_oh,emb)=find_smallest_square_chimera_emb(one_hot_enc['J_core']+one_hot_enc['J_int'],lin_size_oh) # find the smallest linear size for one hot
		qubit_count_oh=count_qubits(emb) # number of qubits involved in one hot encoding
		(lin_size_dw,emb)=find_smallest_square_chimera_emb(domain_wall_enc['J_core']+domain_wall_enc['J_int'],lin_size_dw) # find the smallest linear size for domain wall
		qubit_count_dw=count_qubits(emb) # number of qubits involved in domain wall encoding
		etime=time.time()-startEmb_time # elapsed time for both embeddings
		print([sizeList[iSize],lin_size_oh,lin_size_dw,etime]) # print out different linear sizes and elapsed time
		one_hot_min_sizes[iSize]=lin_size_oh # store one hot linear size
		one_hot_counts[iSize]=qubit_count_oh # store one hot qubit count
		one_hot_counts_noEmb[iSize]=len(one_hot_enc['J_core']) # qubit count before embedding
		domain_wall_min_sizes[iSize]=lin_size_dw # store domain wall linear size
		domain_wall_counts[iSize]=qubit_count_dw # store domain wall qubit count
		domain_wall_counts_noEmb[iSize]=len(domain_wall_enc['J_core']) # qubit count before embedding
	return(one_hot_min_sizes,one_hot_counts,one_hot_counts_noEmb,domain_wall_min_sizes,domain_wall_counts,domain_wall_counts_noEmb) # return minimum sizes and qubit counts


def smallest_chimera_embedding_random_fixed_color_ratio(edgeProb,colorRatio,nColorList):
	one_hot_min_sizes=np.zeros(len(nColorList))
	domain_wall_min_sizes=np.zeros(len(nColorList))
	one_hot_counts=np.zeros(len(nColorList))
	domain_wall_counts=np.zeros(len(nColorList))
	one_hot_counts_noEmb=np.zeros(len(nColorList))
	domain_wall_counts_noEmb=np.zeros(len(nColorList))
	lin_size_oh=1 # start with chimera unit cell
	lin_size_dw=1 # start with chimera unit cell
	for inColor in range(len(nColorList)): # iterate over different graph sizes
		colorGraph_adj=np.array(rand.random([colorRatio*nColorList[inColor],colorRatio*nColorList[inColor]])<edgeProb,dtype=int) # create random graph to try to color
		(one_hot_enc,domain_wall_enc)=generate_color_encoding(colorGraph_adj,nColorList[inColor]) # create domain wall and one hot encodings
		startEmb_time=time.time() # time when embedding started
		(lin_size_oh,emb)=find_smallest_square_chimera_emb(one_hot_enc['J_core']+one_hot_enc['J_int'],lin_size_oh) # find the smallest linear size for one hot
		qubit_count_oh=count_qubits(emb) # number of qubits involved in one hot encoding
		(lin_size_dw,emb)=find_smallest_square_chimera_emb(domain_wall_enc['J_core']+domain_wall_enc['J_int'],lin_size_dw) # find the smallest linear size for domain wall
		qubit_count_dw=count_qubits(emb) # number of qubits involved in domain wall encoding		
		etime=time.time()-startEmb_time # elapsed time for both embeddings
		print([nColorList[inColor],lin_size_oh,lin_size_dw,etime]) # print out different linear sizes and elapsed time
		one_hot_min_sizes[inColor]=lin_size_oh # store one hot linear size
		one_hot_counts[inColor]=qubit_count_oh # store one hot qubit count
		one_hot_counts_noEmb[inColor]=len(one_hot_enc['J_core']) # qubit count before embedding
		domain_wall_min_sizes[inColor]=lin_size_dw # store domain wall linear size
		domain_wall_counts[inColor]=qubit_count_dw # store one hot qubit count
		domain_wall_counts_noEmb[inColor]=len(domain_wall_enc['J_core']) # qubit count before embedding
	return(one_hot_min_sizes,one_hot_counts,one_hot_counts_noEmb,domain_wall_min_sizes,domain_wall_counts,domain_wall_counts_noEmb) # return minimum sizes and qubit counts


def find_smallest_pegasus_emb(J,start_size,*args,**kwargs): # finds the smallest square chimera which the coupling matrix J is embedable into
	edges=J_to_edges(J) # create list of edges 
	if start_size==1: # check if start size is 1
		start_size=2 # produces empty graph when called with one, so call with 2
	lin_size=start_size # start with given linear size
	peg=dnx.pegasus_graph(lin_size) # create pegasus graph
	if 'chainlength_patience' in kwargs: # externally set chain length patience
		emb=mm.find_embedding(edges,peg,chainlength_patience=kwargs['chainlength_patience']) # try to create embedding
	else:
		emb=mm.find_embedding(edges,peg) # try to create embedding
	if emb=={}: # if the embedding has failed increment size until succes
		while emb=={}: # while the embedding has failed
			lin_size=lin_size+1 # increment linear size
			peg=dnx.pegasus_graph(lin_size) # create pegasus graph
			if 'chainlength_patience' in kwargs: # externally set chain length patience
				emb=mm.find_embedding(edges,peg,chainlength_patience=kwargs['chainlength_patience']) # try to create embedding
			else:
				emb=mm.find_embedding(edges,peg) # try to create embedding
	else: # otherwise decrement linear size until failure occurs or single unit cell
		emb1=emb # dummy variable to prevent empty embedding from being output
		if lin_size==2: # let us know if we started with a linear size of 2
			start2=1
		else:
			start2=0		
		while not ((emb1=={}) or (lin_size==2)): # stop if embedding fails or p_1
			emb=emb1 # update emb with non-failed value
			lin_size=lin_size-1 # decrement linear size
			peg=dnx.pegasus_graph(lin_size) # create pegasus graph
			if 'chainlength_patience' in kwargs: # externally set chainlength patience
				emb1=mm.find_embedding(edges,peg,chainlength_patience=kwargs['chainlength_patience']) # try to create embedding
			else:
				emb1=mm.find_embedding(edges,peg) # try to create embedding
		if start2==0: # only undo decrement if we did not start with a size of 2
			lin_size=lin_size+1 # undo last decrement
	return (lin_size,emb) # return linear size and smallest embedding


def smallest_pegasus_embedding_independent_random_color(edgeProb,nColor,sizeList):
	one_hot_min_sizes=np.zeros(len(sizeList))
	domain_wall_min_sizes=np.zeros(len(sizeList))
	one_hot_counts=np.zeros(len(sizeList))
	domain_wall_counts=np.zeros(len(sizeList))
	one_hot_counts_noEmb=np.zeros(len(sizeList))
	domain_wall_counts_noEmb=np.zeros(len(sizeList))
	lin_size_oh=2 # start with p_2
	lin_size_dw=2 # start with p_2
	for iSize in range(len(sizeList)): # iterate over different graph sizes
		colorGraph_adj=np.array(rand.random([sizeList[iSize],sizeList[iSize]])<edgeProb,dtype=int) # create random graph to try to color
		(one_hot_enc,domain_wall_enc)=generate_color_encoding(colorGraph_adj,nColor) # create domain wall and one hot encodings
		startEmb_time=time.time() # time when embedding started
		(lin_size_oh,emb)=find_smallest_pegasus_emb(one_hot_enc['J_core']+one_hot_enc['J_int'],lin_size_oh) # find the smallest linear size for one hot
		qubit_count_oh=count_qubits(emb) # number of qubits involved in one hot encoding
		(lin_size_dw,emb)=find_smallest_pegasus_emb(domain_wall_enc['J_core']+domain_wall_enc['J_int'],lin_size_dw) # find the smallest linear size for domain wall
		qubit_count_dw=count_qubits(emb) # number of qubits involved in domain wall encoding
		etime=time.time()-startEmb_time # elapsed time for both embeddings
		print([sizeList[iSize],lin_size_oh,lin_size_dw,etime]) # print out different linear sizes and elapsed time
		one_hot_min_sizes[iSize]=lin_size_oh # store one hot linear size
		one_hot_counts[iSize]=qubit_count_oh # store one hot qubit count
		one_hot_counts_noEmb[iSize]=len(one_hot_enc['J_core']) # qubit count before embedding
		domain_wall_min_sizes[iSize]=lin_size_dw # store domain wall linear size
		domain_wall_counts[iSize]=qubit_count_dw # store domain wall qubit count
		domain_wall_counts_noEmb[iSize]=len(domain_wall_enc['J_core']) # qubit count before embedding
	return(one_hot_min_sizes,one_hot_counts,one_hot_counts_noEmb,domain_wall_min_sizes,domain_wall_counts,domain_wall_counts_noEmb) # return minimum sizes and qubit counts


def smallest_pegasus_embedding_random_fixed_color_ratio(edgeProb,colorRatio,nColorList):
	one_hot_min_sizes=np.zeros(len(nColorList))
	domain_wall_min_sizes=np.zeros(len(nColorList))
	one_hot_counts=np.zeros(len(nColorList))
	domain_wall_counts=np.zeros(len(nColorList))
	one_hot_counts_noEmb=np.zeros(len(nColorList))
	domain_wall_counts_noEmb=np.zeros(len(nColorList))
	lin_size_oh=2 # start with p_2
	lin_size_dw=2 # start with p_2
	for inColor in range(len(nColorList)): # iterate over different graph sizes
		colorGraph_adj=np.array(rand.random([colorRatio*nColorList[inColor],colorRatio*nColorList[inColor]])<edgeProb,dtype=int) # create random graph to try to color
		(one_hot_enc,domain_wall_enc)=generate_color_encoding(colorGraph_adj,nColorList[inColor]) # create domain wall and one hot encodings
		startEmb_time=time.time() # time when embedding started
		(lin_size_oh,emb)=find_smallest_pegasus_emb(one_hot_enc['J_core']+one_hot_enc['J_int'],lin_size_oh) # find the smallest linear size for one hot
		qubit_count_oh=count_qubits(emb) # number of qubits involved in one hot encoding
		(lin_size_dw,emb)=find_smallest_pegasus_emb(domain_wall_enc['J_core']+domain_wall_enc['J_int'],lin_size_dw) # find the smallest linear size for domain wall
		qubit_count_dw=count_qubits(emb) # number of qubits involved in domain wall encoding		
		etime=time.time()-startEmb_time # elapsed time for both embeddings
		print([nColorList[inColor],lin_size_oh,lin_size_dw,etime]) # print out different linear sizes and elapsed time
		one_hot_min_sizes[inColor]=lin_size_oh # store one hot linear size
		one_hot_counts[inColor]=qubit_count_oh # store one hot qubit count
		one_hot_counts_noEmb[inColor]=len(one_hot_enc['J_core']) # qubit count before embedding
		domain_wall_min_sizes[inColor]=lin_size_dw # store domain wall linear size
		domain_wall_counts[inColor]=qubit_count_dw # store one hot qubit count
		domain_wall_counts_noEmb[inColor]=len(domain_wall_enc['J_core']) # qubit count before embedding
	return(one_hot_min_sizes,one_hot_counts,one_hot_counts_noEmb,domain_wall_min_sizes,domain_wall_counts,domain_wall_counts_noEmb) # return minimum sizes and qubit counts


def generate_schedule_encoding(start_times,duration,conflict_map):
	# start_times is a nEvent by 2 vector of the possible range of start times of each event
	# duration is a vector of length nEvent giving the duration of each event
	# conflict_map is an nEvent by nEvent binary upper triangular matrix which tells whether or not two events conflict
	nEvent=len(duration) # the total number of events
	minTime=min(start_times[:,0]) # the fist time used in the mapping
	maxTime=max(start_times[:,1]) # the maximum time used in the mapping
	qubit_ranges_oh=np.zeros([nEvent,2]) # vector for listing used for each variable in one hot encoding
	qubit_ranges_dw=np.zeros([nEvent,2]) # vector for listing used for each variable in domain wall encoding
	current_oh=0 # current qubit for one hot, initialized at zero
	current_dw=0 # current qubit for domain wall, initialized at zero
	for iFindRanges in range(nEvent): # for loop to find ranges
		qubit_ranges_oh[iFindRanges,0]=current_oh # start at current qubit index
		qubit_ranges_oh[iFindRanges,1]=current_oh+int(start_times[iFindRanges,1]-start_times[iFindRanges,0]) # determine number of qubits based on difference between first and last possible start times
		current_oh=current_oh+int(start_times[iFindRanges,1]-start_times[iFindRanges,0]+1) # update current qubit index
		qubit_ranges_dw[iFindRanges,0]=current_dw # start at current qubit index
		qubit_ranges_dw[iFindRanges,1]=current_dw+int(start_times[iFindRanges,1]-start_times[iFindRanges,0]-1) # determine number of qubits based on difference between first and last possible start times
		current_dw=current_dw+int(start_times[iFindRanges,1]-start_times[iFindRanges,0]) # update current qubit index
	Ntot_oh=current_oh # number of qubits involved in a one hot encoding
	Ntot_dw=current_dw # number of qubits involved in a domain wall encoding
	J_oh_core=np.zeros([Ntot_oh,Ntot_oh]) # J matrix for one hot cores
	h_oh_core=np.zeros(Ntot_oh) # h vector for one hot cores
	J_oh_int=np.zeros([Ntot_oh,Ntot_oh]) # J matrix for one hot interactions
	h_oh_int=np.zeros(Ntot_oh) # h vector for one hot interactions
	J_dw_core=np.zeros([Ntot_dw,Ntot_dw]) # J matrix for domain wall cores
	h_dw_core=np.zeros(Ntot_dw) # h vector for domain wall cores
	J_dw_int=np.zeros([Ntot_dw,Ntot_dw]) # J matrix for domain wall interactions
	h_dw_int=np.zeros(Ntot_dw) # h vector for domain wall interactons
	for iMakeCore in range(nEvent): # for loop to create the cores of each variable
		(J_oh1,h_oh1)=make_var_encode_core(list(range(int(qubit_ranges_oh[iMakeCore,0]),int(qubit_ranges_oh[iMakeCore,1]+1))),Ntot_oh,'one_hot') # core couplers and fields for one hot encoding
		J_oh_core=J_oh_core+J_oh1 # add couplers
		h_oh_core=h_oh_core+h_oh1 # add fields
		(J_dw1,h_dw1)=make_var_encode_core(list(range(int(qubit_ranges_dw[iMakeCore,0]),int(qubit_ranges_dw[iMakeCore,1]+1))),Ntot_dw,'domain_wall') # core couplers and fields for domain wall encoding
		J_dw_core=J_dw_core+J_dw1 # add couplers
		h_dw_core=h_dw_core+h_dw1 # add fields
	for iMakeInts in range(nEvent): # outer for loop to create interactions
		qubits1_oh=list(range(int(qubit_ranges_oh[iMakeInts,0]),int(qubit_ranges_oh[iMakeInts,1]+1))) # qubits in first variable for one hot
		qubits1_dw=list(range(int(qubit_ranges_dw[iMakeInts,0]),int(qubit_ranges_dw[iMakeInts,1]+1))) # qubits in first variable for domain wall
		for jMakeInts in range(iMakeInts+1,nEvent): # inner for loop to create interactions
			if not conflict_map[iMakeInts,jMakeInts]==0: # only encode the interaction if there actually is a potential for conflict
				qubits2_oh=list(range(int(qubit_ranges_oh[jMakeInts,0]),int(qubit_ranges_oh[jMakeInts,1]+1))) # qubits in second variable for one hot
				qubits2_dw=list(range(int(qubit_ranges_dw[jMakeInts,0]),int(qubit_ranges_dw[jMakeInts,1]+1))) # qubits in second variable for domain wall
				nPossibleFirst=int(start_times[iMakeInts,1]-start_times[iMakeInts,0]+1) # number of possible start times for first event
				nPossibleSecond=int(start_times[jMakeInts,1]-start_times[jMakeInts,0]+1) # number of possible start times for second event
				for iFirstStart in range(nPossibleFirst): # loop through all possible start times of first event
					firstStartTime=start_times[iMakeInts,0]+iFirstStart # start time of first event
					for iSecondStart in range(nPossibleSecond): # loop through all possible start times of second event
						secondStartTime=start_times[jMakeInts,0]+iSecondStart # start time of second event
						if (firstStartTime==secondStartTime) or ((firstStartTime>secondStartTime) and ((firstStartTime-secondStartTime)<duration[jMakeInts])) or ((firstStartTime<secondStartTime) and ((secondStartTime-firstStartTime)<duration[iMakeInts])): # conditions for which a conflict happens
							(J_oh1,h_oh1)=two_var_interaction(qubits1_oh,qubits2_oh,iFirstStart,iSecondStart,Ntot_oh,'one_hot') # interaction encoding for one hot
							J_oh_int=J_oh_int+J_oh1 # add couplers
							h_oh_int=h_oh_int+h_oh1 # add fields
							(J_dw1,h_dw1)=two_var_interaction(qubits1_dw,qubits2_dw,iFirstStart,iSecondStart,Ntot_dw,'domain_wall') # interaction encoding for one hot
							J_dw_int=J_dw_int+J_dw1 # add couplers
							h_dw_int=h_dw_int+h_dw1 # add fields
	one_hot_enc={} # dictionary for storing one hot encoding
	one_hot_enc['J_core']=J_oh_core # core coupling
	one_hot_enc['h_core']=h_oh_core # core fields
	one_hot_enc['J_int']=J_oh_int # interaction coupling
	one_hot_enc['h_int']=h_oh_int # interaction fields
	domain_wall_enc={} # dictionary for storing domain wall encoding
	domain_wall_enc['J_core']=J_dw_core # core coupling
	domain_wall_enc['h_core']=h_dw_core # core fields
	domain_wall_enc['J_int']=J_dw_int # interaction coupling
	domain_wall_enc['h_int']=h_dw_int # interaction fields
	return (one_hot_enc,domain_wall_enc)	

def generate_random_scheduling_problem(nEvent,total_time,maxDuration,conflictProbability):
	# generates a random secheduling problem with nEvent events with durations selected uniformly at random from 1 to maxDuration, a total time window of total_time, and a probability that any will conflict of conflictProbability
	# initial possible start times are selected uniformly at random from 0 to total_time-2, the final possible start time is selected uniformly at random from the initial+1 to total_time-1 (zero indexed)
	duration=np.zeros(nEvent) # initialize empty variable for durations
	start_times=np.zeros([nEvent,2]) # initialize empty variable for start times
	conflict_map=np.zeros([nEvent,nEvent]) # initialize empty variable for conflicts
	for iEvent in range(nEvent): # iterate events
		duration[iEvent]=rand.randint(1,maxDuration+1) # assign duration
		start1=rand.randint(0,total_time-2+1) # earliest start time
		start2=rand.randint(start1+1,total_time-1+1) # latest start time
		start_times[iEvent,:]=[start1,start2] # add to list of start times
		for jEvent in range(iEvent+1,nEvent): # iterate events for conflicts
			if rand.uniform()<conflictProbability: # determine if conflict should be added
				conflict_map[iEvent,jEvent]=1 # add conflict
	return(start_times,duration,conflict_map)

def smallest_chimera_embedding_random_fixed_time_nEvent_ratio(conflictProbability,maxDuration,time_nEvent_ratio,nEventList):
	one_hot_min_sizes=np.zeros(len(nEventList))
	domain_wall_min_sizes=np.zeros(len(nEventList))
	one_hot_counts=np.zeros(len(nEventList))
	domain_wall_counts=np.zeros(len(nEventList))
	one_hot_counts_noEmb=np.zeros(len(nEventList))
	domain_wall_counts_noEmb=np.zeros(len(nEventList))
	lin_size_oh=1 # start with chimera unit cell
	lin_size_dw=1 # start with chimera unit cell
	for inEvent in range(len(nEventList)): # iterate over different graph sizes
		(start_times,duration,conflict_map)=generate_random_scheduling_problem(nEventList[inEvent],time_nEvent_ratio*nEventList[inEvent],maxDuration,conflictProbability) # create random problem
		(one_hot_enc,domain_wall_enc)=generate_schedule_encoding(start_times,duration,conflict_map) # create domain wall and one hot encodings
		startEmb_time=time.time() # time when embedding started
		(lin_size_oh,emb)=find_smallest_square_chimera_emb(one_hot_enc['J_core']+one_hot_enc['J_int'],lin_size_oh) # find the smallest linear size for one hot
		qubit_count_oh=count_qubits(emb) # number of qubits involved in one hot encoding
		(lin_size_dw,emb)=find_smallest_square_chimera_emb(domain_wall_enc['J_core']+domain_wall_enc['J_int'],lin_size_dw) # find the smallest linear size for domain wall
		qubit_count_dw=count_qubits(emb) # number of qubits involved in domain wall encoding		
		etime=time.time()-startEmb_time # elapsed time for both embeddings
		print([nEventList[inEvent],lin_size_oh,lin_size_dw,etime]) # print out different linear sizes and elapsed time
		one_hot_min_sizes[inEvent]=lin_size_oh # store one hot linear size
		one_hot_counts[inEvent]=qubit_count_oh # store one hot qubit count
		one_hot_counts_noEmb[inEvent]=len(one_hot_enc['J_core']) # qubit count before embedding
		domain_wall_min_sizes[inEvent]=lin_size_dw # store domain wall linear size
		domain_wall_counts[inEvent]=qubit_count_dw # store one hot qubit count
		domain_wall_counts_noEmb[inEvent]=len(domain_wall_enc['J_core']) # qubit count before embedding
	return(one_hot_min_sizes,one_hot_counts,one_hot_counts_noEmb,domain_wall_min_sizes,domain_wall_counts,domain_wall_counts_noEmb) # return minimum sizes and qubit counts

def smallest_pegasus_embedding_random_fixed_time_nEvent_ratio(conflictProbability,maxDuration,time_nEvent_ratio,nEventList):
	one_hot_min_sizes=np.zeros(len(nEventList))
	domain_wall_min_sizes=np.zeros(len(nEventList))
	one_hot_counts=np.zeros(len(nEventList))
	domain_wall_counts=np.zeros(len(nEventList))
	one_hot_counts_noEmb=np.zeros(len(nEventList))
	domain_wall_counts_noEmb=np.zeros(len(nEventList))
	lin_size_oh=2 # start with p_2
	lin_size_dw=2 # start with p_2
	for inEvent in range(len(nEventList)): # iterate over different graph sizes
		(start_times,duration,conflict_map)=generate_random_scheduling_problem(nEventList[inEvent],time_nEvent_ratio*nEventList[inEvent],maxDuration,conflictProbability) # create random problem
		(one_hot_enc,domain_wall_enc)=generate_schedule_encoding(start_times,duration,conflict_map) # create domain wall and one hot encodings
		startEmb_time=time.time() # time when embedding started
		(lin_size_oh,emb)=find_smallest_pegasus_emb(one_hot_enc['J_core']+one_hot_enc['J_int'],lin_size_oh) # find the smallest linear size for one hot
		qubit_count_oh=count_qubits(emb) # number of qubits involved in one hot encoding
		(lin_size_dw,emb)=find_smallest_pegasus_emb(domain_wall_enc['J_core']+domain_wall_enc['J_int'],lin_size_dw) # find the smallest linear size for domain wall
		qubit_count_dw=count_qubits(emb) # number of qubits involved in domain wall encoding		
		etime=time.time()-startEmb_time # elapsed time for both embeddings
		print([nEventList[inEvent],lin_size_oh,lin_size_dw,etime]) # print out different linear sizes and elapsed time
		one_hot_min_sizes[inEvent]=lin_size_oh # store one hot linear size
		one_hot_counts[inEvent]=qubit_count_oh # store one hot qubit count
		one_hot_counts_noEmb[inEvent]=len(one_hot_enc['J_core']) # qubit count before embedding
		domain_wall_min_sizes[inEvent]=lin_size_dw # store domain wall linear size
		domain_wall_counts[inEvent]=qubit_count_dw # store one hot qubit count
		domain_wall_counts_noEmb[inEvent]=len(domain_wall_enc['J_core']) # qubit count before embedding
	return(one_hot_min_sizes,one_hot_counts,one_hot_counts_noEmb,domain_wall_min_sizes,domain_wall_counts,domain_wall_counts_noEmb) # return minimum sizes and qubit counts
