#-------------------------------------------------------------------------------
# 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, jpype, math, numpy
import jdshot, share
import disco

dirInput = "/home/jeff/duke/donaldLab/dshot/input"
dirWork = "/home/jeff/duke/donaldLab/workspace/jdshot/work"

Subunit = share.f.protein.Subunit
Plotter = jdshot.f.disco.Plotter
OrientationCalculator = jdshot.f.disco.OrientationCalculator
AnnulusCalculator = jdshot.f.disco.AnnulusCalculator
PositionCalculator = jdshot.f.disco.PositionCalculator
AlignmentTensorAxis = share.f.nmr.AlignmentTensorAxis
ResultsWriter = jdshot.f.disco.ResultsWriter
ResultsReader = jdshot.f.disco.ResultsReader
SubunitVariationAnalyzer = jdshot.f.disco.SubunitVariationAnalyzer
PseudoatomBuilder = share.f.pseudoatoms.PseudoatomBuilder
SubunitOrder = share.f.protein.SubunitOrder
Vector2 = share.f.geom.Vector2
ProteinGeometry = share.f.protein.tools.ProteinGeometry
StructureAligner = share.f.analysis.StructureAligner
ReferenceAnalyzer = jdshot.f.disco.ReferenceAnalyzer
FaceInfo = jdshot.f.disco.FaceInfo


def computeSubunitVariation( outDir, pathInEnsemble, subunitId, pathInDistanceRestraints ):
	
	# read inputs and build pseudoatoms
	proteins = share.loadProteins( pathInEnsemble )
	for protein in proteins:
		PseudoatomBuilder.getInstance().build( protein )
	
	mainProtein = proteins.get( 0 )
	distanceRestraintsContext = jdshot.newDistanceRestraintsContext( mainProtein, pathInDistanceRestraints )
	
	# get just the subunits
	mainSubunit = mainProtein.getSubunit( subunitId )
	subunits = jvm.toArrayList( [] )
	for protein in proteins:
		subunits.add( protein.getSubunit( subunitId ) )
	
	# show the ensemble
	share.saveEnsembleKinemage( jvm.newFile( outDir, "subunitEnsemble.kin" ), subunits, mainSubunit )
	
	variation  = SubunitVariationAnalyzer.getVariation( subunits, distanceRestraintsContext )
	ResultsWriter().writeVariation( jvm.newFile( outDir, "variation.dat" ), variation )


def computeSingleOrientation( outDir, pathInSubunit, numSubunits, pathInRdcs, axis ):

	# load the inputs
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	rdcsContext = jdshot.newRdcsContext( symmetryContext, pathInRdcs )
	
	# sample orientations
	tensor = share.newAlignmentTensor( symmetryContext.getMonomer(), rdcsContext.getInternalRdcs() )
	orientations = jvm.toArrayList( [ tensor.getAxis( axis ) ] )
	disco.ResultsWriter().writeOrientations( jvm.newFile( outDir, "orientation.dat" ), orientations )


def computeSampledOrientations( outDir, pathInSubunit, numSubunits, pathInRdcs, axis, percentile, resolution ):
	
	# load the inputs
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	rdcsContext = jdshot.newRdcsContext( symmetryContext, pathInRdcs )
	
	# sample orientations
	orientations = disco.getSampledOrientations( outDir, symmetryContext, rdcsContext, axis, percentile, resolution )
	disco.ResultsWriter().writeOrientations( jvm.newFile( outDir, "sampledOrientations.dat" ), orientations )

	
