class ReferenceAlignment(Step):
	def __init__(self, workspace, merge, ref_file):
		self.workspace=workspace
		self.merge=merge
		self.ref_file=ref_file
		self.quality=None

		file_data=self.ref_file.split('/')
		prefix_data=file_data[len(file_data)-1].split('.')

		self.output_prefix=prefix_data[0]
		self.send_output_to_file=True
		self.send_error_to_file=True
		self.output_veto_regex="_intervals.txt$"
		self.res=2.9
		self.pval=1e-10 ###
		self.fp=0.6 ###
		self.fn=0.06 ###
		self.sf=0.20
		self.sd=0.10
#		self.sd=0.20
#		self.sr=0.03
		self.allow_overhang=True
		self.outlier_pval=0.0001
		self.end_outlier_pval=0.001
		self.max_query_alignment_interval=12
		self.max_reference_alignment_interval=12
		self.min_sites_for_chimera=14
		self.hash_window=5
		self.hash_min_sites=3
		self.hash_sd_max=2.2
		self.hash_sd_rms=1.2
		self.hash_relative_error=0.05
		self.hash_offset_kb=3.0
		self.hash_max_insert_errors=1
		self.hash_max_probe_errors=1
		self.hash_max_unresolved_sites=1
		self.hash_delta=50
		self.target_resolution=1e-3
		self.resolution_reduction=1.2
		self.allow_no_splits=True
		self.allow_infinite_splits=False
		self.scale_bias_wt=0
		self.overwrite_output=True
		self.print_indel_file=True

		self.autoGeneratePrereqs()

	def __hash__(self):
		return hash((self.workspace.input_file, self.workspace.work_dir, self.ref_file))

	def __str__(self):
		return("Comparison of " + self.workspace.input_file + " to " + self.ref_file)

	@staticmethod
	def generateFromStepDir(step_dir, blocks=None):
		match_result=re.match("comparison_([^_]+)_([\d\.]+)_([\d\.]+)_([\d\.e-]+)_([\d]+)_([\d]+)_([^/]+)", step_dir)

		work_dir=os.getcwd()
		input_file=match_result.group(1)
		workspace=Workspace(work_dir, input_file)

		fp=match_result.group(2)
		fn=match_result.group(3)
		pval=match_result.group(4)
		minlen=match_result.group(5)
		minsites=match_result.group(6)
		vital_parameters=VitalParameters(fp, fn, pval, minlen, minsites)
		vital_parameters.blocks=blocks

		merge=Merge(workspace, vital_parameters)

		ref_file=match_result.group(7)

		return ReferenceAlignment(workspace, merge, ref_file)

	def writeCode(self):
		code="cd " + self.workspace.work_dir + "\n"
		code+="mkdir " + self.getStepDir() + "\n"
		code+="cd " + self.getStepDir() + "\n"
		code+="pwd" + "\n"

		param_values=OrderedDict()
		param_values["-ref"]=  "../" + self.anchor.getOutputFile()
		param_values["-i"]=  "../" + self.query.getOutputFile()
		param_values["-o"]=  self.output_prefix
		param_values["-maxthreads"]=  str(self.getThreads())
		param_values["-insertThreads"]=  str(self.getThreads())
		maxmem=self.getMem() / self.getThreads()
		if maxmem < 1:
			maxmem=1
		param_values["-maxmem"]=  str(maxmem)
		param_values["-output-veto-filter"]=  self.output_veto_regex
		param_values["-res"]=  str(self.res)
		param_values["-T"]=  str(self.pval)
		param_values["-FP"]=  str(self.fp)
		param_values["-FN"]=  str(self.fn)
		param_values["-sf"]=  str(self.sf)
		param_values["-sd"]=  str(self.sd)
		param_values["-extend"]=  "1" if self.allow_overhang else "0"
		param_values["-outlier"]=  str(self.outlier_pval)
		param_values["-endoutlier"]=  str(self.end_outlier_pval)
		param_values["-deltaX"]=  str(self.max_query_alignment_interval)
		param_values["-deltaY"]=  str(self.max_reference_alignment_interval)
		param_values["-xmapchim"]=  str(self.min_sites_for_chimera)
		param_values["-hashgen"]=  " ".join([str(self.hash_window), str(self.hash_min_sites), str(self.hash_sd_max), str(self.hash_sd_rms), str(self.hash_relative_error), str(self.hash_offset_kb), str(self.hash_max_insert_errors), str(self.hash_max_probe_errors), str(self.hash_max_unresolved_sites)])
		param_values["-hash"]=  ""
		param_values["-hashdelta"]=  str(self.hash_delta)
		param_values["-mres"]=  str(self.target_resolution)
		param_values["-rres"]=  str(self.resolution_reduction)
		param_values["-nosplit"]=  "2" if self.allow_no_splits else "0" if self.allow_infinite_splits else "1"
		param_values["-biaswt"]=  str(self.scale_bias_wt)

		if self.send_output_to_file:
			param_values["-stdout"]=""
		if self.send_error_to_file:
			param_values["-stderr"]=""
		if self.overwrite_output:
			param_values["-force"]=""
		if self.print_indel_file:
			param_values["-indel"]=""

		param_list=[self.workspace.binaries["bng_ref_aligner"]]
		for key in param_values:
			param_list.append(key)
			param_list.append(param_values[key])
		code+=" ".join(param_list)
		return [code]

	def getStepDir(self):
