# RLM MODEL RESULTS

import numpy as np


def _shift_intercept(arr):
    """
    A convenience function to make the SAS covariance matrix
    compatible with statsmodels.rlm covariance
    """
    arr = np.asarray(arr)
    side = int(np.sqrt(len(arr)))
    return np.roll(np.roll(arr.reshape(side, side), -1, axis=1), -1, axis=0)


class Huber(object):
    """
    """
    def __init__(self):
        self.params = np.array([0.82937387, 0.92610818,
                                -0.12784916, -41.02653105])
        self.bse = np.array([0.11118035, 0.3034081, 0.12885259, 9.8073472])
        self.scale = 2.4407137948148447
        self.weights = np.array([
            1.,  1.,  0.7858871,  0.50494094,  1.,
            1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
            0.36814106])
        self.resid = np.array([
            3.05027584, -2.07757332,  4.17721071,
            6.50163171, -1.64615192, -2.57226011, -1.73127333, -0.73127333,
            -2.25476463,  0.48083217, 1.63147461,  1.42973363, -2.26346951,
            -0.78323693,  2.26646556, 0.88291808, -0.83307835,  0.06186577,
            0.26360675,  1.54306186, -8.91752986])
        self.df_model = 3
        self.df_resid = 17
        self.bcov_unscaled = np.array([
            [1.72887367e-03,  -3.47079127e-03,
             -6.79080082e-04, 2.73387119e-02],
            [-3.47079127e-03,   1.28754242e-02,   9.95952051e-07,
             -6.19611175e-02],
            [-6.79080082e-04,   9.95952051e-07,   2.32216722e-03,
             -1.59355028e-01],
            [2.73387119e-02,  -6.19611175e-02,  -1.59355028e-01,
             1.34527267e+01]])  # From R
        self.fittedvalues = np.array([
            38.94972416,  39.07757332,  32.82278929, 21.49836829,
            19.64615192,  20.57226011,  20.73127333,  20.73127333,
            17.25476463,  13.51916783,  12.36852539,  11.57026637,
            13.26346951,  12.78323693,   5.73353444,   6.11708192,
            8.83307835,   7.93813423,   8.73639325,  13.45693814,
            23.91752986])
        self.tvalues = np.array([7.45971657,  3.0523516,
                                 -0.99221261, -4.18324448])
        # from R this is equivalent to
        # self.res1.params/np.sqrt(np.diag(self.res1.bcov_scaled))

    # The below are taken from SAS
    huber_h1 = [
        95.8813, 0.19485, -0.44161, -1.13577, 0.1949, 0.01232,
        -0.02474, -0.00484, -0.4416, -0.02474, 0.09177, 0.00001, -1.1358,
        -0.00484, 0.00001, 0.01655]
    h1 = _shift_intercept(huber_h1)

    huber_h2 = [
        82.6191, 0.07942, -0.23915, -0.95604, 0.0794, 0.01427,
        -0.03013, -0.00344, -0.2392, -0.03013, 0.10391, -0.00166, -0.9560,
        -0.00344, -0.00166, 0.01392]
    h2 = _shift_intercept(huber_h2)

    huber_h3 = [
        70.1633, -0.04533, -0.00790, -0.78618, -0.0453, 0.01656,
        -0.03608, -0.00203, -0.0079, -0.03608,  0.11610, -0.00333, -0.7862,
        -0.00203, -0.00333,  0.01138]
    h3 = _shift_intercept(huber_h3)


