<commit_before># -*- coding: utf-8 -*-
''' 
Contains a function for restricted maximum likelihood (REML)
hyperparameter estimation.
'''
import numpy as np
import logging
from scipy.optimize import fmin
from pygeons.mp import parmap
from pygeons.filter.gprocs import gpcomp,get_units
logger = logging.getLogger(__name__)


def fmin_pos(func,x0,*args,**kwargs):
  '''fmin with positivity constraint and multiple start points'''
  def pos_func(x,*blargs):
    return func(np.exp(x),*blargs)

  xopt,fopt,_,_,_ = fmin(pos_func,np.log(x0),*args,full_output=True,**kwargs)
  xopt = np.exp(xopt)
  return xopt,fopt


def reml(y,d,s,model,params,
         fix=(), 
         order=1,
         annual=False,
         semiannual=False,
         procs=0):
  ''' 
  Returns the Restricted Maximum Likelihood (REML) estimatates of the
  unknown hyperparameters.

  Parameters
  ----------
  y : (N,D) array
    Observation points.

  d : (...,N) array
    Observed data at *y*.
  
  s : (...,N) array
    Data uncertainty.
  
  model : str, optional
    string indicating the stochastic model being optimized. Can either
    be 'se','fogm', or 'se+fogm'. The required length of *params* will
    change depending on the number of hyperparameters in the model.
    'se' and 'fogm' take 2 hyperparameters, and 'se+fogm' takes 4
    hyperparameters.
    
  params : (P,) array 
    Initial guess for the hyperparameters. If 4 are specified then the
    model consists of SE and FOGM. If 2 are specified then the model
    consists of SE
  
  fix : (L,) int array
    Indices of hyperparameters that will remain fixed
  
  order : int, optional
    Order of the polynomial improper basis functions.
  
  annual : bool, optional  
    Indicates whether to include annual sinusoids in the model.

  semiannual : bool, optional  
    Indicates whether to include semiannual sinusoids in the model.

  procs : int, optional
    Distribute the tasks among this many subprocesses. 
  
  Returns
  -------
  a : (...,P) float array
    hyperparameters for each dataset. The hyperparameters are the
    standard deviation for the SE, the characteristic-length scale for
    the SE, the standard deviation for the FOGM, and the cutoff
    frequency for the FOGM.
  
  b : (...) float array  
    likelihoods associated with each set of hyperparameters

  c : (...) int array
    number of observations used to constrain the hyperparameters
    
  d : (P,) str array
    units of the hyperparameters 
    
  '''
  y = np.asarray(y,dtype=float)
  d = np.asarray(d,dtype=float)
  s = np.asarray(s,dtype=float)
  params = np.asarray(params,dtype=float)
  fix = np.asarray(fix,dtype=int)
  is_free = np.ones(params.shape,dtype=bool)
  is_free[fix] = False

  bcast_shape = d.shape[:-1]
  q = int(np.prod(bcast_shape))
  n = y.shape[0]
  d = d.reshape((q,n))
  s = s.reshape((q,n))
  units = get_units(model)

  def objective(theta,pos,data,sigma):
    test_params = np.copy(params)
    test_params[is_free] = theta 
    gp = gpcomp(model,test_params)
    # return negative log likelihood
    return -gp.likelihood(pos,data,sigma)  

  def task(i):
    logger.debug('Finding REML hyperparameters for dataset %s of %s ...' % (i+1,q))
    # if the uncertainty is inf then the data is considered missing
    is_missing = np.isinf(s[i])
    yi,di,si = y[~is_missing],d[i,~is_missing],s[i,~is_missing]
    opt,l = fmin_pos(objective,params[is_free],args=(yi,di,si),disp=False)<commit_msg>removed unused parameters in reml<commit_after># -*- coding: utf-8 -*-
''' 
Contains a function for restricted maximum likelihood (REML)
hyperparameter estimation.
'''
import numpy as np
import logging
from scipy.optimize import fmin
from pygeons.mp import parmap
from pygeons.filter.gprocs import gpcomp,get_units
logger = logging.getLogger(__name__)


def fmin_pos(func,x0,*args,**kwargs):
  '''fmin with positivity constraint and multiple start points'''
  def pos_func(x,*blargs):
    return func(np.exp(x),*blargs)

  xopt,fopt,_,_,_ = fmin(pos_func,np.log(x0),*args,full_output=True,**kwargs)
  xopt = np.exp(xopt)
  return xopt,fopt


def reml(y,d,s,model,params,
         fix=(),procs=0):
  ''' 
  Returns the Restricted Maximum Likelihood (REML) estimatates of the
  unknown hyperparameters.

  Parameters
  ----------
  y : (N,D) array
    Observation points.

  d : (...,N) array
    Observed data at *y*.
  
  s : (...,N) array
    Data uncertainty.
  
  model : str,
    string indicating the stochastic model being optimized.
    
  params : (P,) array 
    Initial guess for the hyperparameters.
  
  fix : (L,) int array
    Indices of hyperparameters that will remain fixed
  
  procs : int, optional
    Distribute the tasks among this many subprocesses. 
  
  Returns
  -------
  a : (...,P) float array
    hyperparameters for each dataset. The hyperparameters are the
    standard deviation for the SE, the characteristic-length scale for
    the SE, the standard deviation for the FOGM, and the cutoff
    frequency for the FOGM.
  
  b : (...) float array  
    likelihoods associated with each set of hyperparameters

  c : (...) int array
    number of observations used to constrain the hyperparameters
    
  units : (P,) str array
    units of the hyperparameters 
    
  '''
  y = np.asarray(y,dtype=float)
  d = np.asarray(d,dtype=float)
  s = np.asarray(s,dtype=float)
  params = np.asarray(params,dtype=float)
  fix = np.asarray(fix,dtype=int)
  is_free = np.ones(params.shape,dtype=bool)
  is_free[fix] = False

  bcast_shape = d.shape[:-1]
  q = int(np.prod(bcast_shape))
  n = y.shape[0]
  d = d.reshape((q,n))
  s = s.reshape((q,n))
  units = get_units(model)

  def objective(theta,pos,data,sigma):
    test_params = np.copy(params)
    test_params[is_free] = theta 
    gp = gpcomp(model,test_params)
    # return negative log likelihood
    return -gp.likelihood(pos,data,sigma)  

  def task(i):
    logger.debug('Finding REML hyperparameters for dataset %s of %s ...' % (i+1,q))
    # if the uncertainty is inf then the data is considered missing
    is_missing = np.isinf(s[i])
    yi,di,si = y[~is_missing],d[i,~is_missing],s[i,~is_missing]
    opt,l = fmin_pos(objective,params[is_free],args=(yi,di,si),disp=False)
{'fraction_non_alphanumeric': 0.11533505154639176, 'fraction_numerical': 0.0030605670103092783, 'mean_word_length': 3.3148019457956917, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1965482', 'n_tokens_mistral': 2108, 'n_tokens_neox': 2039, 'n_words': 611}