#		return self.step_dir
		return "_".join(["comparison", self.workspace.input_file, str(self.merge.assembly.vital_parameters.fp), str(self.merge.assembly.vital_parameters.fn), str(self.merge.assembly.vital_parameters.pval), str(self.merge.assembly.vital_parameters.min_molecule_len), str(self.merge.assembly.vital_parameters.min_molecule_sites), self.ref_file])

	def getOutputFile(self):
		return self.getStepDir() + "/" + self.output_prefix + "." + self.getOutputFileExtension()

	def getOutputFileExtension(self):
		return "xmap"

	def autoGeneratePrereqs(self):
		work_dir=self.workspace.work_dir
		self.anchor=Input(Workspace(work_dir, self.ref_file))
		self.anchor.prereq=self.merge
		self.query=Input(Workspace(work_dir, self.merge.getOutputFile()))
		self.query.prereq=self.anchor
		self.query.step_dir="input_for_"+self.getStepDir()

	def getPrereq(self):
		return self.query

	def loadQualityReportItems(self):
		if self.quality is None:
			self.loadQualityObjectFromFile()

		report_items=OrderedDict()
		report_items["Num alignments: " + str(self.quality.num_alignments)]=1
		num_query_contigs=self.query.loadQuality_count()
		report_items["Num query contigs total: " + str(num_query_contigs)]=3
		report_items["Num query contigs that don't align: " + str(num_query_contigs-self.quality.aligned_query_contig_num)]=2 
		self.quality.unaligned_query_contig_num=num_query_contigs-self.quality.aligned_query_contig_num
		report_items["Total length: " + str(self.quality.total_length)]=2
		report_items["Query length: " + str(self.quality.total_query_length)]=3
		report_items["Proportion of query with match: "  + str(self.quality.proportion_query)]=1
		report_items["Average proportion of query contig within match: " + str(self.quality.average_proportion_of_query_matching)]=2
		report_items["Anchor length: " + str(self.quality.total_anchor_length)]=3
		report_items["Proportion of anchor with match: " + str(self.quality.proportion_anchor)]=1

		report_items["Total confidence: " + str(self.quality.total_confidence)]=2
		report_items["Max confidence: " + str(self.quality.max_confidence)]=2
		report_items["Min confidence: " + str(self.quality.min_confidence)]=2
		report_items["Average confidence: " + str(self.quality.average_confidence)]=3
		report_items["Weighted average confidence: " + str(self.quality.weighted_average_confidence)]=1
		return report_items

	def createQualityObject(self):
		proportion_of_query_matching=0.0
		total_query_length=0.0
		total_anchor_length=0.0
		num_alignments=0
		aligned_query_contigs=set()
		total_confidence=0.0
		confidence_times_length=0.0
		total_length=0.0
		max_confidence=0.0
		min_confidence=maxint
		
		for alignment in XmapFile(self.getOutputFile()).parse():
			num_alignments+=1
			aligned_query_contigs.add(alignment.query_id)

			conf=alignment.confidence 
			total_confidence+=conf

			if conf > max_confidence:
				max_confidence=conf
			if conf < min_confidence:
				min_confidence=conf

			query_length=abs(alignment.query_end-alignment.query_start)
			total_query_length+=query_length
			proportion_of_query_matching+=query_length/alignment.query_len
			anchor_length=abs(alignment.anchor_end-alignment.anchor_start)
			total_anchor_length+=anchor_length
