import itertools
import numpy as np
from scipy import sparse
from scipy.stats import norm
from scipy.optimize import minimize, minimize_scalar
from scipy.sparse import csc_matrix, linalg as sla
from functools import partial
from collections import deque
from pygfl.solver import TrailSolver

class GaussianKnown:
    '''
    A simple Gaussian distribution with known mean and stdev.
    '''
    def __init__(self, mean, stdev):
        self.mean = mean
        self.stdev = stdev

    def pdf(self, data):
        return norm.pdf(data, loc=self.mean, scale=self.stdev)

    def sample(self):
        return np.random.normal(loc=self.mean, scale=self.stdev)

    def noisy_pdf(self, data):
        return norm.pdf(data, loc=self.mean, scale=np.sqrt(self.stdev**2 + 1))

    def __repr__(self):
        return 'N({:.2f}, {:.2f}^2)'.format(self.mean, self.stdev)


class SmoothedFdr(object):
    def __init__(self, signal_dist, null_dist, penalties_cross_x=None):
        self.signal_dist = signal_dist
        self.null_dist = null_dist

        if penalties_cross_x is None:
            self.penalties_cross_x = np.dot
        else:
            self.penalties_cross_x = penalties_cross_x

        self.w_iters = []
        self.beta_iters = []
        self.c_iters = []
        self.delta_iters = []

        # ''' Load the graph fused lasso library '''
        # graphfl_lib = cdll.LoadLibrary('libgraphfl.so')
        # self.graphfl_weight = graphfl_lib.graph_fused_lasso_weight_warm
        # self.graphfl_weight.restype = c_int
        # self.graphfl_weight.argtypes = [c_int, ndpointer(c_double, flags='C_CONTIGUOUS'), ndpointer(c_double, flags='C_CONTIGUOUS'),
        #                     c_int, ndpointer(c_int, flags='C_CONTIGUOUS'), ndpointer(c_int, flags='C_CONTIGUOUS'),
        #                     c_double, c_double, c_double, c_int, c_double,
        #                     ndpointer(c_double, flags='C_CONTIGUOUS'), ndpointer(c_double, flags='C_CONTIGUOUS'), ndpointer(c_double, flags='C_CONTIGUOUS')]
        self.solver = TrailSolver()

    def add_step(self, w, beta, c, delta):
        self.w_iters.append(w)
        self.beta_iters.append(beta)
        self.c_iters.append(c)
        self.delta_iters.append(delta)

    def finish(self):
        self.w_iters = np.array(self.w_iters)
        self.beta_iters = np.array(self.beta_iters)
        self.c_iters = np.array(self.c_iters)
        self.delta_iters = np.array(self.delta_iters)

    def reset(self):
        self.w_iters = []
        self.beta_iters = []
        self.c_iters = []
        self.delta_iters = []

    def solution_path(self, data, penalties, dof_tolerance=1e-4,
            min_lambda=0.20, max_lambda=1.5, lambda_bins=30,
            converge=0.00001, max_steps=100, m_converge=0.00001,
            m_max_steps=20, cd_converge=0.00001, cd_max_steps=1000, verbose=0, dual_solver='graph',
            admm_alpha=1., admm_inflate=2., admm_adaptive=False, initial_values=None,
            grid_data=None, grid_map=None):
        '''Follows the solution path of the generalized lasso to find the best lambda value.'''
        lambda_grid = np.exp(np.linspace(np.log(max_lambda), np.log(min_lambda), lambda_bins))
        aic_trace = np.zeros(lambda_grid.shape) # The AIC score for each lambda value
        aicc_trace = np.zeros(lambda_grid.shape) # The AICc score for each lambda value (correcting for finite sample size)
        bic_trace = np.zeros(lambda_grid.shape) # The BIC score for each lambda value
        dof_trace = np.zeros(lambda_grid.shape) # The degrees of freedom of each final solution
        log_likelihood_trace = np.zeros(lambda_grid.shape)
        beta_trace = []
        u_trace = []
        w_trace = []
        c_trace = []
        results_trace = []
        best_idx = None
        best_plateaus = None
        flat_data = data.flatten()
        edges = penalties[3] if dual_solver == 'graph' else None
        if grid_data is not None:
                grid_points = np.zeros(grid_data.shape)
                grid_points[:,:] = np.nan
        for i, _lambda in enumerate(lambda_grid):
            if verbose:
                print('#{0} Lambda = {1}'.format(i, _lambda))

            # Clear out all the info from the previous run
            self.reset()

            # Fit to the final values
            results = self.run(flat_data, penalties, _lambda=_lambda, converge=converge, max_steps=max_steps,
                           m_converge=m_converge, m_max_steps=m_max_steps, cd_converge=cd_converge,
                           cd_max_steps=cd_max_steps, verbose=verbose, dual_solver=dual_solver,
                           admm_alpha=admm_alpha, admm_inflate=admm_inflate, admm_adaptive=admm_adaptive,
                           initial_values=initial_values)

            if verbose:
                print('Calculating degrees of freedom')

            # Create a grid structure out of the vector of betas
            if grid_map is not None:
                grid_points[grid_map != -1] = results['beta'][grid_map[grid_map != -1]]
            else:
                grid_points = results['beta'].reshape(data.shape)

            # Count the number of free parameters in the grid (dof)
            plateaus = calc_plateaus(grid_points, dof_tolerance, edges=edges)
            dof_trace[i] = len(plateaus)
            #dof_trace[i] = (np.abs(penalties.dot(results['beta'])) >= dof_tolerance).sum() + 1 # Use the naive DoF

            if verbose:
                print('Calculating AIC')

            # Get the negative log-likelihood
            log_likelihood_trace[i] = -self._data_negative_log_likelihood(flat_data, results['c'])

            # Calculate AIC = 2k - 2ln(L)
            aic_trace[i] = 2. * dof_trace[i] - 2. * log_likelihood_trace[i]
            
            # Calculate AICc = AIC + 2k * (k+1) / (n - k - 1)
            aicc_trace[i] = aic_trace[i] + 2 * dof_trace[i] * (dof_trace[i]+1) / (flat_data.shape[0] - dof_trace[i] - 1.)

            # Calculate BIC = -2ln(L) + k * (ln(n) - ln(2pi))
            bic_trace[i] = -2 * log_likelihood_trace[i] + dof_trace[i] * (np.log(len(flat_data)) - np.log(2 * np.pi))

            # Track the best model thus far
            if best_idx is None or bic_trace[i] < bic_trace[best_idx]:
                best_idx = i
                best_plateaus = plateaus

            # Save the final run parameters to use for warm-starting the next iteration
            initial_values = results

            # Save the trace of all the resulting parameters
            beta_trace.append(results['beta'])
            u_trace.append(results['u'])
            w_trace.append(results['w'])
            c_trace.append(results['c'])

            if verbose:
                print('DoF: {0} AIC: {1} AICc: {2} BIC: {3}'.format(dof_trace[i], aic_trace[i], aicc_trace[i], bic_trace[i]))

        if verbose:
            print('Best setting (by BIC): lambda={0} [DoF: {1}, AIC: {2}, AICc: {3} BIC: {4}]'.format(lambda_grid[best_idx], dof_trace[best_idx], aic_trace[best_idx], aicc_trace[best_idx], bic_trace[best_idx]))

        return {'aic': aic_trace,
                'aicc': aicc_trace,
                'bic': bic_trace,
                'dof': dof_trace,
                'loglikelihood': log_likelihood_trace,
                'beta': np.array(beta_trace),
                'u': np.array(u_trace),
                'w': np.array(w_trace),
                'c': np.array(c_trace),
                'lambda': lambda_grid,
                'best': best_idx,
                'plateaus': best_plateaus}


    def run(self, data, penalties, _lambda=0.1, converge=0.00001, max_steps=100, m_converge=0.00001,
            m_max_steps=100, cd_converge=0.00001, cd_max_steps=100, verbose=0, dual_solver='graph',
            admm_alpha=1., admm_inflate=2., admm_adaptive=False, initial_values=None):
        '''Runs the Expectation-Maximization algorithm for the data with the given penalty matrix.'''
        delta = converge + 1
        
        if initial_values is None:
            beta = np.zeros(data.shape)
            prior_prob = np.exp(beta) / (1 + np.exp(beta))
            u = initial_values
        else:
            beta = initial_values['beta']
            prior_prob = initial_values['c']
            u = initial_values['u']

        prev_nll = 0
        cur_step = 0
        
        while delta > converge and cur_step < max_steps:
            if verbose:
                print('Step #{0}'.format(cur_step))

            if verbose:
                print('\tE-step...')

            # Get the likelihood weights vector (E-step)
            post_prob = self._e_step(data, prior_prob)

            if verbose:
                print('\tM-step...')

            # Find beta using an alternating Taylor approximation and convex optimization (M-step)
            beta, u = self._m_step(beta, prior_prob, post_prob, penalties, _lambda,
                                   m_converge, m_max_steps,
                                   cd_converge, cd_max_steps,
                                   verbose, dual_solver,
                                   admm_adaptive=admm_adaptive,
                                   admm_inflate=admm_inflate,
                                   admm_alpha=admm_alpha,
                                   u0=u)

            # Get the signal probabilities
            prior_prob = ilogit(beta)
            cur_nll = self._data_negative_log_likelihood(data, prior_prob)

            if dual_solver == 'admm':
                # Get the negative log-likelihood of the data given our new parameters
                cur_nll += _lambda * np.abs(u['r']).sum()
            
            # Track the change in log-likelihood to see if we've converged
            delta = np.abs(cur_nll - prev_nll) / (prev_nll + converge)

            if verbose:
                print('\tDelta: {0}'.format(delta))

            # Track the step
            self.add_step(post_prob, beta, prior_prob, delta)

            # Increment the step counter
            cur_step += 1

            # Update the negative log-likelihood tracker
            prev_nll = cur_nll

            # DEBUGGING
            if verbose:
                print('\tbeta: [{0:.4f}, {1:.4f}]'.format(beta.min(), beta.max()))
                print('\tprior_prob:    [{0:.4f}, {1:.4f}]'.format(prior_prob.min(), prior_prob.max()))
                print('\tpost_prob:    [{0:.4f}, {1:.4f}]'.format(post_prob.min(), post_prob.max()))
                if dual_solver != 'graph':
                    print('\tdegrees of freedom: {0}'.format((np.abs(penalties.dot(beta)) >= 1e-4).sum()))

        # Return the results of the run
        return {'beta': beta, 'u': u, 'w': post_prob, 'c': prior_prob}

    def _data_negative_log_likelihood(self, data, prior_prob):
        '''Calculate the negative log-likelihood of the data given the weights.'''
        signal_weight = prior_prob * self.signal_dist.pdf(data)
        null_weight = (1-prior_prob) * self.null_dist.pdf(data)
        return -np.log(signal_weight + null_weight).sum()

    def _e_step(self, data, prior_prob):
        '''Calculate the complete-data sufficient statistics (weights vector).'''
        signal_weight = prior_prob * self.signal_dist.pdf(data)
        null_weight = (1-prior_prob) * self.null_dist.pdf(data)
        post_prob = signal_weight / (signal_weight + null_weight)
        return post_prob

    def _m_step(self, beta, prior_prob, post_prob, penalties,
                _lambda, converge, max_steps,
                cd_converge, cd_max_steps,
                verbose, dual_solver, u0=None,
                admm_alpha=1., admm_inflate=2., admm_adaptive=False):
        '''
        Alternating Second-order Taylor-series expansion about the current iterate
        and coordinate descent to optimize Beta.
        '''
        prev_nll = self._m_log_likelihood(post_prob, beta)
        delta = converge + 1
        u = u0
        cur_step = 0
        while delta > converge and cur_step < max_steps:
            if verbose > 1:
                print('\t\tM-Step iteration #{0}'.format(cur_step))
                print('\t\tTaylor approximation...')

            # Cache the exponentiated beta
            exp_beta = np.exp(beta)

            # Form the parameters for our weighted least squares
            if dual_solver != 'admm' and dual_solver != 'graph':
                # weights is a diagonal matrix, represented as a vector for efficiency
                weights = 0.5 * exp_beta / (1 + exp_beta)**2
                y = (1+exp_beta)**2 * post_prob / exp_beta + beta - (1 + exp_beta)
                if verbose > 1:
                    print('\t\tForming dual...')
                x = np.sqrt(weights) * y
                A = (1. / np.sqrt(weights))[:,np.newaxis] * penalties.T
            else:
                weights = (prior_prob * (1 - prior_prob))
                y = beta - (prior_prob - post_prob) / weights
                print(weights)
                print(y)

            if dual_solver == 'cd':
                # Solve the dual via coordinate descent
                u = self._u_coord_descent(x, A, _lambda, cd_converge, cd_max_steps, verbose > 1, u0=u)
            elif dual_solver == 'sls':
                # Solve the dual via sequential least squares
                u = self._u_slsqp(x, A, _lambda, verbose > 1, u0=u)
            elif dual_solver == 'lbfgs':
                # Solve the dual via L-BFGS-B
                u = self._u_lbfgsb(x, A, _lambda, verbose > 1, u0=u)
            elif dual_solver == 'admm':
                # Solve the dual via alternating direction methods of multipliers
                #u = self._u_admm_1dfusedlasso(y, weights, _lambda, cd_converge, cd_max_steps, verbose > 1, initial_values=u)
                #u = self._u_admm(y, weights, _lambda, penalties, cd_converge, cd_max_steps, verbose > 1, initial_values=u)
                u = self._u_admm_lucache(y, weights, _lambda, penalties, cd_converge, cd_max_steps,
                                        verbose > 1, initial_values=u, inflate=admm_inflate,
                                        adaptive=admm_adaptive, alpha=admm_alpha)
                beta = u['x']
            elif dual_solver == 'graph':
                u = self._graph_fused_lasso(y, weights, _lambda, penalties[0], penalties[1], penalties[2], penalties[3], cd_converge, cd_max_steps, max(0, verbose - 1), admm_alpha, admm_inflate, initial_values=u)
                beta = u['beta']
                # if np.abs(beta).max() > 20:
                #     beta = np.clip(beta, -20, 20)
                #     u = None
            else:
                raise Exception('Unknown solver: {0}'.format(dual_solver))

            if dual_solver != 'admm' and dual_solver != 'graph':
                # Back out beta from the dual solution
                beta = y - (1. / weights) * penalties.T.dot(u)

            # Get the current log-likelihood
            cur_nll = self._m_log_likelihood(post_prob, beta)

            # Track the convergence
            delta = np.abs(prev_nll - cur_nll) / (prev_nll + converge)

            if verbose > 1:
                print('\t\tM-step delta: {0}'.format(delta))

            # Increment the step counter
            cur_step += 1

            # Update the negative log-likelihood tracker
            prev_nll = cur_nll

        return beta, u
    
    def _m_log_likelihood(self, post_prob, beta):
        '''Calculate the log-likelihood of the betas given the weights and data.'''
        return (np.log(1 + np.exp(beta)) - post_prob * beta).sum()

    def _graph_fused_lasso(self, y, weights, _lambda, ntrails, trails, breakpoints, edges, converge, max_steps, verbose, alpha, inflate, initial_values=None):
        '''Solve for u using a super fast graph fused lasso library that has an optimized ADMM routine.'''
        if verbose:
            print('\t\tSolving via Graph Fused Lasso')
        # if initial_values is None:
        #     beta = np.zeros(y.shape, dtype='double')
        #     z = np.zeros(breakpoints[-1], dtype='double')
        #     u = np.zeros(breakpoints[-1], dtype='double')
        # else:
        #     beta = initial_values['beta']
        #     z = initial_values['z']
        #     u = initial_values['u']
        # n = y.shape[0]
        # self.graphfl_weight(n, y, weights, ntrails, trails, breakpoints, _lambda, alpha, inflate, max_steps, converge, beta, z, u)
        # return {'beta': beta, 'z': z, 'u': u }
        self.solver.alpha = alpha
        self.solver.inflate = inflate
        self.solver.maxsteps = max_steps
        self.solver.converge = converge
        self.solver.set_data(y, edges, ntrails, trails, breakpoints, weights=weights)
        if initial_values is not None:
            self.solver.beta = initial_values['beta']
            self.solver.z = initial_values['z']
            self.solver.u = initial_values['u']
        self.solver.solve(_lambda)
        return {'beta': self.solver.beta, 'z': self.solver.z, 'u': self.solver.u }
        

    def _u_admm_lucache(self, y, weights, _lambda, D, converge_threshold, max_steps, verbose, alpha=1.8, initial_values=None, inflate=2., adaptive=False):
        '''Solve for u using alternating direction method of multipliers with a cached LU decomposition.'''
        if verbose:
            print('\t\tSolving u via Alternating Direction Method of Multipliers')

        n = len(y)
        m = D.shape[0]
        a = inflate * _lambda # step-size parameter

        # Initialize primal and dual variables from warm start
        if initial_values is None:
            # Graph Laplacian
            L = csc_matrix(D.T.dot(D) + csc_matrix(np.eye(n)))

            # Cache the LU decomposition
            lu_factor = sla.splu(L, permc_spec='MMD_AT_PLUS_A')
            
            x = np.array([y.mean()] * n) # likelihood term
            z = np.zeros(n) # slack variable for likelihood
            r = np.zeros(m) # penalty term
            s = np.zeros(m) # slack variable for penalty
            u_dual = np.zeros(n) # scaled dual variable for constraint x = z
            t_dual = np.zeros(m) # scaled dual variable for constraint r = s
        else:
            lu_factor = initial_values['lu_factor']
            x = initial_values['x']
            z = initial_values['z']
            r = initial_values['r']
            s = initial_values['s']
            u_dual = initial_values['u_dual']
            t_dual = initial_values['t_dual']

        primal_trace = []
        dual_trace = []
        converged = False
        cur_step = 0
        D_full = D
        while not converged and cur_step < max_steps:
            # Update x
            x = (weights * y + a * (z - u_dual)) / (weights + a)
            x_accel = alpha * x + (1 - alpha) * z # over-relaxation

            # Update constraint term r
            arg = s - t_dual
            local_lambda = (_lambda - np.abs(arg) / 2.).clip(0) if adaptive else _lambda
            r = _soft_threshold(arg, local_lambda / a)
            r_accel = alpha * r + (1 - alpha) * s

            # Projection to constraint set
            arg = x_accel + u_dual + D.T.dot(r_accel + t_dual)
            z_new = lu_factor.solve(arg)
            s_new = D.dot(z_new)
            dual_residual_u = a * (z_new - z)
            dual_residual_t = a * (s_new - s)
            z = z_new
            s = s_new

            # Dual update
            primal_residual_x = x_accel - z
            primal_residual_r = r_accel - s
            u_dual = u_dual + primal_residual_x
            t_dual = t_dual + primal_residual_r

            # Check convergence
            primal_resnorm = np.sqrt((np.array([i for i in primal_residual_x] + [i for i in primal_residual_r])**2).mean())
            dual_resnorm = np.sqrt((np.array([i for i in dual_residual_u] + [i for i in dual_residual_t])**2).mean())
            primal_trace.append(primal_resnorm)
            dual_trace.append(dual_resnorm)
            converged = dual_resnorm < converge_threshold and primal_resnorm < converge_threshold

            if primal_resnorm > 5 * dual_resnorm:
                a *= inflate
                u_dual /= inflate
                t_dual /= inflate
            elif dual_resnorm > 5 * primal_resnorm:
                a /= inflate
                u_dual *= inflate
                t_dual *= inflate

            # Update the step counter
            cur_step += 1

            if verbose and cur_step % 100 == 0:
                print('\t\t\tStep #{0}: dual_resnorm: {1:.6f} primal_resnorm: {2:.6f}'.format(cur_step, dual_resnorm, primal_resnorm))

        return {'x': x, 'r': r, 'z': z, 's': s, 'u_dual': u_dual, 't_dual': t_dual,
                'primal_trace': primal_trace, 'dual_trace': dual_trace, 'steps': cur_step,
                'lu_factor': lu_factor}

    def _u_admm(self, y, weights, _lambda, D, converge_threshold, max_steps, verbose, alpha=1.0, initial_values=None):
        '''Solve for u using alternating direction method of multipliers.'''
        if verbose:
            print('\t\tSolving u via Alternating Direction Method of Multipliers')

        n = len(y)
        m = D.shape[0]

        a = _lambda # step-size parameter

        # Set up system involving graph Laplacian
        L = D.T.dot(D)
        W_over_a = np.diag(weights / a)
        x_denominator = W_over_a + L
        #x_denominator = sparse.linalg.inv(W_over_a + L)

        # Initialize primal and dual variables
        if initial_values is None:
            x = np.array([y.mean()] * n)
            z = np.zeros(m)
            u = np.zeros(m)
        else:
            x = initial_values['x']
            z = initial_values['z']
            u = initial_values['u']

        primal_trace = []
        dual_trace = []
        converged = False
        cur_step = 0
        while not converged and cur_step < max_steps:
            # Update x
            x_numerator = 1.0 / a * weights * y + D.T.dot(a * z - u)
            x = np.linalg.solve(x_denominator, x_numerator)
            Dx = D.dot(x)

            # Update z
            Dx_relaxed = alpha * Dx + (1 - alpha) * z # over-relax Dx
            z_new = _soft_threshold(Dx_relaxed + u / a, _lambda / a)
            dual_residual = a * D.T.dot(z_new - z)
            z = z_new
            primal_residual = Dx_relaxed - z

            # Update u
            u = u + a * primal_residual

            # Check convergence
            primal_resnorm = np.sqrt((primal_residual ** 2).mean())
            dual_resnorm = np.sqrt((dual_residual ** 2).mean())
            primal_trace.append(primal_resnorm)
            dual_trace.append(dual_resnorm)
            converged = dual_resnorm < converge_threshold and primal_resnorm < converge_threshold

            # Update step-size parameter based on norm of primal and dual residuals
            # This is the varying penalty extension to standard ADMM
            a *= 2 if primal_resnorm > 10 * dual_resnorm else 0.5

            # Recalculate the x_denominator since we changed the step-size
            # TODO: is this worth it? We're paying a matrix inverse in exchange for varying the step size
            #W_over_a = sparse.dia_matrix(np.diag(weights / a))
            W_over_a = np.diag(weights / a)
            #x_denominator = sparse.linalg.inv(W_over_a + L)
            
            # Update the step counter
            cur_step += 1

            if verbose and cur_step % 100 == 0:
                print('\t\t\tStep #{0}: dual_resnorm: {1:.6f} primal_resnorm: {2:.6f}'.format(cur_step, dual_resnorm, primal_resnorm))

        dof = np.sum(Dx > converge_threshold) + 1.
        AIC = np.sum((y - x)**2) + 2 * dof

        return {'x': x, 'z': z, 'u': u, 'dof': dof, 'AIC': AIC}

    def _u_admm_1dfusedlasso(self, y, W, _lambda, converge_threshold, max_steps, verbose, alpha=1.0, initial_values=None):
        '''Solve for u using alternating direction method of multipliers. Note that this method only works for the 1-D fused lasso case.'''
        if verbose:
            print('\t\tSolving u via Alternating Direction Method of Multipliers (1-D fused lasso)')

        n = len(y)
        m = n - 1

        a = _lambda

        # The D matrix is the first-difference operator. K is the matrix (W + a D^T D)
        # where W is the diagonal matrix of weights. We use a tridiagonal representation
        # of K.
        Kd = np.array([a] + [2*a] * (n-2) + [a]) + W # diagonal entries
        Kl = np.array([-a] * (n-1)) # below the diagonal
        Ku = np.array([-a] * (n-1)) # above the diagonal

        # Initialize primal and dual variables
        if initial_values is None:
            x = np.array([y.mean()] * n)
            z = np.zeros(m)
            u = np.zeros(m)
        else:
            x = initial_values['x']
            z = initial_values['z']
            u = initial_values['u']

        primal_trace = []
        dual_trace = []
        converged = False
        cur_step = 0
        while not converged and cur_step < max_steps:
            # Update x
            out = _1d_fused_lasso_crossprod(a*z - u)
            x = tridiagonal_solve(Kl, Ku, Kd, W * y + out)
            Dx = np.ediff1d(x)

            # Update z
            Dx_hat = alpha * Dx + (1 - alpha) * z # Over-relaxation
            z_new = _soft_threshold(Dx_hat + u / a, _lambda / a)
            dual_residual = a * _1d_fused_lasso_crossprod(z_new - z)
            z = z_new
            primal_residual = Dx - z
            #primal_residual = Dx_hat - z

            # Update u
            u = (u + a * primal_residual).clip(-_lambda, _lambda)

            # Check convergence
            primal_resnorm = np.sqrt((primal_residual ** 2).mean())
            dual_resnorm = np.sqrt((dual_residual ** 2).mean())
            primal_trace.append(primal_resnorm)
            dual_trace.append(dual_resnorm)
            converged = dual_resnorm < converge_threshold and primal_resnorm < converge_threshold
            
            # Update step-size parameter based on norm of primal and dual residuals
            a *= 2 if primal_resnorm > 10 * dual_resnorm else 0.5
            Kd = np.array([a] + [2*a] * (n-2) + [a]) + W # diagonal entries
            Kl = np.array([-a] * (n-1)) # below the diagonal
            Ku = np.array([-a] * (n-1)) # above the diagonal

            cur_step += 1

            if verbose and cur_step % 100 == 0:
                print('\t\t\tStep #{0}: dual_resnorm: {1:.6f} primal_resnorm: {2:.6f}'.format(cur_step, dual_resnorm, primal_resnorm))

        dof = np.sum(Dx > converge_threshold) + 1.
        AIC = np.sum((y - x)**2) + 2 * dof

        return {'x': x, 'z': z, 'u': u, 'dof': dof, 'AIC': AIC}


    def _u_coord_descent(self, x, A, _lambda, converge, max_steps, verbose, u0=None):
        '''Solve for u using coordinate descent.'''
        if verbose:
            print('\t\tSolving u via Coordinate Descent')
        
        u = u0 if u0 is not None else np.zeros(A.shape[1])

        l2_norm_A = (A * A).sum(axis=0)
        r = x - A.dot(u)
        delta = converge + 1
        prev_objective = _u_objective_func(u, x, A)
        cur_step = 0
        while delta > converge and cur_step < max_steps:
            # Update each coordinate one at a time.
            for coord in range(len(u)):
                prev_u = u[coord]
                next_u = prev_u + A.T[coord].dot(r) / l2_norm_A[coord]
                u[coord] = min(_lambda, max(-_lambda, next_u))
                r += A.T[coord] * prev_u - A.T[coord] * u[coord]

            # Track the change in the objective function value
            cur_objective = _u_objective_func(u, x, A)
            delta = np.abs(prev_objective - cur_objective) / (prev_objective + converge)

            if verbose and cur_step % 100 == 0:
                print('\t\t\tStep #{0}: Objective: {1:.6f} CD Delta: {2:.6f}'.format(cur_step, cur_objective, delta))

            # Increment the step counter and update the previous objective value
            cur_step += 1
            prev_objective = cur_objective

        return u

    def _u_slsqp(self, x, A, _lambda, verbose, u0=None):
        '''Solve for u using sequential least squares.'''
        if verbose:
            print('\t\tSolving u via Sequential Least Squares')

        if u0 is None:
            u0 = np.zeros(A.shape[1])

        # Create our box constraints
        bounds = [(-_lambda, _lambda) for u0_i in u0]

        results = minimize(_u_objective_func, u0,
                           args=(x, A),
                           jac=_u_objective_deriv,
                           bounds=bounds,
                           method='SLSQP',
                           options={'disp': False, 'maxiter': 1000})

        if verbose:
            print('\t\t\t{0}'.format(results.message))
            print('\t\t\tFunction evaluations: {0}'.format(results.nfev))
            print('\t\t\tGradient evaluations: {0}'.format(results.njev))
            print('\t\t\tu: [{0}, {1}]'.format(results.x.min(), results.x.max()))

        return results.x

    def _u_lbfgsb(self, x, A, _lambda, verbose, u0=None):
        '''Solve for u using L-BFGS-B.'''
        if verbose:
            print('\t\tSolving u via L-BFGS-B')

        if u0 is None:
            u0 = np.zeros(A.shape[1])

        # Create our box constraints
        bounds = [(-_lambda, _lambda) for _ in u0]

        # Fit
        results = minimize(_u_objective_func, u0, args=(x, A), method='L-BFGS-B', bounds=bounds, options={'disp': verbose})

        return results.x

    def plateau_regression(self, plateaus, data, grid_map=None, verbose=False):
        '''Perform unpenalized 1-d regression for each of the plateaus.'''
        weights = np.zeros(data.shape)
        for i,(level,p) in enumerate(plateaus):
            if verbose:
                print('\tPlateau #{0}'.format(i+1))
            
            # Get the subset of grid points for this plateau
            if grid_map is not None:
                plateau_data = np.array([data[grid_map[x,y]] for x,y in p])
            else:
                plateau_data = np.array([data[x,y] for x,y in p])

            w = single_plateau_regression(plateau_data, self.signal_dist, self.null_dist)
            for idx in p:
                weights[idx if grid_map is None else grid_map[idx[0], idx[1]]] = w
        posteriors = self._e_step(data, weights)
        weights = weights.flatten()
        return (weights, posteriors)


