def plotting(file, ext):
    ts = pd.read_csv(file, header=[0, 1], index_col=[0])
    ts.index = [dt.datetime.strptime(x, "%H:%M:%S").time() for x in ts.index]
    ts.columns.set_levels(
        [dt.datetime.strptime(x, "%Y-%m-%d").date() for x in ts.columns.levels[0].values], 0, inplace=True
    )

    # gets used in the smoothing
    xtime = [int(x.hour) + int(x.minute) / 60 for x in ts.index]
    timeind = pd.date_range("00:00", "23:59", freq="min").to_pydatetime()  ###Very important!

    # for each month in the data, construct a df of just that month
    for mo in set([x.month for x in ts.columns.levels[0]]):
        monthnum = mo
        motxt = time.strftime("%B", time.strptime(str(monthnum), "%m"))
        mocheck = [x.month == monthnum for x in ts.columns.levels[0]]
        moflag = []
        for ans in mocheck:
            moflag.append(ans)
            moflag.append(ans)

        onemo = ts.loc[:, moflag]

        sleepy = onemo.xs("Sleep", level=1, axis=1).sum(axis=1) / (len(onemo.columns) / 2)
        eaty = onemo.xs("Eat", level=1, axis=1).sum(axis=1) / (len(onemo.columns) / 2)

        # begin plotting
        fig = plt.figure(figsize=(18, 6))
        ax = fig.add_subplot(111)

        ####Plot Sleep
        filtereds = lowess(sleepy, xtime, is_sorted=True, frac=0.025, it=0)
        ax.plot(timeind, filtereds[:, 1], "b", linewidth=2, label="Sleeping")
        ax.fill_between(timeind, 0, filtereds[:, 1], alpha=0.3, facecolor="b")
        # ax.plot(ts.index,sleepy,'b',linewidth=2,label='Sleeping')#raw data, not smoothed
        # ax.fill_between(ts.index, 0, sleepy,alpha=0.3,facecolor='b')

        ####Plot Eat
        filterede = lowess(eaty, xtime, is_sorted=True, frac=0.025, it=0)
        ax.plot(timeind, filterede[:, 1], "orange", linewidth=2, label="Eating")
        ax.fill_between(timeind, 0, filterede[:, 1], alpha=0.3, facecolor="orange")
        # ax.plot(ts.index,eaty,'orange',linewidth=2,label='Eating')
        # ax.fill_between(ts.index, 0, eaty,alpha=0.3,facecolor='orange')

        ####Axis formatting
        xax = ax.get_xaxis()
        xax.set_major_locator(mdates.HourLocator(byhour=range(0, 24, 2)))
        xax.set_major_formatter(mdates.DateFormatter("%H:%M"))
        ax.set_title("Activity Fraction at a Given Time of Day", fontsize="xx-large")
        ax.text("16:00", 0.9, motxt, fontsize="xx-large", color="k", fontweight="bold")
        ax.legend(fontsize="x-large")
        ax.set_ylim(0, 1.1)
        fig.autofmt_xdate()
        filename = "b2_TimeSeries/Activity_" + str(monthnum) + "." + ext
        fig.savefig(filename)
    return