#			length=max(query_length, anchor_length)
			length=query_length
			total_length+=length
			confidence_times_length+=conf*length

		weighted_average_confidence=confidence_times_length/total_length
		proportion_query=query_length / self.query.loadQuality_length()
		proportion_anchor=total_anchor_length / self.anchor.loadQuality_length()
		average_confidence=total_confidence / num_alignments

		self.quality=Quality(
			num_alignments=num_alignments,
			aligned_query_contig_num=len(aligned_query_contigs),
			total_length=total_length,
			total_query_length=total_query_length,
			proportion_query=proportion_query,
			average_proportion_of_query_matching=(proportion_of_query_matching/num_alignments),
			total_anchor_length=total_anchor_length,
			proportion_anchor=proportion_anchor,

			total_confidence=total_confidence,
			min_confidence=min_confidence,
			max_confidence=max_confidence,
			average_confidence=average_confidence,
			weighted_average_confidence=weighted_average_confidence
		)
		self.saveQualityObjectToFile()

	def getQuality_weightedAverageConfidence(self):
		if self.quality is None:
			self.loadQualityObjectFromFile()
		return self.quality.weighted_average_confidence

	def getQuality_totalLength(self):
		if self.quality is None:
			self.loadQualityObjectFromFile()
		return self.quality.total_length

	def getMem(self):
		return self.workspace.resources.getSmallMemory()
	def getTime(self):
		return self.workspace.resources.getSmallTime()
	def getThreads(self):
		return self.workspace.resources.getSmallThreads() class Merge(Step):
	def __init__(self, workspace, assembly):
		if isinstance(assembly,VitalParameters):
			assembly=Assembly(workspace, assembly)
		self.workspace=workspace
		self.assembly=assembly
		self.quality=assembly.quality

		self.output_prefix="merge_of_"+assembly.output_prefix
		self.overwrite_output=True
		self.write_output_to_file=True
		self.write_error_to_file=True

		self.autoGeneratePrereqs()

	def __eq__(self, other):
		if other is None:
			return False
		if self.__class__ != other.__class__:
			return False
		return self.assembly==other.assembly

	def __ne__(self, other):
		return not self==other

	def __hash__(self):
		return hash((self.workspace.input_file, self.workspace.work_dir, self.assembly.vital_parameters.pval, self.assembly.vital_parameters.fp, self.assembly.vital_parameters.fn, self.assembly.vital_parameters.min_molecule_len, self.assembly.vital_parameters.min_molecule_sites, self.__class__.__name__))

	def __str__(self):
		return "merge of " + str(self.assembly)
		
	def writeCode(self):
		code="cd " + self.workspace.work_dir + "\n"
		code+="mkdir -p " + self.getStepDir() + "\n"
		code+="cd " + self.getStepDir() + "\n"
		code+="pwd\n"

		param_values=OrderedDict()
		param_values["-if"]="../" + self.assembly_summary.getOutputFile()
		param_values["-merge"]=""
		param_values["-o"]=self.output_prefix
		if self.overwrite_output:
			param_values["-f"]=""
		if self.write_output_to_file:
			param_values["-stdout"]=""
		if self.write_error_to_file:
			param_values["-stderr"]=""
		
		param_list=[self.workspace.binaries["bng_ref_aligner"]]
		for key in param_values:
			param_list.append(key)
			param_list.append(param_values[key])

		code+=" ".join(param_list) + "\n"

		code += "result=`tail -n 1 ../" + self.getStepDir()  + "/" + self.output_prefix + ".stdout`\n"
		code += "if [[ \"$result\" != \"END of output\" ]]; then exit 1; else touch Complete.status; fi\n"

		return [code]

	def getStepDir(self):
		return "merged_" + self.assembly.getStepDir()
	def getOutputFile(self):
		return self.getStepDir() + "/" + self.output_prefix + "." + self.getOutputFileExtension()
	def getOutputFileExtension(self):
		return "cmap"
	def autoGeneratePrereqs(self):
		self.assembly_summary=Summarize(self.workspace, self.assembly)
	def getPrereq(self):
		return self.assembly_summary

	def getMem(self):
		return self.workspace.resources.getSmallMemory()
	def getTime(self):
		return self.workspace.resources.getSmallTime()
	def getThreads(self):
		return 1 class VitalParameters:
	def __init__(self, fp, fn, pval, min_molecule_len, min_molecule_sites):
		self.fp=fp
		self.fn=fn
		self.pval=pval
		self.min_molecule_len=min_molecule_len
		self.min_molecule_sites=min_molecule_sites
		self.blocks=None
	def __eq__(self, other):
		if other is None:
			return False
		if self.__class__!=other.__class__:
			return False
		return self.__dict__==other.__dict__
	def __ne__(self, other):
		return not self == other class Workspace:
	def __init__(self, work_dir=None, input_file=None):
		self.work_dir=work_dir
		self.input_file=input_file
		self.binaries={}

		self.resources=Resources()

		self.errorNotificationEmail=None

	def __eq__(self, other):
		if other is None:
			return False
		return self.__dict__==other.__dict__

	def __ne__(self, other):
		return not self == other
		
	def addBinary(self, name, path):
		self.binaries[name]=path class CD:
	def __init__(self, newPath):
		self.newPath = os.path.expanduser(newPath)
	def __enter__(self):
		self.savedPath = os.getcwd()
		os.chdir(self.newPath)
	def __exit__(self, etype, value, traceback):
		os.chdir(self.savedPath) class XmapFile(File):
	@staticmethod
	def getExtension():
		return "xmap"

	def parse(self):
		return XmapFile_iter(self.input_file)
	def write(self, alignment, o_file):
		fields=[str(alignment.alignment_id),
			str(alignment.query_id),
			str(alignment.anchor_id),
			str(alignment.query_start),
			str(alignment.query_end),
			str(alignment.anchor_start),
			str(alignment.anchor_end),
			alignment.orientation,
			str(alignment.confidence),
			alignment.hit_enum,
			str(alignment.query_len),
			str(alignment.anchor_len),
			alignment.label_channel,
			alignment.alignment
			]
		o_file.write("\t".join(fields) + "\n") class BedFile(File):
	@staticmethod
	def getExtension():
		return "bed"

	def getHeaders(self):
		headers=[]
		with open(self.input_file, 'r') as i_file:
			for line in i_file:
				fields=line.split("\t")
				if len(fields) < 12:
					headers.append(line)
				else:
					break
		return headers

	def parse(self):
		return BedFile_iter(self.input_file)
	def write(self, feature, o_file):
		fields=[feature.chrom,
			str(feature.start),
			str(feature.stop),
			feature.name,
			str(feature.score),
			feature.strand,
			str(feature.thick_start),
			str(feature.thick_end),
			",".join([str(x) for x in feature.item_rgb]),
			str(feature.block_count),
			",".join([str(x) for x in feature.block_sizes]),
			",".join([str(x) for x in feature.block_starts])]
		o_file.write("\t".join(fields)+"\n")

	def writeDefaultHeaders(self, o_file):
		o_file.write("track name=custom_track description=\"Default description\" useScore=1\n")
		o_file.write("itemRgb=\"On\"\n") class LenFile(File):
	@staticmethod
	def getExtension():
		return "len"

	def parse(self):
		return LenFile_iter(self.input_file)
	def write(self, chromosome, o_file):
		fields=[chromosome.name,
			str(chromosome.length)]
		o_file.write("\t".join(fields)+"\n")

	def writeDefaultHeaders(self, o_file):
		pass class CmapFile(File):
	@staticmethod
	def getExtension():
		return "cmap"
	def parse(self):
		return CmapFile_iter(self.input_file)
	def write(self, label, o_file):
		fields=[str(label.contig_id),
			str(label.contig_len),
			str(label.contig_site_count),
			str(label.label_id),
			label.channel,
			str(label.position),
			str(label.stdev),
			str(label.coverage),
			str(label.occurrences)]
		if hasattr(label, 'snr_mean'):
			fields.extend([str(label.snr_mean),
				str(label.snr_stdev),
				str(label.snr_count)])
		o_file.write("\t".join(fields) + "\n")
	def writeDefaultHeaders(self, o_file):
		o_file.write("""# CMAP File Version:    0.1
# Label Channels:       1
# Nickase Recognition Site 1:   cctcagc
# Enzyme1:      Nt.BbvCI
#h CMapId       ContigLength    NumSites        SiteID  LabelChannel    Positio
#f int  float   int     int     int     float   float   int     int\n""")
	def loadContigs(self):
		contigs={}
		final_labels={}
		for label in self.parse():
			if not label.contig_id in contigs:
				contigs[label.contig_id]=Contig()
			contigs[label.contig_id].addLabel(label, sort=False)

			final_labels[label.contig_id]=label

		for contig_id in final_labels.keys():
			contigs[contig_id].addLabel(final_labels[contig_id], sort=True)

		return contigs class FileConverter:
	def __init__(self, input_file, output_file=None):
		self.input_file=input_file
		self.output_file=output_file

		self.digestor=None

	def convert(self):
		if self.output_file is None:
			raise Exception("Can't convert without target file format info")

		# This is poorly designed
		if self.input_file.getExtension()=="fa":
			if self.output_file.getExtension()=="cmap":
				self.convert_from_fasta_to_cmap()
		elif self.input_file.getExtension()=="xmap":
			if self.output_file.getExtension()=="bed":
				self.convert_from_xmap_to_bed()
			elif self.output_file.getExtension()=="sam":
				self.convert_from_xmap_to_sam()
		elif self.input_file.getExtension()=="cmap":
			if self.output_file.getExtension()=="len":
				self.convert_from_cmap_to_len()

	def convert_from_fasta_to_cmap(self):
		if self.digestor is None:
			raise Exception("Can't convert from fasta to cmap without a motif")
		o_file=open(self.output_file.input_file, 'w')
		self.output_file.writeDefaultHeaders(o_file)
		for i, record in enumerate(self.input_file.parse()):
			try:
				contig_id=int(re.sub('[\D]+', '', record.name))
			except:
				contig_id=i
			length=len(record.seq)
			num_ns=len(re.split('N', str(record.seq)))
			prop_ns='{0:.2f}'.format(float(num_ns)/length)
			if length > 2030:
				if float(prop_ns) < .1:
					labels=self.digestor.digest(record, contig_id)
					for label in labels:
						self.output_file.write(label, o_file)

	def convert_from_xmap_to_bed(self):
		with open(self.output_file.input_file, 'w') as o_file:
			self.output_file.writeDefaultHeaders(o_file)
			for i, alignment in enumerate(self.input_file.parse()):
				chrom="chr"+str(alignment.anchor_id)
				if alignment.orientation=="+":
					start=alignment.anchor_start-alignment.query_start
					stop=alignment.anchor_end+(alignment.query_len-alignment.query_end)
				elif alignment.orientation=="-":
					start=alignment.anchor_start-(alignment.query_len-alignment.query_end)
					stop=alignment.anchor_end+alignment.query_start
				if start < 0:
					start=0
				if stop > alignment.anchor_len:
					stop=alignment.anchor_len
				name="contig"+str(alignment.query_id)
				feature=Feature(chrom, int(start), int(stop), name)
				feature.score=int(alignment.confidence*10)
				feature.strand=alignment.orientation
				feature.thick_start=int(alignment.anchor_start)
				feature.thick_end=int(alignment.anchor_end)
				feature.item_rgb=[0,0,0]
				feature.block_count=0
				feature.block_sizes=[0]
				feature.block_starts=[0]
				
				self.output_file.write(feature,o_file)

	def convert_from_xmap_to_sam(self):
		anchors={}
		confidence_maxima={}
		for align in self.input_file.parse():
			anchors[align.anchor_id]=align.anchor_len

			if not align.query_id in confidence_maxima:
				confidence_maxima[align.query_id]=align.confidence
			elif align.confidence > confidence_maxima[align.query_id]:
				confidence_maxima[align.query_id]=align.confidence

		with open(self.output_file.input_file, 'w') as o_file:
			o_file.write("@HD\tVN:1.0\tSO:unsorted\n")
			for anchor in sorted(anchors.keys()):
				o_file.write("@SQ\tSN:"+str(anchor)+"\tLN:"+str(int(anchors[anchor]))+"\n")

			for align in self.input_file.parse():
				flag=0
				if align.orientation=="-":
					flag+=16
				if align.confidence!=confidence_maxima[align.query_id]:
					flag+=256
				start_pos=align.anchor_start
				if start_pos==0:
					start_pos=1
				output=[
					str(align.query_id),
					str(flag),
					str(align.anchor_id),
					str(int(start_pos)),
					str(int(align.confidence)),
					align.hit_enum.replace('M','='),
					"*","0",
					str(int(align.query_len)),
					"*","*"
				]
				o_file.write("\t".join(output)+"\n")

	def convert_from_cmap_to_len(self):
		with open(self.output_file.input_file, 'w') as o_file:
			ids=set()
			for i, label in enumerate(self.input_file.parse()):
				if label.contig_id in ids:
					continue
				ids.add(label.contig_id)
				name="chr"+str(label.contig_id)
				length=int(label.contig_len)
				chromosome=Chromosome(name, length)
				self.output_file.write(chromosome,o_file)