def computeMsrs( outDir, pathInSubunit, numSubunits, pathInDistanceRestraints, pathInOrientations, orientationIndex, pathInVariation=None, padPercent=None, tag=None ):
	
	# load the inputs
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	distanceRestraintsContext = jdshot.newDistanceRestraintsContext( symmetryContext, pathInDistanceRestraints )
	
	# apply variation if needed
	if pathInVariation != None:
		variation = ResultsReader().readVariation( jvm.newFile( pathInVariation ) )
		distanceRestraintsContext.setInternalRestraints(
			SubunitVariationAnalyzer.applyVariation( variation, distanceRestraintsContext.getInternalRestraints() ),
			symmetryContext.getClonedMonomer()
		)
		print "Applied variation!"
		
	# apply beta if needed
	if padPercent != None:
		share.padRestraints( distanceRestraintsContext.getInternalRestraints(), padPercent / 100.0 )
	
	# read the orientation
	orientations = disco.ResultsReader().readOrientations( jvm.newFile( pathInOrientations ) )
	orientation = orientations.get( orientationIndex )
	print "Analyzing orientation %d of %d..." % ( orientationIndex + 1, orientations.size() )
	
	# normalize a copy of the subunit
	monomer = Subunit( symmetryContext.getPseudoMonomer() )
	OrientationCalculator.normalizeProtein( monomer, orientation )
	
	# compute the annuli
	annulusIndex = AnnulusCalculator.computeAnnuli(
		monomer,
		distanceRestraintsContext.getInternalRestraints(),
		symmetryContext.getNumSubunits()
	)
	
	# build the file tag if needed
	if tag is None:
		tag = "%d" % orientationIndex
	
	# compute the arrangement and MSRs
	arrangement = PositionCalculator.computeArrangementIncremental( annulusIndex )
	msrs = PositionCalculator.computeMsrs( arrangement, annulusIndex )
	print "Orientation satisfies %d of %d distance restraints" % (
		PositionCalculator.getNumSatisfiedRestraints( msrs ),
		distanceRestraintsContext.getInternalRestraints().size()
	)
	ResultsWriter().writeMsrs( jvm.newFile( outDir, "msrs.%s.dat" % tag ), msrs )
	
	# print out the MSRs
	disco.writeChart(
		jvm.newFile( outDir, "arrangement.%s.png" % tag ),
		Plotter.plotArrangement( msrs, annulusIndex )
	)
	disco.writeChart(
		jvm.newFile( outDir, "arrangementZoomed.%s.png" % tag ),
		Plotter.plotArrangement( msrs, annulusIndex, True )
	)
	

def computeSampledStructures( outDir, pathInSubunit, numSubunits, pathInDistanceRestraints, pathInOrientations, samplingResolution ):
	
	# load the inputs
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	distanceRestraintsContext = jdshot.newDistanceRestraintsContext( symmetryContext, pathInDistanceRestraints )
	
	# read the orientations
	orientations = disco.ResultsReader().readOrientations( jvm.newFile( pathInOrientations ) )
	
	# for each orientation...
	allComputedStructures = jvm.toArrayList( [] )
	numSampledStructures = []
	numSatisfiedRestraints = []
	for orientationIndex in range( 0, orientations.size() ):
		
		# do we have msrs for this orientation?
		msrsFile = jvm.newFile( outDir, "msrs.%d.dat" % orientationIndex )
		if msrsFile.exists():
			
			print "Sampling msrs %d of %d..." % ( orientationIndex + 1, orientations.size() )
			orientation = orientations.get( orientationIndex )
			
			# normalize a copy of the subunit
			monomer = Subunit( symmetryContext.getPseudoMonomer() )
			OrientationCalculator.normalizeProtein( monomer, orientation )
			
			# compute the annuli
			annulusIndex = AnnulusCalculator.computeAnnuli(
				monomer,
				distanceRestraintsContext.getInternalRestraints(),
				symmetryContext.getNumSubunits()
			)
			
			msrs = ResultsReader().readMsrs( msrsFile )
			
			# sample discrete structures
			computedStructures = disco.getDiscreteStructures(
				outDir,
				symmetryContext,
				monomer,
				annulusIndex,
				msrs,
				None,
				samplingResolution,
				orientationIndex
			)
			print "Sampled %d computed structures from the MSRs" % len( computedStructures )
			
			allComputedStructures.addAll( jvm.toArrayList( computedStructures ) )
			numSampledStructures.append( len( computedStructures ) )
			numSatisfiedRestraints.append( PositionCalculator.getNumSatisfiedRestraints( msrs ) )
			
		else:
			
			numSampledStructures.append( 0 )
			numSatisfiedRestraints.append( 0 )
	
	print "Sampled %d total structures" % allComputedStructures.size()
	for i in range( 0, orientations.size() ):
		print "Orientation %d:\tnumSatisfied %d\tnumSampled %d" % ( i, numSatisfiedRestraints[i], numSampledStructures[i] )
	
	# write out the computed structures in a PDB file
	share.writeProteins(
		jvm.newFile( outDir, "structures.pdb" ),
		jvm.toArrayList( allComputedStructures )
	)


