##### https://github.com/tansey/smoothfdr
Revision 49cb69c5a18fdb262964fbfeb47ab2099eb32c5c authored by Wesley Tansey on 03 May 2018, 19:46:59 UTC, committed by Wesley Tansey on 03 May 2018, 19:46:59 UTC
2 parent s b713c52 + 63bb329
Tip revision: 49cb69c
smoothed_fdr.py
``````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({0}, {1}^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 '''
# 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',
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 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,
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 - 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',
'''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,
u0=u)

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

# 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

# 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,
'''
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)
# 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,
beta = u['x']
elif dual_solver == 'graph':
u = self._graph_fused_lasso(y, weights, _lambda, penalties, penalties, penalties, penalties, 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
# 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
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

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)

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)

# 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\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)

# 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, idx]] = 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, 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 = 
beta = 
n = len(f)
x =  * 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:
local_check.append(idx - 1) # left
if idx < beta.shape - 1:
local_check.append(idx + 1) # right

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

# 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

`````` Computing file changes ...