# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#
#	LINE ANALYSIS WITH REGULAR GRID OR MARKOV CHAIN MONTE CARLO METHODS
#
# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

# CASSIS scripts
import ScriptEnvironment
import ScriptUtilities
from ScriptRadexLA      import LteRadexLA
from ScriptLineAnalysis import LineAnalysis
from Parameter          import Parameter
from Components         import Components
from Component          import Component
from Constants          import Constants
from FrequencyScale     import FrequencyScale
from PhysicalModel      import PhysicalModel
from cassisStats        import printStats
from FileReader         import FileReader
from LineSpectrumData   import LineSpectrumData
from Range import Range

# CASSIS
from eu.omp.irap.cassis.gui                           import PanelView
from eu.omp.irap.cassis.file                          import FileReaderCassis
from eu.omp.irap.cassis.common.axes                   import XAxisVelocity
from eu.omp.irap.cassis.common                        import UtilArrayList
from eu.omp.irap.cassis.common                        import TypeFrequency
from eu.omp.irap.cassis.common                        import RadexException
from eu.omp.irap.cassis.gui.model.parameter.continuum import ContinuumParameters
from eu.omp.irap.cassis.database.access               import AccessDataBase
from eu.omp.irap.cassis.gui.model.parameter.observing import ObservingMode
from eu.omp.irap.cassis.cassisd.model.synthetic       import RadexConfiguration
from eu.omp.irap.cassis.file                          import FileManagerFits

# Python
import math
import sys
import os
import copy
import time
import decimal

# Java
from java.io import File
from javax.swing import JOptionPane
from java.lang import System
from java.util import Random


# FUNCTIONS {{{1
# =============================================================================

# HELPERS {{{2

def dmap(func, record):
    """A map for a directory"""
    values = (func(v) for k, v in record.items())
    return dict(zip(record, values))


def wrapLines(lines):
	"""Make sure the lines are contained into lists"""
	values = lines.values()
	if all(isinstance(x, (Range, int, float)) for x in values):
		return dmap(lambda x: [x], lines)
	elif all(isinstance(x, (list, tuple)) for x in values):
		return lines
	else:
		msg = 'selectedLines, rmsLines and calLines must have the format:\n'
		msg += '{"1": r1, "2": r2}\n'
		msg += 'or\n'
		msg += '{"1": [r1a, r1b], "2": [r2a, r2b]}'
		raise TypeError(msg)


def dec(value):
	"""Return the decimal value"""
	if value is None:
		return value
	else:
		return decimal.Decimal(str(value))


# TODO Refactor updateParameterValues because it's ugly
def updateParameterValues(parameterValues, logModeForParameters):
	newParameterValues = copy.deepcopy(parameterValues)
	for i in range(len(newParameterValues)):
		for name, useLog in logModeForParameters[i].items():
			if useLog:
				value = newParameterValues[i][name]
				if isinstance(value, decimal.Decimal):
					loggedValue = value.log10()
				else:
					loggedValue = math.log10(value)
				newParameterValues[i][name] = loggedValue
	return newParameterValues


# TODO Use these apply-functions
def applyOnParameter(parameter, fn):
	"""Apply a function on a Parameter"""
	_parameter = copy.deepcopy(parameter)
	_parameter.values = [dec(fn(x)) for x in _parameter.values]
	return _parameter


def applyOnComponent(component, names, fn):
	"""Apply a function on a the parameter with the given name in the component"""
	_component = copy.deepcopy(component)
	for name in names:
		newParameter = applyOnParameter(getattr(_component, name), fn)
		setattr(_component, name, newParameter)
	return _component


def applyOnComponents(components, names, fn):
	"""Apply a function on a the parameter with the given name in the components"""
	_components = copy.deepcopy(components)
	newComponents = map(lambda c: applyOnComponent(c, names, fn), _components.components)
	_components.components = newComponents
	return _components


def printHeader(title, width=80, symbol="="):
	nbOfSymbols = width - len(title) - 2
	remains = width - (2 * (nbOfSymbols / 2) + len(title) + 2)
	print("\n%s %s %s" % (symbol * (nbOfSymbols / 2), title, symbol * ((nbOfSymbols / 2) + remains)))


def warning(msg):
	print("WARNING: %s" % msg)


def uniq(sample):
	return list(set(sample))


def is_uniq(sample):
	return len(uniq(sample)) == 1


def determineTelescope(telescope, numLineUser):
	"""Determine the right telescope for the given line number"""
	if isinstance(telescope, dict):
		# telescope = [v for (k,v) in telescope.items() if numLineUser == k][0]
		telescope = telescope[numLineUser]
	return telescope


def analyzeTypeFreq(vlsrFile):
	"""Determine the right type frequency from the value of vlsrFile"""
	if vlsrFile == 0:
		typeFreq = TypeFrequency.SKY
	else:
		typeFreq = TypeFrequency.REST
	return typeFreq


def analyzeUnitRanges(selectedLines):
	if len(selectedLines.keys()) > 0:
		rangeProvided = selectedLines.values()[0][0]
		return rangeProvided.unit
	else:
		return "GHz"


def analyzeSelectedLines(selectedLines, inputData):
	if len(selectedLines.keys()) == 1 and "*-*" in selectedLines:
		indexes = range(1, len(inputData) + 1)
		selectedRange = selectedLines["*-*"]
		selectedLines = dict((str(x), selectedRange) for x in indexes)
	return selectedLines


def normalizeCollisioner(name):
	"""Rename properly the collisioner
	n_x_mol -> n_x-mol	

	where 
	- x can be o (ortho) or p (para)
	- mol is a molecule name (ex: H2)

	Parameters
	----------
	name : str
		The name of the collisioner

	"""
	return reverse(reverse(name).replace("_", "-", 1))


def parameterIsCollisioner(name, collisioners):
	"""
	
	Parameters
	----------
	name : str
		The name of the parameter

	"""
	return name.startswith("n_") and name[2:] in collisioners


def handleColumnNames(i, cols, col, paramName, collisioners, algorithm):
	# Logarithmic unit
	isLog = False
	if 'Log' in col:
		col = col.strip('Log[').strip(']')
		isLog = True

	s = col.split(" ")
	s.insert(1, "_%d" % (i))
	if len(s) == 2:
		# Case 1: param_i
		# For collisioners replace "_" with "-" when ortho and para
		if parameterIsCollisioner(paramName, collisioners):
			s[0] = normalizeCollisioner(s[0])

		cols.append("%s%s" % (s[0], s[1]))

	# Case 2: n_param_i and isLog == True
	elif len(s) > 2 and not isLog:
		# For collisioners replace "_" with "-" when ortho and para
		if parameterIsCollisioner(paramName, collisioners):
			s[0] = normalizeCollisioner(s[0])

		if algorithm == "RG":
			cols.append("%s%s%s" % (s[0], s[1], s[2]))
		else:
			cols.append("%s%s" % (s[0], s[1]))

	# Case 3: not param_i and (not n_param_i or isLog == False)
	else:
		# For collisioners replace "_" with "-" when ortho and para
		if parameterIsCollisioner(paramName, collisioners):
			s[0] = normalizeCollisioner(s[0])

		if algorithm == "RG":
			cols.append("%s%s %s" % (s[0], s[1], s[2]))
		else:
			cols.append("%s%s" % (s[0], s[1]))

	if isLog and algorithm == "RG":
		cols[-1] = "Log[%s]" % cols[-1]

	return cols


# CLASSES  {{{1

# EXCEPTIONS  {{{2

class LinesValidityError(Exception):
	"""Raised when all indexes of lines are not valid."""
	pass


class LinesEqualityError(Exception):
	"""Raised when all indexes of lines are not equal."""
	pass


class InvalidFormatError(Exception):
	"""Raised when the format of an input is invalid"""
	pass


def makeRange(pattern):
    if not isinstance(pattern, str):
        msg = """The input must be a string with a correct format.
		example: 
			{'1,3-5': f1, '2': f2} returns {1: f1, 2: f2, 3: f1, 4: f1, 5: f1}
			'*-*' returns all the lines in the data
		"""
        raise TypeError(msg)

    singles = pattern.split(",")
    rangeOfNumbers = []
    #
    for single in singles:
        if "-" in single:
            lower, upper  = single.split("-")[0], single.split("-")[-1]
            r = range(int(lower), int(upper)+1)
            rangeOfNumbers.extend(r)
        else:
            rangeOfNumbers.append(int(single))
    return rangeOfNumbers


def analyzeLines(linesAndRanges):
    newLinesAndRanges = {}
    for selectedLines, freqs in linesAndRanges.items():
        for line in makeRange(selectedLines):
            newLinesAndRanges[line] = freqs
    return newLinesAndRanges


def analyzeLinesAndRanges(inputData, linesAndRanges, name):
	newLinesAndRanges = analyzeLines(linesAndRanges)
	checkAllIndexesOfLinesAreValid(inputData, newLinesAndRanges, name)
	return newLinesAndRanges


def selectLinesWithEnoughPoints(nbOfPointsInLines, independenceFactor, nbOfIndependentPointsMin):
	"""Return the lines with enough points for the computation of the reduced chi2Min

	Indeed:
		chi2MinReduced = sum(chi2s*weight/((nbOfPoints/nbOfIndepPoints)-dof)/nbOfLines)

	which means that if nbOfPoints/nbOfIndepPoints = dof, the chi2MinReduced can't be computed
	because of a division by 0 (because there are not enough points).
	This function ensures that the number of points is sufficient for
	a relevant (and a possible) computation

	Usage
	-----
	linesWithEnoughPoints = selectLinesWithEnoughPoints(nbOfPointsInLines,
														independenceFactor,
														nbOfIndependentPointsMin)

	Parameters
	----------
	nbOfPointsInLines : dict
		a dictionary whose keys are the indexes of lines and the values
		the corresponding number of points
	independenceFactor : int
		the factor of independence selected by the user that
		allows to determine the number of independent points
		(nbOfIndependentPoints = nbOfPointsInLines / independenceFactor)
	nbOfIndependentPointsMin : int
		the miminum number of independent points per line
		(this number must be greater than the dof)

	Returns
	-------
	linesWithEnoughPoints : list
		the list of indexes of lines whose number of points satifisfies the criterion
		nbOfIndependentPoints > nbOfIndependentPointsMin

	"""
	linesWithEnoughPoints = []
	for indexOfLine, nbsOfPoints in nbOfPointsInLines.items():
		for nbOfPoints in nbsOfPoints:
			# If the criterion is satisfied and the line is not already in the
			# list, append the index of this line
			if math.floor(nbOfPoints / independenceFactor) > nbOfIndependentPointsMin and\
			   indexOfLine not in linesWithEnoughPoints:
				linesWithEnoughPoints.append(indexOfLine)
	return sorted(linesWithEnoughPoints)