def computeAlignedBackbones( outDir, pathInReference, pathInEnsemble, subunitOrder ):
	
	# read the reference structure
	referenceBackbone = share.loadProtein( pathInReference ).getBackbone()
	
	# read the ensemble and convert to the reference order
	ensemble = share.loadProteins( pathInEnsemble )
	ensembleBackbones = [protein.getBackbone() for protein in ensemble]
	for protein in ensembleBackbones:
		subunitOrder.convertComputedToReference( protein )
	
	# align the ensemble and write out
	share.alignEnsembleOptimally( referenceBackbone, ensembleBackbones )
	share.writeProteins( jvm.newFile( outDir, "backbonesAlignedToReference.pdb" ), jvm.toArrayList( ensembleBackbones ) )
	share.writeProtein( jvm.newFile( outDir, "centeredReference.pdb" ), referenceBackbone )


def loadMinimizationData( pathIn ):
	return numpy.loadtxt( pathIn, usecols=(4,5), unpack=True )


def analyzeMinimizedEnsemble( dirOut, pathInComputed, pathInReference, dirInProteins ):
	
	# read in data
	( computedEnergies, computedRmsds ) = loadMinimizationData( pathInComputed )
	( referenceEnergies, referenceRmsds ) = loadMinimizationData( pathInReference )
	
	# convert the data into points
	reference = Vector2( referenceEnergies, referenceRmsds )
	computed = jvm.toArrayList( [] )
	for i in range( 0, len( computedEnergies ) ):
		computed.add( Vector2( computedEnergies[i], computedRmsds[i] ) )
	
	# plot structure scores
	disco.writeChart(
		jvm.newFile( dirOut, "structureScores" ),
		Plotter.plotStructureScoresWithReference( computed, reference )
	)
	
	# compute variance
	structures = share.loadProteinsFromDir( dirInProteins )
	print "All-atom average variance %f A^2 or %f A RMSD" % (
		share.getAverageVariance( structures ),
		share.getAverageRmsd( structures )
	)
	backbones = jvm.toArrayList( [protein.getBackbone() for protein in structures] )
	print "Backbone atom average variance %f A^2 or %f A RMSD" % (
		share.getAverageVariance( backbones ),
		share.getAverageRmsd( backbones )
	)


