text:
/* -*- c++ -*- */
/*
 * Copyright 2014 Communications Engineering Lab, KIT.
 *
 * This 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, or (at your option)
 * any later version.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this software; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <gnuradio/io_signature.h>
#include <boost/circular_buffer.hpp>
#include "estimator_rcs_impl.h"
#include <numeric>

namespace gr {
  namespace radar {

    estimator_rcs::sptr
    estimator_rcs::make(int num_mean, float center_freq, float antenna_gain_tx, float antenna_gain_rx, float usrp_gain_rx, float power_tx, float corr_factor, float exponent)
    {
      return gnuradio::get_initial_sptr
        (new estimator_rcs_impl(num_mean, center_freq, antenna_gain_tx, antenna_gain_rx, usrp_gain_rx, power_tx, corr_factor, exponent));
    }

    /*
     * The private constructor
     */
    estimator_rcs_impl::estimator_rcs_impl(int num_mean, float center_freq, float antenna_gain_tx, float antenna_gain_rx, float usrp_gain_rx, float power_tx, float corr_factor, float exponent)
      : gr::block("estimator_rcs",
              gr::io_signature::make(0,0,0),
              gr::io_signature::make(0,0,0))
    {
        d_num_mean = num_mean;
        d_center_freq = center_freq;
        d_antenna_gain_tx = antenna_gain_tx;
        d_antenna_gain_rx = antenna_gain_rx;
        d_usrp_gain_rx = usrp_gain_rx;
        d_power_tx = power_tx; // needs to be calibrated for every usage
        d_corr_factor = corr_factor;
        d_exponent = exponent;

        d_rcs_vals.resize(d_num_mean);

        // Register input message port
        d_port_id_in = pmt::mp("Msg in");
        message_port_register_in(d_port_id_in);
        set_msg_handler(d_port_id_in, boost::bind(&estimator_rcs_impl::handle_msg, this, _1));

        // Register output message port
        d_port_id_out = pmt::mp("Msg out");
        message_port_register_out(d_port_id_out);

        d_loop_counter = 0;

        // constant factors in radar equation
        d_antenna_gain_abs_rx = pow(10, d_antenna_gain_rx/10);
        d_antenna_gain_abs_tx = pow(10, d_antenna_gain_tx/10);
        d_lambda = c_light/d_center_freq;
        d_fak = pow(4.0*M_PI, 3) / (d_antenna_gain_abs_rx * d_antenna_gain_abs_tx * pow(d_lambda, 2));

    }

    /*
     * Our virtual destructor.
     */
    estimator_rcs_impl::~estimator_rcs_impl()
    {
    }

    void
    estimator_rcs_impl::set_num_mean(int val){
      d_num_mean = val;
      d_rcs_vals.clear();
      d_rcs_vals.resize(d_num_mean);
      d_loop_counter = 0;
    }

    void
    estimator_rcs_impl::set_center_freq(float val){
      d_center_freq = val;
      d_lambda = c_light/d_center_freq;
      d_fak = pow(4.0*M_PI, 3) / (d_antenna_gain_abs_rx * d_antenna_gain_abs_tx * pow(d_lambda, 2));
    }

    void
    estimator_rcs_impl::set_antenna_gain_tx(float val){
      d_antenna_gain_tx = val;
      d_antenna_gain_abs_tx = pow(10, d_antenna_gain_tx/10);
      d_fak = pow(4.0*M_PI, 3) / (d_antenna_gain_abs_rx * d_antenna_gain_abs_tx * pow(d_lambda, 2));
    }

    void
    estimator_rcs_impl::set_antenna_gain_rx(float val){
      d_antenna_gain_rx = val;
      d_antenna_gain_abs_rx = pow(10, d_antenna_gain_rx/10);
      d_fak = pow(4.0*M_PI, 3) / (d_antenna_gain_abs_rx * d_antenna_gain_abs_tx * pow(d_lambda, 2));
    }

    void
    estimator_rcs_impl::set_usrp_gain_rx(float val){
      d_usrp_gain_rx = val;
    }

    void
    estimator_rcs_impl::set_power_tx(float val){
      d_power_tx = val;
    }

    void
    estimator_rcs_impl::set_corr_factor(float val){
      d_corr_factor = val;
    }

    float
    estimator_rcs_impl::calculate_vector_mean(boost::circular_buffer<float>* rcs_vals) {
      float sum_of_elems = 0;
      for(int k=0; k<rcs_vals->size(); k++){
        sum_of_elems += (*rcs_vals)[k];
      }
      return sum_of_elems/rcs_vals->size();
    }

    float
    estimator_rcs_impl::calculate_rcs() {
      // catch errors
			if(d_range.size() == 0) throw std::runtime_error("range vector has size zero");
			if(d_power.size() == 0) throw std::runtime_error("power vector has size zero");

      // regard usrp gain and signal path
			float power_rx = pow(d_power[0], d_exponent) / d_power_tx / pow(10, d_usrp_gain_rx/10);

			float fak = d_fak * pow(d_range[0], 4);

      // debug output
      // std::cout << "PowerTx: " << d_power_tx << std::endl;
      // std::cout << "PowerRx: " << power_rx << std::endl;
      // std::cout << "Lambda: " << d_lambda << std::endl;
      // std::cout << "GainRx: " << d_antenna_gain_rx << std::endl;
      // std::cout << "GainTx: " << d_antenna_gain_tx << std::endl;
      // std::cout << "fak: " << fak << std::endl;

			return power_rx/d_power_tx * fak * d_corr_factor;
    }

    void
    estimator_rcs_impl::handle_msg(pmt::pmt_t msg)
    {
        // Read msg from peak detector
        d_msg_hold.clear();
        pmt::pmt_t msg_part;
        bool found_range = false;
        bool found_power = false;
        for(int k=0; k<pmt::length(msg); k++){
          msg_part = pmt::nth(k,msg);
          if(pmt::symbol_to_string(pmt::nth(0,msg_part))=="range"){
            d_prange = pmt::nth(1,msg_part);
            d_msg_hold.push_back(msg_part);
            found_range = true;
          }
          else if(pmt::symbol_to_string(pmt::nth(0,msg_part))=="power"){
            d_ppower = pmt::nth(1,msg_part);
            found_power = true;
          }
          else{
            d_msg_hold.push_back(msg_part);
          }
        }

        if(not(found_range&&found_power)) throw std::runtime_error("range or power identifier (symbol) not found");

        d_range = pmt::f32vector_elements(d_prange);
        d_power = pmt::f32vector_elements(d_ppower);
        d_rcs.clear();

        if(d_range.size()!=d_power.size()) throw std::runtime_error("range and power vectors do not have same size");

        // Calculate RCS
        float rcs_mean = 0.0;
        if(d_range.size() == 0 && d_power.size() == 0){
          //std::cout << "ERROR: No target detected for RCS calculation" << std::endl;
        }
        else {
          d_rcs_vals.push_back(calculate_rcs());
        }
        //std::cout << "RCS_TEMP: " << calculate_rcs() << std::endl;

        if(d_loop_counter+1 >= d_num_mean) {
            rcs_mean = calculate_vector_mean(&d_rcs_vals);
        }
        else {
            d_loop_counter++;
        }

        for(int k=0; k<d_range.size(); k++){
            d_rcs.push_back(rcs_mean);
        }

        // Push pmt to output msg port
        d_rcs_key = pmt::string_to_symbol("rcs"); // identifier velocity
        d_rcs_value = pmt::init_f32vector(d_rcs.size(), d_rcs); // vector to pmt
        d_rcs_pack = pmt::list2(d_rcs_key, d_rcs_value); // make list for velocity information

        d_value = pmt::list1(d_rcs_pack);
        for(int k=0; k<d_msg_hold.size(); k++){
            d_value = pmt::list_add(d_value, d_msg_hold[k]);
        }

        message_port_pub(d_port_id_out,d_value); // publish message
    }

  } /* namespace radar */
} /* namespace gr */

