text:
/*
 * Copyright 2009-2020 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
#define BOOST_TEST_MAIN

#define BOOST_TEST_MODULE qmmolecule_test

// Third party includes
#include <boost/test/unit_test.hpp>

// Local VOTCA includes
#include "votca/xtp/qmmolecule.h"

using namespace votca::xtp;
using namespace std;

BOOST_AUTO_TEST_SUITE(qmmolecule_test)

BOOST_AUTO_TEST_CASE(constructors_test) { QMMolecule seg("seg1", 1); }

BOOST_AUTO_TEST_CASE(load_xyz_test) {
  ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << endl;
  xyzfile << " methane" << endl;
  xyzfile << " C            .000000     .000000     .000000" << endl;
  xyzfile << " H            .629118     .629118     .629118" << endl;
  xyzfile << " H           -.629118    -.629118     .629118" << endl;
  xyzfile << " H            .629118    -.629118    -.629118" << endl;
  xyzfile << " H           -.629118     .629118    -.629118" << endl;
  xyzfile.close();

  QMMolecule seg("seg1", 1);
  seg.LoadFromFile("molecule.xyz");

  auto extension = seg.CalcSpatialMinMax();
  Eigen::Vector3d max(0.629118, 0.629118, 0.629118);
  max *= votca::tools::conv::ang2bohr;
  Eigen::Vector3d min = -max;

  bool check_min = extension.first.isApprox(min, 0.00001);

  if (!check_min) {
    std::cout << "ref" << std::endl;
    std::cout << min << std::endl;
    std::cout << "result" << std::endl;
    std::cout << extension.first << std::endl;
  }
  BOOST_CHECK_EQUAL(check_min, true);
  bool check_max = extension.second.isApprox(max, 0.00001);

  if (!check_max) {
    std::cout << "ref" << std::endl;
    std::cout << max << std::endl;
    std::cout << "result" << std::endl;
    std::cout << extension.second << std::endl;
  }
  BOOST_CHECK_EQUAL(check_max, true);

  seg.WriteXYZ("moltest.xyz", "heloo");
}

BOOST_AUTO_TEST_CASE(readwritehdf) {
  ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << endl;
  xyzfile << " methane" << endl;
  xyzfile << " C            .000000     .000000     .000000" << endl;
  xyzfile << " H            .629118     .629118     .629118" << endl;
  xyzfile << " H           -.629118    -.629118     .629118" << endl;
  xyzfile << " H            .629118    -.629118    -.629118" << endl;
  xyzfile << " H           -.629118     .629118    -.629118" << endl;
  xyzfile.close();

  QMMolecule seg("seg1", 1);
  seg.LoadFromFile("molecule.xyz");

  CheckpointFile ff("qmmolecule.hdf5");
  CheckpointWriter ww = ff.getWriter();
  seg.WriteToCpt(ww);

  CheckpointReader rr = ff.getReader();
  QMMolecule seg2(rr);
}

BOOST_AUTO_TEST_SUITE_END()

ext: cc

max_stars_repo_name: fossabot/xtp

max_issues_repo_name: fossabot/xtp

max_forks_repo_name: fossabot/xtp


text:
/*===================================================================

The Medical Imaging Interaction Toolkit (MITK)

Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.

This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.

See LICENSE.txt or http://www.mitk.org for details.

===================================================================*/

#include <mitkCoreObjectFactory.h>
#include "mitkImage.h"
#include <boost/lexical_cast.hpp>
#include "mitkCommandLineParser.h"
#include <mitkIOUtil.h>
#include <itksys/SystemTools.hxx>
#include <mitkITKImageImport.h>
#include <mitkImageCast.h>
#include <mitkProperties.h>

// ITK
#include <itkImageRegionIterator.h>

// MITK
#include <mitkIOUtil.h>

// Classification
#include <mitkCLUtil.h>
#include <mitkVigraRandomForestClassifier.h>

#include <QDir>
#include <QString>
#include <QStringList>

using namespace mitk;

/**
 *
 */