def analyzeAnnuli( outDir, pathInSubunit, numSubunits, pathInDistanceRestraints, pathInOrientations, orientationIndex, pathInVariation=None ):
	
	# load the inputs
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	distanceRestraintsContext = jdshot.newDistanceRestraintsContext( symmetryContext, pathInDistanceRestraints )
	
	# apply variation if needed
	if pathInVariation != None:
		variation = ResultsReader().readVariation( jvm.newFile( pathInVariation ) )
		distanceRestraintsContext.setInternalRestraints(
			SubunitVariationAnalyzer.applyVariation( variation, distanceRestraintsContext.getInternalRestraints() ),
			symmetryContext.getClonedMonomer()
		)
		print "Applied variation!"
	
	# read the orientation
	orientations = disco.ResultsReader().readOrientations( jvm.newFile( pathInOrientations ) )
	orientation = orientations.get( orientationIndex )
	print "Analyzing orientation %d of %d..." % ( orientationIndex + 1, orientations.size() )
	
	# normalize a copy of the subunit
	monomer = Subunit( symmetryContext.getPseudoMonomer() )
	OrientationCalculator.normalizeProtein( monomer, orientation )
	
	# compute the annuli
	annulusIndex = AnnulusCalculator.computeAnnuli(
		monomer,
		distanceRestraintsContext.getInternalRestraints(),
		symmetryContext.getNumSubunits()
	)
	
	# how much noise is there? (assuming one assignment is correct)
	# how many annuli are empty?
	numAnnuli = annulusIndex.getNumAnnuli()
	numEmptyAnnuli = annulusIndex.getNumEmptyAnnuli()
	numNoiseAnnuli = annulusIndex.getNumNoiseAnnuli()
	numRemainingAnnuli = numAnnuli - numEmptyAnnuli
	numNonEmptyNoiseAnnuli = annulusIndex.getNumNonEmptyNoiseAnnuli()
	print "Assignments per restraint: %.1f" % ( 1.0 * numAnnuli / annulusIndex.getNumRestraints() )
	print "Noise: %d/%d (%.1f)%%" % ( numNoiseAnnuli, numAnnuli, ( 100.0 * numNoiseAnnuli / numAnnuli ) )
	print "Empty annuli: %d/%d (%.1f)%%" % ( numEmptyAnnuli, numAnnuli, ( 100.0 * numEmptyAnnuli / numAnnuli ) )
	print "Noise annuli after RDC pruning: %d/%d (%.1f)%%" % ( numNonEmptyNoiseAnnuli, numRemainingAnnuli, ( 100.0 * numNonEmptyNoiseAnnuli / numRemainingAnnuli ) )
	print "Assignments per restraint after RDC pruning: %.1f" % ( 1.0 * numRemainingAnnuli / annulusIndex.getNumRestraints() )


def computeVdwEnergy( pathInMonomer ):
	
	# load the inputs
	print "van der Waals energy: %f" % share.calculateEnergy( share.loadProtein( pathInMonomer ) )


def perturbationStudy( outDir, pathInSubunit, numSubunits, pathInDistanceRestraints, pathInOrientations ):
	
	#padPercents = [f * 0.25 for f in range( 21 ) ]
	#padPercents = [5.25 + f * 0.25 for f in range( 20 ) ]
	padPercents = [10.25 + f * 0.25 for f in range( 20 ) ]
	
	# TEMP
	#print padPercents
	#return

	for padPercent in padPercents:
		computeMsrs(
			outDir,
			pathInSubunit,
			numSubunits,
			pathInDistanceRestraints,
			pathInOrientations,
			0,
			pathInVariation = None,
			padPercent = padPercent,
			tag = "%03d" % ( padPercent * 100 )
		)


def plotMsrsOverlap( outDir, pathInSubunit, numSubunits, pathInReference, pathInOrientation ):

	#padPercents = [ f * 0.25 for f in range( 61 ) ]
	#tags = [ "%03d" % ( padPercent * 100 ) for padPercent in padPercents ]
	# 0,22
	# 22,3
	# 25,4
	# 29,1
	padPercents = [ ( f + 0 ) * 0.25 for f in range( 22 ) ]
	tags = [ "%03d" % ( padPercent * 100 ) for padPercent in padPercents ]
	
	# get the monomer
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	monomer = share.f.protein.Subunit( symmetryContext.getPseudoMonomer() )
	
	# get the reference
	referenceStructure = share.loadProtein( pathInReference )
	referencePosition = ReferenceAnalyzer.getSymmetryAxisPosition( referenceStructure )
	
	# get the orientation
	orientation = disco.ResultsReader().readOrientations( jvm.newFile( pathInOrientation ) ).get( 0 )
	
	# normalize
	OrientationCalculator.normalizeProteinAndAxis(
		monomer,
		referencePosition,
		orientation
	)
	
	referencePosition = Vector2( referencePosition.x, referencePosition.y )
	
	# read the MSRs
	msrsStack = jvm.toArrayList( [] )
	for tag in tags:
		
		print "reading %s MSRs..." % tag
		
		# read the MSRs
		msrsFile = jvm.newFile( outDir, "msrs.%s.dat" % tag )
		msrs = ResultsReader().readMsrs( msrsFile )
		
		msrsStack.add( msrs )
	
	disco.writeChart(
		jvm.newFile( outDir, "msrsOverlap" ),
		Plotter.plotMsrsOverlap( msrsStack, referencePosition )
	)