def analyzeValuesSelected(values, rangesSelected):
	"""Return the values selected in the given range

	Usage:
	------
	analyzeValuesSelected(values, rangeSelected)

	Parameters
	----------
	values : an iterable
        all the values
	rangeSelected : iterable
        a value range [minValue, maxValue]

	Returns
	-------
	valuesSelected

	"""
	if len(rangesSelected) == 0:
		return []

	valuesSelected = []
	for rangeSelected in rangesSelected:
		valuesSelected.append(filter(lambda f: min(rangeSelected) <= f <= max(rangeSelected), values))
	return valuesSelected


# Unshamelessly stolen from
# http://docs.python.org/2/library/itertools.html#itertools.combinations
# because Python 2.5 doesn't have itertools.combinations (actually that's the
# real shame!)
def combinations(iterable, r):
	# combinations('ABCD', 2) --> AB AC AD BC BD CD
	# combinations(range(4), 3) --> 012 013 023 123
	pool = tuple(iterable)
	n = len(pool)
	if r > n:
		return
	indices = range(r)
	yield tuple(pool[i] for i in indices)
	while True:
		for i in reversed(range(r)):
			if indices[i] != i + n - r:
				break
		else:
			return
		indices[i] += 1
		for j in range(i+1, r):
			indices[j] = indices[j-1] + 1
		yield tuple(pool[i] for i in indices)


def analyzeReducePhysicalParam(reducePhysicalParam, parameterName):
	if parameterName == "tex" or parameterName == "tkin":
		reductionFactor = reducePhysicalParam["temp"]
	else:
		reductionFactor = reducePhysicalParam[parameterName]

	if reductionFactor == 0:
		msg = "The reducePhysicalParam must be different from 0. "
		msg += "Please correct %s" % parameterName
		raise ValueError(msg)
	return reductionFactor


# INITIALIZATION/CREATION {{{2

def initializeParameters(component, parameters):
	"""Initialize the parameters

	Usage
	-----
	initializeParameters(component, parameters)

	
	Parameters
	----------
	component : Component

	parameters : dictionary
	{nmol: nmol1, size: size1, ...}


	Returns
	-------
	initializedParameters : dictionary

	"""
	initializedParameters = {}
	#
	for name in component.names:
		if name in parameters.keys():
			initializedParameters[name] = parameters[name]
		else:
			initializedParameters[name] = randomizeParameter(component, name)
	return initializedParameters


def makeRandomParameters(component):
	"""Store parameters in a convenient data structure"""

	temperatureName = Constants.MODEL_TEMP[component.model]
	randomParameters = {}
	randomParameters["nmol"] = randomizeParameter(component, "nmol")
	randomParameters[temperatureName] = randomizeParameter(component, temperatureName)
	randomParameters["fwhm"] = randomizeParameter(component, "fwhm")
	randomParameters["size"] = randomizeParameter(component, "size")
	randomParameters["vlsr"] = randomizeParameter(component, "vlsr")

	# Handle iso when available
	if hasattr(component, 'iso'):
		randomParameters["iso"] = randomizeParameter(component, "iso")

	# Handle collisioner, in RADEX, when available
	if component.model.lower() == "radex":
		for collisioner in Constants.COLLISIONER_PARAMETERS:
			if hasattr(component, collisioner):
				randomParameters[collisioner] = randomizeParameter(component, collisioner)

	return randomParameters


def makeBounds(component):
	names = component.names
	bounds = {}
	for name in names:
		minValue = float(getattr(component, name).get_min())
		maxValue = float(getattr(component, name).get_max())
		bounds[name] = (minValue, maxValue)
	return bounds


def randomizeParameter(component, name):
	r = Random().nextDouble()
	parameter = getattr(component, name)
	delta = parameter.get_max() - parameter.get_min()
	value = float(parameter.get_min()) + r * float(delta)
	return value


def createRowOfData(thetas, logModeForParameters):
	rowOfData = ""
	orderedNames = Constants.ORDERED_PARAMETER_NAMES
	for theta in thetas:
		names  = [name for name in orderedNames if name in theta.keys()]
		values = [theta[value] for value in names]
		for name, value in zip(names, values):
			if name in logModeForParameters and logModeForParameters[name]:
				value = math.pow(10, value)
				rowOfData += "%.4g\t" % value
			elif name.startswith(Constants.COLLISIONER_PREFIX):
				rowOfData += "%.4g\t" % value
			else:
				rowOfData += "%f\t" % value
	return rowOfData


def createParameters(parameterValues, parameterLogModes):
	"""Create a dictionary of parameters from the given Component object
	
	Parameters
	----------
	parameterValues : dict
		parameterValues = {paramName1: value1, paramName2: value2, ...}

	parameterLogModes : dict
		parameterLogModes = {paramName1: TrueOrFalse, paramName2: TrueOrFalse, ...}

	Returns
	-------
	parameters : dict
		parameters = {paramName1: value1', paramName2: value2'}

	"""
	parameters = {}
	names = parameterValues.keys()
	for name in names:
		if parameterLogModes.has_key(name) and parameterLogModes[name]:
			value = 10**parameterValues[name]
		else:
			value = parameterValues[name]
		# value = parameterValues[name]
		parameters[name] = value
	return parameters

	
# TODO Don't use components. 
# Use a dictionary instead:
# interacting, geometry, collisionFile, collisioner
def createPhysicalModel(parameters, component):
	physicalModel = PhysicalModel(parameters)

	# Set the extra attributes necessary for the creation
	# of the component
	physicalModel.interacting   = component.interacting
	physicalModel.geometry      = component.geometry
	physicalModel.collisionFile = component.collisionFile
	physicalModel.collisioner   = component.collisioner

	
	return physicalModel


# TODO Don't use components. 
# Use a dictionary instead:
# interacting, geometry, collisionFile, collisioner
def createPhysicalModels(parameterValues, components):
	physicalModels = []
	for i in range(len(components.components)):
		component  = components.components[i]
		values     = parameterValues[i]
		logModes   = getAttrOfParameters("log_mode", component)
		parameters = createParameters(values, logModes)
		physicalModels.append(createPhysicalModel(parameters, component))
	return physicalModels
	

def createCassisComponent(tags, physicalModel):
	"""Create a component from a set of parameters and an isotopic ratio (optional).

	Usage:
		component = createCassisComponent(tags, physicalModel)

	This will create a molecule 1 with molecular density nmol1 and a molecule 2
	with a molecular density nmol2 = nmol1/iso1.

	If N isotopic ratios are provided, N+1 components are created.

	Inputs:
		* tags         : the list of tags of the molecule for the current
						component
		* physicalModel  : the set of parameters

	"""
	molecules = []
	for i in range(len(tags)):
		if i > 0 and hasattr(physicalModel, 'iso'):
			iso = physicalModel.iso
		else:
			iso = 1

		if physicalModel.model == "Lte":
			molecules.append(ScriptUtilities.molecule(tag=tags[i],
				nmol=float(physicalModel.nmol / iso),
				tex=float(physicalModel.tex),
				fwhm=float(physicalModel.fwhm),
				size=float(physicalModel.size)))

		elif physicalModel.model == "Radex":

			msg = "The number of moltags must be equal to the number of collision files."
			assert len(physicalModel.collisionFile) == len(tags), msg

			collisionFile = physicalModel.collisionFile[i]
			collisionString = "%s" % collisionFile

			for collisioner in physicalModel.collisioner:
				paramValue = filter(lambda x: x[0].endswith(collisioner), physicalModel.parameters.items())[0][1]
				collisioner = collisioner.replace("_", "-")
				collisionString += "&(%s;%g)" % (collisioner, paramValue)

			molecules.append(ScriptUtilities.molecule(tag=tags[i],
							nmol=float(physicalModel.nmol / iso),
							tkin=float(physicalModel.tkin),
							fwhm=float(physicalModel.fwhm),
							size=float(physicalModel.size),
							collision=collisionString))


		# Interacting, geometry
	if physicalModel.interacting is not None:
		interacting = physicalModel.interacting
	else:
		interacting = True

	if physicalModel.geometry != "":
		geometry = physicalModel.geometry
	else:
		geometry = "slab"

	# if physicalModel.tex != "":
	# 	tex = physicalModel.tex
	# else:
	# 	tex = "100.0"
	#
	# if physicalModel.tkin != "":
	# 	tkin = physicalModel.tkin
	# else:
	# 	tkin = "10.0"
	#
	# if physicalModel.size != "":
	# 	size = physicalModel.size
	# else:
	# 	size = "3.0"

	component = ScriptUtilities.component(model=physicalModel.model,
					vlsr=float(physicalModel.vlsr),
					molecules=molecules,
					interacting=interacting,
					geometry=geometry)
					# tex=float(tex),
					# tkin=float(tkin),
					# size=float(size))
	return component


# TODO Don't use components. 
# Use a dictionary instead:
# interacting, geometry, collisionFile, collisioner
def createCassisComponents(moltags, components, parameterValues):
	"""Create all the Cassis components of a moltag
	
	Parameters
	----------
	moltags : list
		The moltags 

	components : Components
		The Components created by the user with Component

	parametereValues: dict
		A dictionary {paramName1: value1, paramName2: value2, ...}

	"""
	physicalModels = createPhysicalModels(parameterValues, components)
	return [createCassisComponent(moltags, pm) for pm in physicalModels]


# COMPUTATIONS  {{{2

def computeChi2Reduced(chi2s, nbOfPoints, dof, nbOfIndependentPoints, weights):
	"""Return the reduced chi2
	Usage:
		chi2Reduced = computeChi2Reduced(chi2s, lineCount, pointCounts, dof, nbOfIndependentPoints)

	Inputs:
	-------
	chi2s : list
		the list of chi2
	pointCounts : list
		the list of number of points
	dof : int
		the degree of freedom
	nbOfIndependentPoints : int
		the number of independent points

	Returns
	-------
	chi2Reduced : float
		the reduced chi2

	"""
	nbOfLines   = len(chi2s)
	indexes     = range(nbOfLines)
	wks         = sorted(weights.keys())
	# sumOfWeights = sum(x[0] for x in weights.values())
	chi2Reduced = sum(chi2s[i] * weights[wks[i]][0]/(nbOfPoints[i]/nbOfIndependentPoints - dof) for i in indexes) / nbOfLines #* sumOfWeights
	return chi2Reduced


def chi2Computation(intensityData, intensityModel, rms, cal, continuum):
	return (intensityData - intensityModel)**2 / (rms**2 + cal**2 * (intensityData - continuum)**2)


