def group_log_means_predict_manual(df_train, df_test, vali):
    start_time = time.time()
    all_mean = df_train['log_Demanda_uni_equil'].mean()
    P_mean = df_train.groupby(by=['short_name'])['log_Demanda_uni_equil'].mean()
    C_mean = df_train.groupby(by=['Cliente_ID'])['log_Demanda_uni_equil'].mean()
    PA_mean = df_train.groupby(by=['short_name', 'Agencia_ID'])['log_Demanda_uni_equil'].mean()
    PR_mean = df_train.groupby(by=['short_name', 'Ruta_SAK'])['log_Demanda_uni_equil'].mean()
    PCA_mean = df_train.groupby(by=['short_name', 'Cliente_ID', 'Agencia_ID'])['log_Demanda_uni_equil'].mean()

    print 'mean calculating time=', time.time()-start_time

    start_time = time.time()
    if not vali:
        df_test['Demanda_uni_equil']=np.apply_along_axis((lambda x:log_means_pred_demand_manual_func(x,\
            P_mean, C_mean, PA_mean, PR_mean, PCA_mean, all_mean)), 1, df_test.values)
        df_test.to_csv('output/'+'manual_group_log_mean_'+ \
                str(datetime.datetime.now().strftime('%Y-%m-%d-%H-%M'))+'.csv', \
                columns=['id','Demanda_uni_equil'], index=False)
    else:
        # global pred_demand, true_demand
        pred_demand = np.apply_along_axis((lambda x:log_means_pred_demand_manual_func(x, \
            P_mean, C_mean, PA_mean, PR_mean, PCA_mean, all_mean)), 1, df_test[labels].values)
        true_demand = df_test['Demanda_uni_equil'].values
        RMSLE = np.sqrt(MSE(np.log1p(pred_demand), np.log1p(true_demand)))
        print 'RMSLE=', RMSLE

    print 'predicting time=', time.time()-start_time

    def run_cv(self, num_round, params):
        '''
        Using FoldTubeID split, loop over CV to get RMSLE for each split
        params is a list of parameters for XGBoost.

        After finishing CV, run score() to get the results
        '''
        self.pred = []
        self.real = []
        if len(params) == 0:
            raise ValueError('Please read in parameters')

        for tr, te in self.cv:
            self.train = self.trainset.loc[tr,:].copy()
            self.test = self.trainset.loc[te,:].copy()

            # Randomize and set seed
            # np.random.permutation(len(trainp1))
            np.random.seed(1)
            self.train = self.train.iloc[np.random.permutation(len(self.train))]
            np.random.seed(2)
            self.test = self.test.iloc[np.random.permutation(len(self.test))]
            y_real = np.array(self.test.iloc[:,-1])

            # Section for training multi-models if you like
            y_pred_xgb = xgboost_model(self.train, self.test, num_round, params)

            y_pred = y_pred_xgb
            self.pred += [y_pred]
            self.real += [y_real]
            self.rmsle_score += [np.sqrt(mean_squared_error(np.log1p(y_real),
                                 np.log1p(y_pred)))]
        print '==========================================================='
        print 'Finished Cross-validation'
        print '==========================================================='

def _gpinv(p, k, sigma):
    """Inverse Generalized Pareto distribution function"""
    x = np.full_like(p, np.nan)
    if sigma <= 0:
        return x
    ok = (p > 0) & (p < 1)
    if np.all(ok):
        if np.abs(k) < np.finfo(float).eps:
            x = - np.log1p(-p)
        else:
            x = np.expm1(-k * np.log1p(-p)) / k
        x *= sigma
    else:
        if np.abs(k) < np.finfo(float).eps:
            x[ok] = - np.log1p(-p[ok])
        else:
            x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
        x *= sigma
        x[p == 0] = 0
        if k >= 0:
            x[p == 1] = np.inf
        else:
            x[p == 1] = - sigma / k

    return x