def plotMsrRanges( outDir, pathInSubunit, numSubunits, pathInDistanceRestraints, pathInReference, pathInOrientation, subunitOrder ):

	#padPercents = [ f * 0.25 for f in range( 61 ) ]
	#tags = [ "%03d" % ( padPercent * 100 ) for padPercent in padPercents ]
	padPercents = [ f * 0.25 for f in range( 10 ) ]
	tags = [ "%03d" % ( padPercent * 100 ) for padPercent in padPercents ]
	
	# TEMP
	#print tags
	#return
	
	# load the inputs
	symmetryContext = jdshot.newSymmetryContext( jdshot.Cn, pathInSubunit, numSubunits )
	distanceRestraintsContext = jdshot.newDistanceRestraintsContext( symmetryContext, pathInDistanceRestraints )
	
	# read the orientations
	orientation = disco.ResultsReader().readOrientations( jvm.newFile( pathInOrientation ) ).get( 0 )
	
	# get the reference structure (just the backbone) and the axis position
	referenceStructure = share.loadProtein( pathInReference )
	referencePosition = ReferenceAnalyzer.getSymmetryAxisPosition( referenceStructure )
	
	# normalize a copy of the subunit
	monomer = Subunit( symmetryContext.getPseudoMonomer() )
	OrientationCalculator.normalizeProteinAndAxis(
		monomer,
		referencePosition,
		orientation
	)
	
	referencePosition = Vector2( referencePosition.x, referencePosition.y )
	
	# prepare the reference backbone for alignment
	referenceStructure = referenceStructure.getBackbone()
	ProteinGeometry.center( referenceStructure )
	
	# TEMP
	sampleResolution = 0.0125 # A
	
	# for each padding...
	closestDists = []
	farthestDists = []
	closestRmsds = []
	farthestRmsds = []
	for tag in tags:
		
		print "analyzing %s MSRs..." % tag
		
		# read the MSRs
		msrsFile = jvm.newFile( outDir, "msrs.%s.dat" % tag )
		msrs = ResultsReader().readMsrs( msrsFile )
		
		# handle the closest
		closestPosition = FaceInfo.getPointClosestTo( msrs, referencePosition )
		if closestPosition is not None:
			closestDists.append( referencePosition.getDistance( closestPosition ) )
#			closestStructure = symmetryContext.getOligomer(
#				jdshot.f.grid.cn.CnGridPoint( closestPosition.x, closestPosition.y, 0.0, 0.0 ),
#				monomer
#			).getBackbone()
#			subunitOrder.convertComputedToReference( closestStructure )
#			StructureAligner.alignOptimally( referenceStructure, closestStructure )
#			closestRmsds.append( share.getRmsd( referenceStructure, closestStructure ) )
		else:
			closestDists.append( 0.0 )
#			closestRmsds.append( 0.0 )
		
		# handle the farthest
		farthestPosition = FaceInfo.getPointFarthestFrom( msrs, referencePosition )
		farthestDists.append( referencePosition.getDistance( farthestPosition ) )
