
	
/******************************************************************
MAIN Function Find PCA component
******************************************************************/
number MSA112_FindNIPALS(image &dataMatrix, image &Pvec, image &Tvec,number Pyth,image &EVs,image &devia, image &spikes, image &spikesXY,number thres)	
{	//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)

//decode spikes paprameters
number spikesCode = spikes.getPixel(0,0)
number TreatOut,OutCorr,ShiftCor,LoadThres
number spikesNful	//number of found spikes over all iterations
image OutliersOld
if(spikesCode != 0)
	{
	TreatOut = MSA112_Decode(spikesCode,10)	//1-treat score outliers, 2-treat loading outliers
	OutCorr = MSA112_Decode(spikesCode,10)	//1-remove by shifting pixels,2-remove ba shifting channles
	ShiftCor = MSA112_Decode(spikesCode,10)
	LoadThres = spikesCode	//threshold in loadings
	if(spikesXY.ImageIsValid()) spikesNful = spikesXY.imageGetDimensionSize(1)
	
	//make OutliersOld mask	to discard previously found spikes
	if(spikesXY.ImageIsValid())
		{
		number SpikesN  = spikesXY.ImageGetDimensionSize(1)
		if(TreatOut==1) OutliersOld = MSA112_SetOldOutliers(Tvec,spikesXY[0,0,spikesN,1])
		if(TreatOut==2) OutliersOld = MSA112_SetOldOutliers(Pvec,spikesXY[0,1,spikesN,2])
		}
	}
//result("\ntreatOut "+treatOut+" OutCorr "+OutCorr+" LoadThres "+LoadThres+ " thres "+thres)


//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
	
image dataMatrixT 
if(Pyth <=1) 
	{
	//0th Tvec-the itteration
	Tvec = matrixMultiply(dataMatrix,Pvec)		
	//to save time for transpose
	dataMatrixT = MatrixTranspose(dataMatrix)
	}
	
//to shortcut for turbo
image dataMatrixCC
if(Pyth ==2) 
	{
	dataMatrixCC = MSA112_matrixMultiply(MatrixTranspose(dataMatrix),dataMatrix)

	if(treatout ==2) dataMatrixT = MatrixTranspose(dataMatrix)
	}
	
//LOOP
number Nitt = 0
number Ccos,Del
number converged

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

	//NIPALS step
	if(Pyth ==0) Pvec = matrixMultiply(dataMatrixT,Tvec) 
	if(Pyth ==1) Pvec = MSA112_MatrixMultiply(dataMatrixT,Tvec) 
	if(Pyth ==2) Pvec = MSA112_MatrixMultiply(dataMatrixCC,Pvec) 
	
	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
	
	
	if(Pyth <=1) 
		{
		if(Pyth ==0) Tvec = matrixMultiply(dataMatrix,Pvec)
		if(Pyth ==1) Tvec = MSA112_MatrixMultiply(dataMatrix,Pvec)
				
		EVcurr = sum(MatrixMultiply(MatrixTranspose(Tvec),Tvec))//MSA112_Variance(Tvec)
		EVs.setPixel(0,Nitt, EVcurr)
		Del = modulus((EVcurr-EVold)/EVold)	
		}
	//if Pyth ==2 (turbo) no calculation of intermediate Tvec, EVs

	Ccos = sum(MatrixMultiply(MatrixTranspose(PvecOld),Pvec)) 
	devia.setPixel(0,Nitt,Acos(Ccos))
	
	//spikes treatment
	if(TreatOut>0)
		{

		image outliers	//clone of score/loading with deviation at outlier position, otherwise ZERO
		if((TreatOut==1)&&(Pyth!=2)) outliers = MSA112_FindScoreOutliers(Tvec,thres)			//score	
		if(TreatOut==2) outliers = MSA112_FindLoadingOutliers(Pvec,LoadThres)	//loading

		if(outliersOld.ImageIsValid()) outliers = MSA112_ReduceOldOutliers(outliers,outliersOld,0.000001)
	
		//get number of outliers
		number spikesN = MSA112_CountNonZeros(outliers,thres)	
//result("\nspikeN "+spikesN+ "      "+EVcurr)
		spikes.setPixel(0,Nitt,spikesN)
//if(Nitt==1) outliers.matrixTranspose().showimage()	
		
		//remove outliers one by one	
		if(spikesN>0)
			{ 		
			image spikesXYitt:= integerimage("loc",2,0,4,spikesN)
			
			number CurSpike
			number pixel,channel
			
			for(number i=0;i<spikesN;i++)
				{				
				number Value,pixel,channel
//result("\nspike "+i)				
				Value = MSA112_RemoveOutlierFromMask(dataMatrix,dataMatrixT,outliers,pixel,channel,thres,treatOut,OutCorr,ShiftCor,LoadThres)	
//result(" pixel "+pixel+" channel "+channel+" "+Value)				
				Value *= 256	//resolution 1/256*sigma 
				if(Value > 65536) Value = 65536//max 256*sigma 				
				Value = round(Value)
//result(" int "+ value)	
			
				//add to outliear coordinate list
				spikesXYitt.setPixel(0,i,pixel)
				spikesXYitt.setPixel(1,i,channel)
				spikesXYitt.setPixel(2,i,Value)	

				//update OutliersOld mask

				if(outliersOld.ImageIsValid()) 
					{
					if(TreatOut==1) OutliersOld.setPixel(0,pixel,1)
					if(TreatOut==2) OutliersOld.setPixel(0,channel,1)
					}

				}
			
		
			//attach fresh found spikes
			spikesXY := MSA112_AugmentImage(spikesXY,spikesXYitt,spikesNful)
			}
			
		}
	
	//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
		}
	}