from Operations.BioNano.Compare.ReferenceAlignment import ReferenceAlignment
from Operations.BioNano.Assemble.Merge import Merge
from Operations.BioNano.Assemble.VitalParameters import VitalParameters
from Utils.Workspace import Workspace
from Utils.CD import CD
from Operations.BioNano.files import XmapFile
from Operations.BioNano.FileConverter import BedFile
from Operations.BioNano.FileConverter import LenFile
from Operations.BioNano.files import CmapFile
from Operations.BioNano.FileConverter import FileConverter
#!/usr/bin/python


work_dir="/path/to/work/dir" ### SET ME
input_file="input.bnx" ### SET ME
workspace=Workspace(work_dir, input_file)

ref_file="reference.cmap" ### SET ME

false_positive=0.5 ### SET ME
false_negative=0.15 ### SET ME
genome_size_mb=900 ### SET ME
p_val=1e-05/genome_size_mb
min_len=100 ### SET ME
min_sites=6 ### SET ME

vital_parameters=VitalParameters(false_positive, false_negative, p_val, min_len, min_sites)
merge=Merge(workspace, vital_parameters)

alignment=ReferenceAlignment(workspace, merge, ref_file)

with CD(work_dir):
	len_file_path=alignment.anchor.getStepDir()+"/"+alignment.output_prefix+".len"
	with open(len_file_path,'w'):
	len_file_converter=FileConverter(CmapFile(alignment.anchor.getOutputFile()), LenFile(len_file_path))
	len_file_converter.convert()

        bed_file_path=alignment.getStepDir()+"/"+alignment.output_prefix+".bed"
	with open(bed_file_path,'w'):
		pass

	bed_file_converter=FileConverter(XmapFile(alignment.getOutputFile()), BedFile(bed_file_path))