<commit_before>/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 * 
 * Written (W) 2013 <NAME>
 * Written (W) 2013 <NAME>
 */

#include <shogun/lib/config.h>
#include <shogun/lib/SGVector.h>
#include <shogun/base/Parameter.h>
#include <shogun/mathematics/Math.h>
#include <shogun/mathematics/JacobiEllipticFunctions.h>
#include <shogun/mathematics/logdet/LinearOperator.h>
#include <shogun/mathematics/logdet/LinearSolver.h>
#include <shogun/mathematics/logdet/EigenSolver.h>
#include <shogun/mathematics/logdet/RationalApproximation.h>
#include <shogun/lib/computation/engine/IndependentComputationEngine.h>

namespace shogun
{

CRationalApproximation::CRationalApproximation()
	: COperatorFunction<float64_t>()
{
	init();

	SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
}

CRationalApproximation::CRationalApproximation(
	CLinearOperator<float64_t>* linear_operator,
	CIndependentComputationEngine* computation_engine,
	CEigenSolver* eigen_solver,
	float64_t desired_accuracy,
	EOperatorFunction function_type)
	: COperatorFunction<float64_t>(linear_operator, computation_engine,
	  function_type)
{
	init();

	m_eigen_solver=eigen_solver;
	SG_REF(m_eigen_solver);

	m_desired_accuracy=desired_accuracy;

	SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
}

CRationalApproximation::~CRationalApproximation()
{
	SG_UNREF(m_eigen_solver);

	SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this)
}

void CRationalApproximation::init()
{
	m_eigen_solver=NULL;
	m_constant_multiplier=0.0;
	m_num_shifts=0;
	m_desired_accuracy=0.0;

	SG_ADD((CSGObject**)&m_eigen_solver, "eigen_solver",
		"Eigen solver for computing extremal eigenvalues", MS_NOT_AVAILABLE);

	SG_ADD(&m_shifts, "complex_shifts", "Complex shifts in the linear system",
		MS_NOT_AVAILABLE);

	SG_ADD(&m_weights, "complex_weights", "Complex weights of the linear system",
		MS_NOT_AVAILABLE);

	SG_ADD(&m_constant_multiplier, "constant_multiplier",
		"Constant multiplier in the rational approximation",
		MS_NOT_AVAILABLE);

	SG_ADD(&m_num_shifts, "num_shifts",
		"Number of shifts in the quadrature rule", MS_NOT_AVAILABLE);

	SG_ADD(&m_desired_accuracy, "desired_accuracy",
		"Desired accuracy of the rational approximation", MS_NOT_AVAILABLE);
}

SGVector<complex64_t> CRationalApproximation::get_shifts() const
{
	return m_shifts;
}

SGVector<complex64_t> CRationalApproximation::get_weights() const
{
	return m_weights;
}

float64_t CRationalApproximation::get_constant_multiplier() const
{
	return m_constant_multiplier;
}

index_t CRationalApproximation::get_num_shifts() const
{
	return m_num_shifts;
}

void CRationalApproximation::set_num_shifts(index_t num_shifts)
{
	m_num_shifts=num_shifts;
}

void CRationalApproximation::precompute()
{
	// compute extremal eigenvalues
	m_eigen_solver->compute();
	SG_DEBUG("max_eig=%.15lf\n", m_eigen_solver->get_max_eigenvalue());
	SG_DEBUG("min_eig=%.15lf\n", m_eigen_solver->get_min_eigenvalue());

	// compute number of shifts from accuracy if shifts are not set yet
	if (m_num_shifts==0)
		m_num_shifts=compute_num_shifts_from_accuracy();

	SG_DEBUG("Computing %d shifts\n", m_num_shifts);
	compute_shifts_weights_const();
}

int32_t CRationalApproximation::compute_num_shifts_from_accuracy()
{
	REQUIRE(m_desired_accuracy>0, "Desired accuracy must be positive but is %f\n",
			m_desired_accuracy);

	float64_t max_eig=m_eigen_solver->get_max_eigenvalue();
	float64_t min_eig=m_eigen_solver->get_min_eigenvalue();

	float64_t log_cond_number=CMath::log(max_eig)-CMath::log(min_eig);
	float64_t two_pi_sq=2.0*M_PI*M_PI;

	int32_t num_shifts=static_cast<index_t>(-1.5*(log_cond_number+6.0)
		  	*CMath::log(m_desired_accuracy)/two_pi_sq);

	return num_shifts;
}

