--GetLineWidthsFromPeakLists
-- Written by F.Damberger Aug 31, 2014 (version 1)
-- See section below which must be edited by user for particular case
--
-- Writes out file in csv format one line per peak: "peaknumber, linewidth" (along the selected dimension)
-- if a BatchList is defined for the selected peaklist (with a set of spectra) then the linewidth for the peak in each spectrum
-- for the peak will be written out on the line in the order specified in the BatchList
-- Script should also work if no peaklist is defined for the peaklist. It will write out one linewidth per line.
--
-- script finds maximum along selected DimTrace and first minima in each direction
-- upfield and downfield, LW=0.5 * ( HalfLineWidth(Upfield) + HalfLineWidth(Downfield) )
-- unless assymetry > MaxAssymetry, then it estimates
-- LW = 2*min(HalfLineWidth(Upfield),HalfLineWidth(Downfield))
-- Peaks with asymmetry larger than MaxAssymetry are reported in the terminal window

--version 1
-- written by F. Damberger, 03-Sept-2014
-- version 2
-- F. Damberger 25-April-2016
-- fixed error caused because table t was not defined
-- added the capability to write out linewidths for a series of spectra defined by a batchlist

-- !!!!!! USER edit this section for your use case !!!!!! -----------------------
-- edit the USER variables below for your application ---------------------------
ProjectName="MyProjectName"
SpecId=15 -- enter SpectrumId of Spectrum with PeakList that has BatchList defined (also works with no BatchList)
PeakListId=1 -- edit to be equal to ID of PeakList with BatchList defined
DimTrace=2 -- 1= F1(direct dim), 2= F2 (indirect dim) -- direction of traces whose linewidths are evaluated
DimOrthogonalToTrace=1 -- see above, should be OTHER dimension if DimTrace=1, DimOrthogonalToTrace=2
DeltaPpm = 0.1 -- this is the range of chemical shift around the selected peak to use for lineshape determination
Unit = "Hz" -- unit of the horizontal scale: Hz or ppm

MaxAssymetry = 0.2 -- this parameter excludes residues if their Lineshape is too nonsymmetric at halfheight
                    -- Assymetry is 2*(LeftHalfWidth-RightHalfWidth/FullWidth)
					-- If a reference residue is excluded then this is reported by the script in the Terminal window

--Formatting 
sp=";" -- spacer symbol used between datapoints in outputfile

FileSuffix=".csv" -- for output file


-- Do not edit below this line --------------------------------------------------
--------------------- end of user editable section ------------------------------

AbsMinBufferSize = 2048
t = {}
-- -----------------------------------------------------------------------------
-- -----------------------------------------------------------------------------
-- function definitions: note Project,SpecId,DimTrace are used as global variables in functions
-- -----------------------------------------------------------------------------
function GetDeltaPpmMax( PpmDim1, Delta, Units, BatchList )
	DeltaPpmMax = 0
	for Index,SpecId in pairs( BatchList ) do
		Spectrum = t.P:getSpectrum( SpecId )
		startPpmDim1,endPpmDim1 = Spectrum:getPpmRange( DimTrace )
		if Units == "Hz" then
			RfFreq1 = Spectrum:getRfFreq(1)
			DeltaPpm = Delta/RfFreq1
		else
			DeltaPpm = Delta
		end
		if DeltaPpm > DeltaPpmMax then DeltaPpmMax = DeltaPpm end
	end
	return DeltaPpmMax
	
end -- function GetDeltaPpmMax
-- -----------------------------------------------------------------------------
function GetMaxIntensityInSlice( SliceBuffer, i_offset, SpecId, Unit )
	local Max = 0 -- assume slices have positive maxima
	local HalfOfMinBufferSize = tonumber( MinBufferSize*0.5 )
	local Xcenter = GetXvalue( SliceBuffer, HalfOfMinBufferSize+i_offset, SpecId, Unit )
	local RfFreq = Spectrum:getRfFreq( 1 )

	for i=1,MinBufferSize do
		Y=GetYvalue( SliceBuffer, i+i_offset, 1.0 )
		X=GetXvalue( SliceBuffer, i+i_offset, SpecId, Unit )
		if Y > Max and math.abs( Xcenter - X ) <  DeltaPpm*RfFreq*0.2 then
			Max = Y
			MaxPos = GetXvalue( SliceBuffer, i, SpecId, Unit )
			MaxIndex = i
		end
	end
	--print("Max intensity in slice is"..Max)
	return Max,MaxPos,MaxIndex
