#include <euredit_sys.h>


int main (int argc, char *argv[])
{
    long          i, j, m = 8, num = 200;
    long          missing_row[1600], missing_column[1600];
    long          max_cycle, count,sx[8],info[2],nrepl;
    double        miss_val;
    double        rep_val[1600];
    double        x[200][8], mean[8], tol, tmp1;
    const char    file[] = "continuous.dat";
    FILE         *fp;

    
    printf ("Results for program em_ex.c \n\n");

    tol = 1.0e-6;
    max_cycle = 20;
	info[0] = 0;
	info[1] = 0;

    if ((fp = fopen (file,"r")) == NULL)
    {
        printf ("Error cannot open data file: %s.\n", file);
	    exit (2);
    }

    /* Read in the data including missing values*/
    for (i = 0; i<num; ++i)
    {
        for (j = 0; j<m; ++j)
        {
            fscanf (fp, "%lf ", &x[i][j]);
        }
    }
    fclose (fp);
    
    miss_val = -999.0;

    for (i=0; i<8; ++i) sx[i] = 1; /* include all the variables */

    for (j = 0; j<m; ++j)
    {
	    tmp1 = 0.0;
		count = 0;
        for (i=0; i<num; ++i)
        {
			if (x[i][j] != miss_val)
            {
                tmp1 += x[i][j];
                ++count;
			}
		}
		mean[j] = tmp1/(double)count;
    }

    printf ("\nThe mean variate values before calling EM \n");
    for (j=0; j<m; ++j)
    {
        printf ("%8.4f \n",mean[j]);
    }

    /* Call multivariate normal E-M component */
    em (num, m, &sx[0], &x[0][0], miss_val, tol, max_cycle, NULL,
        &nrepl, &missing_row[0], &missing_column[0], &rep_val[0], info);

	if (info[0] == -999)
	{
		printf (" Obtain a licence from NAG to run this component.\n");
		exit (2);
	}
	else if (info[0] > 0)
	{
		printf (" Component EM returns the value: %i.\n", info[0]);
		exit (2);
	}

    for (j = 0; j<m; ++j)
    {
	    tmp1 = 0.0;
		count = 0;
        
        for (i=0; i<num; ++i)
        {
            tmp1 += x[i][j];
		}
        
		mean[j] = tmp1/(double)num;
    }

    printf ("\n The mean variate values after calling EM: \n");
    for (j=0; j<m; ++j)
    {
        printf ("%8.4f \n", mean[j]);
    }

    return 0;
} 