void CRationalApproximation::compute_shifts_weights_const()
{
	float64_t PI=M_PI;
	
	// eigenvalues are always real in this case
	float64_t max_eig=m_eigen_solver->get_max_eigenvalue();<commit_msg>min eigenvalue check added (log-det)<commit_after>/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 * 
 * Written (W) 2013 <NAME>
 * Written (W) 2013 <NAME>
 */

#include <shogun/lib/config.h>
#include <shogun/lib/SGVector.h>
#include <shogun/base/Parameter.h>
#include <shogun/mathematics/Math.h>
#include <shogun/mathematics/JacobiEllipticFunctions.h>
#include <shogun/mathematics/logdet/LinearOperator.h>
#include <shogun/mathematics/logdet/LinearSolver.h>
#include <shogun/mathematics/logdet/EigenSolver.h>
#include <shogun/mathematics/logdet/RationalApproximation.h>
#include <shogun/lib/computation/engine/IndependentComputationEngine.h>

namespace shogun
{

CRationalApproximation::CRationalApproximation()
	: COperatorFunction<float64_t>()
{
	init();

	SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
}

CRationalApproximation::CRationalApproximation(
	CLinearOperator<float64_t>* linear_operator,
	CIndependentComputationEngine* computation_engine,
	CEigenSolver* eigen_solver,
	float64_t desired_accuracy,
	EOperatorFunction function_type)
	: COperatorFunction<float64_t>(linear_operator, computation_engine,
	  function_type)
{
	init();

	m_eigen_solver=eigen_solver;
	SG_REF(m_eigen_solver);

	m_desired_accuracy=desired_accuracy;

	SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
}

CRationalApproximation::~CRationalApproximation()
{
	SG_UNREF(m_eigen_solver);

	SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this)
}

void CRationalApproximation::init()
{
	m_eigen_solver=NULL;
	m_constant_multiplier=0.0;
	m_num_shifts=0;
	m_desired_accuracy=0.0;

	SG_ADD((CSGObject**)&m_eigen_solver, "eigen_solver",
		"Eigen solver for computing extremal eigenvalues", MS_NOT_AVAILABLE);

	SG_ADD(&m_shifts, "complex_shifts", "Complex shifts in the linear system",
		MS_NOT_AVAILABLE);

	SG_ADD(&m_weights, "complex_weights", "Complex weights of the linear system",
		MS_NOT_AVAILABLE);

	SG_ADD(&m_constant_multiplier, "constant_multiplier",
		"Constant multiplier in the rational approximation",
		MS_NOT_AVAILABLE);

	SG_ADD(&m_num_shifts, "num_shifts",
		"Number of shifts in the quadrature rule", MS_NOT_AVAILABLE);

	SG_ADD(&m_desired_accuracy, "desired_accuracy",
		"Desired accuracy of the rational approximation", MS_NOT_AVAILABLE);
}

SGVector<complex64_t> CRationalApproximation::get_shifts() const
{
	return m_shifts;
}

SGVector<complex64_t> CRationalApproximation::get_weights() const
{
	return m_weights;
}

float64_t CRationalApproximation::get_constant_multiplier() const
{
	return m_constant_multiplier;
}

index_t CRationalApproximation::get_num_shifts() const
{
	return m_num_shifts;
}

void CRationalApproximation::set_num_shifts(index_t num_shifts)
{
	m_num_shifts=num_shifts;
}

void CRationalApproximation::precompute()
{
	// compute extremal eigenvalues
	m_eigen_solver->compute();
	SG_DEBUG("max_eig=%.15lf\n", m_eigen_solver->get_max_eigenvalue());
	SG_DEBUG("min_eig=%.15lf\n", m_eigen_solver->get_min_eigenvalue());

	REQUIRE(m_eigen_solver->get_min_eigenvalue()>0,
		"Minimum eigenvalue is negative, please provide a Hermitian matrix\n");

	// compute number of shifts from accuracy if shifts are not set yet
	if (m_num_shifts==0)
		m_num_shifts=compute_num_shifts_from_accuracy();

	SG_DEBUG("Computing %d shifts\n", m_num_shifts);
	compute_shifts_weights_const();
}

int32_t CRationalApproximation::compute_num_shifts_from_accuracy()
{
	REQUIRE(m_desired_accuracy>0, "Desired accuracy must be positive but is %f\n",
			m_desired_accuracy);

	float64_t max_eig=m_eigen_solver->get_max_eigenvalue();
	float64_t min_eig=m_eigen_solver->get_min_eigenvalue();

	float64_t log_cond_number=CMath::log(max_eig)-CMath::log(min_eig);
	float64_t two_pi_sq=2.0*M_PI*M_PI;

	int32_t num_shifts=static_cast<index_t>(-1.5*(log_cond_number+6.0)
		  	*CMath::log(m_desired_accuracy)/two_pi_sq);

	return num_shifts;
}

void CRationalApproximation::compute_shifts_weights_const()
{
	float64_t PI=M_PI;
	
	// eigenvalues are always real in this case
	float64_t max_eig=m_eigen_solver->get_max_eigenvalue();
{'fraction_non_alphanumeric': 0.11306743322553599, 'fraction_numerical': 0.013534555036531322, 'mean_word_length': 7.35, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 61, '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': '2492403', 'n_tokens_mistral': 3284, 'n_tokens_neox': 3105, 'n_words': 539}