<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) 2012 <NAME>
 * Copyright (C) 2010-2012 <NAME>, <NAME>
 */

#include <shogun/lib/slep/slep_tree_mt_lsr.h>
#include <shogun/mathematics/Math.h>
#include <shogun/lib/slep/tree/altra.h>
#include <shogun/lib/slep/tree/general_altra.h>
#include <shogun/mathematics/lapack.h>

namespace shogun
{

SGMatrix<double> slep_tree_mt_lsr(
		CDotFeatures* features,
		double* y,
		double z,
		const slep_options& options)
{
	int i,j,t;
	int n_feats = features->get_dim_feature_space();
	int n_vecs = features->get_num_vectors();
	double lambda, lambda_max, beta;
	double funcp = 0.0, func = 0.0;

	int n_tasks = options.n_tasks;

	int iter = 1;
	bool done = false;
	bool gradient_break = false;

	double* ATy = SG_CALLOC(double, n_feats*n_tasks);
	for (t=0; t<n_tasks; t++)
	{
		int task_ind_start = options.ind[t]+1;
		int task_ind_end = options.ind[t+1];
		for (i=task_ind_start; i<task_ind_end; i++)
			features->add_to_dense_vec(y[i],i,ATy+t*n_feats,n_feats);
	}

	if (options.regularization!=0)
	{
		if (options.general)
			lambda_max = general_findLambdaMax_mt(ATy, n_feats, n_tasks, options.G, options.ind_t, options.n_nodes);
		else
			lambda_max = findLambdaMax_mt(ATy, n_feats, n_tasks, options.ind_t, options.n_nodes);
		lambda = z*lambda_max;
	}
	else 
		lambda = z;

	SGMatrix<double> w(n_feats,n_tasks);
	if (options.initial_w)
	{
		for (j=0; j<n_tasks; j++)
			for (i=0; i<n_feats; i++)
				w(i,j) = options.initial_w[j*n_feats+i];
	}
	else
	{
		for (j=0; j<n_tasks*n_feats; j++)
				w[j] = 0.0;
	}

	double* s = SG_CALLOC(double, n_feats*n_tasks);
	double* g = SG_CALLOC(double, n_feats*n_tasks);
	double* v = SG_CALLOC(double, n_feats*n_tasks);

	double* Aw = SG_CALLOC(double, n_vecs);
	for (t=0; t<n_tasks; t++)
	{
		int task_ind_start = options.ind[t]+1;
		int task_ind_end = options.ind[t+1];
		for (i=task_ind_start; i<task_ind_end; i++)
			Aw[i] = features->dense_dot(i,w.matrix+t*n_feats,n_feats);
	}
	double* Av = SG_MALLOC(double, n_vecs);
	double* As = SG_MALLOC(double, n_vecs);
	double* ATAs = SG_MALLOC(double, n_feats*n_tasks);
	double* resid = SG_MALLOC(double, n_vecs*n_tasks);

	double L = 1.0;

	double* wp = SG_CALLOC(double, n_feats*n_tasks);
	for (i=0; i<n_feats*n_tasks; i++)
		wp[i] = w[i];
	double* Awp = SG_MALLOC(double, n_vecs);
	for (i=0; i<n_vecs; i++)
		Awp[i] = Aw[i];
	double* wwp = SG_CALLOC(double, n_feats*n_tasks);

	double* w_row = SG_MALLOC(double, n_tasks);

	double alphap = 0.0;
	double alpha = 1.0;
	
	while (!done && iter < options.max_iter) 
	{
		beta = (alphap-1.0)/alpha;

		for (i=0; i<n_feats*n_tasks; i++)
			s[i] = w[i] + beta*wwp[i];

		for (i=0; i<n_vecs; i++)
			As[i] = Aw[i] + beta*(Aw[i]-Awp[i]);

		// ATAs = A'*As
		for (i=0; i<n_feats*n_tasks; i++)
			ATAs[i] = 0.0;

		for (t=0; t<n_tasks; t++)
		{
			int task_ind_start = options.ind[t]+1;
			int task_ind_end = options.ind[t+1];
			for (i=task_ind_start; i<task_ind_end; i++)
				features->add_to_dense_vec(As[i],i,ATAs+t*n_feats,n_feats);
		}
		
		// g = ATAs - ATy
		for (i=0; i<n_feats*n_tasks; i++)
			g[i] = ATAs[i] - ATy[i];

		for (i=0; i<n_feats*n_tasks; i++)
			wp[i] = w[i];

		for (i=0; i<n_vecs; i++)
			Awp[i] = Aw[i];

		while (true)
		{
			// v = s - g / L
			for (i=0; i<n_feats*n_tasks; i++)
				v[i] = s[i] - g[i]*(1.0/L);

			if (options.general)
				general_altra_mt(w.matrix, v, n_feats, n_tasks, options.G, options.ind_t, options.n_nodes, lambda/L);
			else
				altra_mt(w.matrix, v, n_feats, n_tasks, options.ind_t, options.n_nodes, lambda/L);

			// v = x - s
			for (i=0; i<n_feats*n_tasks; i++)
				v[i] = w[i] - s[i];
	
			for (t=0; t<n_tasks; t++)
			{
				int task_ind_start = options.ind[t]+1;
				int task_ind_end = options.ind[t+1];
				for (i=task_ind_start; i<task_ind_end; i++)
					Aw[i] = features->dense_dot(i,w.matrix+t*n_feats,n_feats);<commit_msg>Fixes for multitask tree guided least squares<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) 2012 <NAME>
 * Copyright (C) 2010-2012 <NAME>, <NAME>
 */

#include <shogun/lib/slep/slep_tree_mt_lsr.h>
#include <shogun/mathematics/Math.h>
#include <shogun/lib/slep/tree/altra.h>
#include <shogun/lib/slep/tree/general_altra.h>
#include <shogun/mathematics/lapack.h>

namespace shogun
{

SGMatrix<double> slep_tree_mt_lsr(
		CDotFeatures* features,
		double* y,
		double z,
		const slep_options& options)
{
	int i,j,t;
	int n_feats = features->get_dim_feature_space();
	int n_vecs = features->get_num_vectors();
	double lambda, lambda_max, beta;
	double funcp = 0.0, func = 0.0;

	int n_tasks = options.n_tasks;

	int iter = 1;
	bool done = false;
	bool gradient_break = false;

	double* ATy = SG_CALLOC(double, n_feats*n_tasks);
	for (t=0; t<n_tasks; t++)
	{
		int task_ind_start = options.ind[t];
		int task_ind_end = options.ind[t+1];
		for (i=task_ind_start; i<task_ind_end; i++)
		{
			features->add_to_dense_vec(y[i],i,ATy+t*n_feats,n_feats);
		}
	}

	if (options.regularization!=0)
	{
		if (options.general)
			lambda_max = general_findLambdaMax_mt(ATy, n_feats, n_tasks, options.G, options.ind_t, options.n_nodes);
		else
			lambda_max = findLambdaMax_mt(ATy, n_feats, n_tasks, options.ind_t, options.n_nodes);
		lambda = z*lambda_max;
	}
	else 
		lambda = z;

	SGMatrix<double> w(n_feats,n_tasks);
	if (options.initial_w)
	{
		for (j=0; j<n_tasks; j++)
			for (i=0; i<n_feats; i++)
				w(i,j) = options.initial_w[j*n_feats+i];
	}
	else
	{
		for (j=0; j<n_tasks*n_feats; j++)
				w[j] = 0.0;
	}

	double* s = SG_CALLOC(double, n_feats*n_tasks);
	double* g = SG_CALLOC(double, n_feats*n_tasks);
	double* v = SG_CALLOC(double, n_feats*n_tasks);

	double* Aw = SG_CALLOC(double, n_vecs);
	for (t=0; t<n_tasks; t++)
	{
		int task_ind_start = options.ind[t];
		int task_ind_end = options.ind[t+1];
		for (i=task_ind_start; i<task_ind_end; i++)
			Aw[i] = features->dense_dot(i,w.matrix+t*n_feats,n_feats);
	}
	double* Av = SG_MALLOC(double, n_vecs);
	double* As = SG_MALLOC(double, n_vecs);
	double* ATAs = SG_MALLOC(double, n_feats*n_tasks);
	double* resid = SG_MALLOC(double, n_vecs*n_tasks);

	double L = 1.0;

	double* wp = SG_CALLOC(double, n_feats*n_tasks);
	for (i=0; i<n_feats*n_tasks; i++)
		wp[i] = w[i];
	double* Awp = SG_MALLOC(double, n_vecs);
	for (i=0; i<n_vecs; i++)
		Awp[i] = Aw[i];
	double* wwp = SG_CALLOC(double, n_feats*n_tasks);

	double* w_row = SG_MALLOC(double, n_tasks);

	double alphap = 0.0;
	double alpha = 1.0;
	
	while (!done && iter <= options.max_iter) 
	{
		beta = (alphap-1.0)/alpha;

		for (i=0; i<n_feats*n_tasks; i++)
			s[i] = w[i] + beta*wwp[i];

		for (i=0; i<n_vecs; i++)
			As[i] = Aw[i] + beta*(Aw[i]-Awp[i]);

		// ATAs = A'*As
		for (i=0; i<n_feats*n_tasks; i++)
			ATAs[i] = 0.0;

		for (t=0; t<n_tasks; t++)
		{
			int task_ind_start = options.ind[t];
			int task_ind_end = options.ind[t+1];
			for (i=task_ind_start; i<task_ind_end; i++)
				features->add_to_dense_vec(As[i],i,ATAs+t*n_feats,n_feats);
		}
		
		// g = ATAs - ATy
		for (i=0; i<n_feats*n_tasks; i++)
			g[i] = ATAs[i] - ATy[i];

		for (i=0; i<n_feats*n_tasks; i++)
			wp[i] = w[i];

		for (i=0; i<n_vecs; i++)
			Awp[i] = Aw[i];

		while (true)
		{
			// v = s - g / L
			for (i=0; i<n_feats*n_tasks; i++)
				v[i] = s[i] - g[i]*(1.0/L);

			if (options.general)
				general_altra_mt(w.matrix, v, n_feats, n_tasks, options.G, options.ind_t, options.n_nodes, lambda/L);
			else
				altra_mt(w.matrix, v, n_feats, n_tasks, options.ind_t, options.n_nodes, lambda/L);

			// v = x - s
			for (i=0; i<n_feats*n_tasks; i++)
				v[i] = w[i] - s[i];
	
			for (t=0; t<n_tasks; t++)
			{
				int task_ind_start = options.ind[t];
				int task_ind_end = options.ind[t+1];
				for (i=task_ind_start; i<task_ind_end; i++)
					Aw[i] = features->dense_dot(i,w.matrix+t*n_feats,n_feats);
{'fraction_non_alphanumeric': 0.1720443199805187, 'fraction_numerical': 0.01363691708267381, 'mean_word_length': 3.981200727713766, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 35, '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': '6633714', 'n_tokens_mistral': 4137, 'n_tokens_neox': 3773, 'n_words': 751}