#		farthestStructure = symmetryContext.getOligomer(
#			jdshot.f.grid.cn.CnGridPoint( farthestPosition.x, farthestPosition.y, 0.0, 0.0 ),
#			monomer
#		).getBackbone()
#		subunitOrder.convertComputedToReference( farthestStructure )
#		StructureAligner.alignOptimally( referenceStructure, farthestStructure )
#		farthestRmsds.append( share.getRmsd( referenceStructure, farthestStructure ) )
		
		# TEMP: get rmsds from sampling
		sampledAxisPositions = PositionCalculator.sampleMsrs( msrs, sampleResolution )
		
		minRmsd = float( 1e3000 ) # ie infinity
		maxRmsd = 0.0
		for axisPosition in sampledAxisPositions:
			sampledStructure = symmetryContext.getOligomer(
				jdshot.f.grid.cn.CnGridPoint( axisPosition.x, axisPosition.y, 0.0, 0.0 ),
				monomer
			).getBackbone()
			subunitOrder.convertComputedToReference( sampledStructure )
			StructureAligner.alignOptimally( referenceStructure, sampledStructure )
			rmsd = share.getRmsd( referenceStructure, sampledStructure )
			minRmsd = min( minRmsd, rmsd )
			maxRmsd = max( maxRmsd, rmsd )
		closestRmsds.append( minRmsd )
		farthestRmsds.append( maxRmsd )
		
	# plot it!
#	disco.writeChart(
#		jvm.newFile( outDir, "msrRangesDist" ),
#		disco.Plotter.plotRanges(
#			jvm.toArrayList( padPercents ),
#			jvm.toArrayList( closestDists ),
#			jvm.toArrayList( farthestDists )
#		)
#	)
#	disco.writeChart(
#		jvm.newFile( outDir, "msrRangesRmsd" ),
#		disco.Plotter.plotRanges(
#			jvm.toArrayList( padPercents ),
#			jvm.toArrayList( closestRmsds ),
#			jvm.toArrayList( farthestRmsds )
#		)
#	)
	
	# TEMP: print numbers
	for i in range( 0, len( closestDists ) ):
		print "%.4f\t%.4f\t%.4f\t%.4f" % ( closestDists[i], farthestDists[i], closestRmsds[i], farthestRmsds[i] )

	print "Done!"


try:
	
	# DAGK
#	computeSubunitVariation(
#		"output/dagk",
#		dirInput + "/raw/2KDC.pdb",
#		0,
#		dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe"
#	)

#	computeSampledOrientations(
#		"output/dagk",
#		dirInput + "/2KDC.monomer.protein",
#		3,
#		dirInput + "/2KDC.experimental.rdc",
#		AlignmentTensorAxis.Z,
#		1.0, # percent/100
#		1.0 # degrees
#	)

#	for i in range( 0, ResultsReader().readOrientations( jvm.newFile( "output/dagk/sampledOrientations.dat" ) ).size() ):
#		computeMsrs(
#			"output/dagk",
#			dirInput + "/2KDC.monomer.protein",
#			3,
#			dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe",
#			"output/dagk/sampledOrientations.dat",
#			i,
#			"output/dagk/variation.dat"
#		)

#	computeMsrs(
#		"output/dagk",
#		dirInput + "/2KDC.monomer.protein",
#		3,
#		dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe",
#		"output/dagk/sampledOrientations.dat",
#		0,
#		"output/dagk/variation.dat"
#	)

#	computeSampledStructures(
#		"output/dagk",
#		dirInput + "/2KDC.monomer.protein",
#		3,
#		dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe",
#		"output/dagk/sampledOrientations.dat",
#		0.75
#	)

#	analyzeMinimizedEnsemble(
#		"output/dagk/minimized",
#		"output/dagk/minimized/scores.txt",
#		"output/dagk/minimized/scoresReference.txt",
#		"output/dagk/minimized/structures"
#	)

#	computeAlignedBackbones(
#		"output/dagk",
#		dirInput + "/2KDC.oligomer.protein",
#		"output/dagk/structures.pdb",
#		SubunitOrder( "ACB" )
#	)
	
#	computeVdwEnergy(
#		dirInput + "/2KDC.monomer.protein"
#	)

	
	# GB1
#	computeSubunitVariation(
#		"output/gb1",
#		dirInput + "/raw/1Q10.pdb",
#		0,
#		dirInput + "/1Q10.experimental.fixed.noe"
#	)
	
#	computeSampledOrientations(
#		"output/gb1",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		dirInput + "/1Q10.experimental.rdc",
#		AlignmentTensorAxis.X,
#		1.0, # percent/100
#		0.75, # degrees
#	)
	
