// version 1.0.0  14Feb2018
// Functions used in mycalcs.cpp
//
#include <stdio.h>
#include "stplugin.h"
#include "mymatrix.h"
#include "mycalcsv.h"

ST_retcode MyAve(MyMatrix &bmat, ST_int first, ST_int last, 
    ST_int nvars, ST_int *nobs)
{
    ST_int     i, j ;   
    ST_double  value ;
    ST_retcode rc ;
    char       msg[80] ;

    *nobs = (ST_int) 0 ;
    for(i=first-1; i<last; i++) { 
        if (!SF_ifobs(i+1)) {
            continue ;
        }
        (*nobs)++ ;
        for(j=0; j<nvars; j++) {
            rc = SF_vdata(j+1, i+1, &value);    
            if(rc>0) {
                sprintf(msg, "Problem accessing Stata data\n") ;
                SF_error(msg) ;
                return(rc) ;
            }
            if (SF_is_missing(value)) {
                sprintf(msg, "missing values encountered\n") ;
                SF_error(msg) ;
                return(static_cast<ST_retcode>(416)) ;
            }
            bmat.IncrementByValue(0, j, value) ;
        }
    }

    bmat.DivideByScalar( static_cast<ST_double>(*nobs) ) ; 
    return(rc) ;
}

ST_retcode MyV(MyMatrix &bmat, MyMatrix &vmat, ST_int first, 
    ST_int last, ST_int nvars, ST_int nobs) 
{

    ST_int      i, j, j2 ;
    ST_double  value ;
    char       msg[80] ;
    ST_retcode rc ;
   
    rc = (ST_retcode) 0 ;
    MyMatrix emat2(1, nvars) ;

    for(i=first-1; i<last; i++) {
        if (!SF_ifobs(i+1)) {
	    continue ;
	}
        for(j=0; j<nvars; j++) {
            rc = SF_vdata(j+1, i+1, &value);    
            if(rc>0) {
                sprintf(msg, "Problem accessing Stata data\n") ;
                SF_error(msg) ;
                return(rc) ;
            }
            if (SF_is_missing(value)) {
                sprintf(msg, "missing values encountered\n") ;
                SF_error(msg) ;
                return(static_cast<ST_retcode>(416)) ;
            }
            emat2.StoreValue(0, j, (bmat.GetValue(0,j)-value) ) ;
        }
        for(j=0; j<nvars; j++) {
// matrix is symmetric only fill in lower diag in loop
//     over observations
            for(j2=0; j2<=j; j2++) {
                vmat.IncrementByValue(j, j2, 
                    (emat2.GetValue(0,j)*emat2.GetValue(0,j2)) ) ;
            }
        }
    }
    

// Fill in above the diagonal
    for(j=0; j<nvars; j++) {
        for(j2=j+1; j2<nvars; j2++) {
            vmat.StoreValue(j, j2, vmat.GetValue(j2, j) ) ;
        }
    }

// Divide V by n*(n-1)
    vmat.DivideByScalar( (static_cast<ST_double> (nobs*(nobs-1))) ) ;

    return(rc) ;
}