int main(int argc, char* argv[])
{
  mitkCommandLineParser parser;
  parser.setArgumentPrefix("--", "-");

  // required params
  parser.addArgument("inputdir", "i", mitkCommandLineParser::InputDirectory, "Input Directory", "Contains input feature files.", us::Any(), false);
  parser.addArgument("outputdir", "o", mitkCommandLineParser::OutputDirectory, "Output Directory", "Destination of output files.", us::Any(), false);
  parser.addArgument("classmask", "m", mitkCommandLineParser::InputFile, "Class mask image", "Contains several classes.", us::Any(), false);

  // optional params
  parser.addArgument("select", "s", mitkCommandLineParser::String, "Item selection", "Using Regular expression, seperated by space e.g.: '*.nrrd *.vtk *test*'",std::string("*.nrrd"),true);
  parser.addArgument("treecount", "tc", mitkCommandLineParser::Int, "Treecount", "Number of trees.",50,true);
  parser.addArgument("treedepth", "td", mitkCommandLineParser::Int, "Treedepth", "Maximal tree depth.",50,true);
  parser.addArgument("minsplitnodesize", "min", mitkCommandLineParser::Int, "Minimum split node size.", "Minimum split node size.",2,true);
  parser.addArgument("precision", "p", mitkCommandLineParser::Float, "Split precision.", "Precision.", mitk::eps,true);
  parser.addArgument("fraction", "f", mitkCommandLineParser::Float, "Fraction of samples per tree.", "Fraction of samples per tree.", 0.6f,true);
  parser.addArgument("replacment", "r", mitkCommandLineParser::Bool, "Sample with replacement.", "Sample with replacement.", true,true);

  // Miniapp Infos
  parser.setCategory("Classification Tools");
  parser.setTitle("Random Forest Training");
  parser.setDescription("Vigra RF impl.");
  parser.setContributor("MBI");

  // Params parsing
  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
  if (parsedArgs.size()==0)
    return EXIT_FAILURE;

  std::string inputdir = us::any_cast<std::string>(parsedArgs["inputdir"]);
  std::string outputdir = us::any_cast<std::string>(parsedArgs["outputdir"]);
  std::string classmask = us::any_cast<std::string>(parsedArgs["classmask"]);

  int treecount = parsedArgs.count("treecount") ? us::any_cast<int>(parsedArgs["treecount"]) : 50;
  int treedepth = parsedArgs.count("treedepth") ? us::any_cast<int>(parsedArgs["treedepth"]) : 50;
  int minsplitnodesize = parsedArgs.count("minsplitnodesize") ? us::any_cast<int>(parsedArgs["minsplitnodesize"]) : 2;
  float precision = parsedArgs.count("precision") ? us::any_cast<float>(parsedArgs["precision"]) : mitk::eps;
  float fraction = parsedArgs.count("fraction") ? us::any_cast<float>(parsedArgs["fraction"]) : 0.6;
  bool withreplacement = parsedArgs.count("replacment") ? us::any_cast<float>(parsedArgs["replacment"]) : true;
  std::string filt_select =/* parsedArgs.count("select") ? us::any_cast<std::string>(parsedArgs["select"]) :*/ "*.nrrd";

  QString filter(filt_select.c_str());

  // **** in principle repeat this block to create a feature matrix X_all for all subjects (in dir)
  // Get nrrd filepath
  QDir dir(inputdir.c_str());
  auto strl = dir.entryList(filter.split(" "),QDir::Files);

  // load class mask
  mitk::Image::Pointer mask = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(classmask)[0].GetPointer());
  unsigned int num_samples = 0;
  mitk::CLUtil::CountVoxel(mask,num_samples);

  // initialize featurematrix [num_samples, num_featureimages]
  Eigen::MatrixXd X(num_samples, strl.size());

  for(int i = 0 ; i < strl.size(); i++)
  {
    // load feature image
    mitk::Image::Pointer img = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(inputdir + strl[i].toStdString())[0].GetPointer());
    // transfom it into a [num_samples, 1] vector depending on the classmask
    Eigen::MatrixXd _x = mitk::CLUtil::Transform<double>(img,mask);
    // replace i-th (empty) col with feature vector in _x
    X.block(0,i,num_samples,1) = _x;
  }
  // ****

  // transform classmask into the label-vector [num_samples, 1]
  Eigen::MatrixXi Y = mitk::CLUtil::Transform<int>(mask,mask);

  mitk::VigraRandomForestClassifier::Pointer classifier = mitk::VigraRandomForestClassifier::New();
  classifier->SetTreeCount(treecount);
  classifier->SetMaximumTreeDepth(treedepth);
  classifier->SetMinimumSplitNodeSize(minsplitnodesize);
  classifier->SetPrecision(precision);
  classifier->SetSamplesPerTree(fraction);
  classifier->UseSampleWithReplacement(withreplacement);

  classifier->PrintParameter();
  classifier->Train(X,Y);

  MITK_INFO << classifier->IsEmpty();

  // no metainformations are saved currently
  // only the raw vigra rf data
  mitk::IOUtil::Save(classifier, outputdir + "RandomForest.hdf5");

  Eigen::MatrixXi Y_pred = classifier->Predict(X);
  Eigen::MatrixXd Probs = classifier->GetPointWiseProbabilities();

  MITK_INFO << Y_pred.rows() << " " << Y_pred.cols();
  MITK_INFO << Probs.rows() << " " << Probs.cols();

  //    mitk::Image::Pointer prediction = mitk::CLUtil::Transform<int>(Y_pred,mask);
  mitk::Image::Pointer probs_1 = mitk::CLUtil::Transform<double>(Probs.col(0),mask);
  mitk::Image::Pointer probs_2 = mitk::CLUtil::Transform<double>(Probs.col(1),mask);
  mitk::Image::Pointer probs_3 = mitk::CLUtil::Transform<double>(Probs.col(2),mask);

  mitk::IOUtil::Save(probs_1, outputdir + "probs_1.nrrd");
  mitk::IOUtil::Save(probs_2, outputdir + "probs_2.nrrd");
  mitk::IOUtil::Save(probs_3, outputdir + "probs_3.nrrd");
  //    mitk::IOUtil::Save(probs_2, outputdir + "test.nrrd");

  return EXIT_SUCCESS;
}