if(Pyth==2)
	{
	//get Tvec and EV in case of turbo
	Tvec = MSA112_MatrixMultiply(dataMatrix,Pvec)			
	EVcurr = sum(MatrixMultiply(MatrixTranspose(Tvec),Tvec))//MSA112_Variance(Tvec)
	EVs.setPixel(0,Nitt, EVcurr)
	}
	
return Nitt
}


/***************** NIPALS with transparent input (full) ******************************/
number MSA112_FindNIPALS(image &dataMatrix, image &Pvec, image &Tvec,number Pyth, image &EVs,image &devia,number converg,number DeviaSpec,number maxItt, image &spikes,image &spikesXY,number outTreat,number OutCorr,number ShiftCor,number LoadThres, number thres)	
{
//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
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)
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

//prepare history of spikes
spikes:= integerimage("outliers",2,0,1,maxItt+1)//the 0-pixel must consist of encoded spikes tretment info
if(OutTreat!=0)
	{
	number SpikesCode = MSA112_Encode(OutTreat,OutCorr,ShiftCor,LoadThres)
	spikes.setPixel(0,0,SpikesCode)
	}

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,Pyth,EVs,devia,spikes,spikesXY,thres)	

return Nitt
}


/***************** compact simplified NIPALS with transparent input with Python ******************************/
number MSA112_FindNIPALS(image &dataMatrix, image &Pvec, image &Tvec, number Pyth, number &EigenValue,number converg,number DeviaSpec,number maxItt,image &spikes,image &spikesXY,number outTreat,number OutCorr,number ShiftCor,number LoadThres, number thres)	
{
image devia,EVs

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,Pyth,EVs,devia,converg,DeviaSpec,maxItt,spikes,spikesXY,outTreat,OutCorr,ShiftCor,LoadThres,thres)

EigenValue = EVs.getPixel(0,Nitt)

return Nitt
}

/***************** compact simplified NIPALS with transparent input no Python (DOCUMENTED ADVANCED ) ******************************/
number MSA112_FindNIPALS(image &dataMatrix, image &Pvec, image &Tvec,number &EigenValue,number converg,number DeviaSpec,number maxItt,image &spikes,image &spikesXY,number outTreat,number OutCorr,number ShiftCor,number LoadThres, number thres)	
{
image devia,EVs
number Pyth =0

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,Pyth,EVs,devia,converg,DeviaSpec,maxItt,spikes,spikesXY,outTreat,OutCorr,ShiftCor,LoadThres,thres)

EigenValue = EVs.getPixel(0,Nitt)

return Nitt
}

/***************** simplified NIPALS without spikes with Python******************************/
number MSA112_FindNIPALS(image dataMatrix, image &Pvec, image &Tvec,number Pyth, image &EVs,image &devia,number converg,number deviaSpec,number maxItt)	
{
//dummy input for outliers
number OutTreat = 0
number OutCorr,ShiftCor,LoadThres,thres
image spikes,spikesXY

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,Pyth,EVs,devia,converg,DeviaSpec,maxItt,spikes,spikesXY,outTreat,OutCorr,ShiftCor,LoadThres,thres)	

return 	Nitt
}

/***************** simplified NIPALS without spikes no Python ******************************/
number MSA112_FindNIPALS(image dataMatrix, image &Pvec, image &Tvec,image &EVs,image &devia,number converg,number deviaSpec,number maxItt)	
{
//dummy input for outliers
number OutTreat = 0
number OutCorr,ShiftCor,LoadThres,thres
image spikes,spikesXY
number Pyth =0

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,Pyth,EVs,devia,converg,DeviaSpec,maxItt,spikes,spikesXY,outTreat,OutCorr,ShiftCor,LoadThres,thres)	

return 	Nitt
}

/***************** simplified NIPALS without spikes, without deviation/eigenValue history with Python ******************************/
number MSA112_FindNIPALS(image dataMatrix, image &Pvec, image &Tvec, number Pyth,number &EigenValue,number converg,number deviaSpec,number maxItt)	
{
//prepare dummy devia/EigenValues history arrays
image devia,EVs

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,Pyth,EVs,devia,converg,deviaSpec,maxItt)	

EigenValue = EVs.getPixel(0,Nitt)

return Nitt
}

/***************** simplified NIPALS without spikes, without deviation/eigenValue history no Python(DOCUMENTED BASIC ) ******************************/
number MSA112_FindNIPALS(image dataMatrix, image &Pvec, image &Tvec,number &EigenValue,number converg,number deviaSpec,number maxItt)	
{
//prepare dummy devia/EigenValues history arrays
image devia,EVs
number Pyth =0

number Nitt = MSA112_FindNIPALS(dataMatrix,Pvec,Tvec,Pyth,EVs,devia,converg,deviaSpec,maxItt)	

EigenValue = EVs.getPixel(0,Nitt)

return Nitt
}

/***************** simplified NIPALS without spikes, without deviation/eigenValue history, with ANGLE convergence only (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 without spikes, without deviation/eigenValue history, with ANGLE convergence only, with default deviaSpec and maxItt with Python*****/
number MSA112_FindNIPALS(image dataMatrix, image &Pvec, image &Tvec,number Pyth, number &EigenValue)	
{
number deviaSpec = 0.02
number maxItt = 50
number converg = 0

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

return Nitt
}

/****** simplified NIPALS without spikes, without deviation/eigenValue history, with ANGLE convergence only, with default deviaSpec and maxItt, no Python(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
}