def loess(x,y,frac=0.2,it=None,scatter=True):
    from statsmodels.nonparametric.smoothers_lowess import lowess 
    y = np.array(y)
    x = np.array(x)
    y = y[x.argsort()] # Sort y according to order of x.  
    x.sort() # Sort x in place.  
    if it is not None: # Helps if you are getting NaN's in the output.  
        d = lowess(y,x,frac=frac,it=it)
    else:
        d = lowess(y,x,frac=frac)
    return d

    def test_frac(self):
        rfile = os.path.join(rpath, "test_lowess_frac.csv")
        test_data = np.genfromtxt(open(rfile, "rb"), delimiter=",", names=True)

        expected_lowess_23 = np.array([test_data["x"], test_data["out_2_3"]]).T
        expected_lowess_15 = np.array([test_data["x"], test_data["out_1_5"]]).T

        actual_lowess_23 = lowess(test_data["y"], test_data["x"], frac=2.0 / 3)
        actual_lowess_15 = lowess(test_data["y"], test_data["x"], frac=1.0 / 5)

        assert_almost_equal(expected_lowess_23, actual_lowess_23, decimal=testdec - 1)
        assert_almost_equal(expected_lowess_15, actual_lowess_15, decimal=testdec)

    def test_iter(self):
        rfile = os.path.join(rpath, "test_lowess_iter.csv")
        test_data = np.genfromtxt(open(rfile, "rb"), delimiter=",", names=True)

        expected_lowess_no_iter = np.array([test_data["x"], test_data["out_0"]]).T
        expected_lowess_3_iter = np.array([test_data["x"], test_data["out_3"]]).T

        actual_lowess_no_iter = lowess(test_data["y"], test_data["x"], it=0)
        actual_lowess_3_iter = lowess(test_data["y"], test_data["x"], it=3)

        assert_almost_equal(expected_lowess_no_iter, actual_lowess_no_iter, decimal=testdec)
        assert_almost_equal(expected_lowess_3_iter, actual_lowess_3_iter, decimal=testdec)

    def test_frac(self):
        rfile = os.path.join(rpath, 'test_lowess_frac.csv')
        test_data = np.genfromtxt(open(rfile, 'rb'),
                                  delimiter = ',', names = True)

        expected_lowess_23 = np.array([test_data['x'], test_data['out_2_3']]).T
        expected_lowess_15 = np.array([test_data['x'], test_data['out_1_5']]).T

        actual_lowess_23 = lowess(test_data['y'], test_data['x'] ,frac = 2./3)
        actual_lowess_15 = lowess(test_data['y'], test_data['x'] ,frac = 1./5)

        assert_almost_equal(expected_lowess_23, actual_lowess_23, decimal = testdec-1)
        assert_almost_equal(expected_lowess_15, actual_lowess_15, decimal = testdec)

    def test_iter(self):
        rfile = os.path.join(rpath, 'test_lowess_iter.csv')
        test_data = np.genfromtxt(open(rfile, 'rb'),
                                  delimiter = ',', names = True)

        expected_lowess_no_iter = np.array([test_data['x'], test_data['out_0']]).T
        expected_lowess_3_iter = np.array([test_data['x'], test_data['out_3']]).T

        actual_lowess_no_iter = lowess(test_data['y'], test_data['x'], it = 0)
        actual_lowess_3_iter = lowess(test_data['y'], test_data['x'], it = 3)

        assert_almost_equal(expected_lowess_no_iter, actual_lowess_no_iter, decimal = testdec)
        assert_almost_equal(expected_lowess_3_iter, actual_lowess_3_iter, decimal = testdec)

    def test_delta(self):
        rfile = os.path.join(rpath, "test_lowess_delta.csv")
        test_data = np.genfromtxt(open(rfile, "rb"), delimiter=",", names=True)

        expected_lowess_del0 = np.array([test_data["x"], test_data["out_0"]]).T
        expected_lowess_delRdef = np.array([test_data["x"], test_data["out_Rdef"]]).T
        expected_lowess_del1 = np.array([test_data["x"], test_data["out_1"]]).T

        actual_lowess_del0 = lowess(test_data["y"], test_data["x"], frac=0.1)
        actual_lowess_delRdef = lowess(test_data["y"], test_data["x"], frac=0.1, delta=0.01 * np.ptp(test_data["x"]))
        actual_lowess_del1 = lowess(test_data["y"], test_data["x"], frac=0.1, delta=1.0 + 1e-10)

        assert_almost_equal(expected_lowess_del0, actual_lowess_del0, decimal=testdec)
        assert_almost_equal(expected_lowess_delRdef, actual_lowess_delRdef, decimal=testdec)
        assert_almost_equal(expected_lowess_del1, actual_lowess_del1, decimal=10)  # testdec)

    def test_simple(self):
        x = np.arange(20, dtype='float32')
        #standard normal noise
        noise = np.array([-0.76741118, -0.30754369,
                            0.39950921, -0.46352422, -1.67081778,
                            0.6595567 ,  0.66367639, -2.04388585,
                            0.8123281 ,  1.45977518,
                            1.21428038,  1.29296866,  0.78028477,
                            -0.2402853 , -0.21721302,
                            0.24549405,  0.25987014, -0.90709034,
                            -1.45688216, -0.31780505])
        y = x + noise

        # R output
        out = [-0.6260344553, 0.565071712, 1.759627189,
                2.9579633258, 4.1560636154, 5.3473396937,
                6.522298218, 7.708159388, 8.8759055519,
                9.9409758603, 10.8981138458, 11.7851424728,
                12.6188717297, 13.4098497374, 14.1516996585,
                14.9180658147, 15.6956600199, 16.4783034134,
                17.2617441531, 18.0459201716]

        expected_lowess = np.array([x, out]).T

        actual_lowess = lowess(y,x)

        assert_almost_equal(expected_lowess, actual_lowess)

