edgemodel.cpp

MS-Clustering is designed to rapidly cluster large MS/MS datasets. The program merges similar spectr

				multi_log_odds_offset[i][j]=multi_log_odds_offset[i][j-1];
		}


		for (j=0; j<num_states; j++)
		{
			log_odds_transfer[i][j]=log(good_trans_counts[i][j]/n_good_trans_counts) - 
									log(bad_trans_counts[i][j]/n_bad_trans_counts);
			if ((i==0 || j==0) && log_odds_transfer[i][j]>0)
				log_odds_transfer[i][j]=0;
		}
	}

	has_values = true;

	if (verbose)
	{
		cout << "Single offset:" << endl;
		for (i=0;  i<num_states; i++)
		{
			cout << setprecision(0) << fixed << i << "\t" << n_good_single_off_counts[i] << "/" << 
			n_bad_single_off_counts[i] << "\t" << setprecision(2) << log_odds_state[i] << " : ";
			for (j=0; j<num_limits; j++)
				cout << "\t" << log_odds_offset[i][j];
			cout << endl ;
		}
		cout << endl;

		cout << "Multi offset:" << endl;
		for (i=0;  i<num_states; i++)
		{
			cout << setprecision(0) << fixed << i << "\t" << n_good_multi_off_counts[i] << "/" << 
			n_bad_multi_off_counts[i] << "\t" << setprecision(2) << multi_log_odds_state[i] << " : ";
			for (j=0; j<num_limits; j++)
				cout << "\t" << multi_log_odds_offset[i][j];
			cout << endl ;
		}
		cout << endl;

		cout << "Transfer odds:" << endl;
		cout << "\t";
		for (i=0; i<num_states; i++)
			cout << "  " << i << "\t";
		cout << endl;
		for (i=0; i<num_states; i++)
		{
			cout << i << "\t";
			for (j=0; j<num_states; j++)
				cout << log_odds_transfer[i][j] << "\t";
			cout << endl;
		}

	}
}


bool RegionalEdgeModel::write_model(const char *path) const
{
	if (! has_values)
		return false;

	int i;

	ofstream ofs(path);
	if (! ofs.is_open() || ! ofs.good())
	{
		cout << "Error: couldn't open edge model file for writing: " << path << endl;
		exit(1);
	}

	ofs << charge << " " << size_idx << " " << region_idx << " " << num_states << endl;
	ofs << fixed << setprecision(3);

	ofs << strong_frag_idxs.size();
	for (i=0; i<strong_frag_idxs.size(); i++)
		ofs << " " << strong_frag_idxs[i];
	ofs << endl;

	ofs << strong_frag_charges.size();
	for (i=0; i<strong_frag_charges.size(); i++)
		ofs << " " << strong_frag_charges[i];
	ofs << endl;

	ofs << log_odds_state.size();
	for (i=0; i<log_odds_state.size(); i++)
		ofs << " " << log_odds_state[i];
	ofs << endl;
	
	ofs << log_odds_offset.size() << endl;
	for (i=0; i<log_odds_offset.size(); i++)
	{
		int j;
		ofs << log_odds_offset[i].size();
		for (j=0; j<log_odds_offset[i].size(); j++)
			ofs << " " << log_odds_offset[i][j];
		ofs << endl;
	}

	ofs << multi_log_odds_state.size();
	for (i=0; i<multi_log_odds_state.size(); i++)
		ofs << " " << multi_log_odds_state[i];
	ofs << endl;
	
	ofs << multi_log_odds_offset.size() << endl;
	for (i=0; i<multi_log_odds_offset.size(); i++)
	{
		int j;
		ofs << multi_log_odds_offset[i].size();
		for (j=0; j<multi_log_odds_offset[i].size(); j++)
			ofs << " " << multi_log_odds_offset[i][j];
		ofs << endl;
	}

	ofs << log_odds_transfer.size() << endl;
	for (i=0; i<log_odds_transfer.size(); i++)
	{
		int j;
		ofs << log_odds_transfer[i].size();
		for (j=0; j<log_odds_transfer[i].size(); j++)
			ofs << " " << log_odds_transfer[i][j];
		ofs << endl;
	}

	return true;
}