#	for i in range( 0, ResultsReader().readOrientations( jvm.newFile( "output/gb1/sampledOrientations.dat" ) ).size() ):
#		computeMsrs(
#			"output/gb1",
#			dirInput + "/1Q10.monomer.protein",
#			2,
#			dirInput + "/1Q10.experimental.fixed.noe",
#			"output/gb1/sampledOrientations.dat",
#			i,
#			"output/gb1/variation.dat"
#		)

#	computeMsrs(
#		"output/gb1",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		dirInput + "/1Q10.experimental.fixed.noe",
#		"output/gb1/sampledOrientations.dat",
#		19,
#		"output/gb1/variation.dat"
#	)

#	computeSampledStructures(
#		"output/gb1",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		dirInput + "/1Q10.experimental.fixed.noe",
#		"output/gb1/sampledOrientations.dat",
#		0.25
#	)
	
#	analyzeMinimizedEnsemble(
#		"output/gb1/minimized",
#		"output/gb1/minimized/scores.txt",
#		"output/gb1/minimized/scoresReference.txt",
#		"output/gb1/minimized/structures"
#	)

#	computeAlignedBackbones(
#		"output/gb1",
#		dirInput + "/1Q10.oligomer.protein",
#		"output/gb1/structures.pdb",
#		SubunitOrder( "AB" )
#	)

#	computeVdwEnergy(
#		dirInput + "/1Q10.monomer.protein"
#	)


	# GB1 reassigned
#	computeSubunitVariation(
#		"output/gb1-reassigned",
#		dirInput + "/raw/1Q10.pdb",
#		0,
#		"output/gb1-reassigned/1Q10.reassigned.noe"
#	)

#	analyzeAnnuli(
#		"output/gb1-reassigned",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		"output/gb1-reassigned/1Q10.reassigned.noe",
#		"output/gb1-reassigned/sampledOrientations.dat",
#		0
#		#"output/gb1-reassigned/variation.dat"
#	)
	
#	computeMsrs(
#		"output/gb1-reassigned",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		"output/gb1-reassigned/1Q10.reassigned.noe",
#		"output/gb1-reassigned/sampledOrientations.dat",
#		0,
#		"output/gb1-reassigned/variation.dat"
#	)
	
#	computeSampledStructures(
#		"output/gb1-reassigned",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		dirInput + "/1Q10.experimental.fixed.noe",
#		"output/gb1-reassigned/sampledOrientations.dat",
#		0.05
#	)


	# DAGK perturbation study

#	computeSingleOrientation(
#		"output/dagk-perturbation",
#		dirInput + "/2KDC.monomer.protein",
#		3,
#		dirInput + "/2KDC.experimental.rdc",
#		AlignmentTensorAxis.Z
#	)

#	perturbationStudy(
#		"output/dagk-perturbation",
#		dirInput + "/2KDC.monomer.protein",
#		3,
#		dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe",
#		"output/dagk-perturbation/orientation.dat"
#	)

#	plotMsrRanges(
#		"output/dagk-perturbation",
#		dirInput + "/2KDC.monomer.protein",
#		3,
#		dirInput + "/2KDC.experimental.onlyDisulfideBonds.noe",
#		dirInput + "/2KDC.oligomer.protein",
#		"output/dagk-perturbation/orientation.dat",
#		SubunitOrder( "ACB" )
#	)
	
#	plotMsrsOverlap(
#		"output/dagk-perturbation",
#		dirInput + "/2KDC.monomer.protein",
#		3,
#		dirInput + "/2KDC.oligomer.protein",
#		"output/dagk-perturbation/orientation.dat"
#	)


	# gb1 no padding
	
#	computeSingleOrientation(
#		"output/gb1-noPadding",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		dirInput + "/1Q10.experimental.rdc",
#		AlignmentTensorAxis.X
#	)
	
#	computeMsrs(
#		"output/gb1-noPadding",
#		dirInput + "/1Q10.monomer.protein",
#		2,
#		dirInput + "/1Q10.experimental.fixed.noe",
#		"output/gb1-noPadding/orientation.dat",
#		0
#	)
	
	
	# RELEASE TESTING
	
	# get the dimer symmetry axis orientation (works flawlessly)