def computeChi2(data, model, rangesSelected, rms, cals, weights, continuums):
	"""Compute the Chi2 value between the data and the model

	The computation takes into account the frequency ranges to be kept.
	"""
	frequencyData  = data.getFrequenciesSignal()
	intensityData  = data.getIntensities()
	intensityModel = model.getIntensities()
	chi2		   = 0
	nbOfPointsOfChi2   = 0

	# Loop through the channels
	for i in range(len(frequencyData)):

		# Case 1 - The range is provided:
		# For each frequency in the data, check that it is within the
		# frequency ranges not to be excluded from the chi2 computation
		if len(rangesSelected) > 0:

			# Loop through the frequency ranges 
			for j in range(len(rangesSelected)):
				if min(rangesSelected[j]) <= frequencyData[i] <= max(rangesSelected[j]):
					chi2 += chi2Computation(intensityData[i], intensityModel[i], rms[j], cals[j], continuums[i])# / weights[j]
					nbOfPointsOfChi2 += 1

		# Case 2 - The range is not provided:
		# Take the whole spectrum for the computation of the chi2 (the
		# spectrum was previously selected for this computation thanks to
		# selectedLines)
		else:
			print("Taking the whole spectrum for the computation of the chi2 ")
			j = 0
			chi2 += chi2Computation(intensityData[i], intensityModel[i], rms[j], cals[j], continuums[i])# / weights[j]
			nbOfPointsOfChi2 += 1

	return (chi2, nbOfPointsOfChi2)


# TODO Refactoring: Simplify this function
def createColumnNames(components, algorithm):
	componentCount = len(components.components)
	cols = []
	cpt_component_species = 0
	for i in range(componentCount):
		component    = components.components[i]
		
		cpt_component_species = cpt_component_species + 1
		collisioners = component.collisioner
		paramNames   = component.names

		for paramName in paramNames:
			# isLog = False
			if paramName == "tex" and not hasattr(component, "tex") and hasattr(component, "tkin"):
				paramName = "tkin"
			elif paramName == "tkin" and not hasattr(component, "tkin") and hasattr(component, "tex"):
				paramName = "tex"

			col = getattr(component, paramName).x_axis_label
			cols = handleColumnNames(cpt_component_species, cols, col, paramName, collisioners, algorithm)

	cols.insert(0, "#")
	if algorithm == "MCMC":
		cols.extend(["chi2", "rate"])
	else:
		cols.extend(["chi2"])

	cols = "\t".join(cols)
	cols += "\n"
	return cols


def calculateFacWithCutOff(i, cutOff, ratioAtCutOff=0.01):
	if i < cutOff:
		fac = (-(1-ratioAtCutOff) * i / cutOff + 1.0)
	else:
		fac = ratioAtCutOff
	return fac


def calculateAlpha(fac, parameterRange, reductionFactor):
	"""Calculate alpha using the range of the given parameter

	reductionFactor - a number > 1

	"""
	alpha = fac * (max(parameterRange) - min(parameterRange)) / reductionFactor
	return alpha


# CHECKS {{{2

def checkRangeValidity(linesAndRanges, data):
	"""Check that the selected ranges are included in the right spectra

	Usage
	-----
	checkRangeValidity(linesAndRanges, data)

	linesAndRanges : a dictionary
					the spectra associated with the ranges
	data :

	Throws
	------
	RangeError if the condition is not satisfied

	"""
	invalidRanges = []
	validRanges   = []
	for line, ranges in linesAndRanges.items():
		idx = line - 1
		frequency = data[idx].getFrequenciesSignal()
		for r in ranges:
			if not isinrange([x*1E3 for x in r], frequency):
				invalidRanges.append({line: r})
				validRanges.append({line: [min(frequency), max(frequency)]})

	if len(invalidRanges) > 0:
		msg = "Invalid ranges. Please check your inputs in selectedLines.\n"
		for i in range(len(invalidRanges)):
			for line, _range in invalidRanges[i].items():
				msg += "Line %d: %r is not in %r\n" % (line, [x*1E3 for x in _range], validRanges[i][line])
		raise RangeError(msg)


def checkAllIndexesOfLinesAreValid(inputData, linesAndRanges, name):
	invalidIndexesOfLines = filter(lambda x: x not in range(1, len(inputData)+1), linesAndRanges.keys())
	if len(invalidIndexesOfLines) > 0:
		msg = "All the lines in %s must be valid lines. Lines: %r are not correct." % (name, invalidIndexesOfLines)
		msg += "\nPlease choose your lines between %d and %d" % (1, len(inputData)+1)
		raise LinesValidityError(msg)


def checkAllIndexesOfLinesAreEqual(allLinesAndRanges):
	# Make combinations of indexes of lines grouping them 2 by 2 to compare them
	allIndexesOfLines = [linesAndRanges.keys() for linesAndRanges in allLinesAndRanges]
	combins = combinations(allIndexesOfLines, 2)
	allIndexesOfLinesAreEqual = all(sorted(c[0]) == sorted(c[1]) for c in combins)

	if not allIndexesOfLinesAreEqual:
		msg = "The lines of selectedLines, rmsLines, calLines and telescope must all be the same.\n"
		msg += "Please check your inputs. "
		raise LinesEqualityError(msg)


def checkInitComponents(initComponent):
	"""Check the initialization of the components in MCMC"""
	msg = "The arguments must be a list of 2 elements."
	if not isinstance(initComponent, (list, tuple)):
		raise TypeError(msg)

	msg = "The %s must contain only 2 elements" % initComponent.__class__.__name__
	if len(initComponent) != 2:
		raise ValueError(msg)

	msg = "The first argument must be a Component."
	if not isinstance(initComponent[0], Component):
		raise TypeError(msg)

	msg = "The second argument must be a dictionary."
	if not isinstance(initComponent[1], dict):
		raise TypeError(msg)


# PREDICATES {{{2

def testThetaIsInRange(thetas, thetaRange):
	componentCount = len(thetas)
	test = True
	for i in range(componentCount):
		names = thetas[i].keys()
		for name in names:
			test = test and (min(thetaRange[i][name]) <= thetas[i][name] <= max(thetaRange[i][name]))
	return test


def testChi2ValIsBetter(previousChi2Val, currentChi2Val):
	# ???
	r = math.log(Random().nextDouble())
	return (previousChi2Val - currentChi2Val) / 2 > r


def fileExists(filename, defaultPath=""):
	if os.path.exists(filename) or\
		os.path.exists(os.path.join(defaultPath, filename)):
		return True
	else:
		return False


def computeAcceptanceRate(count1, count2):
	"""docstring for computeAcceptanceRate"""
	return float(count2) / float(count1)


def writeDataToFile(acceptanceCount, thetas, rate, chi2Val, filename):
	textChi2 = "%f\t" % acceptanceCount
	textChi2 += '\t'.join(str(thetas[i][name]) for i in range(len(thetas)) for name in ['nmol', 'tex', 'fwhm', 'size', 'vlsr'])
	textChi2 += "\t%f\t%d\n" % (rate, chi2Val)

	validModelsFile = open(filename, 'a')
	validModelsFile.write(textChi2)
	validModelsFile.close()


def select_range(initial_range, selected_range, action="exclude"):
	"""Cut a initial_range from a range using the given action as a rule

	Usage
	-----
	final_range = select_range(initial_range, selected_range, action="exclude")

	Parameters
	----------
	initial_range : an iterable
	selected_range : the selected range (a min and a max value for keeping or excluding the
	elements of the initial range)
	action : determines whether the selected range must be kept ("keep") or
	excluded (exclude) from the initial range

	Returns
	-------
	final_range : the sample correctly cut

	Examples:
	>>> print(select_range(range(10), (2, 8), action="exclude"))
	[0, 1, 9]

	>>> print(select_range(range(10), (2, 8), action="keep"))
	[2, 3, 4, 5, 6, 7, 8]

	"""
	final_range = []
	choice = {"exclude": False, "keep": True}
	for x in initial_range:
		if (min(selected_range) <= x <= max(selected_range)) == choice[action]:
			final_range.append(x)
	return final_range


def keep(initial_range, selected_range):
	return select_range(initial_range, selected_range, action="keep")


def exclude(initial_range, selected_range):
	return select_range(initial_range, selected_range, action="exclude")


def isiterable(obj):
	"""Determine whether an object is iterable or not."""
	try:
		iter(obj)
		return True
	except TypeError:
		return False


def isinrange(arange, rangeRef):
	"""Return True if the range is included in the reference range

	Usage:
		a = isinrange(arange, rangeRef)

	Inputs:
	`arange`  : a range of values defined by a min and a max
	`rangeRef`: the reference range defined by a min and a max

	Details:
	This functions is basically a wrapper around:
		a = (min(arange) >= min(rangeRef)) and (max(arange) <= max(rangeRef))

	"""
	# if isinstance(arange, (list, tuple)) and len(arange) == 0:
	if isiterable(arange) and len(arange) == 0:
		return False
	else:
		# if isinstance(arange, (list, tuple)):
		if isiterable(arange):
			return (min(arange) >= min(rangeRef)) and (max(arange) <= max(rangeRef))
		# elif isinstance(arange, (int, float)):
		else:
			return (arange >= min(rangeRef)) and (arange <= max(rangeRef))


def calculateAverage(sample):
	return (sum(sample) / len(sample))


def reverse(iterable):
	"""Reverse an iterable"""
	# With reversed iterator and list comprehension
	# return "".join([c for c in reversed(iterable)])
	# With stride notation
	return iterable[::-1]


def computeNbOfModelsInRG(comps, isoUnique):
	"""Compute the number of models in regular grid

	nbOfModels = computeNbOfModelsInRG(components, isoUnique)

	Inputs:
		cps      - the componentParameter objects
		isoUnique - Take the intersection of all the values of the isos (True)
					or take the union of all the values of the isos (False)
					[OPTIONAL]
					[DEFAULT: False]

	Returns:
	nbOfModels : the number of models to ne computed with the regular grid
				 method

	"""
	components = Components(comps, isoUnique=isoUnique, method="RG")

	c = 0
	for component in components.components:
		print("Component %d:" % (c + 1))
		for name, values in component.parameters.items():

			# Case: 1 value if min = max value
			if values["min"] == values["max"]:
				print("%10s\t: 1 value" % name)

			# Case 2: if nstep is present
			elif values.has_key('nstep'):
				print("%10s\t: %d values" % (name, values['nstep']))

			# Case 3: if nstep is not present, recompute it
			else:
				p = Parameter(name, values['min'], values['max'], values['step'])
				print("%10s\t: %d values" % (name, p.nstep))
		c += 1
		print

	return len(components.models)


def askConfirmation(msg, title="WARNING"):
	answer = JOptionPane.showConfirmDialog(None,
										msg + "\nDo you want to continue?",
										title,
										JOptionPane.YES_NO_OPTION)

	if (answer != JOptionPane.YES_OPTION):
		sys.exit("Exited before the computation.")