max_stars_repo_name: Aekai/Wi-Mind

max_issues_repo_name: Aekai/Wi-Mind

max_forks_repo_name: Aekai/Wi-Mind


text:
#include <assert.h>
#include <stdlib.h>
#include <algorithm>
#include <ctime>
#include <fstream>
#include <boost/algorithm/string.hpp>
#include <boost/lambda/bind.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/random.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/tokenizer.hpp>
#include <boost/tuple/tuple.hpp>
#include "../common/read.h"

using namespace std;
using namespace boost;
using namespace boost::lambda;

int syntax_err() {
	L_(error) << "Usage:\nclassify\n	log_level(0,1,2,...) proj_name out_prefix\n	isoform_format isoforms_path g2i_format g2i_path gene_begin_idx gene_end_idx\n";
	return 1;
}

int lexical_cast_err() {
	L_(error) << "Lexical_cast error when converting arguments to numeric values";
	return 1;
}

int unknown_format_err(string const & format) {
	L_(error) << "Unknown file format error: " << format;
	return 1;
}

int unknown_readtype_err(string const & type) {
	L_(error) << "Unknown read type error: " << type;
	return 1;
}

class FileFormats {
	public:
		static const string UCSC_GENE2ISOFORM;
		static const string LH_GENE_TXT;
};
const string FileFormats::UCSC_GENE2ISOFORM = "UCSC_GENE2ISOFORM";
const string FileFormats::LH_GENE_TXT = "LH_GENE_TXT";

