
//This script investigates the dependence of a eigenvaue on the number of pixels in dataset

/************** This function generates a random mask ******************/
image randomCluster(number N, number fract)
	{
	image Cluster:= realImage("cluster",4,1,N)
	//mask image of total N pixels
	
	number range = 1/(1-fract)
	
	Cluster = random()*range
	Cluster = tert(cluster>1,1,0)
	//fract*N of pixels are 1, (1-fract)*N pixels are 0
	//choice is random
	
	return cluster
	}




//parameters for convergence  
number deviaSpec = 0.02
number maxItt = 50

//get an MSA image where the dataMatriox is kept as an attachment 
image MSAimage := getFrontImage()
image DataMatrix = MSA112_GetMSAMatrix(MSAimage)
							//DataMatrix.showimage()	//you might look at matrix if you want
							
//get the number of pixels					
number m = DataMatrix.ImageGetDimensionSize(1)

//define range and steps for fractioning the dataMatrix
number fractIni = 0.1
number fractFin = 1
number fractStep = 0.1


//fractioning
for(number fract = fractIni;fract <= fractFin;fract = fract + fractStep)
	{
	//make a random mask to select a certain fraction from DataMatrix
	image Cluster:= integerImage("cluster",2,0,1,m)
	Cluster = randomCluster(m,fract) 
	number actualFract = mean(Cluster)	
	//actual fraction can slightly differ from the nominal one
	result("\nfraction "+actualFract)
	
	//apply the mask to DataMatrix
	image DataMatrixSmaller  = MSA112_TakeCluster(dataMatrix,1,cluster,1)	
			//1: do centering; 0: do not
			//mask with index 1 is chosen

	//declare values and images
	number EigenValue	//eigenvalue of a component
	image Pvec,Tvec	//loading and score of a component
	number Nitt		//number of itterations required for convergence
	
	//extract 1st PCA component
             Nitt=MSA112_FindNIPALS(DataMatrixSmaller,Pvec,Tvec,EigenValue)
							//Pvec.MatrixTranspose().showimage()	
	//you might look at loading if you want
							
	EigenValue /= (m*actualFract)
	//now the eigenvalue is the variance within a component
	result(" eigenValue1 "+EigenValue)

	//reduce DataMatrix by the 1stPCA component
	DataMatrixSmaller = MSA112_ReduceDataMatrix(DataMatrixSmaller,Pvec,Tvec)
	
	//extract 2nd PCA component
           Nitt = MSA112_FindNIPALS(DataMatrixSmaller,Pvec,Tvec,EigenValue)
	
	EigenValue /= (m*actualFract)	
	//now the eigenvalue is the variance within a component
	result(" eigenValue2 "+EigenValue)
	}
