Finite mixture model : decompose a surface into "n" subadjacent gaussian bells.

A Mathematica (7.0) application

This script implements the Laplacian filter and the finite mixture approach to decompose a surface into an arbitrary number “n” of bidimensional Gaussian probability distributions. Finite misture are widely used in data mining or pattern recognition. Here, this tool in implemented to decompose an arbitrary surface f(x,y) with some peaks and shoulders. This approach can be useful to decompose any bidimensional data array. For instance a bidimensional chromatogram, an excitation-emssion fluorescence spectra, a bidimensional cytogram...etc...etc.....

The scrip implements three functions:

SURFACE[p_, SizeMatr_]: it generate an irregular smoothed surface (f(x,y)). In addition it estimates the Laplacian surface (∇2𝒻) with the built-in function “LaplacianGaussianFilter”. "SizeMatr" change the size the surface (at larger surface sizes, peaks are more spread). The output are the two contour plots of f(x,y) and  ∇2𝒻, for instance:

MinSearch[iter_, SizeMatr_]: it founds the local minima in the ∇2𝒻surface. The number “n” of local minima and their location provide the information about the number of potential peaks that might decompose the surace f(x,y) and position of maxima of the peaks. The search is performed with the built-in function “NMinimize”. Search method is the Nelder-Mead. “iter” is an arbitrary number that  defines the “RandomSeed”. Higher “RandomSeed” values allow to execute a more exhaustive search of the local minima (the computational time increases). The output is a list that includes the coordinates of each peak maxima (MinSearch[iter, SizeMatr][[1]]) and the contour plot of the   (∇2𝒻 surface with the position of each maxima (that are minima in the  ∇2𝒻surface):

Deconvul[ MinVal_, SizeMatr_]. It executes the deconvolution and visualizes the results. The deconvolution is executed with the   built-in function”NonlinearModelFit”. the model is:

z(x,y)=(∑_(i=1)^n)(((e^(-((((x-(Subscript[μ1, i])))^2)/(2 (Subscript[σ1, i]^2)))-((((y-(Subscript[μ2, i])))^2)/(2 (Subscript[σ2, i]^2))))) (Subscript[ϕ, i]))/(2 π (Subscript[σ1, i]) (Subscript[σ2, i]))).

Where μ1i and μ2iare obtained with MinSearch[iter, SizeMatr] and the function ”NonlinearModelFit” allows to estimate the unknown parameters σ1i, μ2i (the variance) and ϕi (height) of each i peak. “ iter” is the “RandomSeed” defined previously and “MinVal” is MinSearch[iter][[1]]. The only requirement imposed in this model fit is that selected peaks must have positive heights (ϕi>0). The output is a) the contour plot of the model z(x,y); b) the fit between z(x,y) and f(x,y); c) a  table that  summarize the statistics of each peak; d)  a grid of contour plots of each single sub-peak:

Instructions: First, evaluate the cell with the functions SURFACE[p_, SizeMatr_], MinSearch[iter_, SizeMatr_] and Deconvul[ iter_, MinVal_]. (Shift+Enter)

Sucessively the three functions can be executed separately:

Step 1: Execute SURFACE[p, SizeMatr] to generate the surface f(x,y) and its (∇^2)𝒻 (Shift+Enter)

Step 2: Execute MinSearch[iter, SizeMatr]. Search the local minima  (∇^2)𝒻 (Shift+Enter)

Step 3: Execute Deconvul[ iter, MinVal]. Deconvolution and visualize the  results. (Shift+Enter)

Andrea Butturini, Barcelona, Dep. of Ecology, Universitat de Barcelona. Avd. Diagonal 643 08028 Barcelona,  Spain

for any question contact: abutturini@ub.edu

Here start the code:

Needs["MultivariateStatistics`"];
SURFACE[p_,SizeMatr_]:=Module[{RR,Co,Cov,𝒟ist,Sum𝒟ist,m,RS,Wheights,SD,surfacePlots},

Co[]:=RandomReal[{1.25,SizeMatr-0.25}];
Cov[]:=RandomReal[{1,2}];
𝒟ist[n_]:=MultinormalDistribution[{Co[],Co[]},{{((Cov[])^2),0},{0,((Cov[])^2)}}];(*here, with Mathematica 8 the function "BinormalDistribution" can repalce the "MultinormalDistribution" *)
Sum𝒟ist[m_]:=Map[𝒟ist, Range[1,m]];

m=RandomInteger[{4,20}]; (*number of subsurfaces that assemble the surface f(x,y)*)
RS=RandomSample[Range[m],m];
Wheights=RS/Plus@@RS;
SD=Sum𝒟ist[m];
surf=Plus@@Map[ Transpose[Table[Wheights[[#]]PDF[SD[[#]],{x,y}], {x,1,SizeMatr},{y,1,SizeMatr}]]&, Range[1,Length[Wheights]]]//Rescale;(*here, with Mathematica 8 the function "MixtureDistribution" can be used *)
LAPLACIAN=ReplaceAll[LaplacianGaussianFilter[surf,1], x_?(#>0&)->0]; (*Ti is the Laplacian. The Gaussian filter is used. Radius is set 1. It can be changed*)

LaplacianPlot=ListContourPlot[( Rescale[ LAPLACIAN]-1),InterpolationOrder->1,Contours->10,ContourStyle->White,ColorFunction->"GrayTones", PlotRange->All , Frame->True, FrameStyle->Directive[Black,0],PlotLabel->Style["(∇^2)𝒻",12, Red]];


LPSurf=ListContourPlot[-surf, PlotRange->All,ColorFunction->"GrayTones",ContourStyle->White,InterpolationOrder->1,
PlotLabel->Style["f(x,y)",12, Red] ,Frame->True, FrameStyle->Directive[Black,0]];

surfacePlots=GraphicsGrid[{{LPSurf,LaplacianPlot}}, ImageSize->400];

Grid[{{Style["The Surface and its Laplacian (∇^2)𝒻 ", Red, 30,Bold]},{surfacePlots}}]

] (*here I create a ficticius surface..with some peaks and shoulders *)

MinSearch[iter_,SizeMatr_]:= Module[{ps0,MinSearchPlot},ps0=Module[{MATRINTERP,MI,PEAKSPOS, row, column},
{ row, column}=Dimensions[LAPLACIAN];
 
MATRINTERP=ListInterpolation[Transpose[LAPLACIAN], InterpolationOrder->2];
MI[{a_,b_}]:=MATRINTERP[a,b];
Off[NMinimize::"cvmit"];
PEAKSPOS=Module[{yyy0,iii,PPP,PPP2,eeee},
yyy0=Table[ NMinimize[{MATRINTERP[x,y],2≤y≤SizeMatr-1&&2≤x≤SizeMatr-1},{{y,1, column},{x,1, row}},Method->{"NelderMead",
"RandomSeed"->i}, WorkingPrecision->10],{i, iter}][[All,2]];


iii=Transpose[{yyy0[[All,1,2]],yyy0[[All,2,2]]}];
PPP=Union[Round[iii]];
PPP2=Table[Nearest[iii, PPP[[n]]], {n,Length[PPP]}];
eeee=Table[{Mean[PPP2[[n,All,1]]], Mean[PPP2[[n,All,2]]]}, {n, 1, Length[PPP2]}]
];
Map[Reverse,PEAKSPOS]
];

MinSearchPlot=Show[LaplacianPlot,ListPlot[ps0, PlotStyle->{Red, PointSize[0.04]}]];
{ps0, Grid[{{Style["The Laplacian (∇^2)𝒻 and Minima location", Red, 30,Bold]},{MinSearchPlot}}]}
](*here I search the local minima in the laplacian (∇^2)𝒻 surface with the NelderMead optitimization Algorithim*)

Deconvul[ MinVal_,SizeMatr_]:=Module[{AsymBidnormalDistrib,ABD,variables,PARAMINITIALS,MODEL,DS,SURFACEXY,fs, SOLCPpeaks,sim,SinglPeaks,ps,pw,COST,COSTNEGATIVESPOSITION,peak,CPpeaks,lpsim,lppoints,Volum,lm,SummaryResults,LP,NewCoord},

AsymBidnormalDistrib[{μ1_,μ2_, σ1_,σ2_,ϕ_}]:=ϕ PDF[MultinormalDistribution[{μ1,μ2},{{((σ1)^2),0},{0,((σ2)^2)}}],{x,y}];

ABD[n_]:=Module[{iu},
iu=ToString[n];
AsymBidnormalDistrib[{
Symbol["μ1"<>iu],Symbol["μ2"<>iu], Symbol["σ1"<>iu],Symbol["σ2"<>iu],  Symbol["ϕ"<>iu]}
]
];

variables[VALINIT_ ]:=Module[{ParDistrNorm },
ParDistrNorm=Module[{pars,ParsL, rrr},
rrr=Length[Transpose[VALINIT][[1]]];
pars[x_]:={Symbol["μ1"<>ToString[x] ],Symbol["μ2"<>ToString[x] ], Symbol["σ1"<>ToString[x] ], Symbol["σ2"<>ToString[x] ],Symbol["ϕ"<>ToString[x]]
};
ParsL = Sort[Flatten[Map[pars,  Range[rrr]]]];
Join[Transpose[{Take[ParsL , 2 rrr],Flatten[Transpose[VALINIT]]}],Take[ParsL , -3rrr]]
]];

COST[nPeaks_ ]:=Module[{pars,ParsL,rrr},
rrr=nPeaks;
pars[x_]:=Symbol["ϕ"<>ToString[x]];
ParsL = Sort[Flatten[Map[pars,  Range[rrr]]]]
];
 

ps=Map[Reverse,MinVal];
NewCoord=ps;

fs=Module[{Deconv},
Label[begin];
PARAMINITIALS=Union[Table[Symbol["μ1"<>ToString[n]]->NewCoord[[n,1]], {n,Length[NewCoord] }],
Table[Symbol["μ2"<>ToString[n]]->NewCoord[[n,2]], {n,Length[NewCoord] }]] ;


MODEL=Plus@@Map[ABD,Range[1,Length[NewCoord ]]]/. PARAMINITIALS;

Off[NonlinearModelFit::"sszero"];

DS=Dimensions[surf];
SURFACEXY=Flatten[Table[{x,y,surf[[x,y]]}, {x,1,DS[[1]],1},{y,1,DS[[2]],1}],1];

Deconv= NonlinearModelFit[SURFACEXY,MODEL,DeleteDuplicates[Drop[variables[NewCoord], Length[NewCoord]*2 ]],{x,y},MaxIterations->Automatic,PrecisionGoal->∞,ConfidenceLevel->.999,AccuracyGoal->1000] ;

(*Here start the loop to remove potential peaks with negative height*)
COSTNEGATIVESPOSITION=Module[
{yty,val,elem, rr},
yty=Position[Take[Deconv["BestFitParameters"],  -Length[NewCoord]][[All,2]], x_?(#<0&)] ;



val= COST[Length[NewCoord]][[Flatten[yty]]];

rr=Table[
elem= Drop[Characters[ToString[val[[v]]]] ,1] ;
If[Length[elem]<2, ToExpression[elem[[1]]],10ToExpression[elem[[-2]]] +ToExpression[elem[[-1]]] ], {v,1,Length[val]}];
Partition[rr,1]
];

 

If [Length[COSTNEGATIVESPOSITION]≥  1, Goto[follow], Goto[end]];
Label[follow];
NewCoord=Delete[ NewCoord,COSTNEGATIVESPOSITION] ;
 

Goto[begin];
Label[end];
Deconv
];

peak[n_]:=Table[ABD[n]/. Join[fs[ "BestFitParameters"],  PARAMINITIALS], {x,1,SizeMatr},{y,1,SizeMatr}];
CPpeaks[n_]:=ListContourPlot[peak[n], ContourShading->None, ContourStyle->Black,Frame->True, FrameStyle->Directive[Black,0], PlotRange->{0,surf//Max}, InterpolationOrder->1];
sim=Partition[fs[ "PredictedResponse"],SizeMatr];

pw=Show[ ListContourPlot[sim,ContourShading->None,ContourStyle->White],ListPlot[Map[Reverse,NewCoord], PlotStyle->{Red, PointSize[0.04]}],Frame->True, FrameStyle->Directive[Black,0]];
{lpsim,lppoints}={ ListContourPlot[-sim,InterpolationOrder->1, PlotRange->All, ColorFunction->"GrayTones",ContourStyle->White,PlotLabel->Style["Model z(x,y)",12, Red],Frame->True, FrameStyle->Directive[Black,0]],
 ListPlot[Map[Reverse,NewCoord] ,PlotMarkers->Style["+",Yellow, Bold,30],Frame->True, FrameStyle->Directive[Black,0] ]};


GridSinglPEaks=Module[{},
SinglPeaks[n_]:=Show[pw,CPpeaks[n],Frame->True, FrameStyle->Directive[Black,0]];
Grid[{{Style["THE SINGLE GAUSSIAN BELLS", Red, 30,Bold],SpanFromLeft},{GraphicsGrid[Partition[Table[SinglPeaks[n], {n,1,Length[ps]}],3,3,1,{}]],SpanFromLeft}}, Frame->All]
];


Volum=Table[Plus@@Flatten[peak[n]], {n,1,Length[NewCoord]}]/(Plus@@Flatten[sim]) 100;

lm=LinearModelFit[ Transpose[fs[{"Response","PredictedResponse"}]],x,x];


LP=ListPlot[Transpose[fs[{"Response","PredictedResponse"}]], AspectRatio->1, Frame->True, Frame->True, FrameLabel->{Style["Model z(x,y) ",16,Bold, Red],Style["Surface f(x,y) ",16, Bold,Red]},FrameStyle->Directive[Black,12], AxesOrigin->{0,0}, PlotStyle->{Black, PointSize[0.02]}];

SummaryResults=Grid[Module[{Phi,Labels,bf,RoY,RoX,MuX,MuY,SR,PHight},
Labels={"Peaks#", "PEak Volume (%)", "Peak High","μx","μy","σx","σy" };

bf=fs[ "BestFitParameters"];
RoY=Partition[bf,NewCoord//Length
][[-2,All,2]];
RoX=Partition[bf,NewCoord//Length
][[-3,All,2]];

Phi=Partition[bf,NewCoord//Length
][[-1,All,2]];

PHight=Phi/(2 Pi RoY RoX);

MuX=NewCoord[[All,1]];
MuY=NewCoord[[All,2]];
SR={Range[1,Length[Volum]], Volum,PHight,MuX,MuY,RoY,RoX}//Transpose;

Prepend[SR,Labels]
],
Frame->All,Background->{None,{Red}}];
 Grid[{{Style["RESULTS OF DECONVOLUTION", Red, 30,Bold],SpanFromLeft},{ LPSurf,Show[lpsim,lppoints],LP},{SummaryResults,SpanFromLeft},
{GridSinglPEaks,SpanFromLeft}}, Frame->All] 
 
] (*I execute the deconvolution an generate some plots..........*)

Here start the simulations:

SizeMatr=20;(*This parameter set the dimension of the matrix that geneate the f(x,y). It can be changed by the user*)
 SURFACE[p,SizeMatr]

Step 2: Search the local minima in the Laplacian

(*This parameter set the persistence of the search of local minima in the Laplacian (∇^2)𝒻. It can be changed by the user. Higher values best search and longer computational times. To get god fit Randomseed>50 is recomended*)

Randomseed=20;
ms=MinSearch[Randomseed,SizeMatr];
MinVal=ms[[1]];
ms[[2]]

Step 3 :execute the deconvolution and visualize the results

Deconvul[ MinVal,SizeMatr]