
	
/******************************************************************
MAIN Function Find PCA component
******************************************************************/
number MSA112_FindNIPALS(image &dataMatrix, image &Pvec, image &Tvec, image &EVs,image &devia)	
{	//NIPALS algorithm	
//Pvec loadings (m,Depth)
//Tvec scorea   (m,N)
//Pvec is important as an entrance, Tvec is not important

//EVs (1,maxItt+1)matrix		0th-pixel to keep the specified precision of EV
//devia (1,maxItt+1)matrix		0th-pixel to keep the specified precision of angle
//spikes (1,maxItt+1)matrix		0th-pixel to keep the spikes code where all parameters are
//spikesXY(spikesNful,2)	keeps x,y coordinates of spikes,spikesNful is predefined (adjusted in situ)

	
//get max number of itterations
number maxItt = EVs.ImageGetDimensionSize(1) - 1
	
//set convergence parameters
//if delSpec > 10(-5) VALUE convergence, otherwise ANGLE convergence
number DelSpec = EVs.getPixel(0,0)	//percentage		
number CcosSpec = devia.getPixel(0,0)
//result(" DelSpec "+DelSpec+" CcosSpec "+CcosSpec)
if(DelSpec>0.00001) CcosSpec = 2//never converge on angle
else DelSpec = -1				//never converges on value
//result(" DelSpec "+DelSpec+" CcosSpec "+CcosSpec)


//prepare
number EVcurr
number EVold= -1 //to avoid divergence
	//ensure normalization
number Plength = sqrt(sum(MatrixMultiply(MatrixTranspose(Pvec),Pvec)))
if(Plength==0) throw("Zero P-vector at entrance !!")
Pvec /= Plength
image PvecOld = Pvec
	

//0th Tvec-the itteration
Tvec = MSA112_MatrixMultiply(dataMatrix,Pvec)		
//to save time for transpose
image dataMatrixT = MatrixTranspose(dataMatrix)

	

//LOOP
number Nitt = 0
number Ccos,Del
number converged

while(!converged)	//itteration on vector, converge at Ccos==1
	{
	Nitt++

	//NIPALS step
	Pvec = matrixMultiply(dataMatrixT,Tvec) 

	
	Plength = sqrt(sum(MatrixMultiply(MatrixTranspose(Pvec),Pvec)))
	if(Plength==0) 
		{
		OKdialog("Zero P-vector at "+Nitt+"-itteration !!")
		throw("")
		}
	if(Plength==infinity()) 
		{
		number dim =Pvec.imagegetDimensionSize(1)
		number overflow =mean(Pvec)/sqrt(1e38/dim)
		OKdialog("overfloating !!\nreduce intensity scale\nat least in "+overflow +" times")
		throw("")
		}
		
	Pvec /= Plength
	

	Tvec = MSA112_MatrixMultiply(dataMatrix,Pvec)
	EVcurr = sum(MatrixMultiply(MatrixTranspose(Tvec),Tvec))//MSA112_Variance(Tvec)
	EVs.setPixel(0,Nitt, EVcurr)
	Del = modulus((EVcurr-EVold)/EVold)	


	Ccos = sum(MatrixMultiply(MatrixTranspose(PvecOld),Pvec)) 
	devia.setPixel(0,Nitt,Acos(Ccos))
	
	
	//check if converged
	if((Ccos>=CcosSpec)||(Del<DelSpec))
			{
			//result(" |Ccos "+Ccos+" |del "+del)
			converged = 1
			}
	EVold = EVcurr
	PvecOld = Pvec
//result("\nNitt "+Nitt)
		
	//check if max iterations is reached 
	if(Nitt >=MaxItt) 
		{
		result(" PRECISION? ")
		break
		}
	}
	
return Nitt
}



/***************** compact simplified NIPALS with transparent input (DOCUMENTED BASIC ) ******************************/
number MSA112_FindNIPALS(image &dataMatrix, image &Pvec, image &Tvec,number &EigenValue,number converg,number DeviaSpec,number maxItt)	
{
//get size of matrix
number kN,Depth
DataMatrix.getSize(Depth,kN)

//define Pvec,Tvec
Tvec:= realimage("T-vector",4,1,kN)
Pvec:= realimage("P-vector",4,1,Depth)
Pvec.setPixel(0,0,1) 

//prepare history of EigenValues
image EVs := realimage("Variance",4,1,maxItt+1)
//the 0-pixel in EVs array must consist the maxRelDiffInEigenValues
if(converg = 0) EVs.setPixel(0,0,0)	//angular convergence; 0th-pixel signals about that
else EVs.setPixel(0,0,deviaSpec/100) //value convergence, 0th-pixel in EVs is deviaSpec interpreted in % 

//prepare history of deviations
number minCos = cos(deviaSpec)
image devia:= realimage("Precision",4,1,maxItt+1)
//the 0-pixel in devia array must consist the minCosBetweenLoadingVectors
devia.setPixel(0,0,minCos) //rest is unuased

number Pyth =0

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,EVs,devia)

EigenValue = EVs.getPixel(0,Nitt)

return Nitt
}





/***************** simplified NIPALS  witho default convergence(DOCUMENTED BASIC )******************************/
								//CANNOT INCLUDE PYTHON
number MSA112_FindNIPALS(image dataMatrix, image &Pvec, image &Tvec,number &EigenValue,number deviaSpec,number maxItt)	
{
number converg = 0
number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,EigenValue,converg,deviaSpec,maxItt)	

return Nitt
}





/****** simplified NIPALS with default deviationSpec amd max iterations (DOCUMENTED BASIC )*****/
number MSA112_FindNIPALS(image dataMatrix, image &Pvec, image &Tvec,number &EigenValue)	
{
number deviaSpec = 0.02
number maxItt = 50

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,EigenValue,deviaSpec,maxItt)	

return Nitt
}