#-------------------------------------------------------------------------------
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
# 
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
# USA
# 
# Contact Info:
# 	Bruce Donald
# 	Duke University
# 	Department of Computer Science
# 	Levine Science Research Center (LSRC)
# 	Durham
# 	NC 27708-0129 
# 	USA
# 	brd@cs.duke.edu
# 
# Copyright (C) 2011 Jeffrey W. Martin and Bruce R. Donald
# 
# <signature of Bruce Donald>, April 2011
# Bruce Donald, Professor of Computer Science
#-------------------------------------------------------------------------------


import share

share.getJvm().start()


def writeCoords( file, addresses, protein ):
	
	file.write( "%i\t" % addresses.size() )
	for address in addresses:
		
		# always use the first subunit for the atom position
		address = share.f.protein.AtomAddressInternal( 0, address.getResidueId(), address.getAtomId() )
		
		pos = protein.getAtom( address ).getPosition()
		file.write( "%0.3f\t%0.3f\t%0.3f\t" % ( pos.x, pos.y, pos.z ) )


def generate( pathOut, pathPdb, pathDistanceRestraints, printComment, padPercent ):

	# load the protein and noes
	protein = share.loadProtein( pathPdb )
	restraintsReadable = share.loadDistanceRestraints( pathDistanceRestraints )
	share.interpretDistanceRestraintsWithPseudoatoms( restraintsReadable, protein )
	restraints = share.mapDistanceRestraintsToProtein( restraintsReadable, protein )
	share.reportMismatchedRestraints( restraints, protein.getSubunits().size() )
	restraints = share.filterRestraints( restraints )
	share.printRestraintsStats( restraints )
	
	# pad the distances if needed
	if padPercent > 0:
		share.padRestraints( restraints, padPercent )
	
	# sort the restraints by distance
	#restraints = sorted( restraints, key=lambda restraint: restraint.getMaxDistance() )
	
	# write out the restraints
	file = open( pathOut, "w" )
	for restraint in restraints:
		
		# render the comment if needed before we map to a subunit
		comment = ""
		if printComment:
			restraintReadable = share.f.nmr.DistanceRestraintMapper.mapInternalToReadable( restraint, protein )
			comment = share.f.nmr.DistanceRestraintWriter().writeToString( restraintReadable )
		
		share.f.protein.AtomAddressInternal.mapAddressesToSubunit( restraint.getLefts(), 0 )
		share.f.protein.AtomAddressInternal.mapAddressesToSubunit( restraint.getRights(), 0 )
		
		writeCoords( file, restraint.getLefts(), protein )
		writeCoords( file, restraint.getRights(), protein )
		file.write( "%0.3f\t%0.3f" % ( restraint.getMaxDistance(), restraint.getMinDistance() ) );
		if printComment:
			file.write( "\t#%s" % comment )
		file.write( "\n" )
	
	print "Restraints written to", pathOut


if __name__ == "__main__":
	
	# set up paths
	dirInput = "/home/jeff/duke/donaldLab/dshot/input"
	dirOut = "/home/jeff/workspace/jdshot/output"
	dirResults = "/home/jeff/duke/donaldLab/noe+rdc/results"
	
	# DAGK
#	dirResults += "/DAGK/annuli/padding=0%"
#	generate(
#		dirResults + "/2KDC-restraints.txt",
#		dirInput + "/2KDC.oligomer.protein",
#		dirInput + "/2KDC.experimental.noe",
#		#dirResults + "/consistent.noe",
#		True,
#		0.0
#	)
	
	# DAGK test
	dirResults += "/DAGK/temp/allRestraints"
	generate(
		dirResults + "/2KDC-restraints.txt",
		dirInput + "/2KDC.oligomer.protein",
		dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe",
		True,
		0.0
	)

#	# GB1
#	dirResultsGB1 = dirResults + "/GB1/normal"
#	generate(
#		dirResultsGB1 + "/1Q10-restraints.txt",
#		dirInput + "/1Q10.oligomer.protein",
#		dirResultsGB1 + "/1Q10.original.noe",
#		True,
#		0.0
#	)
#
#	# GB1
#	dirResultsGB1 = dirResults + "/GB1/reassigned"
#	generate(
#		dirResultsGB1 + "/1Q10-restraints.txt",
#		dirInput + "/1Q10.oligomer.protein",
#		dirResultsGB1 + "/1Q10.relaxed.noe",
#		True,
#		0.0
#	)
	
	# BM2
	#generate(
	#	dirOut + "/2KJ1-restraints.txt",
	#	dirData + "/2KJ1.oligomer.protein",
	#	dirData + "/2KJ1.experimental.noe",
	#	False
	#)