def transform_log(series, robust=True):
    """Perform element-wise logarithm transformation in a numerical series.
    
    Parameters
    ----------
    series : pandas.Series
        series to transform
    robust : bool
        True - handle negative and zero values properly
        False - transform negative value to nan, zero to -inf

    Returns
    -------
    log_series : pandas.Series
        ANOTHER series consisting of the transformed values
    """
    # TODO: support log10
    # TODO: separate log1p and log explicitly
    if not isinstance(series, pd.Series):
        raise TypeError("argument 'series' is NOT 'pandas.Series' type")
    if not is_numerical_type(series):
        raise ValueError("value type of argument 'series' is NOT numerical")
    
    if robust:
        return series.apply(lambda x: np.log1p(x) if x>=0 else -np.log1p(-x))
    else:
        return series.apply(np.log)

    def pdf(self, x: Array, log=False):

        n, d = x.shape

        theta = self.params
        ok = valid_rows_in_u(x)
        log_pdf = np.repeat(np.nan, n)
        if not ok.any():
            return log_pdf
        elif theta == 0:
            log_pdf[ok] = 0
            return log_pdf

        lu = np.log(x).sum(1)
        t = self.ipsi(x).sum(1)

        if theta < 0:  # dim == 2
            pos_t = t < 1
            log_pdf = np.log1p(theta) - (1 + theta) * lu - (d + 1 / theta) * np.log1p(-t)
            log_pdf[~ok] = np.nan
            log_pdf[ok & ~pos_t] = -np.inf
        else:
            p = np.log1p(theta * np.arange(1, d)).sum()
            log_pdf = p - (1 + theta) * lu - (d + 1 / theta) * np.log1p(t)

        return log_pdf if log else np.exp(log_pdf)

 def transform(self, X):
   if self.columns:
     for column in self.columns:
       X[column] = np.log1p(X[column])
       return X
   else:
     return np.log1p(X)

def compute_weights(data, Nlive):
    """Returns log_ev, log_wts for the log-likelihood samples in data,
    assumed to be a result of nested sampling with Nlive live points."""

    start_data=concatenate(([float('-inf')], data[:-Nlive]))
    end_data=data[-Nlive:]

    log_wts=zeros(data.shape[0])

    log_vols_start=cumsum(ones(len(start_data)+1)*log1p(-1./Nlive))-log1p(-1./Nlive)
    log_vols_end=np.zeros(len(end_data))
    log_vols_end[-1]=np.NINF
    log_vols_end[0]=log_vols_start[-1]+np.log1p(-1.0/Nlive)
    for i in range(len(end_data)-1):
        log_vols_end[i+1]=log_vols_end[i]+np.log1p(-1.0/(Nlive-i))

    log_likes = concatenate((start_data,end_data,[end_data[-1]]))

    log_vols=concatenate((log_vols_start,log_vols_end))

    log_ev = log_integrate_log_trap(log_likes, log_vols)

    log_dXs = logsubexp(log_vols[:-1], log_vols[1:])
    log_wts = log_likes[1:-1] + log_dXs[:-1]

    log_wts -= log_ev

    return log_ev, log_wts

def _logistic(X, y, w):
    """Compute the logistic function of the data: sum(sigmoid(yXw))

    Parameters
    ----------
    X : ndarray, shape (n_samples, n_features)
        Design matrix.

    y : ndarray, shape (n_samples,)
        Target / response vector. Each entry must be +1 or -1.

    w : ndarray, shape (n_features,)
        Unmasked, ravelized input map.

    Returns
    -------
    energy : float
        Energy contribution due to logistic data-fit term.
    """

    z = np.dot(X, w[:-1]) + w[-1]
    yz = y * z
    idx = yz > 0
    out = np.empty_like(yz)
    out[idx] = np.log1p(np.exp(-yz[idx]))
    out[~idx] = -yz[~idx] + np.log1p(np.exp(yz[~idx]))
    out = out.sum()
    return out

    def __init__(self, past, future, features = None):
        """Create a training pattern.

        Parameters:
        past -- past feature vectors as a tensor of shape [P, V]
            where P is past days and V is the vectors/day
        future -- future feature vectors as a tensor of [F, V]
            where F is future days and V is the vectors/day
        features -- a sequence of feature names to use
            where None means use all features
        """

        # calculate training input from past features
        past_subfeatures = [[self._subfeatures(vector, features)
            for vector in vectors]
                for vectors in past]
        self._input = numpy.array(
            [list(util.flatten(vectors)) for vectors in past_subfeatures])

        # calculate training output from future volatility
        future_returns = numpy.log1p(
            [[vector.ret for vector in vectors] for vectors in future])
        self._output = numpy.std(future_returns, axis = 0, ddof = 1)\
            * numpy.sqrt(252)

        # calculate past returns for forecasts
        self._past_returns = numpy.log1p(
            [[vector.ret for vector in vectors] for vectors in past])

