#pragma rtGlobals=1		// Use modern global access method.
#pragma version =1.1

//Verified release.
//	Verification statement: 	This code has been tested against original implementation of this code which is part of Irena package ("Sizes"). 
//							It has been found to perform as well as the original code.
// 
//Jan Ilavsky, ilavsky@aps.anl.gov
//		date: August 30, 2009
// 			update JIL ver 1.1 August 30, 2009... Fixed two  errors, see comments in the code. 
//   note: There is still this line related to conversion of results into binned distribution, which I do nto know how useful or appropritae is for general purpose code. 

// Total non-negative least square method 
// reference:   http://www.caam.rice.edu/~zhang/reports/tr0408.ps 
// (Michael Merritt and Yin Zhang, Technical Report TR04-08, Department of Computational and Applied Mathematics, Rice University, Houston, Texas 77005, U.S.A., May, 2004)
// this method solves problems, which can be described as :
//  B = A x Model, where B is measured data (with uncertainities), A is n x m matrix and Model is what we are looking for.
// Note, B has associated x-data in Wave BdataBinning
//  Model has associated X-data in Wave ModeldataBinning
// none of the binning needs to be linear. 
// ApproachParameter is "step" - needs to be smaller than 1, usually 0.6 is good. Reasonable range seems to be 0.3 - 0.99
// limit MaxNumIterations to sensible number... Depends on complexity of problem... 


 Function TNNLS(AmatrixInput,ModelWaveOutput, ModelDataBinning,BvectorInput, BdataBinning, UncertainitiesInput, ApproachParameter, MaxNumIterations)
	wave AmatrixInput,ModelWaveOutput,ModelDataBinning,BvectorInput, BdataBinning, UncertainitiesInput
	variable ApproachParameter, MaxNumIterations		
	//Amatrix is n x m matrix relating ModelData to BVector(measured data).
	 
	 //create working space
	 string OldDf=GetDataFolder(1)
	 NewDataFolder/O/S root:Packages
	 NewDataFolder/O/S root:Packages:NNLS
	 //create local copies of input waves
 	 Duplicate/O AmatrixInput, AmatrixInputE, AmatrixOrig
	 Duplicate/O BvectorInput, BvectorInputE
	//note this code defaults to use uncertainities. If these are not available, set following parameter to 0 and modify the code
	 variable useUncertainitiesInput=1
	 if(useUncertainitiesInput)
	 	AmatrixInputE = AmatrixInput[p][q] / UncertainitiesInput[p]
	 	MatrixOp/O BvectorInputE = BvectorInput / UncertainitiesInput
	 endif
	 //create local meaningful waves to work with which account for the errors used.
	Duplicate/O AmatrixInputE, Amatrix
	Duplicate/O ModelWaveOutput, ModelWave
	Duplicate/O BvectorInputE, BVector, BvectorOrig, CurrentResultB
	Duplicate/O UncertainitiesInput, Uncertainities
	//create some parameters, which the user may actually use after the run to find what happened. 
	variable/g NumberIterations, Chisquare
	//first we will need some variables
	variable i,j
	//we will need these multiple times, precalcualting may save time.
	//	Amatrix[][] = AmatrixInputE[p][q] * BdataBinning[p]...error in versxion 1. Removed JIL 8/30/09
	MatrixOp/O AmatrixT=Amatrix^t
	MatrixOp/O BVector = BvectorOrig //* BdataBinning... error in version 1. Removed JIL 8/30/09
	//starting conditions
	redimension/D ModelWave
	//we will start with reliably small positive number. Assume 1e-32 is such, but this may fail in some cases... 
	ModelWave = 1e-32
	//working waves & variables
	variable numIter=0
	variable lasterr=1e3
	variable err=0	
	variable alphaStar, temp1, temp2
	//here the iterations start...
	Do
	//start of NNLS Interior point gradient method itself. Following steps designation relates to the original paper... 
		//step 1
		MatrixOp/O Qk = AmatrixT x Amatrix x ModelWave - AmatrixT x BVector
		MatrixOp/O Dk = ModelWave / (AmatrixT x Amatrix x ModelWave)
		MatrixOp/O Pk = - Dk * Qk	
		//step 2	
		MatrixOp/O AkSTAR= (pk^t x AmatrixT x Amatrix x Pk)
		temp1 = AkSTAR[0]
		MatrixOp/O AkSTAR= - (pk^t x Qk) 
		temp2 = AkSTAR[0]
		alphaStar = temp2 / temp1
		redimension/N=(numpnts(ModelWave)) AkSTAR
		AkSTAR = alphaStar
		//above is ideal step to make//below is limiting the step so we do not get negative values...		
		MatrixOP/O AlphaWv =  - ModelWave/Pk		//this is maximum alpha, which we can make, if pk is negative
		For(i=0;i<numpnts(PK);i+=1)
			if(Pk[i]<0)		//if pk is negative, we may have to limit the step
				AkSTAR[i]=min(ApproachParameter*AlphaWv[i],AkSTAR[i])		//the step is limited to smaller from the two values
			endif
		endfor
		//step 3
		MatrixOp/O ModelWave = ModelWave +(AkSTAR * Pk ) 	//new model, go back after some calculations below
	//end of NNLS Interior point gradient method itself	
	//figure out chisquare so we know if we should bail out... 
		MatrixOp/O CurrentResultB = AmatrixOrig x ModelWave		//thsi is calculated data from our model
		duplicate/O BvectorInput, tempWv
		tempWv = (BvectorInput - CurrentResultB)/UncertainitiesInput
		tempWv = tempWv^2
		Chisquare = sum(tempWv)
		err=Chisquare/numpnts(BvectorInput)
		numIter+=1
		NumberIterations=numIter
	while((err>1 && numIter<MaxNumIterations))	
	//this bails out from TNNLS if :
		//1. Chiaquare/number of point is less than 1 (fit within uncertainities
		//2. Number of iterations is above user input value.
	
	KillWaves/Z Amatrix, tempWv, AkSTAR, Qk, Dk, Pk, AmatrixT, Amatrix, BVector	//cleanup.
	 ModelWaveOutput = ModelWave //I think this should be the right output, but in my code I need to do following line... Why?
	ModelWaveOutput = ModelWave[p]/ (2*TNNLS_BinWidthInRadia(ModeldataBinning,p))		//I am really not sure if this is of any use for general purpose code. 
	// this was part of my original size distribution code and it converts the results to binned distribution 
	// I think it is related to the Amatrix containing the BdataBinning...	
	//returns the result into the user provided wave... ModelWaveOutput
	setDataFolder OldDf
	
end


Function TNNLS_BinWidthInRadia(ModelDataBinning,i)			//calculates the width in radia by taking half distance to point before and after
	Wave ModelDataBinning
	variable i								//returns number in A

	variable width
	variable Imax=numpnts(ModelDataBinning)
	
	if (i==0)
		width=ModelDataBinning[1]-ModelDataBinning[0]
	elseif (i==Imax-1)
		width=ModelDataBinning[i]-ModelDataBinning[i-1]
	else
		width=((ModelDataBinning[i]-ModelDataBinning[i-1])/2)+((ModelDataBinning[i+1]-ModelDataBinning[i])/2)
	endif
	return width
end