#	disco.computeDimerAlignmentTensorAxis(
#		"output/testRelease",
#		"/home/jeff/duke/donaldLab/dshot/input/raw/1Q10.pdb",
#		2,
#		"/home/jeff/duke/donaldLab/dshot/input/1Q10.experimental.rdc",
#		"/home/jeff/duke/donaldLab/dshot/input/1Q10.experimental.fixed.noe",
#		numMsrsProcesses = 3
#	)

	# Protein Science DAGK (works flawlessly)
#	disco.run(
#		"output/testRelease",
#		"/home/jeff/duke/donaldLab/dshot/input/raw/2KDC.pdb",
#		3,
#		"/home/jeff/duke/donaldLab/dshot/input/2KDC.experimental.rdc",
#		"/home/jeff/duke/donaldLab/dshot/input/2KDC.experimental.onlyDisulfideBonds.noe",
#		2.0, #1.0,
#		2.0, #0.75,
#		numMinimizationProcesses = 6,
#		numMsrsProcesses = 6
#	)

	# Protein Science GB1 normal (works flawlessly)
#	disco.run(
#		"output/testRelease",
#		"/home/jeff/duke/donaldLab/dshot/input/raw/1Q10.pdb",
#		2,
#		"/home/jeff/duke/donaldLab/dshot/input/1Q10.experimental.rdc",
#		"/home/jeff/duke/donaldLab/dshot/input/1Q10.experimental.fixed.noe",
#		2.0, #0.75,
#		2.0, #0.25,
#		alignmentTensorAxis = AlignmentTensorAxis.X,
#		numMinimizationProcesses = 6,
#		numMsrsProcesses = 6
#	)

	# Protein Science GB1 reassigned (works flawlessly, takes about an hour)
#	disco.run(
#		"output/testRelease",
#		"/home/jeff/duke/donaldLab/dshot/input/raw/1Q10.pdb",
#		2,
#		"/home/jeff/duke/donaldLab/dshot/input/1Q10.experimental.rdc",
#		"/home/jeff/duke/donaldLab/workspace/jdshot/output/gb1-reassigned/1Q10.reassigned.noe",
#		0.0,
#		2.0, #0.25,
#		alignmentTensorAxis = AlignmentTensorAxis.X,
#		numMinimizationProcesses = 6,
#		numMsrsProcesses = 6
#	)
	
	# RECOMB DAGK (works flawlessly)
#	disco.run(
#		"output/testRelease",
#		"/home/jeff/duke/donaldLab/dshot/input/raw/2KDC.pdb",
#		3,
#		"/home/jeff/duke/donaldLab/dshot/input/2KDC.experimental.rdc",
#		"/home/jeff/duke/donaldLab/dshot/input/2KDC.experimental.onlyDisulfideBonds.noe",
#		0,
#		2.0, #0.75,
#		pathInUntrustedDistanceRestraints = "/home/jeff/duke/donaldLab/dshot/input/2KDC.experimental.onlyPres.noe",
#		trustedPaddingPercent = 0.03,
#		untrustedPaddingPercent = 0.03,
#		numMinimizationProcesses = 6,
#		numMsrsProcesses = 6
#	)

	# RECOMB GB1 reassigned
#	disco.run(
#		"output/testRelease",
#		"/home/jeff/duke/donaldLab/dshot/input/raw/1Q10.pdb",
#		2,
#		"/home/jeff/duke/donaldLab/dshot/input/1Q10.experimental.rdc",
#		"/home/jeff/duke/donaldLab/workspace/jdshot/output/gb1-reassigned/1Q10.reassigned.noe",
#		0.0,
#		0.005,
#		alignmentTensorAxis = AlignmentTensorAxis.X,
#		trustedPaddingPercent = 0.0,
#		numMinimizationProcesses = 6,
#		numMsrsProcesses = 6
#	)
	

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