def my_logaddexp(a, b):
    tmp = a - b
    return np.select([a == b, tmp > 0, tmp <= 0], [
        a + 0.69314718055994529,
        a + np.log1p(np.exp(-tmp)),
        b + np.log1p(np.exp(tmp))
    ], default=tmp)

def log_likelihood_state(params,sender,time):
    #params = [theta,A,alpha,delta,epsilon,sigma]
    tol = 1e-24
    theta = float(params[0])
    alpha = float(params[2])
    if min(theta,alpha)<0:
        ll = -float('inf')
    else:
        (S,X,SX,m1,m2,N)= sufficient_statistics(sender,time)
        
        if theta < 0:
            theta = 0

        if 1 - theta + tol < 0:
            theta = 1   

        puu = alpha*np.log(theta+tol)
        pvv = alpha*np.log(1-theta+tol)
        puv = np.log1p(-np.exp(alpha*np.log(theta+tol)))
        pvu = np.log1p(-np.exp(alpha*np.log(1-theta+tol)))
        try:
            ll = (N[0]*puu+N[1]*puv+
                  N[2]*pvu+N[3]*pvv)
        except:
            print 'll error: theta = %s, alpha = %s'%(theta,alpha)
            ll=0
    return -ll #take negative for minimization

    def fit ( self , X , y ):
        N = len( y )

        # num of happy tweets
        N_1  = np.sum( y )
        # num of sad tweets
        N_0  = N - N_1

        # ratio of happy/sad tweet
        Pi_0 = ( N_0 + 2 / N )
        Pi_1 = ( N_1 + 2 / N )

        #output is an array, N_jc[0] is the count
        #of how many 'obamas' when happy/sad
        N_j0 = (1-y)*X
        N_j1 = y*X

        Theta_j0 = ( ( N_j0 + 1 ) / ( N_0 + 2 ) )
        Theta_j1 = ( ( N_j1 + 1 ) / ( N_1 + 2 ) )

        logpi = [ np.log( Pi_0 ), np.log( Pi_1 ) ]
        self.logpi = np.array( logpi ) 
        self.logtheta = np.array([ np.log( Theta_j0 ), np.log( Theta_j1 ) ])
        self.log1theta = np.array( [ np.log1p( -1*Theta_j0 ), np.log1p( -1*Theta_j0 ) ] )

        save_params( self, 'params' )

