<commit_before>// -*- mode: c++; coding: utf-8 -*-
/*! @file exact.cpp
    @brief Inplementation of ExactModel class
*/
#include "exact.hpp"
#include "util.hpp"

#include <functional>
#include <chrono>

#include <json.hpp>

#include <wtl/debug.hpp>
#include <wtl/exception.hpp>
#include <wtl/iostr.hpp>
#include <wtl/zfstream.hpp>
#include <wtl/algorithm.hpp>
#include <wtl/numeric.hpp>
#include <wtl/math.hpp>
#include <wtl/os.hpp>
#include <wtl/concurrent.hpp>

namespace likeligrid {

const std::vector<double> ExactModel::STEPS_ = {0.32, 0.16, 0.08, 0.04, 0.02, 0.01};
const std::vector<size_t> ExactModel::BREAKS_ = {5, 5, 5, 5, 5, 5};
bool ExactModel::SIGINT_RAISED_ = false;

ExactModel::ExactModel(const std::string& infile,
    const size_t max_sites,
    const unsigned int concurrency)
    : ExactModel(wtl::izfstream(infile), max_sites, concurrency) {HERE;}

ExactModel::ExactModel(
    std::istream& ist,
    const size_t max_sites,
    const unsigned int concurrency)
    : concurrency_(concurrency) {HERE;

    nlohmann::json jso;
    ist >> jso;
    names_ = jso["pathway"].get<std::vector<std::string>>();
    const size_t npath = names_.size();
    annot_.reserve(npath);
    for (const std::string& s: jso["annotation"]) {
        annot_.emplace_back(s);
    }
    std::cerr << "annot_: " << annot_ << std::endl;

    const size_t nsam = jso["sample"].size();
    std::vector<bits_t> all_genotypes;
    all_genotypes.reserve(nsam);
    for (const std::string& s: jso["sample"]) {
        all_genotypes.emplace_back(s);
    }

    genot_.reserve(nsam);
    const size_t ngene = jso["sample"].at(0).get<std::string>().size();
    nsam_with_s_.assign(ngene + 1, 0);  // at most
    std::valarray<double> s_gene(ngene);
    for (const auto& bits: all_genotypes) {
        const size_t s = bits.count();
        ++nsam_with_s_[s];
        if (s > max_sites) continue;
        genot_.push_back(bits);
        for (size_t j=0; j<ngene; ++j) {
            if (bits[j]) ++s_gene[j];
        }
    }
    wtl::rstrip(&nsam_with_s_);
    std::cerr << "Original N_s: " << nsam_with_s_ << std::endl;
    if (max_sites + 1 < nsam_with_s_.size()) {
        nsam_with_s_.resize(max_sites + 1);
        std::cerr << "Using N_s: " << nsam_with_s_ << std::endl;
    } else {
        std::cerr << "Note: -s is too large" << std::endl;
    }
    w_gene_ = s_gene / s_gene.sum();
    std::cerr << "s_gene : " << s_gene << std::endl;
    std::cerr << "w_gene_: " << w_gene_ << std::endl;

    mle_params_.resize(annot_.size());
    mle_params_ = 1.0;
}

void ExactModel::run(const bool writing) {HERE;
    while (stage_ < STEPS_.size()) {
        if (writing) {run_fout();} else {run_cout();}
    }
    --stage_;
    search_limits();
}

void ExactModel::run_fout() {HERE;
    const std::string outfile = init_meta();
    std::cerr << "mle_params_: " << mle_params_ << std::endl;
    if (outfile.empty()) return;
    const auto axes = make_vicinity(mle_params_, breaks(), radius());
    for (size_t j=0; j<names_.size(); ++j) {
        std::cerr << names_[j] << ": " << axes[j] << std::endl;
    }
    {
        wtl::ozfstream fout(outfile, std::ios::out | std::ios::app);
        std::cerr << "Writing: " << outfile << std::endl;
        run_impl(fout, wtl::itertools::product(axes));
    }
}

void ExactModel::run_cout() {HERE;
    const auto axes = make_vicinity(mle_params_, breaks(), radius());
    for (size_t j=0; j<names_.size(); ++j) {
        std::cerr << names_[j] << ": " << axes[j] << std::endl;
    }
    {
        std::stringstream sst;
        run_impl(sst, wtl::itertools::product(axes));
        std::cout << sst.str();
        read_results(sst);
    }
    std::cerr << "mle_params_: " << mle_params_ << std::endl;
    ++stage_;
}

void ExactModel::search_limits() const {HERE;
    {
        const std::vector<std::valarray<double>> axes(names_.size(), wtl::lin_spaced(200, 2.0, 0.01));
        wtl::ozfstream fout("uniaxis.tsv.gz");
        run_impl(fout, wtl::itertools::uniaxis(axes, mle_params_));
    }
    for (const auto& p: find_intersections(*this)) {
        std::cerr << p.first << ": " << p.second << std::endl;
        const auto axes = make_vicinity(p.second, 5, 0.02);
        wtl::ozfstream fout("limit-" + p.first + ".tsv.gz");
        //TODO: if exists
        run_impl(fout, wtl::itertools::product(axes));
    }
}

void ExactModel::run_impl(std::ostream& ost, wtl::itertools::Generator<std::valarray<double>>&& gen) const {HERE;
    std::cerr << skip_ << " to " << gen.max_count() << std::endl;
    if (skip_ == 0) {
        ost << "##max_count=" << gen.max_count() << "\n";
        ost << "##max_sites=" << nsam_with_s_.size() - 1 << "\n";
        ost << "##step=" << STEPS_.at(stage_) << "\n";
        ost << "loglik\t" << wtl::join(names_, "\t") << "\n";
        ost.flush();
    }

    wtl::Semaphore semaphore(concurrency_);
    auto task = [this,&semaphore](const std::valarray<double> th_path) {
        // argument is copied for thread
        auto buffer = wtl::make_oss();
        buffer << calc_loglik(th_path) << "\t"<commit_msg>Round uniaxis<commit_after>// -*- mode: c++; coding: utf-8 -*-
/*! @file exact.cpp
    @brief Inplementation of ExactModel class
*/
#include "exact.hpp"
#include "util.hpp"

#include <functional>
#include <chrono>

#include <json.hpp>

#include <wtl/debug.hpp>
#include <wtl/exception.hpp>
#include <wtl/iostr.hpp>
#include <wtl/zfstream.hpp>
#include <wtl/algorithm.hpp>
#include <wtl/numeric.hpp>
#include <wtl/math.hpp>
#include <wtl/os.hpp>
#include <wtl/concurrent.hpp>

namespace likeligrid {

const std::vector<double> ExactModel::STEPS_ = {0.32, 0.16, 0.08, 0.04, 0.02, 0.01};
const std::vector<size_t> ExactModel::BREAKS_ = {5, 5, 5, 5, 5, 5};
bool ExactModel::SIGINT_RAISED_ = false;

ExactModel::ExactModel(const std::string& infile,
    const size_t max_sites,
    const unsigned int concurrency)
    : ExactModel(wtl::izfstream(infile), max_sites, concurrency) {HERE;}

ExactModel::ExactModel(
    std::istream& ist,
    const size_t max_sites,
    const unsigned int concurrency)
    : concurrency_(concurrency) {HERE;

    nlohmann::json jso;
    ist >> jso;
    names_ = jso["pathway"].get<std::vector<std::string>>();
    const size_t npath = names_.size();
    annot_.reserve(npath);
    for (const std::string& s: jso["annotation"]) {
        annot_.emplace_back(s);
    }
    std::cerr << "annot_: " << annot_ << std::endl;

    const size_t nsam = jso["sample"].size();
    std::vector<bits_t> all_genotypes;
    all_genotypes.reserve(nsam);
    for (const std::string& s: jso["sample"]) {
        all_genotypes.emplace_back(s);
    }

    genot_.reserve(nsam);
    const size_t ngene = jso["sample"].at(0).get<std::string>().size();
    nsam_with_s_.assign(ngene + 1, 0);  // at most
    std::valarray<double> s_gene(ngene);
    for (const auto& bits: all_genotypes) {
        const size_t s = bits.count();
        ++nsam_with_s_[s];
        if (s > max_sites) continue;
        genot_.push_back(bits);
        for (size_t j=0; j<ngene; ++j) {
            if (bits[j]) ++s_gene[j];
        }
    }
    wtl::rstrip(&nsam_with_s_);
    std::cerr << "Original N_s: " << nsam_with_s_ << std::endl;
    if (max_sites + 1 < nsam_with_s_.size()) {
        nsam_with_s_.resize(max_sites + 1);
        std::cerr << "Using N_s: " << nsam_with_s_ << std::endl;
    } else {
        std::cerr << "Note: -s is too large" << std::endl;
    }
    w_gene_ = s_gene / s_gene.sum();
    std::cerr << "s_gene : " << s_gene << std::endl;
    std::cerr << "w_gene_: " << w_gene_ << std::endl;

    mle_params_.resize(annot_.size());
    mle_params_ = 1.0;
}

void ExactModel::run(const bool writing) {HERE;
    while (stage_ < STEPS_.size()) {
        if (writing) {run_fout();} else {run_cout();}
    }
    --stage_;
    search_limits();
}

void ExactModel::run_fout() {HERE;
    const std::string outfile = init_meta();
    std::cerr << "mle_params_: " << mle_params_ << std::endl;
    if (outfile.empty()) return;
    const auto axes = make_vicinity(mle_params_, breaks(), radius());
    for (size_t j=0; j<names_.size(); ++j) {
        std::cerr << names_[j] << ": " << axes[j] << std::endl;
    }
    {
        wtl::ozfstream fout(outfile, std::ios::out | std::ios::app);
        std::cerr << "Writing: " << outfile << std::endl;
        run_impl(fout, wtl::itertools::product(axes));
    }
}

void ExactModel::run_cout() {HERE;
    const auto axes = make_vicinity(mle_params_, breaks(), radius());
    for (size_t j=0; j<names_.size(); ++j) {
        std::cerr << names_[j] << ": " << axes[j] << std::endl;
    }
    {
        std::stringstream sst;
        run_impl(sst, wtl::itertools::product(axes));
        std::cout << sst.str();
        read_results(sst);
    }
    std::cerr << "mle_params_: " << mle_params_ << std::endl;
    ++stage_;
}

void ExactModel::search_limits() const {HERE;
    {
        auto axis = wtl::round(wtl::lin_spaced(200, 2.0, 0.01), 100);
        axis = (axis * 100.0).apply(std::round) / 100.0;
        const std::vector<std::valarray<double>> axes(names_.size(), axis);
        wtl::ozfstream fout("uniaxis.tsv.gz");
        run_impl(fout, wtl::itertools::uniaxis(axes, mle_params_));
    }
    for (const auto& p: find_intersections(*this)) {
        std::cerr << p.first << ": " << p.second << std::endl;
        const auto axes = make_vicinity(p.second, 5, 0.02);
        wtl::ozfstream fout("limit-" + p.first + ".tsv.gz");
        //TODO: if exists
        run_impl(fout, wtl::itertools::product(axes));
    }
}

void ExactModel::run_impl(std::ostream& ost, wtl::itertools::Generator<std::valarray<double>>&& gen) const {HERE;
    std::cerr << skip_ << " to " << gen.max_count() << std::endl;
    if (skip_ == 0) {
        ost << "##max_count=" << gen.max_count() << "\n";
        ost << "##max_sites=" << nsam_with_s_.size() - 1 << "\n";
        ost << "##step=" << STEPS_.at(stage_) << "\n";
        ost << "loglik\t" << wtl::join(names_, "\t") << "\n";
        ost.flush();
    }

    wtl::Semaphore semaphore(concurrency_);
    auto task = [this,&semaphore](const std::valarray<double> th_path) {
        // argument is copied for thread
        auto buffer = wtl::make_oss();
        buffer << calc_loglik(th_path) << "\t"
{'fraction_non_alphanumeric': 0.1907085691977357, 'fraction_numerical': 0.010638297872340425, 'mean_word_length': 3.2064860426929394, 'pattern_counts': {'":': 6, '<': 289, '<?xml version=': 0, '>': 57, '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': '7168897', 'n_tokens_mistral': 4013, 'n_tokens_neox': 3862, 'n_words': 752}