ext: cpp

max_stars_repo_name: ZP-Hust/MITK

max_issues_repo_name: ZP-Hust/MITK

max_forks_repo_name: ZP-Hust/MITK


text:
#ifndef BURST__ITERATOR__BUFFERED_CHUNK_ITERATOR_HPP
#define BURST__ITERATOR__BUFFERED_CHUNK_ITERATOR_HPP

#include <burst/algorithm/copy_at_most_n.hpp>
#include <burst/iterator/end_tag.hpp>
#include <burst/type_traits/iterator_difference.hpp>
#include <burst/type_traits/iterator_value.hpp>

#include <boost/iterator/iterator_facade.hpp>
#include <boost/range/iterator_range.hpp>

#include <cassert>
#include <cstddef>
#include <iterator>
#include <memory>
#include <tuple>

namespace burst
{
    /*!
        \brief
            Буферизованный кусочный итератор

        \details
            Разбивает входящую последовательность на равные куски (последний кусок может быть
            меньше, если элементов не хватает до целого) заданного размера.
            Текущий кусок хранится в буфере внутри итератора.

            Моделирует однопроходный итератор (InputIterator).

        \tparam Iterator
            Тип исходного итератора.
        \tparam Integral
            Тип числа, с помощью которого будет задаваться размер одного куска.

        \see make_buffered_chunk_iterator
     */
    template <typename Iterator, typename Integral>
    class buffered_chunk_iterator:
        public boost::iterator_facade
        <
            buffered_chunk_iterator<Iterator, Integral>,
            boost::iterator_range<iterator_value_t<Iterator> *>,
            std::input_iterator_tag,
            boost::iterator_range<iterator_value_t<Iterator> *>,
            iterator_difference_t<Iterator>
        >
    {
    private:
        using buffer_value_type = iterator_value_t<Iterator>;
        using base_type =
            boost::iterator_facade
            <
                buffered_chunk_iterator,
                boost::iterator_range<buffer_value_type *>,
                std::input_iterator_tag,
                boost::iterator_range<buffer_value_type *>,
                iterator_difference_t<Iterator>
            >;

    public:
        using reference = typename base_type::reference;
        using difference_type = typename base_type::difference_type;

        buffered_chunk_iterator (Iterator first, Iterator last, Integral chunk_size):
            m_current(first),
            m_end(last),
            m_chunk_size(static_cast<difference_type>(chunk_size)),
            m_chunk
            (
                std::shared_ptr<buffer_value_type>
                (
                    new buffer_value_type[static_cast<std::size_t>(m_chunk_size)],
                    std::default_delete<buffer_value_type[]>{}
                )
            ),
            m_chunk_end(m_chunk.get())
        {
            assert(m_chunk_size > 0);
            increment();
        }

        buffered_chunk_iterator (iterator::end_tag_t, const buffered_chunk_iterator & begin):
            m_current(begin.m_end),
            m_end(begin.m_end),
            m_chunk_size(begin.m_chunk_size),
            m_chunk{begin.m_chunk},
            m_chunk_end(m_chunk.get())
        {
        }

    private:
        friend class boost::iterator_core_access;

        bool empty () const
        {
            return m_chunk.get() == m_chunk_end;
        }

        reference dereference () const
        {
            return boost::make_iterator_range(m_chunk.get(), m_chunk_end);
        }

        bool equal (const buffered_chunk_iterator & that) const
        {
            assert(this->m_end == that.m_end);
            assert(this->m_chunk == that.m_chunk);
            assert(this->m_chunk_size == that.m_chunk_size);
            return this->m_current == that.m_current && empty();
        }

        void increment ()
        {
            std::tie(m_current, std::ignore, m_chunk_end) =
                copy_at_most_n(m_current, m_end, m_chunk_size, m_chunk.get());
        }

        Iterator m_current;
        Iterator m_end;
        difference_type m_chunk_size;
        std::shared_ptr<buffer_value_type> m_chunk;
        buffer_value_type * m_chunk_end;
    };