bool RegionalEdgeModel::read_model(const char *path)
{
	
	ifstream ifs(path);
	if (! ifs.is_open() || ! ifs.good())
		return false;

	ifs >> charge >> size_idx >> region_idx >> num_states;

	int i;
	int ns;
	ifs >> ns;
	strong_frag_idxs.resize(ns);
	for (i=0; i<ns; i++)
		ifs >> strong_frag_idxs[i];

	ifs >> ns;
	strong_frag_charges.resize(ns);
	for (i=0; i<ns; i++)
		ifs >> strong_frag_charges[i];

	ifs >> ns;
	log_odds_state.resize(ns);
	for (i=0; i<ns; i++)
		ifs >> log_odds_state[i];
	
	ifs >> ns;
	log_odds_offset.resize(ns);
	for (i=0; i<ns; i++)
	{
		int ls;
		int j;
		ifs >> ls;
		log_odds_offset[i].resize(ls);
		for (j=0; j<ls; j++)
			ifs >> log_odds_offset[i][j];
	}

	ifs >> ns;
	multi_log_odds_state.resize(ns);
	for (i=0; i<ns; i++)
		ifs >> multi_log_odds_state[i];
	
	ifs >> ns;
	multi_log_odds_offset.resize(ns);
	for (i=0; i<ns; i++)
	{
		int ls;
		int j;
		ifs >> ls;
		multi_log_odds_offset[i].resize(ls);
		for (j=0; j<ls; j++)
			ifs >> multi_log_odds_offset[i][j];
	}

	ifs >> ns;
	log_odds_transfer.resize(ns);

	for (i=0; i<ns; i++)
	{
		int ls;
		int j;
		ifs >> ls;
		log_odds_transfer[i].resize(ls);

		for (j=0; j<ls; j++)
			ifs >> log_odds_transfer[i][j];
	}

	has_values = true;

	return true;
}


string EdgeModel::make_regional_name(const char *model_name, int c, int s, int r) const
{
	ostringstream oss;
	oss << model_name << "_edge_" << c << "_" << s << "_" << r << ".txt";

	return oss.str();
}


void EdgeModel::write_edge_models(const char *model_name)
{
	if (! this->ind_was_initialized)
		return;
			
	const string dir_name = config->get_resource_dir() + "/" + model_name + "_EDGE" + "/";
	const string file_path = dir_name + string(model_name) + "_edge_model.txt";

	ofstream ofs(file_path.c_str());
	if (! ofs.is_open() || ! ofs.good())
	{
		cout << "Error: couldn't open edge model for writing: " << file_path << endl;
		exit(1);
	}

	ofs << "#WEIGHT_SINGLE_STATE_SCORE " << setprecision(3) << fixed << weight_single_state_score << endl;
	ofs << "#WEIGHT_SINGLE_OFFSET_SCORE " << weight_single_offset_score << endl;
	ofs << "#WEIGHT_MULTI_STATE_SCORE " << weight_multi_state_score << endl;
	ofs << "#WEIGHT_MULTI_OFFSET_SCORE " << weight_multi_offset_score << endl;
	ofs << "#WEIGHT_TRANSFER_SCORE " << weight_transfer_score << endl;
	ofs << "#WEIGHT_COMBO_SCORE " << weight_combo_score << endl;
	ofs << "#MULTY_AA_PENALTY " << multi_aa_penalty << endl;
	ofs << "#BAD_DOUBLE_PENALTY " << bad_pair_penalty << endl;
	ofs << "#PROBLEMATIC_PAIR_PENALTY " << problematic_pair_penalty << endl;
	ofs << "#TOLS " << tol_limits.size();
	int i;
	for (i=0; i<tol_limits.size(); i++)
		ofs << " " << tol_limits[i];
	ofs << endl;

	if (double_combo_scores.size()<Val)
	{
		cout << "Error: double combo table not filled!" << endl;
		exit(1);
	}
	for (i=Ala; i<=Val; i++)
	{
		int j;
		for (j=Ala; j<=Val; j++)
			ofs << double_combo_scores[i][j] << " ";
		ofs << endl;
	}

	ofs.close();

	int num_written=0;
	for (i=0; i<this->regional_edge_models.size(); i++)
	{
		int j;
		for (j=0; j<regional_edge_models[i].size(); j++)
		{
			int k;
			for (k=0; k<regional_edge_models[i][j].size(); k++)
			{
				if (regional_edge_models[i][j][k] && regional_edge_models[i][j][k]->has_values)
				{
					string path = dir_name + make_regional_name(model_name,i,j,k);
					if (regional_edge_models[i][j][k]->write_model(path.c_str()))
						num_written++;
				}
			}
		}
	}

	cout << endl << "Wrote " << num_written << " edge models" << endl;
}