def getModelIndexes(nbOfModels):
	"""
	Return the values of the model every modulo N
	where N < nbOfModels < 10*N

	"""
	values = []
	for i in range(0, 10):
		min = math.pow(10, i)
		max = math.pow(10, i + 1)
		if min == 10:
			mod = 100
		else:
			mod = min / 10
		if min < nbOfModels+1 <= max:
			values = filter(lambda x: x % min == 0 or (x - min) % mod == 0, range(nbOfModels))
	return values


def getMolName(moltag):
	db = AccessDataBase.getDataBaseConnection()
	return db.getMolName(str(moltag))


def getDatabaseName(tags):
	db = AccessDataBase.getDataBaseConnection()
	return [db.getDataBaseOf(str(tag)) for tag in tags]


def getAttrOfParameters(attr, component):
	"""Return the attribute of all the parameters in a compononent
	
	Parameters
	----------
	attr - str
		Any valid attribute of a Parameter object
		Ex: "log_mode" or "min" or "max" or "step" or "nstep" 

	component - Component
		The component

	Returns
	-------
	{paramName: valueOfAttr, ...}

	Example
	-------
	>>> logModes = getAttrOfParameters("log_mode", component)
	>>> print(logModes)
	{"nmol": True, "vlsr": False, ...}

	"""
	names = component.names
	attrs = [getattr(getattr(component, name), attr) for name in names]
	return dict( zip(names, attrs))


# CLASSES  {{{1
# =============================================================================

class RangeError(Exception):
	"""Raised when the selected range is not within the initial range"""
	pass

def analyzeTelescope(telescope):
	def message(name, value):
		m = "WARNING. The %s file: %s does not exist.\n" % (name, value)
		m += "Please change your input."
		return m

	# A dictionary of telescopes
	if isinstance(telescope, dict):
		for t in telescope.values():
			if not fileExists(t, Constants.TELESCOPE_DIR):
				sys.exit(message("telescope", t))
			else:
				telescope = dict((str(k), v) for k, v in telescope.items())
				telescope = analyzeLines(telescope)

	# A single string telescope
	elif isinstance(telescope, str) or isinstance(telescope, unicode):
		if not fileExists(telescope, Constants.TELESCOPE_DIR):
			sys.exit(message("telescope", telescope))

	# No telescope? Be kind and explain patiently of the script works...
	else:
		msg = """You must provide either a single telescope or the lines on which different telescopes are applied. 
- For a single telescope: telescope=telescope1
- Example for different telescopes: telescope={'1-3': telescope1, '4,5': telescope2}

"""
		sys.exit(msg)
	return telescope

def get_continuum_param(continuum):
		#if isinstance(continuum, float) or isinstance(continuum, int):
		#	 continuumParameters = ContinuumParameters(float(continuum))
		#else:
		continuumParameters = ContinuumParameters(continuum)
		return continuumParameters