def add_lowess(ax, lines_idx=0, frac=.2, **lowess_kwargs):
    """
    Add Lowess line to a plot.

    Parameters
    ----------
    ax : matplotlib Axes instance
        The Axes to which to add the plot
    lines_idx : int
        This is the line on the existing plot to which you want to add
        a smoothed lowess line.
    frac : float
        The fraction of the points to use when doing the lowess fit.
    lowess_kwargs
        Additional keyword arguments are passes to lowess.

    Returns
    -------
    fig : matplotlib Figure instance
        The figure that holds the instance.
    """
    y0 = ax.get_lines()[lines_idx]._y
    x0 = ax.get_lines()[lines_idx]._x
    lres = lowess(y0, x0, frac=frac, **lowess_kwargs)
    ax.plot(lres[:, 0], lres[:, 1], 'r', lw=1.5)
    return ax.figure

def regression_plot(Z,X,band_names=None,visible_only=True,figsize=(12,7)):
    """
    Produce a figure with a plot for each image band that displays the
    relationship between depth and radiance and gives a visual representation
    of the regression carried out in the `slopes` and `regressions` methods.

    Notes
    -----
    This method doesn't come directly from Lyzenga 1978 but the author of this
    code found it helpful.

    Parameters
    ----------
    Z : np.ma.MaskedArray
        Array of depth values repeated for each band so that Z.shape==X.shape.
        The mask needs to be the same too so that Z.mask==X.mask for all the
        bands.
    X : np.ma.MaskedArray
        The array of log transformed radiance values from equation B1 of
        Lyzenga 1978.

    Returns
    -------
    figure
        A matplotlib figure.
    """
    if band_names is None:
        band_names = ['Band'+str(i+1) for i in range(X.shape[-1])]
    nbands = X.shape[-1]
    if np.atleast_3d(Z).shape[-1] == 1:
        Z = np.repeat(np.atleast_3d(Z), nbands, 2)
    if visible_only:
        fig, axs = plt.subplots( 2, 3, figsize=figsize)
    else:
        fig, axs = plt.subplots( 2, 4, figsize=figsize )
    regs = regressions(Z,X)
    for i, ax in enumerate(axs.flatten()):
        if i > nbands-1:
            continue
        slp, incpt, rval = regs[:,i]
        # print X.shape, Z.shape
        x, y = equalize_array_masks(Z[...,i], X[...,i])
        if x.count() < 2:
            continue
        x, y = x.compressed(), y.compressed()
        # print "i = {}, x.shape = {}, y.shape = {}".format(i, x.shape, y.shape)
        ax.scatter( x, y, alpha=0.1, edgecolor='none', c='gold' )
        smth = lowess(y,x,frac=0.2)
        # ax.plot(smth.T[0],smth.T[1],c='black',alpha=0.5)
        ax.plot(smth.T[0],smth.T[1],c='black',alpha=0.5,linestyle='--')
        reglabel = "m=%.2f, r=%.2f" % (slp,rval)
        f = lambda x: incpt + slp * x
        ax.plot( x, f(x), c='brown', label=reglabel, alpha=1.0 )
        ax.set_title( band_names[i] )
        ax.set_xlabel( r'Depth (m)' )
        ax.set_ylabel( r'$X_i$' )
        ax.legend(fancybox=True, framealpha=0.5)
    plt.tight_layout()
    return fig

    def test_delta(self):
        rfile = os.path.join(rpath, 'test_lowess_delta.csv')
        test_data = np.genfromtxt(open(rfile, 'rb'),
                                  delimiter = ',', names = True)

        expected_lowess_del0 = np.array([test_data['x'], test_data['out_0']]).T
        expected_lowess_delRdef = np.array([test_data['x'], test_data['out_Rdef']]).T
        expected_lowess_del1 = np.array([test_data['x'], test_data['out_1']]).T

        actual_lowess_del0    = lowess(test_data['y'], test_data['x'], frac=0.1)
        actual_lowess_delRdef = lowess(test_data['y'], test_data['x'], frac=0.1,
                       delta = 0.01 * np.ptp(test_data['x']))
        actual_lowess_del1    = lowess(test_data['y'], test_data['x'], frac = 0.1, delta = 1.0 + 1e-10)

        assert_almost_equal(expected_lowess_del0, actual_lowess_del0, decimal = testdec)
        assert_almost_equal(expected_lowess_delRdef, actual_lowess_delRdef, decimal = testdec)
        assert_almost_equal(expected_lowess_del1, actual_lowess_del1, decimal = 10) #testdec)

    def test_simple(self):
        rfile = os.path.join(rpath, "test_lowess_simple.csv")
        test_data = np.genfromtxt(open(rfile, "rb"), delimiter=",", names=True)
        expected_lowess = np.array([test_data["x"], test_data["out"]]).T

        actual_lowess = lowess(test_data["y"], test_data["x"])

        assert_almost_equal(expected_lowess, actual_lowess, decimal=testdec)

    def test_simple(self):
        rfile = os.path.join(rpath, 'test_lowess_simple.csv')
        test_data = np.genfromtxt(open(rfile, 'rb'),
                                  delimiter = ',', names = True)
        expected_lowess = np.array([test_data['x'], test_data['out']]).T

        actual_lowess = lowess(test_data['y'], test_data['x'])

        assert_almost_equal(expected_lowess, actual_lowess, decimal = testdec)

 def generate(name, fname, x='x', y='y', out='out', kwargs=None, decimal=7):
     kwargs = {} if kwargs is None else kwargs
     data = np.genfromtxt(os.path.join(rpath, fname), delimiter=',', names=True)
     assert_almost_equal.description = name
     if callable(kwargs):
         kwargs = kwargs(data)
     result = lowess(data[y], data[x], **kwargs)
     expect = np.array([data[x], data[out]]).T
     assert_almost_equal(result, expect, decimal)

 def generate(name, fname,
              x='x', y='y', out='out', kwargs={}, decimal=7):
     data = np.genfromtxt(
         os.path.join(rpath, fname), delimiter=',', names=True)
     assert_equal_at_testdec = partial(
         assert_almost_equal, decimal=decimal)
     assert_equal_at_testdec.description = name
     if callable(kwargs):
         kwargs = kwargs(data)
     result = lowess(data[y], data[x], **kwargs)
     expect = np.array([data[x], data[out]]).T
     return assert_equal_at_testdec, result, expect