class Hampel(object):
    """
    """
    def __init__(self):
        self.params = np.array([0.74108304, 1.22507934,
                                -0.14552506, -40.47473236])
        self.bse = np.array([0.13482596, 0.36793632, 0.1562567, 11.89315426])
        self.scale = 3.0882646556217064
        self.weights = np.array([
            1.,  1.,  1.,  1.,  1.,   1.,  1.,  1.,  1.,
            1., 1.,  1.,  1.,  1.,  1., 1.,  1.,  1.,  1.,  1., 0.80629719])
        self.resid = np.array([
            3.06267708, -2.08284798,  4.36377602,
            5.78635972, -1.7634816,
            -2.98856094, -2.34048993, -1.34048993, -3.02422878,  1.08249252,
            2.39221804,  2.47177232, -1.62645737, -0.25076107,  2.32088237,
            0.88430719, -1.37812296, -0.35944755, -0.43900184,  1.40555003,
            -7.65988702])
        self.df_model = 3
        self.df_resid = 17
        self.bcov_unscaled = np.array([
            [1.72887367e-03,  -3.47079127e-03,
             -6.79080082e-04, 2.73387119e-02],
            [-3.47079127e-03,   1.28754242e-02,   9.95952051e-07,
             -6.19611175e-02],
            [-6.79080082e-04,   9.95952051e-07,   2.32216722e-03,
             -1.59355028e-01],
            [2.73387119e-02,  -6.19611175e-02,  -1.59355028e-01,
             1.34527267e+01]])
        self.fittedvalues = np.array([
            38.93732292,  39.08284798,
            32.63622398,  22.21364028,
            19.7634816,  20.98856094,  21.34048993,  21.34048993,
            18.02422878,  12.91750748,  11.60778196,  10.52822768,
            12.62645737,  12.25076107,   5.67911763,   6.11569281,
            9.37812296,   8.35944755,   9.43900184,  13.59444997,
            22.65988702])
        self.tvalues = np.array([5.49659011,  3.32959607,
                                 -0.93132046, -3.40319578])

    hampel_h1 = [
        141.309,  0.28717, -0.65085, -1.67388, 0.287,  0.01816,
        -0.03646, -0.00713, -0.651, -0.03646,  0.13524,  0.00001, -1.674,
        -0.00713, 0.00001,  0.02439]
    h1 = _shift_intercept(hampel_h1)

    hampel_h2 = [
        135.248,  0.18207, -0.36884, -1.60217, 0.182, 0.02120,
        -0.04563, -0.00567, -0.369, -0.04563,  0.15860, -0.00290, -1.602,
        -0.00567, -0.00290, 0.02329]
    h2 = _shift_intercept(hampel_h2)

    hampel_h3 = [
        128.921,  0.05409, -0.02445, -1.52732, 0.054,  0.02514,
        -0.05732, -0.00392, -0.024, -0.05732,  0.18871, -0.00652, -1.527,
        -0.00392, -0.00652,  0.02212]
    h3 = _shift_intercept(hampel_h3)


class BiSquare(object):
    def __init__(self):
        self.params = np.array([0.9275471,   0.65073222,
                                -0.11233103, -42.28525369])
        self.bse = np.array([0.10805398,  0.29487634,  0.12522928, 9.5315672])
        self.scale = 2.2818858795649497
        self.weights = np.array([
            0.89283149,  0.88496132, 0.79040651,
            0.3358111,   0.94617358, 0.90040725,  0.96630596,  0.99729171,
            0.94968061,  0.99900087, 0.98959903,  0.99831448,  0.84731833,
            0.96455873,  0.91767906, 0.98724523,  0.99762848,  0.99694419,
            0.98650731,  0.95897484, 0.00222999])
        self.resid = np.array([
            2.50917802,   -2.60315301,   3.56070896,
            6.93256033,   -1.76597524,   -2.41670746,  -1.39345348, -.39345348,
            -1.70651907,  -0.23917521,   0.77180408,   0.31020526,
            -3.01451315,  -1.42960401,   2.19218084,   0.85518774,
            -0.36817892,   0.4181383,    0.87973712,   1.53911661,
            -10.43556344])
        self.df_model = 3
        self.df_resid = 17
        self.bcov_unscaled = np.array([
            [1.72887367e-03,  -3.47079127e-03,
             -6.79080082e-04, 2.73387119e-02],
            [-3.47079127e-03,   1.28754242e-02,   9.95952051e-07,
             -6.19611175e-02],
            [-6.79080082e-04,   9.95952051e-07,   2.32216722e-03,
             -1.59355028e-01],
            [2.73387119e-02,  -6.19611175e-02,  -1.59355028e-01,
             1.34527267e+01]])
        self.fittedvalues = np.array([
            39.49082198,  39.60315301,  33.43929104,
            21.06743967,
            19.76597524,  20.41670746,  20.39345348,  20.39345348,
            16.70651907,  14.23917521,  13.22819592,  12.68979474,
            14.01451315,  13.42960401,   5.80781916,   6.14481226,
            8.36817892,   7.5818617,    8.12026288,  13.46088339,
            25.43556344])
        self.tvalues = np.array([8.58410823,  2.20679698,
                                 -0.8970029, -4.43633799])

    bisquare_h1 = [
        90.3354,  0.18358, -0.41607, -1.07007, 0.1836, 0.01161,
        -0.02331, -0.00456, -0.4161, -0.02331,  0.08646, 0.00001, -1.0701,
        -0.00456, 0.00001,  0.01559]
    h1 = _shift_intercept(bisquare_h1)

    bisquare_h2 = [
        67.82521, 0.091288, -0.29038, -0.78124, 0.091288,
        0.013849, -0.02914, -0.00352, -0.29038, -0.02914, 0.101088, -0.001,
        -0.78124, -0.00352,   -0.001, 0.011766]
    h2 = _shift_intercept(bisquare_h2)

    bisquare_h3 = [
        48.8983, 0.000442, -0.15919, -0.53523, 0.000442,
        0.016113, -0.03461, -0.00259, -0.15919, -0.03461, 0.112728,
        -0.00164, -0.53523, -0.00259, -0.00164, 0.008414]
    h3 = _shift_intercept(bisquare_h3)


