#=====================================================================================================================================
#Copyright
#=====================================================================================================================================

#Copyright (C) 2014 Alexander Blaessle, Patrick Mueller and the Friedrich Miescher Laboratory of the Max Planck Society
#This software is distributed under the terms of the GNU General Public License.

#This file is part of PyFRAP.

#PyFRAP is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.

#This program 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 General Public License for more details.

#You should have received a copy of the GNU General Public License
#along with this program.  If not, see <http://www.gnu.org/licenses/>.

#===========================================================================================================================================================================
#Module Description
#===========================================================================================================================================================================

"""Optimization module for PyFRAP toolbox.

Currently contains all functions necessary to transform a constrained FRAP optimization problem into
a unconstrained one, making it suitable to Nelder-Mead optimization algorithm. 

"""

#===========================================================================================================================================================================
#Importing necessary modules
#===========================================================================================================================================================================

#Numpy/Scipy
import numpy as np

#PyFRAP
import pyfrp_fit_module 

from pyfrp_term_module import *

#===========================================================================================================================================================================
#Module Functions
#===========================================================================================================================================================================

def constrObjFunc(x,fit,debug,ax,returnFit):
	
	"""Objective function when using Constrained Nelder-Mead.
	
	Calls :py:func:`pyfrp.modules.pyfrp_optimization_module.xTransform` to transform x into
	constrained version, then uses :py:func:`pyfrp.modules.pyfrp_fit_module.FRAPObjFunc` to 
	find SSD.
	
	Args:
		x (list): Input vector, consiting of [D,(prod),(degr)].
		fit (pyfrp.subclasses.pyfrp_fit): Fit object.
		debug (bool): Display debugging output and plots.
		ax (matplotlib.axes): Axes to display plots in.
		returnFit (bool): Return fit instead of SSD.
	
	Returns:
		 float: SSD of fit. Except ``returnFit==True``, then will return fit itself. 
	"""
	
	
	LBs, UBs = buildBoundLists(fit)
	
	x=xTransform(x,LBs,UBs)

	ssd=pyfrp_fit_module.FRAPObjFunc(x,fit,debug,ax,returnFit)
	
	return ssd

def xTransform(x,LB,UB):
	
	"""Transforms ``x`` into constrained form, obeying upper 
	bounds ``UB`` and lower bounds ``LB``.
	
	.. note:: Will add tiny offset to LB(D), to avoid singularities.
	
	Idea taken from http://www.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd--fminsearchcon
	
	Args:
		x (list): Input vector, consiting of [D,(prod),(degr)].
		LB (list): List of lower bounds for ``D,prod,degr``.
		UB (list): List of upper bounds for ``D,prod,degr``.
	
	Returns:
		list: Transformed x-values. 
	"""
	
	#Make sure everything is float
	x=np.asarray(x,dtype=np.float64)
	LB=np.asarray(LB,dtype=np.float64)
	UB=np.asarray(UB,dtype=np.float64)
	
	#Check if LB_D==0, then add a little noise to it so we do not end up with xtrans[D]==0 and later have singularities when scaling tvec
	if LB[0]==0:
		LB[0]=1E-10
	
	#Determine number of parameters to be fitted
	nparams=len(x)

	#Make empty vector
	xtrans = np.zeros(np.shape(x))
	
	# k allows some variables to be fixed, thus dropped from the
	# optimization.
	k=0

	for i in range(nparams):

		#Upper bound only
		if UB[i]!=None and LB[i]==None:
		
			xtrans[i]=UB[i]-x[k]**2
			k=k+1
			
		#Lower bound only	
		elif UB[i]==None and LB[i]!=None:
			
			xtrans[i]=LB[i]+x[k]**2
			k=k+1
		
		#Both bounds
		elif UB[i]!=None and LB[i]!=None:
			
			xtrans[i] = (np.sin(x[k])+1.)/2.*(UB[i] - LB[i]) + LB[i]
			xtrans[i] = max([LB[i],min([UB[i],xtrans[i]])])
			k=k+1
		
		#No bounds
		elif UB[i]==None and LB[i]==None:
		
			xtrans[i] = x[k]
			k=k+1
			
		#Note: The original file has here another case for fixed variable, but since we made the decision earlier which when we call frap_fitting, we don't need this here.
	
	return xtrans	
		
def transformX0(x0,LB,UB):
	
	"""Transforms ``x0`` into constrained form, obeying upper 
	bounds ``UB`` and lower bounds ``LB``.
	
	Idea taken from http://www.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd--fminsearchcon
	
	Args:
		x0 (list): Input initial vector, consiting of [D,(prod),(degr)].
		LB (list): List of lower bounds for ``D,prod,degr``.
		UB (list): List of upper bounds for ``D,prod,degr``.
	
	Returns:
		list: Transformed x-values. 
	"""
	
	x0u = list(x0)
	
	nparams=len(x0)
	
	k=0
	for i in range(nparams):
		
		#Upper bound only
		if UB[i]!=None and LB[i]==None:
			if UB[i]<=x0[i]:
				x0u[k]=0
			else:
				x0u[k]=sqrt(UB[i]-x0[i])	
			k=k+1
			
		#Lower bound only
		elif UB[i]==None and LB[i]!=None:
			if LB[i]>=x0[i]:
				x0u[k]=0
			else:
				x0u[k]=np.sqrt(x0[i]-LB[i])	
			k=k+1
		
		
		#Both bounds
		elif UB[i]!=None and LB[i]!=None:
			if UB[i]<=x0[i]:
				x0u[k]=np.pi/2
			elif LB[i]>=x0[i]:
				x0u[k]=-np.pi/2
			else:
				x0u[k] = 2*(x0[i] - LB[i])/(UB[i]-LB[i]) - 1;
				#shift by 2*pi to avoid problems at zero in fminsearch otherwise, the initial simplex is vanishingly small
				x0u[k] = 2*np.pi+np.arcsin(max([-1,min(1,x0u[k])]));
			k=k+1
		
		#No bounds
		elif UB[i]==None and LB[i]==None:
			x0u[k] = x[i]
			k=k+1
	
	return x0u

def buildBoundLists(fit):
	
	"""Builds list of lower bounds and upper bounds.
	
	Args:
		fit (pyfrp.subclasses.pyfrp_fit): Fit object.
		
	Returns:
		tuple: Tuple containing: 
		
			* LBs (list): List of lower bounds.
			* UBs (list): List of upper bounds.
			
	
	
	"""
	
	LBs=[fit.LBD]+int(fit.fitProd)*[fit.LBProd]+int(fit.fitDegr)*[fit.LBDegr]+len(fit.ROIsFitted)*[fit.LBEqu]
	UBs=[fit.UBD]+int(fit.fitProd)*[fit.UBProd]+int(fit.fitDegr)*[fit.UBDegr]+len(fit.ROIsFitted)*[fit.UBEqu]
	
	return LBs,UBs