def stl_loess(ts, loess_frac=0.2):
	ts.OriginalReading.interpolate(inplace=True) # if there are NAs
	x_stl = sm.tsa.seasonal_decompose(ts.OriginalReading.values, freq=95)

	# Apply loess filter to get the trend
	trend_loess = lowess(ts.OriginalReading.values, ts.index, frac=loess_frac)[:,1]
	x_ts = ts.OriginalReading.values
	seasonal = x_stl.seasonal # the seasonality is fine
	x_ts -= seasonal # remove seasonality
	remainder = x_ts - trend_loess # remove trend, similar to x_stl.resid but complete

	return remainder

def DeTrend(TheCandData,TheMDI):
    '''' Fit a lowess to the data and remove the curve '''
    ''' default smoothing looks ok on test '''
    ''' Works on grids as well as vectors '''
    
    sizee=np.shape(TheCandData)
    if (len(sizee) < 2):
        gots=np.where(TheCandData > TheMDI)[0]
        los=lowess(TheCandData[gots],range(len(TheCandData[gots])))[:,1]
        TheCandData[gots]=TheCandData[gots]-los
    else:
        for ltt in range(len(TheCandData[0,:,0])):
	    for lnn in range(len(TheCandData[0,0,:])):
                gots=np.where(TheCandData[:,ltt,lnn] > TheMDI)[0]
		if (len(gots) > 12):
                    los=lowess(TheCandData[gots,ltt,lnn],range(len(TheCandData[gots,ltt,lnn])))[:,1]
                    TheCandData[gots,ltt,lnn]=TheCandData[gots,ltt,lnn]-los	 
		else:
		    TheCandData[:,ltt,lnn]=TheMDI           
    
    return TheCandData # DETREND

    def test_iter(self):

        x = np.arange(20, dtype='float32')
        #cauchy noise
        noise = np.array([ 1.86299605, -0.10816866,  1.87761229,
                        -3.63442237,  0.30249022,
                        1.03560416,  0.21163349,  1.14167809,
                        -0.00368175, -2.08808987,
                        0.13065417, -1.8052207 ,  0.60404596,
                        -2.30908204,  1.7081412 ,
                        -0.54633243, -0.93107948,  1.79023999,
                        1.05822445, -1.04530564])

        y = x + noise

        # R output
        out = [0.6264479483, 1.5008396363, 2.3861761926, 3.2716390242,
               4.1397266375, 4.9926614002, 5.9062225, 6.8541464784,
               7.8163358136, 8.6684661827, 9.5321215273, 10.4655376106,
               11.469691774, 12.612670578, 13.8080457514, 14.9355218409,
               16.0491183613, 17.1604998952, 18.2739171976, 19.3834268539]


        expected_lowess_no_iter = np.array([x, out]).T

        out = [1.1091939965, 1.9662338415, 2.8223436958, 3.6741660675,
               4.5153163696, 5.3483205165, 6.2127611584, 7.0371035909,
               7.8823844068, 8.7036783127, 9.5698728732, 10.5011237563,
               11.4924301926, 12.6180333554, 13.8056705213, 14.9280791108,
               16.0363681325, 17.1426206341, 18.2516511313, 19.3581200948]

        expected_lowess_3_iter = np.array([x, out]).T

        actual_lowess_no_iter = lowess(y,x,it=0)
        actual_lowess_3_iter = lowess(y,x,it=3)

        assert_almost_equal(expected_lowess_no_iter, actual_lowess_no_iter)
        assert_almost_equal(expected_lowess_3_iter, actual_lowess_3_iter)