void EdgeModel::read_edge_models(Config *_config, const char *model_name, bool silent_ind)
{
	config = _config;

	const string dir_name = config->get_resource_dir() + "/" + model_name + "_EDGE" + "/";
	const string file_path = dir_name + string(model_name) + "_edge_model.txt";

	

	ifstream ifs(file_path.c_str());
	if (! ifs.is_open() || ! ifs.good())
	{
		cout << "Error: couldn't open edge model for reading: " << file_path << endl;
		exit(1);
	}

	int i,num_tols = 0;
	string dummy;
	ifs >> dummy >> weight_single_state_score;
	ifs >> dummy >> weight_single_offset_score;
	ifs >> dummy >> weight_multi_state_score;
	ifs >> dummy >> weight_multi_offset_score;
	ifs >> dummy >> weight_transfer_score;
	ifs >> dummy >> weight_combo_score;
	ifs >> dummy >> multi_aa_penalty;
	ifs >> dummy >> bad_pair_penalty;
	ifs >> dummy >> problematic_pair_penalty;
	ifs >> dummy >> num_tols;

	tol_limits.resize(num_tols);
	for (i=0; i<num_tols; i++)
		ifs >> tol_limits[i];

	// read double combo table
	double_combo_scores.clear();
	double_combo_scores.resize(Val+1);
	for (i=0; i<=Val; i++)
		double_combo_scores[i].resize(Val+1,0);

	char buff[1024];
	for (i=Ala; i<=Val; i++)
	{
		ifs.getline(buff,1024);
		istringstream iss(buff);
		int j;
		for (j=Ala; j<=Val; j++)
			iss >> double_combo_scores[i][j];
	}

	const vector< vector< mass_t > >& size_thresholds = config->get_size_thresholds();

	// resize according to threhsolds, assume that the last mass is POS_INF
	regional_edge_models.resize(size_thresholds.size());
	
	int num_read=0;
	int c;
	for (c=0; c<size_thresholds.size(); c++)
	{
		if (size_thresholds[c].size()>0)
			regional_edge_models[c].resize(size_thresholds[c].size());

		int s;
		for (s=0; s<regional_edge_models[c].size(); s++)
		{
			const int num_regions = config->get_num_regions(c,s);
			regional_edge_models[c][s].resize(num_regions,NULL);

			int r;
			for (r=0; r<regional_edge_models[c][s].size(); r++)
			{
				string path = dir_name + make_regional_name(model_name,c,s,r);
				
				ifstream ifs(path.c_str());
				if (ifs.is_open() && ifs.good())
				{
					ifs.close();
			
					regional_edge_models[c][s][r] = new RegionalEdgeModel;
					if (regional_edge_models[c][s][r]->read_model(path.c_str()))
						num_read++;
				}
				else
					ifs.close();
			}
		}
	}
	if (! silent_ind)
		cout << "Read " << num_read << " edge models." << endl;

	ind_was_initialized = true;
}




void EdgeModel::train_all_edge_models(const FileManager& fm, void *model_ptr, int specific_charge)
{
	Model *model = (Model *)model_ptr;
	config = model->get_config();
	
	const vector< vector< mass_t > >& size_thresholds = config->get_size_thresholds();