class Andrews(object):
    def __init__(self):
        self.params = [0.9282, 0.6492, -.1123, -42.2930]
        self.bse = [.1061, .2894, .1229, 9.3561]
        self.scale = 2.2801
        self.df_model = 3.
        self.df_resid = 17.
        # bcov_unscaled not defined because it is not given as part of SAS
        self.resid = [
            2.503338458, -2.608934536, 3.5548678338, 6.9333705014,
            -1.768179527, -2.417404513, -1.392991531, -0.392991531,
            -1.704759385, -0.244545418, 0.7659115325, 0.3028635237,
            -3.019999429, -1.434221475, 2.1912017882, 0.8543828047,
            -0.366664104, 0.4192468573, 0.8822948661, 1.5378731634,
            -10.44592783]

        self.sresids = [
            1.0979293816, -1.144242351, 1.5591155202, 3.040879735,
            -0.775498914, -1.06023995, -0.610946684, -0.172360612,
            -0.747683723, -0.107254214, 0.3359181307, 0.1328317233,
            -1.324529688, -0.629029563, 0.9610305856, 0.3747203984,
            -0.160813769, 0.1838758324, 0.3869622398, 0.6744897502,
            -4.581438458]

        self.weights = [
            0.8916509101, 0.8826581922, 0.7888664106, 0.3367252734,
            0.9450252405, 0.8987321912, 0.9656622, 0.9972406688,
            0.948837669, 0.9989310017, 0.9895434667, 0.998360628,
            0.8447116551, 0.9636222149, 0.916330067, 0.9869982597,
            0.9975977354, 0.9968600162, 0.9861384742, 0.9582432444, 0]

    def conf_int(self):  # method to be consistent with sm
        return [(0.7203, 1.1360), (.0819, 1.2165), (-.3532, .1287),
                (-60.6305, -23.9555)]

    andrews_h1 = [
        87.5357, 0.177891, -0.40318, -1.03691, 0.177891,  0.01125,
        -0.02258, -0.00442, -0.40318, -0.02258, 0.083779, 6.481E-6,
        -1.03691, -0.00442, 6.481E-6,  0.01511]
    h1 = _shift_intercept(andrews_h1)

    andrews_h2 = [
        66.50472,  0.10489,  -0.3246, -0.76664, 0.10489, 0.012786,
        -0.02651,  -0.0036, -0.3246, -0.02651,  0.09406, -0.00065,
        -0.76664,  -0.0036, -0.00065, 0.011567]
    h2 = _shift_intercept(andrews_h2)

    andrews_h3 = [
        48.62157, 0.034949, -0.24633, -0.53394, 0.034949, 0.014088,
        -0.02956, -0.00287, -0.24633, -0.02956, 0.100628, -0.00104,
        -0.53394, -0.00287, -0.00104, 0.008441]
    h3 = _shift_intercept(andrews_h3)


# RLM Results with Huber's Proposal 2
# Obtained from SAS