def _u_objective_func(u, x, A):
    return np.linalg.norm(x - A.dot(u))**2

def _u_objective_deriv(u, x, A):
    return 2*A.T.dot(A.dot(u) - x)

def _u_slsqp_constraint_func(idx, _lambda, u):
    '''Constraint function for the i'th value of u.'''
    return np.array([_lambda - np.abs(u[idx])])

def _u_slsqp_constraint_deriv(idx, u):
    jac = np.zeros(len(u))
    jac[idx] = -np.sign(u[idx])
    return jac

def _1d_fused_lasso_crossprod(x):
    '''Efficiently compute the cross-product D^T x, where D is the first-differences matrix.''' 
    return -np.ediff1d(x, to_begin=x[0], to_end=-x[-1])

def _soft_threshold(x, _lambda):
    return np.sign(x) * (np.abs(x) - _lambda).clip(0)

## Tri-Diagonal Matrix Algorithm (a.k.a Thomas algorithm) solver
## Source: http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
def tridiagonal_solve(a,b,c,f):
    alpha = [0]
    beta = [0]
    n = len(f)
    x = [0] * n
 
    for i in range(n-1):
        alpha.append(-b[i]/(a[i]*alpha[i] + c[i]))
        beta.append((f[i] - a[i]*beta[i])/(a[i]*alpha[i] + c[i]))
 
    x[n-1] = (f[n-1] - a[n-2]*beta[n-1])/(c[n-1] + a[n-2]*alpha[n-1])
 
    for i in reversed(range(n-1)):
        x[i] = alpha[i+1]*x[i+1] + beta[i+1]
 
    return np.array(x)