def stitch(record1, record2):
    seq1 = array([record1.seq.tostring()])
    seq2 = array([reverse_complement(record2.seq.tostring())])
    seq1.dtype = '|S1'
    seq2.dtype = '|S1'
    quals1 = array(record1.letter_annotations['phred_quality'])
    quals2 = array(record2.letter_annotations['phred_quality'][::-1])
    
    log10p_consensus_1 = log1p(-power(10, -quals1 / 10.)) / log(10)
    log10p_consensus_2 = log1p(-power(10, -quals2 / 10.)) / log(10)
    log10p_error_1 = -log10(3) - (quals1 / 10.)
    log10p_error_2 = -log10(3) - (quals2 / 10.)
    
    min_overlap = 1
    max_overlap = max(len(record1), len(record2))
    overlaps = {}
    for overlap in range(1, max_overlap):
        s1 = seq1[-overlap:]
        s2 = seq2[:overlap]
        q1 = quals1[-overlap:]
        q2 = quals2[:overlap]
        lpc1 = log10p_consensus_1[-overlap:]
        lpc2 = log10p_consensus_2[:overlap]
        lpe1 = log10p_error_1[-overlap:]
        lpe2 = log10p_error_2[:overlap]
        
        consensus = choose(q1 < q2, [s1, s2])
        score = sum(choose(consensus == s1, [lpe1, lpc1])) + sum(choose(consensus == s2, [lpe2, lpc2])) + len(consensus) * log10(4) * 2    # last term is null hypothesis, p=1/4
        consensus.dtype = '|S%i' % len(consensus)
        overlaps[overlap] = (consensus[0],score)
    
    return overlaps

    def __init__(self, daily_returns, benchmark_daily_returns, risk_free_rate, days, period=DAILY):
        assert(len(daily_returns) == len(benchmark_daily_returns))

        self._portfolio = daily_returns
        self._benchmark = benchmark_daily_returns
        self._risk_free_rate = risk_free_rate
        self._annual_factor = _annual_factor(period)
        self._daily_risk_free_rate = self._risk_free_rate / self._annual_factor

        self._alpha = None
        self._beta = None
        self._sharpe = None
        self._return = np.expm1(np.log1p(self._portfolio).sum())
        self._annual_return = (1 + self._return) ** (365 / days) - 1
        self._benchmark_return = np.expm1(np.log1p(self._benchmark).sum())
        self._benchmark_annual_return = (1 + self._benchmark_return) ** (365 / days) - 1
        self._max_drawdown = None
        self._volatility = None
        self._annual_volatility = None
        self._benchmark_volatility = None
        self._benchmark_annual_volatility = None
        self._information_ratio = None
        self._sortino = None
        self._tracking_error = None
        self._annual_tracking_error = None
        self._downside_risk = None
        self._annual_downside_risk = None
        self._calmar = None
        self._avg_excess_return = None

def logsum_pair(logx, logy):
    """
    Return log(x+y), avoiding arithmetic underflow/overflow.

    logx: log(x)
    logy: log(y)

    Rationale:

    x + y    = e^logx + e^logy
             = e^logx (1 + e^(logy-logx))
    log(x+y) = logx + log(1 + e^(logy-logx)) (1)

    Likewise,
    log(x+y) = logy + log(1 + e^(logx-logy)) (2)

    The computation of the exponential overflows earlier and is less precise
    for big values than for small values. Due to the presence of logy-logx
    (resp. logx-logy), (1) is preferred when logx > logy and (2) is preferred
    otherwise.
    """
    if logx == logzero():
        return logy
    elif logx > logy:
        return logx + np.log1p(np.exp(logy - logx))
    else:
        return logy + np.log1p(np.exp(logx - logy))

    def keras_cv(self, params):
        """
        Using FoldTubeID split, loop over CV to get RMSLE for each split
        params is a list of parameters for Keras Neural Networks.

        After finishing CV, run score() to get the results
        """
        self.pred = []
        self.real = []
        if len(params) == 0:
            raise ValueError("Please read in parameters")

        for tr, te in self.cv:
            self.train = self.trainset.loc[tr, :].copy()
            self.test = self.trainset.loc[te, :].copy()

            # Randomize and set seed
            # np.random.permutation(len(trainp1))
            np.random.seed(1)
            self.train = self.train.iloc[np.random.permutation(len(self.train))]
            np.random.seed(2)
            self.test = self.test.iloc[np.random.permutation(len(self.test))]
            y_real = np.array(self.test.iloc[:, -1])

            # Section for training multi-models if you like
            y_pred = keras_model(self.train, self.test, params)

            self.pred += [y_pred]
            self.real += [y_real]
            self.rmsle_score += [np.sqrt(mean_squared_error(np.log1p(y_real), np.log1p(y_pred)))]
        print "==========================================================="
        print "Finished Keras Cross-validation"
        print "==========================================================="

    def getval(self, keys):
        array = tuple(keys)

        array = np.unique(array)
        array_length = len(array)
        rsmleValues_array = []
        count = 0
        maxValue = np.amax(array)
        minValue = np.amin(array)



        for i in array:
            count = 0
            for j in array:
                count = count + (np.log1p(i) - np.log1p(j))**2
            rsmleValues_array.append(np.sqrt(count / array_length))

        count = -1;
        index = 0
        min_error_value = rsmleValues_array[0]
        for val in rsmleValues_array:
            count += 1
            if val < min_error_value:
                min_error_value = val
                index = count

        demand = array[index]
        return demand

