matrix algebra functions with ImageJ macro

matrixalgebra.ijm

// matrix algebra functions
//Kota Miura (miura@embl.de)
//050131 2D matrix to 1D arary and vice versa
//	matrix transpose, multiplication
//050201 Pivot, Gauss-Jordan (half done)

//  matrix format: m[0] contains number of rows in matrix. In the order of a11(=m[1]), a21(=m[2]), a31, a21,a22.....
// vector format: no info on length. starts from v[0]

//****************************** Global variables ***************************************** 
var ijA=newArray(2);
var jiA=newArray(2);

//****************************** 1D index to 2D indice  ***************************************** 

function Kmatval(mat,i,j){
	matvalue=0.000000;
	matvalue=mat[K_getK2(mat[0],i,j)];
	return matvalue;
}

// known k, uinknown ij
function K_getij(rows,ijA,k) {
	ijA[1]=floor((k-1)/rows);
	ijA[0]=k-1-ijA[1]*rows;
}

//known ij, unknown k
function K_getK(rows,ijA) {
	return ijA[1]*rows+ijA[0]+1;
}
function K_getK2(rows,i,j) {
	return j*rows+i+1;
}

function KviewMat(viewMat) {
	print("-----------");
	colnum=(viewMat.length-1)/viewMat[0];
	for(i=0;i<viewMat[0];i++) {
		strnum="";
		for(j=0;j<colnum;j++) {
			temp="			"+Kmatval(viewMat,i,j);
			strnum=strnum+temp;
		}
		print(strnum);
	}
	print("-----------");
}
//****************************** element controls  ***************************************** 

// swapping of ath and bth value in a 3-points array
function swapArray3(a,b,Array3) {
	tempsort=Array3[a];
	Array3[a]=Array3[b];
	Array3[b]=tempsort;	
}

// matrix augumentation
function Kaugumentation(matLeft,vecRight,augMat) {
	for(i=0;i<matLeft.length;i++) 
		augMat[i]=matLeft[i];
	for(i=0;i<vecRight.length;i++) 
		augMat[matLeft.length+i]=vecRight[i]; // vecRight[i+1]?
	return 0;
}

//****************************** matrix controls  ***************************************** 
//matrix transpose
function K_transposeMat(srcmatA,destmatA) {
//	ijA=newArray(2);
	destmatA[0]=(srcmatA.length-1)/srcmatA[0];
	for(i=1;i<srcmatA.length;i++) {
		K_getij(srcmatA[0],ijA,i);
		swapArray3(0,1,ijA);
		destmatA[i]=srcmatA[K_getK(srcmatA[0],ijA)];
		print(i+": original "+srcmatA[i]+" transposed "+destmatA[i]);
	}
}
// ------------ matrix multiplication ------------------------------
function K_multiplyMat(matC,matD,resmat) {
	if (((matC.length-1)/matC[0]!=matD[0])||((matD.length-1)/matD[0]!=matC[0])) return -1;
	rows=matC[0];
	cols=matD[0];
	for (j=0;j<cols;j++) {
		for (i=0;i<rows;i++) {
			resmat[K_getK2(rows,i,j)]=Kmultiplycore(matC,matD,i,j);
		}
	}
	return 0;
}

function Kmultiplycore(matC,matD,i,j) {
	res=0.0;
	for(k=0;k<matD[0];k++) {
		res=res+Kmatval(matC,i,k)*Kmatval(matD,k,j);
		print(res);
	}
	return res;
}

// ---------- Gauss-Jordan --------------------------------------------------

//partial pivot for the Gauss-Jordan nxn matrix
//	need to implement some way to record the pivotting
function K_pivot(matPivot,startn,pivotrecA) {
	n=matPivot[0]; //row number = columnnumber
	max=0.0;
	temp=0.0;
	maxi=n;
	if ((startn+1)<n) {
		for(i=startn+1;i<n;i++) {		//find maximum position within "startn"th column
//		for(i=0;i<n;i++) {		//find maximum position within "startn"th column //050202

			temp=abs(Kmatval(matPivot,i,startn));
			print(temp);
			if (temp>max) {
				max=temp;
				maxi=i;
			}
		}
	} else {
		max=(Kmatval(matPivot,startn,startn));
		maxi=startn;
		
	}
	print("max:"+max+" index:"+maxi);
	if (maxi==startn) return -1;
	else {
//		for (i=startn;i<=n;i++) {
		for (i=0;i<=n;i++) {
			//temp=Kmatval(matPivot,maxi,i);
			//matPivot[K_getK2(n,maxi,i)]=matPivot[K_getK2(n,startn,i)];
			//matPivot[K_getK2(n,startn,i)]=temp;
			swapArray3(K_getK2(n,maxi,i),K_getK2(n,startn,i),matPivot);
		}
		print("pivot between "+startn+" and "+maxi);
		swapArray3(maxi,startn,pivotrecA);		//pivot indexing
		return 0;
	}
}