def smoothEstimates(points):
    #First we take the points that are passed and turn them into x and y
    [x, y] = points
    #Since our x's are dates, we turn them into Epoch timestamps
    x = [time.mktime(p.timetuple()) for p in x]
    #We want the lowess model to take into account the closest 10 points when smoothing, meaning we pass it frac.
    frac = 10.0/len(x)
    #We then pass the points to a lowess smoother, which considers frac% of points at a time.
    smoothed = smooth.lowess(y, x, frac=frac, is_sorted=True)
    #We change the epoch timestamps back into datetime objects and assign that to x
    x = [datetime.datetime.fromtimestamp(p[0]) for p in smoothed]
    #We then assign the smoothed estimates to y
    y = [p[1] for p in smoothed]
    #Now we recombine the x and y and pass it back as one list.
    smoothed = [x, y]
    return smoothed

 def test_lowess(self):
     if skip_lowess:
         raise SkipTest
     frac = 0.5
     it = 1
     data = self.s.data.copy()
     for i in range(data.shape[0]):
         data[i, :] = lowess(
             endog=data[i, :],
             exog=self.s.axes_manager[-1].axis,
             frac=frac,
             it=it,
             is_sorted=True,
             return_sorted=False,)
     self.s.smooth_lowess(smoothing_parameter=frac,
                          number_of_iterations=it,)
     nose.tools.assert_true(np.allclose(data, self.s.data))