#-------------------------------------------------------------------------------
# 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, util, math, jpype
from Struct import Struct
from figure import Figure

share.getJvm().start( "-Xmx2g" ) #, debugAddress=8787 )


SubunitOrder = share.f.protein.SubunitOrder


def run( pathOutRmsdsPlot, pathOutKinemage, pathInReference, subunitOrder, pathInEnsemble, pathInDistanceRestraints ):

	# read inputs
	ensemble = share.loadProteins( pathInEnsemble )
	reference = share.loadProtein( pathInReference )
	
	# get backbones
	referenceBackbone = reference.getBackbone()
	ensembleBackbones = [protein.getBackbone() for protein in ensemble]
	
	# convert the subunit order to match the reference structure
	for ensembleBackbone in ensembleBackbones:
		subunitOrder.convertComputedToReference( ensembleBackbone )
	
	# load the distance restraints
	restraintsReadable = share.loadDistanceRestraints( pathInDistanceRestraints )
	share.interpretDistanceRestraintsWithPseudoatoms( restraintsReadable, reference )
	restraints = share.mapDistanceRestraintsToProtein( restraintsReadable, reference )
	share.printRestraintsStats( restraints )
	restraints = share.filterRestraintsIntersubunit( restraints )
	restraintsReadable = share.unmapDistanceRestraintsFromProtein( restraints, reference )
	
	# perform all the backbone alignments
	share.alignEnsembleOptimally( referenceBackbone, ensembleBackbones )
	
	# do the analysis
	analyzeRmsd( pathOutRmsdsPlot, referenceBackbone, ensembleBackbones )
	analyzeRestraints( reference, ensemble, restraintsReadable )
	
	share.saveEnsembleKinemage( pathOutKinemage, ensembleBackbones, referenceBackbone )


def analyzeRmsd( pathOut, referenceBackbone, ensembleBackbones ):
	
	exactBounds = False
	
	rmsds = [share.getRmsd( referenceBackbone, backbone ) for backbone in ensembleBackbones]
	
	# print out individual RMSDs
	for rmsd in rmsds:
		print "Backbone RMSD: %f" % rmsd
	print "Min: %f" % min( rmsds )
	print "Max: %f" % max( rmsds )
	
	# render the plot
	fig = Figure()
	( _, bins, _ ) = fig.plot.hist( rmsds, bins=10, facecolor="#004586" )
	fig.plot.set_xlabel( u"RMSD to reference structure in \u212B" )
	fig.plot.set_ylabel( "Number of structures" )
	if exactBounds:
		fig.plot.set_xticks( bins )
		fig.plot.set_xticklabels( ["%.2f" % s for s in bins], rotation=-45, horizontalalignment="left" )
	else:
		fig.plot.set_xbound( lower=0, upper=math.ceil( max( rmsds ) ) )
	
	fig.save( pathOut )
	print "Write RMSDs figure to:\n\t%s" % pathOut


def analyzeRestraintsForStructure( name, structure, restraintsReadable ):
	
	util.setSilence( True )
	restraints = share.mapDistanceRestraintsToProtein( restraintsReadable, structure )
	util.setSilence( False );
	
	rmsd = share.getDistanceRestraintRmsd( structure, restraints )
	numSatisfied = share.getDistanceRestraintNumSatisfied( structure, restraints )
	numRestraints = len( restraints )
	
	print "%s NOE RMSD: %f\tNum Satisfied: %d/%d %.2f" % (
		name,
		rmsd,
		numSatisfied,
		numRestraints,
		numSatisfied * 100.0 / numRestraints
	)
	
	
def analyzeRestraints( referenceStructure, computedStructures, restraintsReadable ):
	
	analyzeRestraintsForStructure( "Reference", referenceStructure, restraintsReadable )
	for computedStructure in computedStructures:
		analyzeRestraintsForStructure( "Computed", computedStructure, restraintsReadable )


try:

	dirResults = "/home/jeff/duke/donaldLab/workspace/jdshot/output"
	dirInput = "/home/jeff/duke/donaldLab/dshot/input"
	
	# GB1
#	dirResultsGb1 = dirResults + "/gb1"
#	run(
#		dirResultsGb1 + "/backboneRmsds.png",
#		dirResultsGb1 + "/ensemble.kin",
#		dirInput + "/1Q10.oligomer.protein",
#		SubunitOrder( "AB" ),
#		dirResultsGb1 + "/structures.pdb",
#		dirInput + "/1Q10.experimental.fixed.noe"
#	)
	
	# GB1 reassigned
	dirResultsGb1 = dirResults + "/gb1-reassigned"
	run(
		dirResultsGb1 + "/backboneRmsds.png",
		dirResultsGb1 + "/ensemble.kin",
		dirInput + "/1Q10.oligomer.protein",
		SubunitOrder( "AB" ),
		dirResultsGb1 + "/structures.pdb",
		dirInput + "/1Q10.experimental.fixed.noe"
	)
	
	# DAGK
#	dirResultsDagk = dirResults + "/dagk"
#	run(
#		dirResultsDagk + "/backboneRmsds.png",
#		dirResultsDagk + "/ensemble.kin",
#		dirInput + "/2KDC.oligomer.protein",
#		SubunitOrder( "ACB" ),
#		dirResultsDagk + "/structures.pdb",
#		dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe"
#	)

except jpype.JavaException as ex:
	import sys
	print >> sys.stderr, "Java Exception: %s\n" % ex
	print >> sys.stderr, ex.stacktrace()
	raise
