import numpy as np

def deseason_harmonic(dat, K, L):
    '''
    deseasoned_data, season, beta = deseason_harmonic(dat, K, L)

    Subtracts the seasonal cycle (season) from the data (data). Season
    is calculated by fitting K harmonics of the annual cycle (as well as the
    mean) to the data. Assumes on year has L elements (i.e., 365 for daily data,
    73 for pentad data, 52 for weekly data, etc.).
    Outputs the deseasonalized data, the season, and the fitting parameters (beta)

    Written by Eric Oliver, Dalhousie University, 2007-2011
    Adapted from original MATLAB script on 28 November 2012
    '''

#   ensure that mat is a matrix and a column vector
    dat = np.mat(dat)
    if dat.shape[1]!=0:
        dat = dat.T

#   length of time series and generate time vector
    n = len(dat)
    time = np.mat(np.arange(1,n+1)/(1.*L))

#   set up mean and harmonics to fit data
    P = np.mat(np.ones((n,1)))
    for k in range(1,K+1):
        P = np.concatenate((P, np.cos(k*2*np.pi*time.T)), 1)
        P = np.concatenate((P, np.sin(k*2*np.pi*time.T)), 1)

#   Remove seasonal cycle by harmonic regression
    beta = (np.linalg.inv(P.T*P)*P.T)*dat
    season = P*beta
    dat_ds = dat - season

    return dat_ds, season, beta