    /*!
        \brief
            Функция для создания буферизованного кусочного итератора

        \param [first, last)
            Последовательность, которая будет побита на куски.
        \param n
            Размер куска.

        \returns
            Экземпляр буферизованного кусочного итератора.

        \see buffered_chunk_iterator
     */
    template <typename Iterator, typename Integral>
    auto make_buffered_chunk_iterator (Iterator first, Iterator last, Integral n)
    {
        return buffered_chunk_iterator<Iterator, Integral>(first, last, n);
    }

    /*!
        \brief
            Функция для создания буферизованного кусочного итератора

        \param range
            Последовательность, которая будет побита на куски.
        \param n
            Размер куска.

        \warning
            Время жизни диапазона `range` должно быть не меньше времени жизни полученного
            буферизованного кусочного итератора, в противном случае поведение неопределено.

        \returns
            Экземпляр буферизованного кусочного итератора.

        \see buffered_chunk_iterator
     */
    template <typename Range, typename Integral>
    auto make_buffered_chunk_iterator (const Range & range, Integral n)
    {
        using std::begin;
        using std::end;
        return make_buffered_chunk_iterator(begin(range), end(range), n);
    }

    /*!
        \brief
            Функция для создания конца буферизованного кусочного диапазона

        \details
            Принимает на вход итератор на начало кусочного диапазона и специальную метку,
            обозначающую, что нужно создать итератор на конец диапазона.

        \returns
            Итератор на конец буферизованного кусочного диапазона.

        \see buffered_chunk_iterator
     */
    template <typename Iterator, typename Integral>
    auto
        make_buffered_chunk_iterator
        (
            iterator::end_tag_t,
            const buffered_chunk_iterator<Iterator, Integral> & begin
        )
    {
        return buffered_chunk_iterator<Iterator, Integral>(iterator::end_tag, begin);
    }
} // namespace burst

#endif // BURST__ITERATOR__BUFFERED_CHUNK_ITERATOR_HPP

ext: hpp

max_stars_repo_name: izvolov/thrust

max_issues_repo_name: izvolov/burst

max_forks_repo_name: izvolov/thrust


text:
#include "Commandline.hxx"

#include <boost/program_options.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <iostream>

Command parse_options(int argc, const char **argv) {
    using namespace boost::program_options;

    // The command line variables
    std::string config, single_page, output, info;
    bool human;
    int start = -1, end = -1, page = -1;

    // All available options
    options_description desc("Allowed options");
    desc.add_options()
        ("help", "produce help message")
        ("config,c", value<std::string>(&config))
        ("human,h", bool_switch(&human)->default_value(false))
        ("single_page,S", value<std::string>(&single_page))
        ("output,o", value<std::string>(&output))
        ("info,i", value<std::string>(&info))
        ("start,s", value<int>(&start))
        ("end,e", value<int>(&end))
        ("page,p", value<int>(&page));

    variables_map vm;
    store(parse_command_line(argc, argv, desc), vm);
    notify(vm);

    // Only print help text if --help or -h is specified
    if (vm.count("help") || argc == 1) {
        return HelpCommand{};
    }

    // We need a valid config, so we exit if there is none
    if (config.empty()) {
        throw std::runtime_error(
            "no configuration file set\n\n" + potoo_description);
    }

    // config is set, but output and human are not
    if (!config.empty() && !(!output.empty() || human || !single_page.empty() || !info.empty())) {
        throw std::runtime_error(
            "please either specify an output file, the human flag, the info or the single_page parameter"
        );
    }

    if (
        (vm.count("start") && start <= -1)
        || (vm.count("end") && end <= -1)
        || (vm.count("page") && page <= -1)
        ) {
        throw std::runtime_error("invalid start or end range or page");
    }

    boost::algorithm::trim(config);
    boost::algorithm::trim(single_page);
    boost::algorithm::trim(output);

    const static auto integer_to_optional = [](int param) -> boost::optional<int> {
        if (param == -1) return boost::none;
        else return boost::make_optional(param);
    };

    if (!single_page.empty()) { // single_page parameter was supplied, so we just save the image

        if (vm.count("start") || vm.count("end")) {
            throw std::runtime_error("start and end are invalid parameters for single_page");
        }

        PageCommand fp;
        fp._config = config;
        fp._path = single_page;
        fp._page = integer_to_optional(page);
        return fp;
    } else if (!info.empty()) {
        InfoCommand ic;
        ic._config = config;
        ic._path = info;
        return ic;
    }
    else { // single_page parameter was not supplied, run the main routine

        if ((vm.count("start") || vm.count("end")) && vm.count("page")) {
            throw std::runtime_error("start and/or end and page cannot be used at the same time");
        }

        if (human) { // We want to print for humans, so let's config that
            HumanCommand hc;
            hc._config = config;
            hc._start = integer_to_optional(start);
            hc._end = integer_to_optional(end);
            hc._page = integer_to_optional(page);
            return hc;
        } else { // Return the path for saving
            OutputCommand oc;
            oc._config = config;
            oc._path = output;
            oc._start = integer_to_optional(start);
            oc._end = integer_to_optional(end);
            oc._page = integer_to_optional(page);
            return oc;
        }
    }
}