typedef tuple<string, string, interval_list<long> > gene_info;

int main(int argc, char * argv[])
{
	typedef tokenizer<char_separator<char> > tok;
	// initialize the logging system
	Log::ReportingLevel = LogLevel::info;

	// get parameters
	if (argc < 10) {
		return syntax_err();
	}
	int argi = 1;
	long log_level = lexical_cast<long>(argv[argi++]);
	Log::ReportingLevel = log_level;
	string proj_name = argv[argi++];
	string out_prefix = argv[argi++];

	// load isoform parameters
	string isoform_format = argv[argi++];
	string isoforms_path = argv[argi++];
	string g2i_format = argv[argi++];
	string g2i_path = argv[argi++];
	unsigned long gene_begin_idx, gene_end_idx;
	try {
		gene_begin_idx = lexical_cast<unsigned long>(argv[argi++]);
		gene_end_idx = lexical_cast<unsigned long>(argv[argi++]);

	} catch (bad_lexical_cast &) {
		return lexical_cast_err();
	}

	// global variables
	ifstream ifs;
	string line;

	// load isoforms
	L_(info) << "Loading isoforms...";
	ifs.clear();
	ifs.open(isoforms_path.c_str());
	assert(ifs.is_open() && !ifs.eof());
	vec_gap iso_gaps;
	if (isoform_format == FileFormats::LH_GENE_TXT) {
		while (getline(ifs, line) && !ifs.eof()) {
			istringstream iss(line);
			ga_ptr gap(new ga_entry());

			string strExonStarts, strExonEnds;
			iss >> gap->gname >> gap->chrom >> gap->strand
				>> gap->txStart >> gap->txEnd
				>> gap->exonCount >> strExonStarts >> strExonEnds;
			gap->cdsStart = gap->txStart;
			gap->cdsEnd = gap->txEnd;

			/* deal with exonStarts and exonEnds */
			char_separator<char> sep(",");
			tok tExonStarts(strExonStarts, sep);
			tok tExonEnds(strExonEnds, sep);

			for_each(tExonStarts.begin(), tExonStarts.end(),
					bind(&vector<long>::push_back,
						ref(gap->exonStarts),
						bind(atol,
							bind(&string::c_str, cref(_1)))));

			for_each(tExonEnds.begin(), tExonEnds.end(),
					bind(&vector<long>::push_back,
						ref(gap->exonEnds),
						bind(atol,
							bind(&string::c_str, cref(_1)))));

			iso_gaps.push_back(gap);
		}
	} else {
		return unknown_format_err(isoform_format);
	}
	ifs.close();
	map<string, ga_ptr> iname2gap;
	BOOST_FOREACH(ga_ptr const & gap, iso_gaps) {
		iname2gap[gap->gname] = gap;
	}
	L_(info) << "Loaded " << iso_gaps.size() << " isoforms";

	L_(info) << "Loading gene->isoform mappings ...";
	ifs.clear();
	ifs.open(g2i_path.c_str());
	assert(ifs.is_open() && !ifs.eof());
	map<string, vector<ga_ptr> > gname2isogaps;
	set<string> gname_set;
	if (g2i_format == FileFormats::UCSC_GENE2ISOFORM) {
		while (getline(ifs, line) && !ifs.eof()) {
			istringstream iss(line);
			string gname, iname;
			iss >> gname >> iname;
			gname_set.insert(gname);
			gname2isogaps[gname].push_back(iname2gap[iname]);
		}
	} else {
		return unknown_format_err(g2i_format);
	}
	ifs.close();
	L_(info) << "Loaded " << gname_set.size() << " genes";

	L_(info) << "Selecting genes: from #" << gene_begin_idx
		<< " to #" << gene_end_idx << ", ignoring single-isoform ones";
	vector<string> gnames;
	unsigned long count = 0;
	BOOST_FOREACH(string const & gname, gname_set) {
		if (count >= gene_begin_idx
				&& count < gene_end_idx) {
			if (gname2isogaps[gname].size() >= 2) 	// remove to report rpkm values for genes with single isoform as well
			{
				gnames.push_back(gname);
			}
		}
		count++;
	}
	L_(info) << "Selected " << gnames.size() << " gene(s)";
	L_(debug) << "Gene names:";
	BOOST_FOREACH(string const & gname, gnames) {
		L_(debug) << gname;
		BOOST_FOREACH(ga_ptr const gap, gname2isogaps[gname]) {
			L_(debug) << "  " << gap->gname;
		}
	}

	// For each gene, build its isoform structure...
	L_(info) << "Building isoform structures for the "
		<< gnames.size() << " selected gene(s)";
	map<string, interval_list<long> > gname2il;
	map<string, ExonSet> gname2exons;
	map<string, IsoformsPtr> gname2isop;
	BOOST_FOREACH(string const & gname, gnames) {
		// get isoforms
		vector<ga_ptr> const & iso_gaps = gname2isogaps[gname];
		// generate exon set
		ExonSet exons;
		IsoformsPtr isop(new Isoforms);
		BOOST_FOREACH(ga_ptr const & gap, iso_gaps) {
			for (unsigned long i = 0; i < gap->exonCount; ++i) {
				exons.insert(gap->exonStarts[i], gap->exonEnds[i]);
				gname2il[gname].add_interval(
						gap->exonStarts[i], gap->exonEnds[i]);
			}
		}
		gname2exons[gname] = exons;
		// build isoforms
		isop->build(exons, iso_gaps);
		gname2isop[gname] = isop;

        ofstream tmp_ofs;
        tmp_ofs.clear();
        string tmp_path = out_prefix + gname + ".matrix";
        tmp_ofs.open(tmp_path.c_str());
		assert(tmp_ofs.is_open() && !tmp_ofs.eof());

    	tmp_ofs << isop->chrom << "\t" << isop->strand << "\t" << exons.exonps << "\n";

		ublas::matrix<int> tmp_iso_array = ublas::zero_matrix<int>(isop->num_known_isoforms, isop->num_exons);
        unsigned long tmp_iso_idx = 0;
		vector<string> tmp_known_iso_names;
        BOOST_FOREACH (ga_ptr const & gap, iso_gaps) {
		  	tmp_known_iso_names.push_back(gap->gname);
		  	unsigned long tmp_exon_idx = 0;
		  	unsigned long tmp_iso_exon_idx = 0;
		  	BOOST_FOREACH (ExonPtr const & exonp, exons.exonps) {
		    	for (unsigned long i = tmp_iso_exon_idx;
			 					i < gap->exonCount; ++i) {
		      		if (exonp->is_inside(gap->exonStarts[i], gap->exonEnds[i])) {
						tmp_iso_array(tmp_iso_idx, tmp_exon_idx) = 1;
						tmp_iso_exon_idx = i;
						break;
		      		}
		    	}
                tmp_ofs << tmp_iso_array(tmp_iso_idx, tmp_exon_idx) << "\t";
		   		tmp_exon_idx++;
		 	}
        	tmp_ofs << "\n";
		 	tmp_iso_idx++;
    	}

		// no need to enumerate all possible isoforms for
		// the purpose of quantifying known isoforms
		// isop->enumerate_possible_isoforms();
		// isop->rebuild_possible_iso_probs();
	}
	L_(info) << "Built isoform structures for the "
		<< gnames.size() << " selected gene(s)";

	return 0;
}

ext: cpp

max_stars_repo_name: gersteinlab/LESSeq

max_issues_repo_name: gersteinlab/LESSeq

max_forks_repo_name: gersteinlab/LESSeq