end -- GetMaxIntensityInSlice
-- ---------------------core function for script ---------------------------------------------------------------------------------------
function GetHalfMaximaAndWidth( SliceBuffer, i_offset, SpecId, Unit )
	local FullWidth=0	-- = LeftHalfWidth + RiteHalfWidth
	local Assymetry=0	-- = 2*math.abs(LeftHalfWidth - RiteHalfWidth)/FullWidth
	local MinFullWidth=0 -- = 2*Smaller(LeftHalfWidth,RiteHalfWidth)
	local MaxFullWidth=0 -- = 2*Larger(LeftHalfWidth,RiteHalfWidth)
	local LeftHalfHeight=0 -- Intensity at point closest to 0.5*Maximum Intensity (left side of maximum)
	local RiteHalfHeight=0 -- Intensity at point closest to 0.5*Maximum Intensity (right side of maximum)
	local LeftHalfMax = 0 -- Position of point at LeftHalfHeight (in units given by global variable 'Unit')
	local LeftHalfWidth = 0 -- math.abs( LeftHalfMax - MaxPos ) (MaxPos is position of maximum)
	local RiteHalfMax = 0 -- Position of point at RiteHalfHeight (in units given by global variable 'Unit')
	local RiteHalfWidth = 0 -- math.abs( RiteHalfMax - MaxPos )

	local MaxSoFar = 0
	Maximum,MaxPos,MaxIndex = GetMaxIntensityInSlice( SliceBuffer, i_offset, SpecId, Unit )
	--[[
	-- The old approach: scan across the whole SliceBuffer
	-- storing the LeftHalfMax when i < MinBufferSize and the RiteHalfMax when i >= MinBufferSize
	-- The problem with this approach is that it ignores the possibility that there are two maxima in the SliceBuffer
	-- Therefore I switched to the approach below
	]]
	-- THIS IS THE NEW APPROACH - scanning for the first HalfMax to the left and
	-- then scanning for the first HalfMax to the right
	-- find the Left HalfMax
	CurrentY=Maximum*2 -- ensure that while loop can begin
	LeftHalfMax=9999
	--LastMinusHalfMax = 1 -- ensure that product CurrentMinusHalfMax*LastMinusHalfMax is positive when the scan starts
	--for j=1,MaxIndex do -- Left Half of SliceBuffer
	j=0
	i=MaxIndex+1-j
	--NoHalfMaxFound = true
	while j<MinBufferSize and CurrentY>Maximum*0.45 and i+i_offset - 1 > 0 do -- search until intensity is less than Maximum*0.5
		j=j+1
		i=MaxIndex+1-j -- from MaxIndex to 1
		if i+i_offset <1 then print("i+i_offset = "..i+i_offset.." decrementing i+i_offset below 1!!!!!!") end
		CurrentY = GetYvalue( SliceBuffer, i+i_offset, 1.0 )
		CurrentX = GetXvalue( SliceBuffer, i+i_offset, SpecId, Unit )
		CurrentMinusHalfMax=CurrentY-Maximum*0.5
		--[[
		if CurrentMinusHalfMax*LastMinusHalfMax < 0 then
			NoHalfMaxFound = false
		end
		]]
		--if NumHalfMax < 2 then -- look for only first HalfMax (i.e. ignore other maxima in the SliceBuffer)
			if math.abs( Maximum*0.5 - CurrentY  ) < math.abs( Maximum*0.5 - LeftHalfMax ) then
				LeftHalfMax = CurrentY
				LeftHalfMaxPos = CurrentX
				LeftHalfWidth = math.abs( LeftHalfMaxPos - MaxPos )
				LeftHalfHeight = LeftHalfMax/Maximum
			end
		--end
		--LastMinusHalfMax = CurrentMinusHalfMax
	end
	-- find the Rite HalfMax
	CurrentY=Maximum*2 -- ensure that while loop can begin
	RiteHalfMax=9999
	--LastMinusHalfMax = 1 -- ensure that product CurrentMinusHalfMax*LastMinusHalfMax is positive when the scan starts
	--for i=MaxIndex,MinBufferSize do -- Rite Half of SliceBuffer
	i=MaxIndex-1
	--NoHalfMaxFound = true
	while i<MinBufferSize and CurrentY>Maximum*0.45 do
		i=i+1
		CurrentY = GetYvalue( SliceBuffer, i+i_offset, 1.0 )
		CurrentX = GetXvalue( SliceBuffer, i+i_offset, SpecId, Unit )
		CurrentMinusHalfMax=CurrentY-Maximum*0.5
		--[[
		if CurrentMinusHalfMax*LastMinusHalfMax < 0 then
			NoHalfMaxFound = false
		end
		]]
		--if NumHalfMax < 2 then -- look for only first HalfMax (i.e. ignore other maxima in the SliceBuffer)
			if math.abs( Maximum*0.5 - CurrentY  ) < math.abs( Maximum*0.5 - RiteHalfMax ) then
				RiteHalfMax = CurrentY
				RiteHalfMaxPos = CurrentX
				RiteHalfWidth = math.abs( RiteHalfMaxPos - MaxPos )
				RiteHalfHeight = RiteHalfMax/Maximum
			end
		--end
		--LastMinusHalfMax = CurrentMinusHalfMax
	end
	FullWidth = LeftHalfWidth+RiteHalfWidth
	Assymetry = math.abs(LeftHalfWidth - RiteHalfWidth)/FullWidth
	LeftHalfHeight = LeftHalfMax/Maximum
	RiteHalfHeight = RiteHalfMax/Maximum
	if LeftHalfWidth < RiteHalfWidth then
		MinFullWidth = LeftHalfWidth*2.0
		MaxFullWidth = RiteHalfWidth*2.0
	else
		MinFullWidth = RiteHalfWidth*2.0
		MaxFullWidth = LeftHalfWidth*2.0
	end
	
	return FullWidth, Assymetry, MinFullWidth, MaxFullWidth,LeftHalfHeight,RiteHalfHeight