ext: cxx

max_stars_repo_name: casavi/potoo

max_issues_repo_name: casavi/potoo

max_forks_repo_name: casavi/potoo


text:
#include "ForestParser.h"

#include <istream>
#include <string>

#include <boost/make_shared.hpp>

#include "util/tokenize_piece.hh"

#include "syntax-common/exception.h"

namespace MosesTraining
{
namespace Syntax
{
namespace PostprocessEgretForests
{

ForestParser::ForestParser()
  : m_input(0)
{
}

ForestParser::ForestParser(std::istream &input)
  : m_input(&input)
{
  ++(*this);
}

ForestParser &ForestParser::operator++()
{
  if (!m_input) {
    return *this;
  }
  m_vertexSet.clear();
  m_entry.forest.vertices.clear();
  if (!std::getline(*m_input, m_tmpLine)) {
    m_input = 0;
    return *this;
  }
  // The first line contains the sentence number.
  ParseSentenceNumLine(m_tmpLine, m_entry.sentNum);
  // The second line contains the sentence string.
  std::getline(*m_input, m_entry.sentence);
  // Subsequent lines contain hyperedges -- or a blank line if there was a
  // parse failure -- terminated by a blank line.
  std::getline(*m_input, m_tmpLine);
  if (m_tmpLine == "") {  // Parse failure
    std::getline(*m_input, m_tmpLine);
    assert(m_tmpLine == "");
    return *this;
  }
  while (m_tmpLine != "") {
    ParseHyperedgeLine(m_tmpLine, m_entry.forest);
    std::getline(*m_input, m_tmpLine);
  }
  return *this;
}

boost::shared_ptr<Forest::Vertex> ForestParser::AddVertex(const VertexSP &v)
{
  std::pair<VertexSet::iterator, bool> ret = m_vertexSet.insert(v);
  if (ret.second) {
    m_entry.forest.vertices.push_back(*ret.first);
  }
  return *ret.first;
}

void ForestParser::ParseSentenceNumLine(const std::string &line,
                                        std::size_t &sentNum)
{
  const util::AnyCharacter delimiter(" \t");
  util::TokenIter<util::AnyCharacter, true> p(line, delimiter);
  if (*p != "sentence") {
    // FIXME
    throw Exception("");
  }
  ++p;
  std::string tmp;
  p->CopyToString(&tmp);
  sentNum = std::atoi(tmp.c_str());
}

void ForestParser::ParseHyperedgeLine(const std::string &line, Forest &forest)
{
  const util::AnyCharacter delimiter(" \t");
  util::TokenIter<util::AnyCharacter, true> p(line, delimiter);
  VertexSP v = AddVertex(ParseVertex(*p));
  HyperedgeSP e = boost::make_shared<Forest::Hyperedge>();
  e->head = v.get();
  ++p;
  if (*p != "=>") {
    // FIXME
    throw Exception("");
  }
  for (++p; *p != "|||"; ++p) {
    v = ParseVertex(*p);
    if (v->start == -1) {
      // Egret does not give start/end for terminals.
      v->start = v->end = e->head->start;
    }
    e->tail.push_back(AddVertex(v).get());
  }
  ++p;
  std::string tmp;
  p->CopyToString(&tmp);
  e->weight = std::atof(tmp.c_str());
  e->head->incoming.push_back(e);
}

boost::shared_ptr<Forest::Vertex> ForestParser::ParseVertex(
  const StringPiece &s)
{
  VertexSP v = boost::make_shared<Forest::Vertex>();
  std::size_t pos = s.rfind('[');
  if (pos == std::string::npos) {
    s.CopyToString(&v->symbol.value);
    v->symbol.isNonTerminal = false;
    v->start = v->end = -1;
    return v;
  }
  if (pos > 2 && s[pos-2] == '^' && s[pos-1] == 'g') {
    s.substr(0, pos-2).CopyToString(&v->symbol.value);
  } else {
    s.substr(0, pos).CopyToString(&v->symbol.value);
  }
  v->symbol.isNonTerminal = true;
  std::size_t begin = pos + 1;
  pos = s.find(',', begin+1);
  std::string tmp;
  s.substr(begin, pos-begin).CopyToString(&tmp);
  v->start = std::atoi(tmp.c_str());
  s.substr(pos+1, s.size()-pos-2).CopyToString(&tmp);
  v->end = std::atoi(tmp.c_str());
  return v;
}

bool operator==(const ForestParser &lhs, const ForestParser &rhs)
{
  // TODO Is this right?  Compare values of istreams if non-zero?
  return lhs.m_input == rhs.m_input;
}

bool operator!=(const ForestParser &lhs, const ForestParser &rhs)
{
  return !(lhs == rhs);
}

}  // namespace PostprocessEgretForests
}  // namespace Syntax
}  // namespace MosesTraining