class HuberHuber(object):
    def __init__(self):
        self.h1 = [
            114.4936, 0.232675, -0.52734, -1.35624, 0.232675, 0.014714,
            -0.02954, -0.00578, -0.52734, -0.02954, 0.10958, 8.476E-6,
            -1.35624, -0.00578, 8.476E-6, 0.019764]
        self.h1 = _shift_intercept(self.h1)
        self.h2 = [
            103.2876, 0.152602, -0.33476, -1.22084, 0.152602, 0.016904,
            -0.03766, -0.00434, -0.33476, -0.03766, 0.132043, -0.00214,
            -1.22084, -0.00434, -0.00214, 0.017739]
        self.h2 = _shift_intercept(self.h2)
        self.h3 = [
            91.7544, 0.064027, -0.11379, -1.08249, 0.064027, 0.019509,
            -0.04702, -0.00278, -0.11379, -0.04702, 0.157872, -0.00462,
            -1.08249, -0.00278, -0.00462, 0.015677]
        self.h3 = _shift_intercept(self.h3)
        self.resid = [
            2.909155172, -2.225912162, 4.134132661, 6.163172632,
            -1.741815737, -2.789321552, -2.02642336, -1.02642336,
            -2.593402734, 0.698655, 1.914261011, 1.826699492, -2.031210331,
            -0.592975466, 2.306098648, 0.900896645, -1.037551854,
            -0.092080512, -0.004518993, 1.471737448, -8.498372406]
        self.sresids = [
            0.883018497, -0.675633129, 1.25483702, 1.870713355,
            -0.528694904, -0.84664529, -0.615082113, -0.311551209,
            -0.787177874, 0.212063383, 0.581037374, 0.554459746,
            -0.616535106, -0.179986379, 0.699972205, 0.273449972,
            -0.314929051, -0.027949281, -0.001371654, 0.446717797,
            -2.579518651]
        self.weights = [
            1, 1, 1, 0.718977066, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
            1, 1, 1, 1, 1, 0.52141511]

        self.params = (.7990, 1.0475, -0.1351, -41.0892)
        self.bse = (.1213, .3310, .1406, 10.7002)
        self.scale = 3.2946
        self.df_model = 3
        self.df_resid = 17

    def conf_int(self):  # method for consistency with sm
        return [(0.5612, 1.0367), (.3987, 1.6963),
                (-.4106, .1405), (-62.0611, -20.1172)]


class HampelHuber(object):
    def __init__(self):
        self.h1 = [
            147.4727, 0.299695, -0.67924, -1.7469, 0.299695, 0.018952,
            -0.03805, -0.00744, -0.67924, -0.03805, 0.141144, 0.000011,
            -1.7469, -0.00744, 0.000011, 0.025456]
        self.h1 = _shift_intercept(self.h1)
        self.h2 = [
            141.148, 0.190007, -0.38493, -1.67206, 0.190007, 0.02213,
            -0.04762, -0.00592, -0.38493, -0.04762, 0.165518, -0.00303,
            -1.67206, -0.00592, -0.00303, 0.024301]
        self.h2 = _shift_intercept(self.h2)
        self.h3 = [
            134.5444, 0.05645, -0.02552, -1.59394, 0.05645, 0.026232,
            -0.05982, -0.00409, -0.02552, -0.05982, 0.196946, -0.0068,
            -1.59394, -0.00409, -0.0068, 0.023083]
        self.h3 = _shift_intercept(self.h3)
        self.resid = [
            3.125725599, -2.022218392, 4.434082972, 5.753880172,
            -1.744479058, -2.995299443, -2.358455878, -1.358455878,
            -3.068281354, 1.150212629, 2.481708553, 2.584584946,
            -1.553899388, -0.177335865, 2.335744732, 0.891912757,
            -1.43012351, -0.394515569, -0.497391962, 1.407968887,
            -7.505098501]
        self.sresids = [
            0.952186413, -0.616026205, 1.350749906, 1.752798302,
            -0.531418771, -0.912454834, -0.718453867, -0.413824947,
            -0.934687235, 0.350388031, 0.756000196, 0.787339321,
            -0.473362692, -0.054021633, 0.711535395, 0.27170242,
            -0.43565698, -0.120180852, -0.151519976, 0.428908041,
            -2.28627005]
        self.weights = [
            1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
            1, 1, 1, 0.874787298]

        self.params = (.7318, 1.2508, -0.1479, -40.2712)
        self.bse = (.1377, .3757, .1596, 12.1438)
        self.scale = 3.2827
        self.df_model = 3
        self.df_resid = 17

    def conf_int(self):
        return [(0.4619, 1.0016), (.5145, 1.9872),
                (-.4607, .1648), (-64.0727, -16.4697)]


class BisquareHuber(object):
    def __init__(self):