class UserInputs(object):
	"""A class storing all the inputs from the user"""
	def __init__(self, selectedLines={}, telescope=None, tuningRange=[],
			tuningBand=60, aijMin=0, moltags=[],
			eup=[0.0, 150.0], jup=["*", "*"], kup=["*", "*"], lup=["*", "*"], mup=["*", "*"],
			template="ISM",
			inputFile="", outputFile="",
		    plotTitle="",
			nbIndPoints=1,
			tmb2ta=False, isoUnique=False,
			tbg=2.73,
			warning=True,
			lineNumber=None,
			rmsLines={}, calLines={}, observing_mode="PSw/DBSw", 
			continuum=0, indice_hdu=0):

		super(UserInputs, self).__init__()
		self.telescope   = analyzeTelescope(telescope)
		self.cal         = 0.0
		self.tmb2ta      = tmb2ta
		self.tbg         = tbg
		self.tuningRange = tuningRange
		self.tuningBand  = tuningBand
		self.aijMin      = aijMin
		self.moltags     = moltags
		self.components  = None

		# Quantum numbers
		self.eup      = eup
		self.jup      = jup
		self.kup      = kup
		self.lup      = lup
		self.mup      = mup

		self.rmsModel = 0.0
		self.chi2Min  = 1e16
		self.chi2Max  = 1e20
		self.text     = ""
		self.thetaOptimal = None

		printHeader("INPUT")

		# Input data file
		self.inputFile = inputFile

		if plotTitle == "":
			self.plotTitle = os.path.basename(self.inputFile)
			print("No plotTitle specified. plotTitle is set to %s" % self.plotTitle)
		else:
			self.plotTitle = plotTitle

		if outputFile == "":
			self.outputFile = self.inputFile + "_out.dat"
			print("No outputFile specified. outputFile is set to %s" % self.outputFile)
		else:
			self.outputFile = outputFile

		# Make sure selectedLines, rmsLines and calLines are lists
		selectedLines = wrapLines(selectedLines)
		rmsLines      = wrapLines(rmsLines)
		calLines      = wrapLines(calLines)

		self.template    = template
		self.nbIndPoints = nbIndPoints
		self.isoUnique   = isoUnique
		self.warning     = warning
		self.unit        = analyzeUnitRanges(selectedLines)
		self.bestPhysicalModels = []

		# Create species if it exists
		self.species = []
		self.lineAnalysis =  LineAnalysis()
		for moltag in moltags:
			spec = self.lineAnalysis.species(template=template, tag=moltag)
			if spec is None:
				raise ValueError("WARNING. species or template does not exist. Please change your input.")
			self.species.append(spec)

		# Check the inputs before reading the data
		self.indice_hdu = indice_hdu
		self.checkInputFile(self.inputFile)
		self.inputData = self.readInputData(self.inputFile, self.indice_hdu)

		# Analyze the selectedLines
		selectedLines = analyzeSelectedLines(selectedLines, self.inputData)

		# ANALYZE AND CHECK THE LINES INDEXES AND RANGES PROVIDED
		# =======================================================
		# TODO Rename notExcludedLines -> selectedLinesAndRanges
		self.selectedLines = analyzeLinesAndRanges(self.inputData, selectedLines, "selectedLines")
		self.rmsLines      = analyzeLinesAndRanges(self.inputData, rmsLines, "rmsLines")
		self.calLines      = analyzeLinesAndRanges(self.inputData, calLines, "calLines")

		if observing_mode == ObservingMode.FSW.toString():
			self.observing_mode = ObservingMode.FSW
		else:
			self.observing_mode = ObservingMode.PSDBS
		self.continuum=continuum
		self.checkAndFixRmsIfNeeded()

		allLines = [self.selectedLines, self.rmsLines, self.calLines]
		if isinstance(self.telescope, dict):
			allLines.append(self.telescope)
		checkAllIndexesOfLinesAreEqual(allLines)

		# Temporarily set the weights to 1.0 instead of letting the user define them
		self.weightLines = dict((x, [1.0]) for x in self.selectedLines)

		# Lines that can't be computed by LteRadexLA
		self.linesNotComputed = []

		# CHECK THE RANGES PROVIDED ARE IN THE SELECTED SPECTRA
		# =====================================================
		if self.unit == "GHz":
			checkRangeValidity(self.selectedLines, self.inputData)

		self.selectedNumLines = sorted(self.selectedLines.keys())
		if lineNumber is not None and\
		   lineNumber not in self.selectedNumLines:
			msg = "%d is not a valid lineNumber.\n" % lineNumber
			msg += "lineNumber must be the number of a line amongst the selected ones %r" % self.selectedNumLines
			raise ValueError(msg)
		else:
			self.lineNumber = lineNumber

		print("Number of lines selected: %d" % len(self.selectedLines))

	
	def checkAndFixRmsIfNeeded(self):
		rmsChanged = False
		for key, val in self.rmsLines.iteritems():
			for i, v in enumerate(val):
				if v == 0.0:
					self.rmsLines[key][i] = 1e-6
					rmsChanged = True
		if rmsChanged:
			print("WARNING: At least one rms value was set to 0.0 which is not allowed. This value was replaced by 1e-6 K.")

	def warnNoValueSelected(self, i):
		msg = "Line %d - No range selected\n" % i
		if self.warning:
			askConfirmation(msg)
		else:
			warning(msg)

	def checkData(self, data):
		try:
			intensityData = data.getIntensities()
			return intensityData
		except AttributeError:
			msg = "Intensity can't be found in the provided data."
			raise AttributeError(msg)

	def checkModel(self, model):
		try:
			intensityModel = model.getIntensities()
			return intensityModel
		except AttributeError:
			msg = "The model can't be computed. Check your inputs.\n"
			msg += "See the output in the terminal for more details."
			raise AttributeError(msg)

	def checkInputFile(self, inputFile):
		def message(name, value):
			m = "WARNING. The %s file: %s does not exist.\n" % (name, value)
			m += "Please change your input."
			return m

		# Input data file
		if not os.path.exists(inputFile):
			sys.exit(message("input data", inputFile))

	def rescaleModel(self, intensityData, intensityModel):
		if len(intensityData) != len(intensityModel):
			print("Data and model don't have the same length.")
			print("Rescaling model...")
			self.lineAnalysis.rescale(data, modelSpectrum)
			intensityData  = data.getIntensities()
			intensityModel = modelSpectrum.getIntensities()
		return intensityData, intensityModel

	def convertRangeToMHz(self, data, numLineUser, silent=1):
		"""Convert the range provided to MHz"""
		rangeSelected = []
		frequencySignal = data.getFrequenciesSignal()
		fScale			= FrequencyScale(data.getVlsr(), data.getFreqRef())
		for selectedRange in self.selectedLines[numLineUser]:

			# Handle case where the range is ]-inf; inf[
			if min(selectedRange) == float("-inf") and max(selectedRange) == float("inf"):
				selectedRange = frequencySignal[:]

			# GHz -> MHz
			elif self.unit == "GHz":
				selectedRange = [f * 1e3 for f in selectedRange]

			# (km/s|micrometer|cm-1) -> MHz
			elif self.unit in ("km/s", "micrometer", "cm-1"):
				if selectedRange.unit == self.unit:
					selectedRange = [fScale.toMHz(v, self.unit) for v in selectedRange]

			# Unknown unit
			else:
				msg = "%s is an unknown unit." % self.unit
				raise ValueError(msg)

			selectedRangeScaled = keep(frequencySignal, selectedRange)
			selectedRangeScaled.sort()

			if len(selectedRangeScaled) > 0:
				rangeSelected.append([min(selectedRangeScaled), max(selectedRangeScaled)])

			elif min(selectedRange) == max(selectedRange):
				selectedRange = []
			else:
				if silent == 0:
					print("Line %d - %r is not in [%f, %f]" % (numLineUser, selectedRange, min(frequencySignal), max(frequencySignal)))

		return rangeSelected


	def readInputData(self, inputFile, indice_hdu=0):
		"""Read the input data given by the user."""
		fname, extension = os.path.splitext(inputFile)
		self.commentedData = self.lineAnalysis.readData("",
										   fname,
										   [""],
										   extension, indice_hdu=0)
		speciesCount   = len(self.species)
		pointsCount    = 0
		numEl = 0

		tmpArray = [[None for x in range(len(self.commentedData))] for i in range(speciesCount)]

		for i in range(speciesCount):
			for cpt, spec in enumerate(self.commentedData):
				sd = self.lineAnalysis.split2(data=spec,
								tuningRange=self.tuningRange,
								tuningBand=self.tuningBand,
								plottingBand="Signal",
								plottingSort="Frequency",
								species=self.species[i],
								eup=self.eup,
								jup=self.jup,
								kup=self.kup,
								lup=self.lup,
								mup=self.mup,
								aijMin=self.aijMin)
				numEl += len(sd)
				tmpArray[i][cpt] = sd
			for spectrum in sd:
				pointsCount += spectrum.getSize()
			print("Number of lines found in moltag %d: %d" % (self.moltags[i], len(sd)))

		pointsTotal = pointsCount / self.nbIndPoints
		print("Number of points: %d; Number of independent points: %d" % (pointsCount, pointsTotal))
		print("Total number of lines in the data: %d" % numEl)

		self.inputData = []
		idSpecCom = [[0 for x in range(len(self.commentedData))] for i in range(speciesCount)]
		for i in range(numEl):
			fMin = float("inf")
			idSpec = 0
			idCom = 0
			for numSpec in range(speciesCount):
				for numCom in range(len(self.commentedData)):
					if len(tmpArray[numSpec][numCom]) > idSpecCom[numSpec][numCom]:
						line = tmpArray[numSpec][numCom][idSpecCom[numSpec][numCom]].getListOfLines().get(0)
						if line.getObsFrequency() < fMin:
							fMin = line.getObsFrequency()
							idSpec = numSpec
							idCom = numCom
			self.inputData.append(tmpArray[idSpec][idCom][idSpecCom[idSpec][idCom]])
			idSpecCom[idSpec][idCom] += 1
		return self.inputData
	
	def initComponentsForMCMC(self, *initComponents):
		"""Initialize the component with the values of the given parameters
		
		Usage
		-----
		initComponentsForMCMC(*initComponents)

		Parameters
		----------
		initComponents : lists of 2 elements
			A pair of elements: {param1: value1, ...}
			The first element is a Component object The second element is
			a dictionary whose keys are the parameter names and the values their
			corresponding initial values.

		Returns
		-------
		None

		"""
		msg = "The arguments must be at least one list of a Component and a tuple."
		if not(len(initComponents) > 0 and isinstance(initComponents[0], (list, tuple))):
			raise TypeError(msg)

		self.initComponents = []
		for initComponent in initComponents:
			checkInitComponents(initComponent)
			
			if isinstance(initComponent[0], Component):
				initComponent[1] = dict((Constants.renameTemp(name, initComponent[0].model), value) for name, value in initComponent[1].items())
			else:
				initComponent[1] = dict((name, value) for name, value in initComponent[1].items())
			
			msg = "The parameters must be valid parameter names."
			invalidParameters = [x for x in initComponent[1].keys() if x not in initComponent[0].names]
			if len(invalidParameters) > 0:
				msg += " Please rename: %r" % ", ".join("%s" % x for x in invalidParameters)
				raise ValueError(msg)

			msg = "The parameter %s must be within the selected bounds [%g; %g]."
			for name, value in initComponent[1].items():
				value = decimal.Decimal(str(value))
				minSelectedBound = getattr(initComponent[0], name).min
				maxSelectedBound = getattr(initComponent[0], name).max

				if getattr(initComponent[0], name).log_mode:
					minSelectedBound = 10**getattr(initComponent[0], name).min
					maxSelectedBound = 10**getattr(initComponent[0], name).max

				if (value < minSelectedBound) or (value > maxSelectedBound):
					raise ValueError(msg % (name, minSelectedBound, maxSelectedBound))

				if getattr(initComponent[0], name).log_mode:
					value = value.log10()
					initComponent[1].update({name: float(value)})

			self.initComponents.append(initComponent)
			
	def computeLteRadexLAModel(self, data, telescope, coms, freqs=None):
		"""docstring for fnacomputeLteRadexLa"""
		vlsrFile      = data.getVlsr()
		typeFreq      = analyzeTypeFreq(vlsrFile)
		lineFrequency = data.getListOfLines().get(0).getObsFrequency()
		if not freqs:
			freqs     = data.getFrequenciesSignal()
		try:
			lteRadexLA    = LteRadexLA()
			modelSpectrum = lteRadexLA.run(freqs=freqs,
										   lineFrequency=lineFrequency,
										   telescope=telescope,
										   tmb2ta=self.tmb2ta,
										   tbg=self.tbg,
										   aijMin=self.aijMin,
										   eup=self.eup,
										   jup=self.jup,
										   kup=self.kup,
										   lup=self.lup,
										   mup=self.mup,
										   rms=self.rmsModel,
										   typeFreq=typeFreq,
										   vlsrOffset=vlsrFile,
										   components=coms,
										   continuum=self.continuum,
										   observing_mode=self.observing_mode)
			self.is_last_check_parameters_model_OK = lteRadexLA.check_parameters
		except RadexException, re:
			sys.exit(re.getMessage())
		return modelSpectrum

		
	def computeChi2MinUsingMCMC(self, iterationCount, cutOff=None, ratioAtCutOff=None):
		"""Compute the minimum chi2 using the Markov chain Monte Carlo method.

		Usage
		-----

		computeChi2MinUsingMCMC(self, iterationCount, cutOff, ratioAtCutOff)

		Parameters
		----------

		iterationCount: the number of iterations
		cutOff : the iteration index for which the reduction factor equals 0.01
		ratioAtCutOff : the ratio of fac at cutOff
		comps: the components to be used

		"""
		printHeader("COMPUTATION OF THE CHI2 WITH MARKOV CHAIN MONTE CARLO METHOD")

		# if not isinstance(comps, (list, tuple)):
		# 	msg = "%s is not a valid list of Components. Please check your inputs.\n" % comps
		# 	msg += "Usage: computeChi2MinUsingMCMC(iterationCount, cutOff, ratioAtCutOff, *comps)"
		# 	raise TypeError(msg)

		if not isinstance(cutOff, (int, float)) or not isinstance(ratioAtCutOff, (int, float)):
			msg = "cutOff and ratioAtCutOff must be set correctly.\n"
			msg += "Usage: computeChi2MinUsingMCMC(iterationCount, cutOff, ratioAtCutOff)"
			raise TypeError(msg)

		if ratioAtCutOff > 1:
			msg = "ratioAtCutOff must be < 1"
			raise ValueError(msg)
			

		comps           = [c[0] for c in self.initComponents]
		params          = [c[1] for c in self.initComponents]
		self.components = Components(comps, isoUnique=self.isoUnique, method="MCMC")
		componentCount  = len(self.components.components)

		# Check contains a radex model and radex is installed.
		if self.components.isRadex:
			self.ensureRadexInstalled()

		acceptanceCount1 = 0 # The number of models accepted in TEST 1
		acceptanceCount2 = 0 # The number of models accepted in TEST 2
		rejectionCount1  = 0 # The number of models rejected in TEST 1
		minVal           = self.chi2Max
		lineAnalysis     = LineAnalysis()
		chi2Totals       = []
		thetas           = []
		self.rates       = []
		self.iterIndexes = []
		lines            = ""
		text             = ""
		nbPoints         = 0
		modelIndex       = 0
		validModelsFile  = open(self.outputFile, "w")
		validModelsFile.close()
		chi2Total = 0
		computationIsInitialized = False
		dof = self.components.dof

		self.continuumValues = {}
		continuumParam = get_continuum_param(self.continuum)
		self.nbOfPointsInLines = {}
		for numLineUser in self.selectedNumLines:
			data = self.inputData[numLineUser-1]
			frequencyData = data.getFrequenciesSignal()

			self.continuumValues[numLineUser] = [continuumParam.getContinuumValue(f) for f in frequencyData]
			# DETERMINE THE LINES TO BE USED
			rangesSelected = self.convertRangeToMHz(data, numLineUser, silent=0)

			self.frequenciesSelected = analyzeValuesSelected(frequencyData, rangesSelected)

			if len(self.frequenciesSelected) == 0:
				self.warnNoValueSelected(numLineUser)

			self.nbOfPointsInLines[numLineUser] = [len(self.frequenciesSelected[i]) for i in range(len(self.frequenciesSelected))]

		self.selectedLinesWithEnoughPoints = selectLinesWithEnoughPoints(self.nbOfPointsInLines, self.nbIndPoints, dof + 2)

		# WARN THE USER IF THE SET OF SELECTED LINES IS NOT WHAT IS REALLY
		# GOING TO BE USED
		# ================================================================
		if len(self.selectedLinesWithEnoughPoints) != len(self.selectedNumLines):
			rejectedLines = sorted(set(self.selectedLinesWithEnoughPoints).symmetric_difference(self.selectedNumLines))
			msg = "The lines %r don't have enough points and have been rejected." % rejectedLines
			if self.warning:
				askConfirmation(msg)
			else:
				warning(msg)

		if len(self.selectedLinesWithEnoughPoints) == 0:
			msg = "No ranges in your selection of lines has enough points"
			msg += " for a relevant computation of the minimal chi2.\n"
			msg += "Please correct your inputs (the number of parameters varying, the frequency ranges...)"
			raise ValueError(msg)

		print("Number of lines selected with enough points: %d\n" % len(self.selectedLinesWithEnoughPoints))

		# TODO Implement computation of the number of models taking into
		# account the iso

		# Display the number of models
		# ----------------------------
		if self.warning:
			msg = "Computation of %d models using %d variable parameters." % (iterationCount, dof)
			askConfirmation(msg, "Markov Chain Monte Carlo method")

		# TODO Refactor

		# Create the initial positions and the bounds of theta
		thetaBounds = map(makeBounds, comps)
		thetas_n = map(lambda c, p: initializeParameters(c, p), comps, params)
		self.is_last_check_parameters_model_OK = True
		for k in range(iterationCount):

			if cutOff and ratioAtCutOff:
				fac = calculateFacWithCutOff(k, cutOff, ratioAtCutOff)
			else:
				msg = "You must provide both cutOff and ratioAtCutOff."
				raise ValueError(msg)

			# DISPLAY THE NUMBER OF THE MODEL
			# ===============================
			modelIndexes = getModelIndexes(iterationCount)
			if modelIndex in modelIndexes:
				print("* model %d/%d |" % (modelIndex + 1, iterationCount)),

			if computationIsInitialized:
				# print(thetas_r[0]["nmol")]
				# TODO Refactor
				logModeForParameters = {}
				for i in range(componentCount):
					component = self.components.components[i]
					names = component.names

					for name in names:
						logModeForParameters[name] = getattr(component, name).log_mode

						reductionFactor = analyzeReducePhysicalParam(component.reducePhysicalParam, name)
						alpha = calculateAlpha(fac, thetaBounds[i][name], reductionFactor)

						v = Random().nextDouble()
						thetas_n[i][name] = thetas_r[i][name] + alpha * v - alpha / 2

					# If isoUnique, take only the iso paramater of the first component into account
					if self.isoUnique and i != 0 and hasattr(component, "iso"):
						thetas_n[i]["iso"] = thetas_n[0]["iso"]
			else:
				if self.isoUnique and hasattr(self.components.components[0], "iso"):
					for i in range(componentCount):
						thetas_n[i]["iso"] = thetas_n[0]["iso"]

			# TEST 1
			# =========================================================================

			if testThetaIsInRange(thetas_n, thetaBounds):
				acceptanceCount1 += 1
				chi2Val           = 0

				# LOOPS OVER THE LINES {{{2
				# =====================================================================

				self.chi2Vals = []
				self.nbOfPointsOfChi2s = []

				for numLineUser in self.selectedLinesWithEnoughPoints:
					numLine = numLineUser - 1
					data = self.inputData[numLine]
					frequencyData = data.getFrequenciesSignal()
					rangesSelected = self.convertRangeToMHz(data, numLineUser, silent=1)

					# Determine the telescope for the numLine
					telescope = determineTelescope(self.telescope, numLineUser)
					coms      = createCassisComponents(self.moltags, self.components, thetas_n)

					#######################################################

					# COMPUTATION OF THE MODEL
					# ========================
					modelSpectrum = self.computeLteRadexLAModel(data, telescope, coms)
					if not self.is_last_check_parameters_model_OK:
						return 
					# Check that the spectrum computed by Radex is not empty
					# If it is empty, continue to the next numLineUser
					if modelSpectrum == None:
						if numLineUser not in self.linesNotComputed:
							self.linesNotComputed.append(numLineUser)
						continue

					# COMPARE DATA AND MODEL
					# ======================
					intensityData  = self.checkData(data)
					intensityModel = self.checkModel(modelSpectrum)
					intensityData, intensityModel = self.rescaleModel(intensityData, intensityModel)
					frequencyModel = modelSpectrum.getFrequenciesSignal()

					# COMPUTE THE CHI2
					# ================
					chi2Val, nbOfPointsOfChi2 = computeChi2(data,
												  modelSpectrum,
												  rangesSelected,
												  self.rmsLines[numLineUser],
												  self.calLines[numLineUser],
												  self.weightLines[numLineUser],
												  self.continuumValues[numLineUser])
					# 22-APR-2013
					# Emmanuel: Weight the chi2 values by divising the chi2 by the number of
					# points used (in MCMC only)
					# chi2Val = chi2Val / nbOfPointsOfChi2
					self.chi2Vals.append(chi2Val)
					self.nbOfPointsOfChi2s.append(nbOfPointsOfChi2)

					# TODO Update the selected lines so they take into account the
					# lines that can't be computed by LteRadexLA
					# if len(self.linesNotComputed) != 0:
					# 	self.selectedLinesWithEnoughPoints = [l for l in
					# 			self.selectedLinesWithEnoughPoints if l not in
					# 			self.linesNotComputed]
					# TODO Implement updateLines
					# 	self.selectedLinesWithEnoughPoints = updateLines(self.selectedLinesWithEnoughPoints, self.linesNotComputed)

				# Sum all the chi2 values and the number of points of chi2 for all the lines
				chi2Total = sum(self.chi2Vals)
				# Check that Radex computed a valid chi2Min:
				# - chi2Min must be a positive floating number
				# if chi2Val == 0:
				if not isinstance(chi2Total, float) or chi2Total < 0:
					msg = "chi2 value = %r is not a valid value. Exiting..." % chi2Total
					raise ValueError(msg)

				if chi2Total < self.chi2Min:
					self.chi2Min = chi2Total
					self.chi2Mins = self.chi2Vals[:]
					if computationIsInitialized:
						self.thetaOptimal = copy.deepcopy(thetas_n)

				nbPoints += 1

				# Initialization of the computation when TEST 1 passed for the
				# first time
				if not computationIsInitialized:
					computationIsInitialized = True
					chi2Val_r = chi2Total
					thetas_r = copy.deepcopy(thetas_n)
					print("\nInitial thetas:")
					cpt_compo_species = 0
					for i in range(componentCount):
						names = self.components.components[i].names
						
						cpt_compo_species = cpt_compo_species + 1
						print("- Component %d:" % (cpt_compo_species)),
						initialParameters = ["nmol = %.2g" % thetas_r[i]["nmol"]]
						initialParameters += ["%s = %.2f" % (name, thetas_r[i][name]) for name in names if name != "nmol"]
						initialParameters = " %s" % (", ".join(initialParameters))
						print(initialParameters)

					print("Initial chi2  : %r" % chi2Val_r)
				else:
					# TEST 2
					# =================================================================
					# Test d'acceptation nouveau parametre non optimum
					if testChi2ValIsBetter(chi2Val_r, chi2Total):
						chi2Val_r = chi2Total
						thetas_r = copy.deepcopy(thetas_n)
						acceptanceCount2 += 1
						acceptanceRate = float(acceptanceCount2)/float(acceptanceCount1)
						self.rates.append(acceptanceRate)
						self.iterIndexes.append(k)

						# WRITE THE THETAS
						# ================
						self.algorithm = "MCMC"
						columnNames = createColumnNames(self.components, self.algorithm)
						if acceptanceCount2 == 1:
							validModelsFile  = open(self.outputFile, "a")
							listmodels = ""
							for cpm in self.components.components:
								listmodels = listmodels + cpm.model + "\t"
							#header = "chi2Min=%g\tChi2MinReduced=%g\tDof=%g\tNbPoints=%d\tPlotTitle%s=%s\n" % (self.chi2Min, 0, 1, acceptanceCount2, testjm, "Density of models with MCMC")
							header = listmodels + "\n"
							validModelsFile.write(header + columnNames)
						else:
							validModelsFile  = open(self.outputFile, "a")

						rowOfData = "%d\t" % k
						rowOfData += createRowOfData(thetas_n, logModeForParameters)
						rowOfData += "%f\t" % chi2Val_r
						rowOfData += "%f\n" % acceptanceRate
						validModelsFile.write(rowOfData)
						validModelsFile.close()
			else:
				rejectionCount1 += 1


			acceptanceRate = float(acceptanceCount2)/float(acceptanceCount1)

			# The minimal chi2 determines the best parameters
			if chi2Val_r <= minVal:
				minVal = chi2Val_r

			if modelIndex in modelIndexes:
				print("acceptance rate = %.2f |" % acceptanceRate),
				print("rejected models = %d |" % rejectionCount1),
				print("Nb of points = %d |" % sum(self.nbOfPointsOfChi2s)),
				print("chi2Min = %.2f |" % self.chi2Min),
				print("reduced chi2Min = %.2f" % computeChi2Reduced(self.chi2Mins, self.nbOfPointsOfChi2s, dof, self.nbIndPoints, self.weightLines))

			modelIndex += 1

		nbPointsUsed = sum(self.nbOfPointsOfChi2s) / self.nbIndPoints

		# Prevent a division by 0
		if self.components.dof == nbPointsUsed:
			msg = "Can't compute the reduced minimum chi2.\n"
			msg += "Please use a wider range and/or give more lines for the chi2 computation."
			raise ZeroDivisionError(msg)

		self.chi2MinReduced = computeChi2Reduced(self.chi2Mins, self.nbOfPointsOfChi2s, dof, self.nbIndPoints, self.weightLines)


		# DETERMINATION OF THE BEST PARAMETERS
		# ====================================

		if self.thetaOptimal:
			self.bestPhysicalModels = createPhysicalModels(self.thetaOptimal, self.components)
		else:
			msg = "MCMC couldn't converge on a solution because:\n"
			msg += "chi2Total > chi2Min    (%e > %e)\n" % (chi2Total, self.chi2Min)
			msg += "Please run the script again.\n"
			sys.exit(msg)

		# PRINT REPORT
		# ============
		print("\nCompleted.")
		linesNotTakenIntoAccount = set(self.selectedLinesWithEnoughPoints).symmetric_difference(set(self.selectedNumLines))
		linesNotTakenIntoAccount = list(linesNotTakenIntoAccount)

		if len(linesNotTakenIntoAccount) > 0:
			printHeader("OUTPUT")
			print("!! WARNING !!: The following lines were not taken into account in the computation of the chi2:"),
			print("%r" % linesNotTakenIntoAccount)
			print()
		else:
			printHeader("OUTPUT")
			print()

		moleculeNames = []
		for tag in self.moltags:
			moleculeNames.append(AccessDataBase.getDataBaseConnection().getMolName(str(tag)))
		#moleculeNames = map(getMolName, self.moltags)
		print("Species: %s" % "; ".join(str(molName) for molName in moleculeNames))
		print("Initial number of models: %d" % iterationCount)
		print("Number of lines used: %d" % len(self.nbOfPointsOfChi2s))
		print("Acceptance rate         : %.2f (%d/%d)" % (acceptanceRate, acceptanceCount2, acceptanceCount1))
		print("Minimal chi2            : %.2f" % self.chi2Min)
		print("Reduced minimal chi2    : %.2f" % self.chi2MinReduced)
		print("Degree of freedom       : %d" % dof)
		print("Best parameters:")

		bp = 1
		for physicalModel in self.bestPhysicalModels:
			parameters = ["nmol = %.2g" % physicalModel.nmol]
			parameters += ["%s = %.2f" % (name, getattr(physicalModel, name)) for name in physicalModel.names if name != "nmol"]
			print("- Component %d: " % bp,)
			bp += 1
			print("%s" % ", ".join(parameters))
			

		print("\nExtra information:")
		print("-" * 20)
		print("* Lines used: %r" % self.selectedLinesWithEnoughPoints)
		
		# Warn the user that some lines couldn't be computed
		if len(self.linesNotComputed) != 0:
			msg = "Warning. Lines %s skipped because they are not in the collisional file."
			print(msg % ", ".join(str(x) for x in self.linesNotComputed))

		print("* Minimal chi2s per line: %r" % [float("%.2f" % chi2Min) for chi2Min in self.chi2Mins])
		print("* List of number of points for chi2s: %r" % self.nbOfPointsOfChi2s)

		# STATISTICS
		# ==========
		reader = FileReader(self.outputFile)
		reader.read_file(lines_to_skip=1)
		chi2s = reader.columns["chi2"]
		print("=" * 80)
		print()
		printStats(["chi2"], self.outputFile, sigmaClip=3)
		print()

	def ensureRadexInstalled(self):
		radexConf = RadexConfiguration.getInstance()
		if radexConf.isDummy():
			raise ValueError("RADEX is not available in this version of CASSIS.")
		if not radexConf.haveFortran():
			raise ValueError("RADEX is not installed. Please install it first.")

	# TODO Improve the collisioners ("nH2", ...)
	def computeChi2MinUsingRG(self, *comps):
		"""Compute the chi2 min and write the data into a file

		Inputs:
			* parameters

		Returns:
			None

		"""
		printHeader("COMPUTATION OF THE CHI2 WITH A REGULAR GRID")

		lineAnalysis = LineAnalysis()
		nbPoints     = 0
		text         = ""
		modelIndex   = 0

		# Create the components from the parameters
		# -----------------------------------------
		self.components = Components(comps, isoUnique=self.isoUnique, method="RG")
		dof             = self.components.dof
		componentCount  = len(self.components.components)
		modelTotal      = len(self.components.models)

		# Check contains a radex model and radex is installed.
		if self.components.isRadex:
			self.ensureRadexInstalled()

		# Display the number of models
		# ----------------------------
		if self.warning:
			nbOfModels = computeNbOfModelsInRG(comps, self.isoUnique)
			msg = "Computation of %d models using %d variable parameters." % (nbOfModels, dof)
			askConfirmation(msg, "Markov Chain Monte Carlo method")

		vlsrParameter = self.components.components[0].vlsr
		vlsrFile      = self.inputData[0].getVlsr()
		freqRef       = self.inputData[0].getFreqRef()
		vlsrAverage   = float(calculateAverage(vlsrParameter.values))
		fScale		  = FrequencyScale(vlsrFile, freqRef)

		print("Average vlsr: %.2f km/s\n" % vlsrAverage)

		#continuum = ContinuumParameters(self.components.components[0].continuum)
		self.continuumValues = {}
		continuumParam = get_continuum_param(self.continuum)
		self.nbOfPointsInLines = {}
		for numLineUser in self.selectedNumLines:
			numLine = numLineUser - 1
			data = self.inputData[numLine]
			frequencyData = data.getFrequenciesSignal()

            # COMPUTE THE CONTINUUM
            # =====================
			self.continuumValues[numLineUser] = [continuumParam.getContinuumValue(f) for f in frequencyData]

			# DETERMINE THE LINES TO BE USED
			# ==============================
			rangesSelected = self.convertRangeToMHz(data, numLineUser, silent=1)
			self.frequenciesSelected = analyzeValuesSelected(frequencyData, rangesSelected)
			self.nbOfPointsInLines[numLineUser] = [len(self.frequenciesSelected[i]) for i in range(len(self.frequenciesSelected))]

		self.selectedLinesWithEnoughPoints = selectLinesWithEnoughPoints(self.nbOfPointsInLines, self.nbIndPoints, dof + 2)

		if len(self.selectedLinesWithEnoughPoints) == 0:
			msg = "No ranges in your selection of lines has enough points"
			msg += " for a relevant computation of the minimal chi2.\n"
			msg += "Please correct your inputs (the number of parameters varying, the frequency ranges...)"
			raise ValueError(msg)

		self.is_last_check_parameters_model_OK = True
		# Loop through the models of components
		# -------------------------------------
		for cpm in self.components.models:
			chi2Val	 = 0
			self.chi2Vals = []
			self.nbOfPointsOfChi2s= []

			modelIndexes = getModelIndexes(len(self.components.models))
			if modelIndex in modelIndexes:
				print("* model %d/%d |" % (modelIndex+1, modelTotal))

			# Loop through the indexes of the lines
			# -------------------------------------

			for numLineUser in self.selectedLinesWithEnoughPoints:
				numLine = numLineUser - 1
				data = self.inputData[numLine]
				vlsrFile = data.getVlsr()

				# Case 1: Ranges provided for the current line
				# --------------------------------------------
				#
				# Use the range provided by the user
				rangesSelected = self.convertRangeToMHz(data, numLineUser, silent=modelIndex)

				# Determine the telescope for the numLine
				telescope = determineTelescope(self.telescope, numLineUser)
				# Create all the components
				coms = map(lambda cp: createCassisComponent(self.moltags, cp), cpm)

				# Compute the model
				modelSpectrum = self.computeLteRadexLAModel(data, telescope, coms)
				if not self.is_last_check_parameters_model_OK:
						return 

				# Check that the spectrum computed by Radex is not empty
				# If it is empty, continue to the next numLineUser
				if modelSpectrum == None:
					if numLineUser not in self.linesNotComputed:
						self.linesNotComputed.append(numLineUser)
					continue

				# COMPARE DATA AND MODEL
				# ======================
				intensityData  = self.checkData(data)
				intensityModel = self.checkModel(modelSpectrum)
				intensityData, intensityModel = self.rescaleModel(intensityData, intensityModel)

				# COMPUTE THE CHI2 MIN
				# ====================
				chi2Val, nbOfPointsOfChi2 = computeChi2(data,
											  modelSpectrum,
											  rangesSelected,
											  self.rmsLines[numLineUser],
											  self.calLines[numLineUser],
											  self.weightLines[numLineUser],
											  self.continuumValues[numLineUser])
				self.chi2Vals.append(chi2Val)
				self.nbOfPointsOfChi2s.append(nbOfPointsOfChi2)

			# Sum all the chi2 values and the number of points of chi2 for all the lines
			chi2Total = sum(self.chi2Vals)

			# Check that Radex computed a valid chi2Min:
			# - chi2Min must be a positive floating number
			# if chi2Val == 0:
			# 	print("Radex computed a chi2 value = 0.")
			if not isinstance(chi2Total, float) or chi2Total < 0:
				msg = "chi2 value = %r is not a valid value. Exiting..." % chi2Total
				raise ValueError(msg)

			# Take the minimum chi2 of sum(chi2_i) where chi2_i is the the chi2 for
			# a spectrum
			if chi2Total < self.chi2Min:
				self.bestPhysicalModels = cpm
				self.chi2Min = chi2Total
				self.chi2Mins = self.chi2Vals[:]
				print("chi2Min = %f\tmodel %d/%d\t%r" % (self.chi2Min, modelIndex+1, modelTotal, cpm))

			# PREPARE THE DATA FOR WRITING
			# ============================

			paramValues = []
			for i in range(componentCount):
				paramNames = self.components.components[i].names
				paramValues.extend(self.components.extractParametersValues(i, nbPoints, paramNames))

			# if model == 0:
			paramValues.insert(0, nbPoints)
			paramValues.append(chi2Total)
			lines = ["%g" % x for x in paramValues]
			text += "\t".join(lines)
			text += "\n"

			if modelIndex in modelIndexes:
				print("Nb of points = %d |" % sum(self.nbOfPointsOfChi2s)),
				print("chi2Min = %.2f" % self.chi2Min)

			nbPoints += 1
			modelIndex += 1


		# TODO Refactor: Write a function
		# Write data
		# ----------

		self.algorithm  = "RG"
		columnNames     = createColumnNames(self.components, self.algorithm)
		validModelsFile = open(self.outputFile, 'w')
		dof             = self.components.dof
		nbPointsUsed    = sum(self.nbOfPointsOfChi2s) / self.nbIndPoints

		# Prevent a division by 0
		if dof == nbPointsUsed:
			msg = "Can't compute the reduced minimum chi2.\n"
			msg += "Please use a wider range and/or give more lines for the chi2 computation."
			raise ZeroDivisionError(msg)

		self.chi2MinReduced = computeChi2Reduced(self.chi2Mins, self.nbOfPointsOfChi2s, dof, self.nbIndPoints, self.weightLines)
		header = "chi2Min=%g\tChi2MinReduced=%g\tDof=%g\tNbPoints=%d\tPlotTitle=%s\n" % (self.chi2Min, self.chi2MinReduced, dof, nbPointsUsed, self.plotTitle)
		validModelsFile.write(header + columnNames)
		validModelsFile.write(text)
		validModelsFile.close()

		# PRINT REPORT
		# ============
		print("\nCompleted.")
		linesNotTakenIntoAccount = set(self.selectedLinesWithEnoughPoints).symmetric_difference(set(self.selectedNumLines))
		linesNotTakenIntoAccount = list(linesNotTakenIntoAccount)

		if len(linesNotTakenIntoAccount) > 0:
			printHeader("OUTPUT")
			print("!! WARNING !!: The following lines were not taken into account in the computation of the chi2:",)
			print("%r" % linesNotTakenIntoAccount)
			print()
		else:
			printHeader("OUTPUT")
			print()

		moleculeNames = []
		for tag in self.moltags:
			moleculeNames.append(AccessDataBase.getDataBaseConnection().getMolName(str(tag)))
		#moleculeNames = map(getMolName, self.moltags)
		print("Species: %s" % "; ".join(str(molName) for molName in moleculeNames))
		print("Nb of models = %d" % modelTotal)
		print("Minimal chi2 = %f" % self.chi2Min)
		print("Reduced chi2 min = %f" % self.chi2MinReduced)
		print("Degree of freedom       : %d" % dof)

		if len(self.bestPhysicalModels) == 0:
			msg = "No chi2 < %g was found. Check your inputs." % self.chi2Min
			raise ValueError(msg)

		i = 0
		for bestParameter in self.bestPhysicalModels:
			print("Best parameters %d: %r" % (i+1, bestParameter))
			i += 1

		totalNbOfPointsInData = sum(map(lambda i: self.inputData[i].getSize(), range(len(self.inputData))))
		print("Total nb of points used in the computation = %d / %d" % (nbPointsUsed, totalNbOfPointsInData))
		print("-" * 80)

	def _handleBestModel(self, moltag, overSampling, tuningBand, telescope="", doesPlot=True):
		"""_handleBestModel is a wrapper for easily dealing with either

		creation of the best model (in the plugin script) or
		creation and plot of the best model (in the standalone script)

		The creation and plot the model is done with doesPlot=True
		The creation of the model only is done with doesPlot=False

		"""
		print("Plotting the best model..."),

		if isinstance(self.telescope, dict):
			msg = "Please choose the telescope between:\n%s" % "\n".join(self.telescope.values())
			msg += "\nand relaunch plotBestModel."

		# telescope is not provided
		if telescope == "":
			# Case 1: Many telescopes
			if isinstance(self.telescope, dict):

				# Take the only telescope if it is provided in the user inputs
				if len(self.telescope) == 1 or is_uniq(self.telescope.values()):
					telescope = self.telescope.values()[0]
				else:
					raise ValueError("[telescope] " + msg)

			# Case 2: A single telescope
			if isinstance(self.telescope, str) or isinstance(self.telescope, unicode):
				telescope = self.telescope

		# telescope is provided
		else:
			# Case 1: Many telescopes
			# Check the telescope provided is a telescope provide in user inputs
			if isinstance(self.telescope, dict):
				if telescope not in self.telescope.values():
					msg = "[telescope] Wrong telescope! " + msg
					raise ValueError(msg)

			# Case 2: A single telescope
			if isinstance(self.telescope, str) or isinstance(self.telescope, unicode):
				if telescope != self.telescope:
					msg = "[telescope] You used %s in the chi2 computation not %s.\n" % (self.telescope, telescope)
					msg += "Please use the right telescope or don't enter the telescope in plotBestModel."
					raise ValueError(msg)

		coms = []
		for bestParameter in self.bestPhysicalModels:
			coms.append(createCassisComponent(self.moltags, bestParameter))

		# TODO Fix memory leak in Java part
		fname, extensionName = os.path.splitext(self.inputFile)
		if (extensionName == '.fits' or extensionName == '.fts') and (self.indice_hdu != -1):
			fileManager = FileManagerFits(File(self.inputFile), False)
			print("indice_hdu = " + str(self.indice_hdu))
			cassisSpectrum = fileManager.read(self.indice_hdu)
		else:
		 	cassisSpectrum = FileReaderCassis.createCassisSpectrumFromFile(File(self.inputFile))
		
		cassisSpectrum.setTitle(self.inputFile)
		res, model = self.lineAnalysis.run(cassisSpectrum = cassisSpectrum,
										   tuningRange=self.tuningRange,
										   tuningBand=tuningBand,
										   plottingBand="Signal",
										   plottingSort="Frequency",
										   species=[self.template, moltag],
										   eup=self.eup,
										   jup=self.jup,
										   kup=self.kup,
										   lup=self.lup,
										   mup=self.mup,
										   aijMin=self.aijMin,
										   telescope=telescope,
										   tmb2ta=self.tmb2ta,
										   tbg=self.tbg,
										   overSampling=overSampling,
										   rms=self.rmsModel,
										   components = coms,
										   continuum=self.continuum,
										   observing_mode = self.observing_mode,
										   warning=False)

		print("OK")
		print("OverSampling: %d" % overSampling)
		print("Telescope: %s" % telescope)

		# Trying to free some memory
		del(cassisSpectrum)

		System.gc()

		if doesPlot:
			self.lineAnalysis.plot()
			return model
		else:
			return model

	def displayLines(self, indexes):
		"""Display lines with given indexes"""
		if not isinstance(indexes, (list, tuple)):
			indexes = [indexes-1]
		else:
			indexes = [x-1 for x in indexes]

		p = PanelView.getInstance().getLineSpectrumView().getControl().getModel().getGalleryNavigationModel()
		p.setFilter(indexes, True)

	def plotBestModel(self, moltag=None, tuningBand=None, overSampling=1, telescope="", lines=None):
		"""Plot the best model and return the LineAnalysisModel

		Usage
		-----
		lineAnalysisModel = plotBestModel(moltag, tuningBand=tuningBand, overSampling=overSampling, telescope=telescope, lines=lines)

		Examples
		--------
		# Plot the best model for the lines 1
		lineAnalysisModel = plotBestModel(moltag, lines="1")
		# or
		lineAnalysisModel = plotBestModel(moltag, lines=1)

		# Plot the best model for the lines 1 and 3
		lineAnalysisModel = plotBestModel(moltag, lines="1,3")

		# Plot the best model for the lines 1 to 3
		lineAnalysisModel = plotBestModel(moltag, lines="1-3")

		Parameters
		----------
		moltag : int or list of int
			the molecular(s) tag

		tuningBand : int (optional)
			the tuning band

		overSampling : int (optional)
			the oversampling factor

		telescope : str (optional)
			the telescope

		lines : int or str (optional)
			The lines to be displayed in LineAnalysis (start at 1)

		Returns
		-------
		lineAnalysisModel : LineAnalysisModel

		"""
		# If only moltag is provided, the script automatically takes its value
		if len(self.moltags) == 1:
			moltag = self.moltags[0]
		else:
			# If there are no moltag, the user MUST choose one
			if (not (isinstance(moltag, int)  or (isinstance(moltag, list) and len(moltag) >= 1))):
				msg = "You must provide one or more moltag. Please choose between % r." % self.moltags
				raise ValueError(msg)

		if isinstance(moltag, int):
			if moltag not in self.moltags:
				raise ValueError("You must provide a moltag corresponding to your input moltags %s" % self.moltags)
		elif isinstance(moltag, list):
			for tag in moltag:
				if tag not in self.moltags:
					raise ValueError("You must provide a moltag corresponding to your input moltags %s" % self.moltags)

		if tuningBand is None:
			tuningBand = self.tuningBand

		# Create, plot and trigger
		lineAnalysisModel = self._handleBestModel(moltag, overSampling, tuningBand, telescope, doesPlot=True)

		# Select the lines to be displayed
		if lines:
			if isinstance(lines, int):
				lines = str(lines)
			self.displayLines(makeRange(lines))
		return lineAnalysisModel

	def plotChi2Contour(self):
		print("The Chi2Contour functionality was removed.")

	def _createLisFile(self, filename, columns, vlsr):
		"""Create a .lis file from the columns list

		Inputs:
			columns - a dictionary corresponding to
			column name: column value

		"""
		pointCount = len(columns["FreqLsb"])

		if os.path.exists(filename):
			f = open(filename, "w")

		f = open(filename, "a")
		f.write("// number of line : %d\n" % pointCount)
		f.write("// vlsr : %.2f\n" % vlsr)
		f.write("FreqLsb VeloLsb FreqUsb VeloUsb Int DeltaF DeltaV\n")

		for i in range(pointCount):
			f.write("%f\t%f\t%3.2f\t%3.2f\t%f\t%f\t%f\t\n" % (columns["FreqLsb"][i],
														columns["VeloLsb"][i],
														columns["FreqUsb"][i],
														columns["VeloUsb"][i],
														columns["Int"][i],
														columns["DeltaF"][i],
														columns["DeltaV"][i]))
		f.close()

	def _createColumnsForInputData(self):
		frequencySignal = []
		frequencyImage  = []
		intensity       = []
		velocitySignal  = []
		velocityImage   = []
		deltaF          = []
		deltaV          = []

		for i in range(len(self.inputData)):
			data = self.inputData[i]
			print("Concatenating spectrum %d..." % (i + 1)),
			frequencySignal.extend(data.getFrequenciesSignal())
			frequencyImage.extend([0.0] * len(data.getFrequenciesSignal()))
			intensity.extend(data.getIntensities())
			velocitySignal.extend(data.getVelocitySignal())
			velocitySignal.reverse()
			print("(%d points)" % data.getSize())
			velocityImage.extend([0.0] * len(data.getVelocitySignal()))
			deltaF.extend(data.getDeltaF())
			deltaV.extend(data.getDeltaV())

		columns = {"FreqLsb": frequencySignal,
				"VeloLsb": velocitySignal,
				"FreqUsb": frequencyImage,
				"VeloUsb": velocityImage,
				"Int"    : intensity,
				"DeltaF" : deltaF,
				"DeltaV" : deltaV}
		return columns

	def concatenateLisFiles(self, filename):
		vlsr = self.inputData[0].vlsr
		columns = self._createColumnsForInputData()
		self._createLisFile(filename, columns, vlsr)

	def getFrequencyForSpectrum(self, numLineUser, overSampling=1):
		"""docstring for getFrequencyForSpectrum"""
		frequency = self.inputData[numLineUser-1].getFrequenciesSignal()
		if overSampling != 1:
			frequency = UtilArrayList.oversample(frequency, overSampling)
		return frequency

	def getVelocityForSpectrum(self, numLineUser, overSampling=1):
		spectrum = self.inputData[numLineUser-1]

		# Convert frequency to velocity
		lineRestFrequency = spectrum.getListOfLines().get(0).getObsFrequency()
		xaxis = XAxisVelocity.getXAxisVelocity()
		xaxis.setCoeff(spectrum.getVlsr(), lineRestFrequency)

		velocity = map(xaxis.convertFromMhzFreq, spectrum.getFrequenciesSignal())
		if overSampling != 1:
			velocity = UtilArrayList.oversample(velocity, overSampling)

		return velocity

	def overSampleModel(self, overSampling):
		resampledFrequencies = []
		resampledVelocities  = []
		resampledIntensities = []

		# The values of the best parameters are at the power of 10
		# we must apply log10
		parameterValues = [x.parameters for x in self.bestPhysicalModels]
		logModes = [getAttrOfParameters("log_mode", c) for c in self.components.components]
		newParameterValues = updateParameterValues(parameterValues, logModes)
		#
		for numLineUser in self.selectedLinesWithEnoughPoints:
			data = self.inputData[numLineUser-1]
			telescope = determineTelescope(self.telescope, numLineUser)
			#
			coms = createCassisComponents(self.moltags, self.components, newParameterValues)
			resampledFrequency = self.getFrequencyForSpectrum(numLineUser, overSampling)
			resampledVelocity  = self.getVelocityForSpectrum(numLineUser, overSampling)
			modelSpectrum      = self.computeLteRadexLAModel(data, telescope, coms, freqs=resampledFrequency)
			resampledIntensity = modelSpectrum.getIntensities()

			resampledFrequencies.append(resampledFrequency)
			resampledVelocities.append(resampledVelocity)
			resampledIntensities.append(resampledIntensity)
		return (resampledFrequencies, resampledVelocities, resampledIntensities)
		
	def writeLineSpectrumFile(self, filename, frequency, velocity, intensity):
		# Assumption: All the vlsrs of the spectra are identical
		vlsr = self.inputData[0].getVlsr()
		lineSpectrumData = LineSpectrumData(frequency, velocity, intensity, vlsr)
		lineSpectrumData.saveToFile(filename)

	def getBestPhysicalModels(self):
		"""Return the best physical models for the input components

		"""
		return self.bestPhysicalModels

	def saveBestPhysicalModels(self, filename, overSampling=1):
		"""Save the spectra corresponding to the best model to a file.
		
		Usage
		----
		saveBestPhysicalModels(filename, overSampling)

		Parameters
		----------
		filename : str
			The filename 

		overSampling : int
			The oversampling factor

		Returns
		-------
		None

		"""
		basename, extension = os.path.splitext(filename)

		if extension not in (".lis", ".fus"):
			filename += ".lis"

		print("Oversampling = %d" % overSampling)
		frequency, velocity, intensity = self.overSampleModel(overSampling)
		self.writeLineSpectrumFile(filename, frequency, velocity, intensity)


# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

# vim: set noet fdm=marker:
