#include <math.h>
#include "mex.h"


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
	/* Variables */
    int n,n1,n2,s,s1,s2,nodeNum,e,Vind,edgeNum,
            nNodes,nFree,nEdges,maxState,sizeEdgePot[3],nInducedEdges,
            *nStates,*edgeEnds,*V,*E,*clamped;
    double *nodePot, *edgePot,*newNodePot,*nodeMap,*newNstates,*newEdgePot,*newEdgeEnds,*edgeMap;
	
	/* Inputs */
    nodePot = mxGetPr(prhs[0]);
    edgePot = mxGetPr(prhs[1]);
    nStates = (int*)mxGetPr(prhs[2]);
    edgeEnds = (int*)mxGetPr(prhs[3]);
    V = (int*)mxGetPr(prhs[4]);
    E = (int*)mxGetPr(prhs[5]);
    clamped = (int*)mxGetPr(prhs[6]);
    
    /* Sizes */
    nNodes = mxGetDimensions(prhs[0])[0];
    maxState = mxGetDimensions(prhs[0])[1];
    nEdges = mxGetDimensions(prhs[3])[0];
    
    /* Compute size of induced sub-graph */
    nFree = 0;
    for(n=0;n < nNodes;n++) {
        if(clamped[n]==-1)
            nFree++;
    }
    nInducedEdges = 0;
    for(e=0;e < nEdges; e++) {
        n1 = edgeEnds[e];
        n2 = edgeEnds[e+nEdges];
        if(clamped[n1]==-1 && clamped[n2]==-1)
            nInducedEdges++;
    }
    sizeEdgePot[0] = maxState;
    sizeEdgePot[1] = maxState;
    sizeEdgePot[2] = nInducedEdges;
        
    /* Output */
    plhs[0] = mxCreateDoubleMatrix(nFree,maxState,mxREAL);
    plhs[1] = mxCreateDoubleMatrix(nNodes,1,mxREAL);
    plhs[2] = mxCreateDoubleMatrix(nFree,1,mxREAL);
    plhs[3] = mxCreateNumericArray(3,sizeEdgePot,mxDOUBLE_CLASS,mxREAL);
	plhs[4] = mxCreateDoubleMatrix(nInducedEdges,2,mxREAL);
	plhs[5] = mxCreateDoubleMatrix(nEdges,1,mxREAL);
    newNodePot = mxGetPr(plhs[0]);
    nodeMap = mxGetPr(plhs[1]);
    newNstates = mxGetPr(plhs[2]);
    newEdgePot = mxGetPr(plhs[3]);
	newEdgeEnds = mxGetPr(plhs[4]);
    edgeMap = mxGetPr(plhs[5]);
	
    nodeNum = 0;
    for(n=0;n < nNodes;n++) {
        if(clamped[n]==-1) {
            
            /* Grab node potential */
            for(s=0;s < nStates[n];s++) {
                newNodePot[nodeNum + nFree*s] = nodePot[n + nNodes*s];
            }
            
            /* Absorb edge potentials from clamped neighbors */
            for(Vind = V[n]; Vind < V[n+1]; Vind++) {
                e = E[Vind];
                n1 = edgeEnds[e];
                n2 = edgeEnds[e+nEdges];
                
                if(n==n1) {
                    if(clamped[n2]!=-1) {
                        for(s=0;s<nStates[n];s++) {
                            newNodePot[nodeNum + nFree*s] *= edgePot[s+maxState*(clamped[n2] + maxState*e)];
                        }
                    }
                }
                else {
                    if(clamped[n1]!=-1) {
                        for(s=0;s<nStates[n];s++) {
                            newNodePot[nodeNum + nFree*s] *= edgePot[clamped[n1]+maxState*(s + maxState*e)];
                        }
                    }
                }
            }
            nodeMap[n] = nodeNum+1;
            newNstates[nodeNum] = nStates[n];
            nodeNum++;
        }
    }
    
    /* Grab edges in subgraph */
    edgeNum = 0;
    for(e=0;e < nEdges; e++) {
        n1 = edgeEnds[e];
        n2 = edgeEnds[e+nEdges];
        if(clamped[n1]==-1 && clamped[n2]==-1) {
			for(s1=0;s1<nStates[n1];s1++) {
                for(s2=0;s2<nStates[n2];s2++) {
                    newEdgePot[s1+maxState*(s2+maxState*edgeNum)] = edgePot[s1+maxState*(s2+maxState*e)];
                }
            }
			newEdgeEnds[edgeNum] = nodeMap[n1];
			newEdgeEnds[edgeNum+nInducedEdges] = nodeMap[n2];
			edgeMap[e] = edgeNum+1;
            edgeNum++;
        }
    }
}