//	for (i=1;i<matPivot.length;i++) {
//		K_getij(n,ijA,i);
//		print(ijA[0]+","+ijA[1]+": "+matPivot[i]+"  Original Position: "+pivotrecA[ijA[0]]);
//	}

//********************c codes ******************************************/
//static int swap( matrix , n )
//double matrix[MATRIX][MATRIX+1] ;
//int n ;
//{
//int    i,mi;
//double max ;
//
//	for ( max = .0 , mi = n , i = n + 1 ; i < MATRIX ; ++i ) {    /*here, mi=n=j and i=j+1*/ /*how does the MATRIX being defined?*/
//		double tmp = fabs( matrix[i][n] ) ;						/*fabs returns the absolute value of floating points argument*/
//		if ( tmp > max )										/* check if the fabs worked and also not 0*/
//			max = tmp , mi = i ;								/* max defined as value at (j, j)? and mi being defined as j*/
//	}
//	if ( mi == n )
//		return -1 ;
//	for ( i = n ; i < MATRIX + 1 ; ++i ) {					/*horizontal flip between column mi and n, starting from nth row*/
//		double tmp = matrix[mi][i] ;
//		matrix[mi][i] = matrix[n][i] ;++++
//		matrix[n][i] = tmp ;
//	}
//	return 0 ;
//}/
/************************************************************************/

macro "test viematt"{
	mattestA=newArray(3,2,4,8,5,13,29,7,20,50);
	KviewMat(mattestA);
}
macro "test calc"{
	print(6.25-(2.444*3.625));
	print(-5-(2.444*-1.5));
}
macro "test Gauss Jordan elimination" {
//	mattestA=newArray(3,2,4,8,5,13,29,7,20,50);
	mattestA=newArray(3,5,8,-10,3,2,5,10,-5,1);
//	vectortestA=newArray(23,58,132);
	vectortestA=newArray(410,-30,30);
	completion=KGaussJordan(mattestA,vectortestA);
	if (completion==0) 
		for (i=0;i<vectortestA.length;i++)
			print(vectortestA[i]);
	else
		print("error");
}

macro "test pivot" {
	mattestA=newArray(3,2,4,8,5,13,29,7,20,50);
	vectortestA=newArray(23,58,132);
	KviewMat(mattestA);
	augumentMat=newArray(mattestA.length+vectortestA.length);
	n=mattestA[0];
	tmpm=0.0;
	pivotrecA=newArray(n);
	for (i=0;i<n;i++) 
		pivotrecA[i]=i;
	completion=Kaugumentation(mattestA,vectortestA,augumentMat);
	for(k=0;k<n;k++) {
		currentelement=Kmatval(augumentMat,k,k); //print("current element "+k+":"+currentelement);
		completion=K_pivot(augumentMat,k,pivotrecA); //print("pivot return "+ completion);
	}	
	for (i=1;i<augumentMat.length;i++) {
		K_getij(n,ijA,i);
		print(ijA[0]+","+ijA[1]+": "+augumentMat[i]);
	}
}