def exp1():
	train,y,test,idx = get_data_1()
	train = np.log1p(train.astype(float))
	test = np.log1p(test.astype(float))
	scaler = StandardScaler().fit(train)
	train = scaler.transform(train)
	test = scaler.transform(test)
	mtrain = pd.read_csv('meta_features_train.csv')
	mtest = pd.read_csv('meta_features_test.csv')
	scaler2 = StandardScaler().fit(mtrain)
	mtrain = scaler2.transform(mtrain)
	mtest = scaler2.transform(mtest)
	train = np.column_stack((train,mtrain))
	test = np.column_stack((test,mtest))
	rtrain_nn,rtest_nn = nn_features(train,y,test,model=build_nn2,random_state=1,n_folds=5,early_stop=50)
	rtrain_nn_total = rtrain_nn
	rtest_nn_total = rtest_nn
	for i in range(9):
		rand_seed = i*113+9201
		rtrain_nn,rtest_nn = nn_features(train,y,test,model=build_nn2,random_state=rand_seed,n_folds=5,early_stop=50)
		rtrain_nn_total += rtrain_nn
		rtest_nn_total += rtest_nn
		pd.DataFrame(data=rtrain_nn_total).to_csv('rtrain_nn_last.csv',index=False)
		pd.DataFrame(data=rtest_nn_total).to_csv('rtest_nn_last.csv',index=False)
	
	pd.DataFrame(data=rtrain_nn_total/10).to_csv('rtrain_nn_final.csv',index=False)
	pd.DataFrame(data=rtest_nn_total/10).to_csv('rtest_nn_final.csv',index=False)

def ndcg_at_k(
    rating_true,
    rating_pred,
    col_user=DEFAULT_USER_COL,
    col_item=DEFAULT_ITEM_COL,
    col_rating=DEFAULT_RATING_COL,
    col_prediction=DEFAULT_PREDICTION_COL,
    relevancy_method="top_k",
    k=DEFAULT_K,
    threshold=DEFAULT_THRESHOLD,
):
    """Normalized Discounted Cumulative Gain (nDCG).
    
    Info: https://en.wikipedia.org/wiki/Discounted_cumulative_gain
    
    Args:
        rating_true (pd.DataFrame): True DataFrame
        rating_pred (pd.DataFrame): Predicted DataFrame
        col_user (str): column name for user
        col_item (str): column name for item
        col_rating (str): column name for rating
        col_prediction (str): column name for prediction
        relevancy_method (str): method for determining relevancy ['top_k', 'by_threshold']
        k (int): number of top k items per user
        threshold (float): threshold of top items per user (optional)

    Returns:
        float: nDCG at k (min=0, max=1).
    """

    df_hit, df_hit_count, n_users = merge_ranking_true_pred(
        rating_true=rating_true,
        rating_pred=rating_pred,
        col_user=col_user,
        col_item=col_item,
        col_rating=col_rating,
        col_prediction=col_prediction,
        relevancy_method=relevancy_method,
        k=k,
        threshold=threshold,
    )

    if df_hit.shape[0] == 0:
        return 0.0

    # calculate discounted gain for hit items
    df_dcg = df_hit.copy()
    # relevance in this case is always 1
    df_dcg["dcg"] = 1 / np.log1p(df_dcg["rank"])
    # sum up discount gained to get discount cumulative gain
    df_dcg = df_dcg.groupby(col_user, as_index=False).agg({"dcg": "sum"})
    # calculate ideal discounted cumulative gain
    df_ndcg = pd.merge(df_dcg, df_hit_count, on=[col_user])
    df_ndcg["idcg"] = df_ndcg["actual"].apply(
        lambda x: sum(1 / np.log1p(range(1, min(x, k) + 1)))
    )

    # DCG over IDCG is the normalized DCG
    return (df_ndcg["dcg"] / df_ndcg["idcg"]).sum() / n_users

def preprocess(df):
    df = create_datetime_features(df)
    df['period'] = df.datetime.map(calculate_period)
    if 'count' in df.columns:
        df['log_count'] = np.log1p(df['count'])
        df['log_registered'] = np.log1p(df['registered'])
        df['log_casual'] = np.log1p(df['casual'])
    return df