--	return FullWidth, Assymetry, LeftHalfWidth, RiteHalfWidth,LeftHalfHeight,RiteHalfHeight
end -- function GetHalfMaximaAndWidth
-- ---------------------------------------------------------------------------------------------------------------------
function GetMinBufferSize( MinBufferSize, Slice )
	result = AbsMinBufferSize
	BufferSize = Slice:getSampleCount( 1 )
	if BufferSize < MinBufferSize then result = BufferSize end
	return result
end --GetMinBufferSize
-- -----------------------------------------------------------------------------
function GetPos( Peak, Spectrum)
	PpmDim1,PpmDim2=Peak:getPos( Spectrum )
	FormattedPpm1 = string.format( "%7.3f", PpmDim1 )
	FormattedPpm2 = string.format( "%7.3f", PpmDim2 )
	return FormattedPpm1,FormattedPpm2
	--return PpmDim1,PpmDim2
end
-- -----------------------------------------------------------------------------
function GetPpmRange( DimTrace, PpmDim1, PpmDim2, DeltaPpm, Spectrum )
	startPpmSourceDim,endPpmSourceDim = Spectrum:getPpmRange( DimTrace )
	if DimTrace == 2 then -- slice is along indirect dimension
		PpmSourceDim = PpmDim2
	else					-- slice is along direct dimension
		PpmSourceDim = PpmDim1
	end
	--print( string.format( "%7.3f", startPpmDim1 ), string.format( "%7.3f", endPpmDim1 ) )
	if PpmSourceDim + DeltaPpm < startPpmSourceDim then startPpmSourceDim = PpmSourceDim + DeltaPpm end
	if PpmSourceDim - DeltaPpm > endPpmSourceDim then endPpmSourceDim = PpmSourceDim - DeltaPpm end
	return startPpmSourceDim,endPpmSourceDim
end
-- -----------------------------------------------------------------------------

function Float( Value )
	return string.format("%3.3f", Value )
end
-- -----------------------------------------------------------------------------
function GetSlice( Spectrum, StartPpm, EndPpm, PpmOrthoDim )
	--print(Spectrum:getDimCount())
	if Verbose then
		print("function GetSlice:"..Spectrum:getId(),Float(StartPpm),Float(EndPpm),Float(PpmOrthoDim) )
	end
	if DimTrace == 2 then
		Slice = Spectrum:getSlicePpm( 2, PpmOrthoDim, PpmOrthoDim, StartPpm, EndPpm )
	else
		Slice = Spectrum:getSlicePpm( 1, StartPpm, EndPpm, PpmOrthoDim, PpmOrthoDim )
	end
	return Slice
end -- function GetSlice
-- -----------------------------------------------------------------------------
function GetXvalue( SliceBuffer, PointIndex, SpecId, Unit )
	--uses global variables Project and DimTrace
	RfFreq = Project:getSpectrum(SpecId):getRfFreq( DimTrace )
	PpmMax,PpmMin = SliceBuffer:getPpmRange( 1 )
	PpmCenter = (PpmMax + PpmMin )/2
	Value = SliceBuffer:getFreq( 1, PointIndex )
	if Unit == "Hz" then 
		Value = Value - PpmCenter
		Value = Value * RfFreq
	end
	return Value
end

function GetYvalue( SliceBuffer, PointIndex, Yscaling )
	SliceMaxPoint=SliceBuffer:getSampleCount(1)
	if PointIndex>SliceMaxPoint then error("GetYvalue is calling PointIndex="..PointIndex.." which is above maximum value of buffer :"..SliceMaxPoint) end
	if PointIndex<1 then print("function GetYvalue: WARNING: GetYvalue is calling PointIndex="..PointIndex.." which is below 1") end
	Value = SliceBuffer:getAt( PointIndex )/Yscaling
	--Value = Value + BatchIndex * OffsetStepY
	--return string.format( "%7.3f", Value )
	return Value
end

