#! /usr/bin/env python
import os,sys,math,numpy
sys.path.append(os.environ['FIASCO'])
from fiasco_numpy import *

options= [("GLM_COMPLEX",GLM_COMPLEX),
          ("GLM_RESIDUALS",GLM_RESIDUALS),
          ("GLM_VARIANCES",GLM_VARIANCES),
          ("GLM_COVARIANCES",GLM_COVARIANCES),
          ("GLM_SSQR",GLM_SSQR),
          ("GLM_ORTHO",GLM_ORTHO),
          ("GLM_DEVIANCE",GLM_DEVIANCE)]

def testRegressor( r, obs, counts, factors,
                   rBeta, rStdErr, rDeviance, rDevianceResiduals ):
    nfactors= factors.shape[1]
    nobs= len(obs)

    print "No options set: nparam is %d"%r.n_params(nfactors)
    for (optName,opt) in options:
        if r.is_settable(opt):
            print "Setting %s"%optName
            r.set(opt,1)
        else:
            print "Option %s is not settable"%optName
    print "All available options set: nparam is %d"%r.n_params(nfactors)

    params= numpy.zeros(r.n_params(nfactors))
    r.fit(obs, factors, counts, params)
    betas= params[0:nfactors]
    print "betas out: %s"%str(betas)
    if r.get(GLM_COVARIANCES):
        covElements= []
        offset= nfactors
        for i in xrange(nfactors):
            covElements.append(params[offset:offset+nfactors])
            offset += nfactors
        covariance= numpy.matrix(covElements)
        print "covariance:"
        print covariance
        print "standard errors:"
        for i in xrange(nfactors):
            print "beta= %g has std err %g"%(betas[i],math.sqrt(covariance[i,i]))
            print "vs. R beta %g with std err %g"%(rBeta[i],rStdErr[i])
    if r.get(GLM_DEVIANCE):
        deviance= params[offset]
        print "deviance: %g vs. %g reported by R"%(deviance,rDeviance)
    if r.get(GLM_RESIDUALS):
        print "deviance residuals: %s"%str(obs)
        print "vs. R deviance residuals: %s"%str(rDevianceResiduals)

# Blindness-by-age example from 'Intro to R'
# Notes: 
print """
###########################
# Blindness-by-age example from Intro to R
###########################
"""
obs= numpy.array([6.,17.,26.,37.,44.],dtype=numpy.float)
counts= numpy.array([50.,50.,50.,50.,50.],dtype=numpy.float)
factorsT= numpy.matrix([[1.,1.,1.,1.,1.],
                       [20.,35.,45.,55.,70.]],dtype=numpy.float)
factors= factorsT.transpose()
rBeta= [-3.53778, 0.08114]
rStdErr= [0.50232, 0.01082]
rDeviance= 0.317066
rDevianceResiduals= [-0.1797,   0.1157,  -0.1182,   0.3791,  -0.3372]


r= create_logistic_regressor();
#r.set(GLM_DEBUG,1)

testRegressor(r,obs,counts,factors, rBeta, rStdErr, rDeviance,
              rDevianceResiduals)


# G. Rodrigues' birth control regression dataset
# Notes: this example treats age group as continuous- should be catagorical
#        Need correct R results
print """
###########################
# G. Rodrigues' birth control regression dataset
###########################
"""
obs= numpy.array([6., 4., 52., 10.,
                  14., 10., 54., 27.,
                  33., 80., 46., 78.,
                  6., 48., 8., 31.], # yes answers by catagory
                 dtype=numpy.float)
counts= numpy.array([59., 14., 264., 60.,
                     74., 29., 209., 92.,
                     145., 157., 164., 146.,
                     41., 94., 16., 43.], # counts by catagory
                    dtype=numpy.float)
# Factors are: const, member_age_group_1, member_age_group_2,
# member_age_group_3, education_upper, contraceptive_use
factors= numpy.matrix([[1.,0.,0.,0.],
                       [1.,0.,0.,1.],
                       [1.,0.,1.,0.],
                       [1.,0.,1.,1.],
                       [1.,1.,0.,0.],
                       [1.,1.,0.,1.],
                       [1.,1.,1.,0.],
                       [1.,1.,1.,1.],
                       [1.,2.,0.,0.],
                       [1.,2.,0.,1.],
                       [1.,2.,1.,0.],
                       [1.,2.,1.,1.],
                       [1.,3.,0.,0.],
                       [1.,3.,0.,1.],
                       [1.,3.,1.,0.],
                       [1.,3.,1.,1.]],dtype=numpy.float )
rBeta= [-1.95714,0.42272,0.32409,0.83182]
rStdErr= [0.15776, 0.06558, 0.12367, 0.11742 ]
rDeviance= 30.822 
rDevianceResiduals= [-0.5293361, 0.3478760, 1.4374130, -2.5384442,
                     0.2644181, 0.1549846, 0.9738833,  -2.2512210,
                     -0.5623718, 1.9903473, -0.8982826, 0.5609409,
                     -2.7372558, -0.4857580, 0.7277568, 1.4640499]


r= create_logistic_regressor();
#r.set(GLM_DEBUG,1)

testRegressor(r,obs,counts,factors, rBeta, rStdErr, rDeviance,
              rDevianceResiduals)


# From R 'example(glm)'
# Notes- why doesn't counts=None work?
print """
###########################
# Poisson example from R's 'example(glm)'
###########################
"""

obs= numpy.array([18., 17., 15., 20., 10., 20., 25., 13., 12.],
                 dtype= numpy.double)
counts= None
counts= numpy.array(9*[0.0],dtype=numpy.double)
# factors are treatment and outcome, treated as catagorical
# So the X matrix which follows is: constant, outcome2, outcome3,
# treatment2, treatment3
factorsT= numpy.matrix([[1.,1.,1.,1.,1.,1.,1.,1.,1.],
                        [0.,1.,0.,0.,1.,0.,0.,1.,0.],
                        [0.,0.,1.,0.,0.,1.,0.,0.,1.],
                        [0.,0.,0.,1.,1.,1.,0.,0.,0.],
                        [0.,0.,0.,0.,0.,0.,1.,1.,1.]])
factors= factorsT.transpose()
rBeta= [ 3.045e+00, -4.543e-01, -2.930e-01, 8.717e-16, 4.557e-16]
rStdErr= [1.709e-01, 2.022e-01, 1.927e-01, 2.000e-01, 2.000e-01]
rDeviance= 5.1291
rDevianceResiduals= [-0.67125, 0.96272,-0.16965,-0.21999,-0.95552, 1.04939,
                     0.84715,-0.09167,-0.96656]
r= create_poisson_regressor();
r.set(GLM_DEBUG,1)

testRegressor(r,obs,counts,factors, rBeta, rStdErr, rDeviance,
              rDevianceResiduals)