function KGaussJordan(matDev,vectorGJ) {
	augumentMat=newArray(matDev.length+vectorGJ.length);
	//for (i=0;i<augumentMat.length;i++)
	//	augumentMat[i]=0.000;
	n=matDev[0];
	//temp1=0.000;
	pivotrecA=newArray(n);
	for (i=0;i<n;i++) 
		pivotrecA[i]=i;
	completion=Kaugumentation(matDev,vectorGJ,augumentMat);
	if (completion!=0) return -1;
//	for (i=1;i<augumentMat.length;i++) {
//		K_getij(n,ijA,i);
//		print(ijA[0]+","+ijA[1]+": "+augumentMat[i]);
//	}
	for(k=0;k<n;k++) {
		currentelement=Kmatval(augumentMat,k,k); print("current element "+k+":"+currentelement);
		KviewMat(augumentMat);
		completion=K_pivot(augumentMat,k,pivotrecA); 
		KviewMat(augumentMat);
		//completion=0;
		if ( (currentelement==0) && (completion==-1) ) {
			return -1;
		}
//		for (j=k+1;j<=n;j++) {	//050202 check del
		division=Kmatval(augumentMat,k,k);
		for (j=0;j<=n;j++) {
//			temp1=(Kmatval(augumentMat,k,j)/Kmatval(augumentMat,k,k));
			temp1=(Kmatval(augumentMat,k,j)/division);

			augumentMat[K_getK2(n,k,j)]=temp1;
		}
		print("division");
		KviewMat(augumentMat);
		
		for (i=0;i <n; i++) {
			if (i != k) {
//				for (j=k+1;j<=n;j++) {
				multiply=Kmatval(augumentMat,i,k);
				for (j=0;j<=n;j++) {
//					temp1=(Kmatval(augumentMat,i,k)*Kmatval(augumentMat,k,j));
					temp1=(multiply*Kmatval(augumentMat,k,j));
					temp2=(Kmatval(augumentMat,i,j)-temp1);
					//print(Kmatval(augumentMat,i,k)+"*"+Kmatval(augumentMat,k,j)+"="+temp1+" : "+temp2);
					augumentMat[K_getK2(n,i,j)]=temp2;
				}
			}
		}
		print("multiply & subtract");
		KviewMat(augumentMat);

		//for(p=1;p<augumentMat.length;p++) {
		//	K_getij(n,ijA,p);
		//	print(ijA[0]+","+ijA[1]+": "+augumentMat[p]);
		//}
	}
	for(i=0;i<n;i++) 
		vectorGJ[i]=Kmatval(augumentMat,pivotrecA[i],n);
	return 0;
}

// java source from "Java algorithm dictionary"
// a part of Gauss Jordan source code from its web site
//http://oku.edu.mie-u.ac.jp/~okumura/java-algo/archive/java-algo.zip

//    public static void solve(double[][] a) {
//        int n = a.length; 
//       for (int k = 0; k < n; k++) {
//            for (int j = k + 1; j <= n; j++) a[k][j] /= a[k][k];		/ j <= n+1???
//            for (int i = 0; i < n; i++)
//                if (i != k)
//                    for (int j = k + 1; j <= n; j++)
//                        a[i][j] -= a[i][k] * a[k][j];
//        }
//    }


//************************************************************************/
//int gauss_jordan( d0 , d1 )
//double d0[MATRIX][MATRIX] ;
//double d1[MATRIX] ;
//{
//double matrix[MATRIX][MATRIX+1] ;
//int i , j , k ;
//int swap() ;
//		/*combining "matrix" and "dim" into a singel matrix, with dim at the bottom row*//
//	for ( j = 0 ; j < MATRIX ; ++j ) {
//		for ( i = 0 ; i < MATRIX ; ++i )
//			matrix[j][i] = d0[j][i] ;
//		matrix[j][MATRIX] = d1[j] ;
//	}
//
//	for ( j = 0 ; j < MATRIX ; ++j ) {
//		double tmpm ;
//		if ( matrix[j][j] == .0 && swap( matrix , j ) == -1 )
//			return -1 ;
//		for ( tmpm=matrix[j][j] , i = 0 ; i < MATRIX + 1 ; ++i )
//			matrix[j][i] /= tmpm ;
//
//		for ( i = 0 ; i < MATRIX ; ++i ) {
//			double tmp = matrix[i][j] ;
//			if ( i != j )
//				for ( k = 0 ; k < MATRIX + 1 ; ++k )
//					matrix[i][k] -= matrix[j][k] * tmp ;
//		}
//	}
//	for ( i = 0 ; i < MATRIX ; ++i )
//		d1[i] = matrix[i][MATRIX] ;
//	return 0 ;
//}/

//***************************************************************************
macro "test pos convet" {
	matA=newArray(4,0,1,2,3,0,1,2,3,0,1,2,3,0,1,2,3);
	ijA=newArray(2);
	for (i=1;i<matA.length;i++) {
		K_getij(matA[0],ijA,i);
		print(i+":"+" i:"+ijA[0]+" j:"+ijA[1]+" k:"+K_getK(matA[0],ijA));  
	}
}

macro "test transpose" {
	matA=newArray(4,0,1,2,3,0,1,2,3,0,1,2,3,0,1,2,3);
	tmatA=newArray(matA.length);
	resmat=newArray(matA.length);
	matX=newArray(4,0,1,2,3,0,4,5,6,0,7,8,9,0,10,11,12);	
	K_transposeMat(matA,tmatA);
	//print("go back");
	//K_transposeMat(tmatA,matA);
	K_multiplyMat(matX,tmatA,resmat);
	for(i=1;i<resmat.length;i++) {
		K_getij(4,ijA,i);
		print("results "+ijA[0]+","+ijA[1]+" :"+resmat[i]);
	}	
}