max_issues_repo_name: Pangeamt/nectm

max_forks_repo_name: Pangeamt/nectm


text:
/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8
// archive_pointer_oserializer.ipp:

// (C) Copyright 2002 Robert Ramey - http://www.rrsd.com .
// Use, modification and distribution is subject to the Boost Software
// License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

//  See http://www.boost.org for updates, documentation, and revision history.

#include <utility>
#include <cassert>

#include <boost/config.hpp> // msvc 6.0 needs this for warning suppression

#include <boost/archive/detail/archive_pointer_oserializer.hpp>
#include <boost/archive/detail/basic_serializer_map.hpp>
#include <boost/serialization/singleton.hpp>

namespace boost {
namespace archive {
namespace detail {

namespace { // anon
    template<class Archive>
    class oserializer_map : public basic_serializer_map
    {
    };
}

template<class Archive>
BOOST_ARCHIVE_OR_WARCHIVE_DECL(BOOST_PP_EMPTY())
archive_pointer_oserializer<Archive>::archive_pointer_oserializer(
    const boost::serialization::extended_type_info & eti
) :
    basic_pointer_oserializer(eti)
{
    std::pair<
        BOOST_DEDUCED_TYPENAME oserializer_map<Archive>::iterator,
        bool
    > result;
    result = serialization::singleton<
            oserializer_map<Archive>
        >::get_mutable_instance().insert(this);
    assert(result.second);
}

template<class Archive>
BOOST_ARCHIVE_OR_WARCHIVE_DECL(const basic_pointer_oserializer *)
archive_pointer_oserializer<Archive>::find(
    const boost::serialization::extended_type_info & eti
){
    const basic_serializer_arg bs(eti);
    basic_serializer_map::const_iterator it;
    it =  boost::serialization::singleton<
            oserializer_map<Archive>
        >::get_const_instance().find(& bs);
    assert(
        it
        !=
        boost::serialization::singleton<
                oserializer_map<Archive>
            >::get_const_instance().end()
    );
    return static_cast<const basic_pointer_oserializer *>(*it);
}

template<class Archive>
BOOST_ARCHIVE_OR_WARCHIVE_DECL(BOOST_PP_EMPTY())
archive_pointer_oserializer<Archive>::~archive_pointer_oserializer(){
    // note: we need to check that the map still exists as we can't depend
    // on static variables being constructed in a specific sequence
    if(! serialization::singleton<
            oserializer_map<Archive>
        >::is_destroyed()
    ){
        unsigned int count;
        count = serialization::singleton<
                oserializer_map<Archive>
            >::get_mutable_instance().erase(this);
        assert(count);
    }
}

} // namespace detail
} // namespace archive
} // namespace boost

ext: ipp


text:
#ifndef _CELERITE2_INTERFACE_HPP_DEFINED_
#define _CELERITE2_INTERFACE_HPP_DEFINED_

#include <Eigen/Core>
#include "forward.hpp"

