#-------------------------------------------------------------------------------
# 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 jvm, jdshot, share
import jpype
import os.path

jdshot.getJvm().start()

import disco

Math = jvm.f.java.lang.Math
AlignmentTensorAxis = share.f.nmr.AlignmentTensorAxis
Plotter = jdshot.f.disco.Plotter
OrientationCalculator = jdshot.f.disco.OrientationCalculator
Vector3 = share.f.geom.Vector3


def run( outDir, pathInSubunit, numSubunits, pathInRdcs, axis, pathInReference, pathInOrientations=None ):
	
	# load the inputs
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	rdcsContext = jdshot.newRdcsContext( symmetryContext, pathInRdcs )
	reference = share.loadProtein( pathInReference )
	
	# compute the symmetry axis orientation
	tensor = share.newAlignmentTensor( symmetryContext.getMonomer(), rdcsContext.getInternalRdcs() )
	computedOrientation = tensor.getAxis( axis )
	
	# compute the symmetry axis of the reference oligomer
	referenceOrientation = jdshot.f.disco.ReferenceAnalyzer.getSymmetryAxisOrientation( reference )
	
	# find the flip of the reference axis that minimizes the angle to the best tensor axis
	negatedReferenceOrientation = Vector3( referenceOrientation )
	negatedReferenceOrientation.negate()
	angleRad = Math.acos( computedOrientation.getDot( referenceOrientation ) / computedOrientation.getLength() / referenceOrientation.getLength() )
	angleRadNegated = Math.acos( computedOrientation.getDot( negatedReferenceOrientation ) / computedOrientation.getLength() / negatedReferenceOrientation.getLength() )
	if angleRadNegated < angleRad:
		referenceOrientation = negatedReferenceOrientation
		angleRad = angleRadNegated
	print "Angle to reference: %f" % Math.toDegrees( angleRad )
	
	# sample orientations
	tensors = share.newAlignmentTensors( symmetryContext.getMonomer(), rdcsContext.getInternalRdcs(), 10000 )
	#tensors = OrientationCalculator.pickTensorCluster( tensors, tensor.getAxis( axis ), axis )
	
	# load the sampled orientations
	sampledOrientations = jdshot.f.disco.ResultsReader().readOrientations( jvm.newFile( pathInOrientations ) )
	
	# finally, plot it
	disco.writeChart(
		jvm.newFile( outDir, "sampledTensorAxesWithReference.png" ),
		Plotter.plotAlignmentTensorsAxisWithReferenceAndSampling( tensors, tensor, axis, referenceOrientation, sampledOrientations )
	)
	
	print "Done!"


dirInput = "/home/jeff/duke/donaldLab/dshot/input"

try:

	# DAGK
	run(
		"output/dagk",
		dirInput + "/2KDC.monomer.protein",
		3,
		dirInput + "/2KDC.experimental.rdc",
		AlignmentTensorAxis.Z,
		dirInput + "/2KDC.oligomer.protein",
		"output/dagk/sampledOrientations.dat"
	)

	# GB1
#	run(
#		"output/gb1",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		dirInput + "/1Q10.experimental.rdc",
#		AlignmentTensorAxis.X,
#		dirInput + "/1Q10.oligomer.protein",
#		"output/gb1/orientations.dat"
#	)
	
	# M2
#	run(
#		dirInput + "/2RLF.monomer.protein",
#		4,
#		dirInput + "/2RLF.experimental.rdc",
#		dirInput + "/2RLF.oligomer.protein",
#	)

	# BM2
#	run(
#		dirInput + "/2KJ1.monomer.protein",
#		4,
#		dirInput + "/2KJ1.experimental.rdc",
#		dirInput + "/2KJ1.oligomer.protein",
#	)
	
	# PLB
#	run(
#		dirInput + "/1ZLL.monomer.protein",
#		4,
#		dirInput + "/1ZLL.experimental.rdc",
#		dirInput + "/1ZLL.oligomer.protein",
#	)

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