function TableSize( table )
	local i=0
	for a,b in pairs( table ) do
		i = i + 1
	end
	return i
end
-- -----------------------------------------------------------------------------
--              END OF FUNCTION DEFINITIONS                                   --
-- -----------------------------------------------------------------------------

-- Main program

-- Check input
Project=cara:getProject( ProjectName )
if not Project then
	error("No project found with name "..ProjectName.." Please edit top of script so ProjectName = a project name existing in the repository")
end
PeakList=Project:getPeakList( PeakListId )
if not PeakList then error("No peaklist found in this project\n with the Id "..PeakListId.. "\nPlease edit line with PeakListId at the top of the script\n so it refers to the PeakListId of the \npeaklist with a BatchList defined") end
print("Selected PeakList "..PeakList:getName())
BatchList=PeakList:getBatchList() -- if batchlist exists,
--script will loop through all spectra in list and
--write out one pair of peak intensity,linewidth columns per spectrum in the batchlist
Spectrum=Project:getSpectrum( SpecId )
if not Spectrum then
	error("No spectrum found with Id="..SpecId.." in the selected project.\nPlease edit the top of script so SpecId= the Id of a spectrum in the selected project")
end
FilePrefix=Project:getName()
ProposedFileName = FilePrefix.."_".."LWf"..DimTrace..FileSuffix
Filename = dlg.getText("Enter the output filename", "", ProposedFileName )

-- b. get SliceBuffers
t.SliceBuffer = {}

outfile = io.output( Filename )
--print("Writing out "..TableSize(BatchList).." rows of spectra in  columns of residues")

Line = {}
Header="PeakNumber"..sp.."PeakWidth"
p=0
Line[p]=Header
NumPeaksWithAssymetry=0
NumPeaks=0
SumFullWidths=0
SumOfSqFullWidth=0
if not BatchList then
	BatchList[1] = SpecId
	BatchTotNum = 1
else
	BatchTotNum = TableSize(BatchList)
end

for PeakId,Peak in pairs( PeakList:getPeaks() ) do
	p=p+1
	Line[p]=PeakId
	if Verbose then print(PeakId) end

	for BatchNum=1,BatchTotNum do
	Spectrum = Project:getSpectrum( BatchList[ BatchNum ] )
	PpmDim1,PpmDim2 = GetPos( Peak, Spectrum ) -- get position of reference spins
	StartPpmSourceDim,EndPpmSourceDim = GetPpmRange( DimTrace, PpmDim1, PpmDim2, DeltaPpm, Spectrum ) -- get ppm range including RefSpinDim1
	if DimTrace == 2 then --OrthoDim is the dimension orthogonal to the slice
		PpmOrthoDim = PpmDim1 --indirect dim is orthogonal to slice
	else
		PpmOrthoDim = PpmDim2 -- direct dimension is orthogonal to slice
	end
	Slice = GetSlice( Spectrum, StartPpmSourceDim, EndPpmSourceDim, PpmOrthoDim ) -- get the 1D slice including the peak of RefSpinDim1
	MinBufferSize = GetMinBufferSize( AbsMinBufferSize, Slice )
	i_offset = tonumber( ( Slice:getSampleCount(1) - MinBufferSize )/2 )
	FullWidth, Assymetry, MinFullWidth, MaxFullWidth,LeftHalfHeight,RiteHalfHeight=GetHalfMaximaAndWidth( Slice, i_offset, SpecId, Unit )
	ForAssymetry = string.format("%3.3f",Assymetry)
	ForFullWidth = string.format("%3.2f",FullWidth)
	if Assymetry > MaxAssymetry then
		print("WARNING: In spectrum with ID"..BatchList[BatchNum].." peak "..PeakId.."has Assymetry= "..ForAssymetry.."!!!".."MinFullWidth="..Float(MinFullWidth).." MaxFullWidth="..Float(MaxFullWidth).." FullWidth="..Float(FullWidth) )
		--Line[p] = Line[p]..sp -- leave empty value if Linewidth is not reliable
		Line[p] = Line[p]..sp..Float(MinFullWidth) -- use the lower FulllWidth (assume the other side has overlap)
		NumPeaksWithAssymetry = NumPeaksWithAssymetry + 1
	else								
		NumPeaks = NumPeaks + 1
		SumFullWidths = SumFullWidths + FullWidth
		SumOfSqFullWidth = SumOfSqFullWidth + FullWidth*FullWidth
		--print(RefResNum,FullWidth) -- TEST TEST TEST
		Line[p] = Line[p]..sp..ForFullWidth
	end -- if Asymmetry is too high
	
	end -- for all entries in BatchList
	if Verbose then print(Line[p]) end
end

-- write lines to terminal window
for m=0,p do print(Line[m]) end

for n=0,p do
	outfile:write( Line[n].."\n" )
end
outfile:close()