namespace celerite2 {
namespace core {

#define MakeEmptyWork(BaseType)                                                                                \
  typedef typename BaseType::Scalar Scalar;                                                                                \
  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Empty;                                                              \
  Empty

/**
 * \brief Compute the Cholesky factorization of the system
 *
 * This computes `d` and `W` such that:
 *
 * `diag(a) + tril(U*V^T) + triu(V*U^T) = L*diag(d)*L^T`
 *
 * where
 *
 * `L = 1 + tril(U*W^T)`
 *
 * This can be safely applied in place: `d_out` can point to `a` and `W_out` can
 * point to `V`, and the memory will be reused.
 *
 * @param t     (N,): The input coordinates (must be sorted)
 * @param c     (J,): The transport coefficients
 * @param a     (N,): The diagonal component
 * @param U     (N, J): The first low rank matrix
 * @param V     (N, J): The second low rank matrix
 * @param d_out (N,): The diagonal component of the Cholesky factor
 * @param W_out (N, J): The second low rank component of the Cholesky factor
 */
template <typename Input, typename Coeffs, typename Diag, typename LowRank, typename DiagOut, typename LowRankOut>
Eigen::Index factor(const Eigen::MatrixBase<Input> &t,         // (N,)
                    const Eigen::MatrixBase<Coeffs> &c,        // (J,)
                    const Eigen::MatrixBase<Diag> &a,          // (N,)
                    const Eigen::MatrixBase<LowRank> &U,       // (N, J)
                    const Eigen::MatrixBase<LowRank> &V,       // (N, J)
                    Eigen::MatrixBase<DiagOut> const &d_out,   // (N,)
                    Eigen::MatrixBase<LowRankOut> const &W_out // (N, J)
) {
  MakeEmptyWork(Diag) S;
  return factor<false>(t, c, a, U, V, d_out, W_out, S);
}

/**
 * \brief Compute the solution of a lower triangular linear equation
 *
 * This computes `Z` such that:
 *
 * `Y = L * Y`
 *
 * where
 *
 * `L = 1 + tril(U*W^T)`
 *
 * This can be safely applied in place.
 *
 * @param t     (N,): The input coordinates (must be sorted)
 * @param c     (J,): The transport coefficients
 * @param U     (N, J): The first low rank matrix
 * @param W     (N, J): The second low rank matrix
 * @param Y     (N, Nrhs): The right hand side
 * @param Z_out (N, Nrhs): The solution of this equation
 */
template <typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void solve_lower(const Eigen::MatrixBase<Input> &t,               // (N,)
                 const Eigen::MatrixBase<Coeffs> &c,              // (J,)
                 const Eigen::MatrixBase<LowRank> &U,             // (N, J)
                 const Eigen::MatrixBase<LowRank> &W,             // (N, J)
                 const Eigen::MatrixBase<RightHandSide> &Y,       // (N, nrhs)
                 Eigen::MatrixBase<RightHandSideOut> const &Z_out // (N, nrhs)
) {
  MakeEmptyWork(Input) F;
  solve_lower<false>(t, c, U, W, Y, Z_out, F);
}

/**
 * \brief Compute the solution of a upper triangular linear equation
 *
 * This computes `Z` such that:
 *
 * `Y = L^T * Y`
 *
 * where
 *
 * `L = 1 + tril(U*W^T)`
 *
 * This can be safely applied in place.
 *
 * @param t     (N,): The input coordinates (must be sorted)
 * @param c     (J,): The transport coefficients
 * @param U     (N, J): The first low rank matrix
 * @param W     (N, J): The second low rank matrix
 * @param Y     (N, Nrhs): The right hand side
 * @param Z_out (N, Nrhs): The solution of this equation
 */
template <typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void solve_upper(const Eigen::MatrixBase<Input> &t,               // (N,)
                 const Eigen::MatrixBase<Coeffs> &c,              // (J,)
                 const Eigen::MatrixBase<LowRank> &U,             // (N, J)
                 const Eigen::MatrixBase<LowRank> &W,             // (N, J)
                 const Eigen::MatrixBase<RightHandSide> &Y,       // (N, nrhs)
                 Eigen::MatrixBase<RightHandSideOut> const &Z_out // (N, nrhs)
) {
  MakeEmptyWork(Input) F;
  solve_upper<false>(t, c, U, W, Y, Z_out, F);
}

/**
 * \brief Apply a strictly lower matrix multiply
 *
 * This computes:
 *
 * `Z += tril(U * V^T) * Y`
 *
 * where `tril` is the strictly lower triangular function.
 *
 * Note that this will *update* the value of `Z`.
 *
 * @param t     (N,): The input coordinates (must be sorted)
 * @param c     (J,): The transport coefficients
 * @param U     (N, J): The first low rank matrix
 * @param V     (N, J): The second low rank matrix
 * @param Y     (N, Nrhs): The matrix to be multiplied
 * @param Z_out (N, Nrhs): The matrix to be updated
 */
template <typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void matmul_lower(const Eigen::MatrixBase<Input> &t,               // (N,)
                  const Eigen::MatrixBase<Coeffs> &c,              // (J,)
                  const Eigen::MatrixBase<LowRank> &U,             // (N, J)
                  const Eigen::MatrixBase<LowRank> &V,             // (N, J)
                  const Eigen::MatrixBase<RightHandSide> &Y,       // (N, nrhs)
                  Eigen::MatrixBase<RightHandSideOut> const &Z_out // (N, nrhs)
) {
  MakeEmptyWork(Input) F;
  matmul_lower<false>(t, c, U, V, Y, Z_out, F);
}

/**
 * \brief Apply a strictly upper matrix multiply
 *
 * This computes:
 *
 * `Z += triu(V * U^T) * Y`
 *
 * where `triu` is the strictly lower triangular function.
 *
 * Note that this will *update* the value of `Z`.
 *
 * @param t     (N,): The input coordinates (must be sorted)
 * @param c     (J,): The transport coefficients
 * @param U     (N, J): The first low rank matrix
 * @param V     (N, J): The second low rank matrix
 * @param Y     (N, Nrhs): The matrix to be multiplied
 * @param Z_out (N, Nrhs): The matrix to be updated
 */
template <typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void matmul_upper(const Eigen::MatrixBase<Input> &t,               // (N,)
                  const Eigen::MatrixBase<Coeffs> &c,              // (J,)
                  const Eigen::MatrixBase<LowRank> &U,             // (N, J)
                  const Eigen::MatrixBase<LowRank> &V,             // (N, J)
                  const Eigen::MatrixBase<RightHandSide> &Y,       // (N, nrhs)
                  Eigen::MatrixBase<RightHandSideOut> const &Z_out // (N, nrhs)
) {
  MakeEmptyWork(Input) F;
  matmul_upper<false>(t, c, U, V, Y, Z_out, F);
}

/**
 * \brief The general lower-triangular dot product of a rectangular celerite system
 *
 * @param t1     (N,): The left input coordinates (must be sorted)
 * @param t2     (M,): The right input coordinates (must be sorted)
 * @param c      (J,): The transport coefficients
 * @param U      (N, J): The first low rank matrix
 * @param V      (M, J): The second low rank matrix
 * @param Y      (M, Nrhs): The matrix that will be multiplied
 * @param Z_out  (N, Nrhs): The result of the operation
 */
template <typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void general_matmul_lower(const Eigen::MatrixBase<Input> &t1,              // (N,)
                          const Eigen::MatrixBase<Input> &t2,              // (M,)
                          const Eigen::MatrixBase<Coeffs> &c,              // (J,)
                          const Eigen::MatrixBase<LowRank> &U,             // (N, J)
                          const Eigen::MatrixBase<LowRank> &V,             // (M, J)
                          const Eigen::MatrixBase<RightHandSide> &Y,       // (M, nrhs)
                          Eigen::MatrixBase<RightHandSideOut> const &Z_out // (N, nrhs)
) {
  MakeEmptyWork(Input) F;
  general_matmul_lower<false>(t1, t2, c, U, V, Y, Z_out, F);
}

/**
 * \brief The general upper-triangular dot product of a rectangular celerite system
 *
 * @param t1     (N,): The left input coordinates (must be sorted)
 * @param t2     (M,): The right input coordinates (must be sorted)
 * @param c      (J,): The transport coefficients
 * @param U      (N, J): The first low rank matrix
 * @param V      (M, J): The second low rank matrix
 * @param Y      (M, Nrhs): The matrix that will be multiplied
 * @param Z_out  (N, Nrhs): The result of the operation
 */
template <typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void general_matmul_upper(const Eigen::MatrixBase<Input> &t1,              // (N,)
                          const Eigen::MatrixBase<Input> &t2,              // (M,)
                          const Eigen::MatrixBase<Coeffs> &c,              // (J,)
                          const Eigen::MatrixBase<LowRank> &U,             // (N, J)
                          const Eigen::MatrixBase<LowRank> &V,             // (M, J)
                          const Eigen::MatrixBase<RightHandSide> &Y,       // (M, nrhs)
                          Eigen::MatrixBase<RightHandSideOut> const &Z_out // (N, nrhs)
) {
  MakeEmptyWork(Input) F;
  general_matmul_upper<false>(t1, t2, c, U, V, Y, Z_out, F);
}

} // namespace core
} // namespace celerite2

#endif // _CELERITE2_INTERFACE_HPP_DEFINED_