def ilogit(x):
    return 1. / (1. + np.exp(-x))

def calc_plateaus(beta, rel_tol=1e-4, edges=None, verbose=0):
    '''Calculate the plateaus (degrees of freedom) of a 1d or 2d grid of beta values in linear time.'''
    to_check = deque(itertools.product(*[range(x) for x in beta.shape])) if edges is None else deque(range(len(beta)))
    check_map = np.zeros(beta.shape, dtype=bool)
    check_map[np.isnan(beta)] = True
    plateaus = []

    if verbose:
        print('\tCalculating plateaus...')

    if verbose > 1:
        print('\tIndices to check {0} {1}'.format(len(to_check), check_map.shape))

    # Loop until every beta index has been checked
    while to_check:
        if verbose > 1:
            print('\t\tPlateau #{0}'.format(len(plateaus) + 1))

        # Get the next unchecked point on the grid
        idx = to_check.popleft()

        # If we already have checked this one, just pop it off
        while to_check and check_map[idx]:
            try:
                idx = to_check.popleft()
            except:
                break

        # Edge case -- If we went through all the indices without reaching an unchecked one.
        if check_map[idx]:
            break

        # Create the plateau and calculate the inclusion conditions
        cur_plateau = set([idx])
        cur_unchecked = deque([idx])
        val = beta[idx]
        min_member = val - rel_tol
        max_member = val + rel_tol

        # Check every possible boundary of the plateau
        while cur_unchecked:
            idx = cur_unchecked.popleft()
            
            # neighbors to check
            local_check = []

            # Generic graph case
            if edges is not None:
                local_check.extend(edges[idx])

            # 1d case -- check left and right
            elif len(beta.shape) == 1:
                if idx[0] > 0:
                    local_check.append(idx[0] - 1) # left
                if idx[0] < beta.shape[0] - 1:
                    local_check.append(idx[0] + 1) # right

            # 2d case -- check left, right, up, and down
            elif len(beta.shape) == 2:
                if idx[0] > 0:
                    local_check.append((idx[0] - 1, idx[1])) # left
                if idx[0] < beta.shape[0] - 1:
                    local_check.append((idx[0] + 1, idx[1])) # right
                if idx[1] > 0:
                    local_check.append((idx[0], idx[1] - 1)) # down
                if idx[1] < beta.shape[1] - 1:
                    local_check.append((idx[0], idx[1] + 1)) # up

            # Only supports 1d and 2d cases for now
            else:
                raise Exception('Degrees of freedom calculation does not currently support more than 2 dimensions unless edges are specified explicitly. ({0} given)'.format(len(beta.shape)))

            # Check the index's unchecked neighbors
            for local_idx in local_check:
                if not check_map[local_idx] \
                    and beta[local_idx] >= min_member \
                    and beta[local_idx] <= max_member:
                        # Label this index as being checked so it's not re-checked unnecessarily
                        check_map[local_idx] = True

                        # Add it to the plateau and the list of local unchecked locations
                        cur_unchecked.append(local_idx)
                        cur_plateau.add(local_idx)

        # Track each plateau's indices
        plateaus.append((val, cur_plateau))

    # Returns the list of plateaus and their values
    return plateaus

def plateau_loss_func(c, data, signal_dist, null_dist):
    '''The negative log-likelihood function for a plateau.'''
    return -np.log(c * signal_dist.pdf(data) + (1. - c) * null_dist.pdf(data)).sum()

def single_plateau_regression(data, signal_dist, null_dist):
    '''Perform unpenalized 1-d regression on all of the points in a plateau.'''
    return minimize_scalar(plateau_loss_func, args=(data, signal_dist, null_dist), bounds=(0,1), method='Bounded').x