text:
/*
 Copyright (C) 2016 Quaternion Risk Management Ltd
 All rights reserved.

 This file is part of ORE, a free-software/open-source library
 for transparent pricing and risk analysis - http://opensourcerisk.org

 ORE is free software: you can redistribute it and/or modify it
 under the terms of the Modified BSD License.  You should have received a
 copy of the license along with this program.
 The license is also available online at <http://opensourcerisk.org>

 This program is distributed on the basis that it will form a useful
 contribution to risk analytics and model standardisation, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 FITNESS FOR A PARTICULAR PURPOSE. See the license for more details.
*/

/*! \file ored/utilities/xmlutils.cpp
    \brief
    \ingroup utilities
*/

#include <boost/algorithm/string/join.hpp>
#include <boost/lexical_cast.hpp>
#include <ored/utilities/log.hpp>
#include <ored/utilities/parsers.hpp>
#include <ored/utilities/to_string.hpp>
#include <ored/utilities/xmlutils.hpp>

// we only want to include these here.
#include <rapidxml.hpp>
#include <rapidxml_print.hpp>

#include <algorithm>
#include <fstream>

using namespace std;
using namespace rapidxml;
using QuantLib::Size;

namespace ore {
namespace data {

namespace {

// helper routine to convert a value of an arbitrary type to string

string convertToString(const Real value) {
    // We want to write out a double that conforms to xs:double, this means no
    // scientific notation, so we check for really small numbers here and explicitly set
    // to 16 decimal places
    string result;
    if (std::abs(value) < 1.0e-6) {
        std::ostringstream obj1;
        obj1.precision(16);
        obj1 << std::fixed << value;
        result = obj1.str();
    } else {
        // And here we just use boost::lexical_cast which is better
        // at precision than std::to_string()
        result = boost::lexical_cast<std::string>(value);
    }
    return result;
}

template <class T> string convertToString(const T& value) { return ore::data::to_string(value); }

// FIXME another variant of the above routine, just for backwards compatibility, we should get
// rid of one of these two
string convertToString2(const std::string& s) { return s; }
template <class T> string convertToString2(const T& value) { return std::to_string(value); }

} // namespace

XMLDocument::XMLDocument() : _doc(new rapidxml::xml_document<char>()), _buffer(NULL) {}

XMLDocument::XMLDocument(const string& fileName) : _doc(new rapidxml::xml_document<char>()), _buffer(NULL) {
    // Need to load the entire file into memory to pass to doc.parse().
    ifstream t(fileName.c_str());
    QL_REQUIRE(t.is_open(), "Failed to open file " << fileName);
    t.seekg(0, std::ios::end);                  // go to the end
    Size length = static_cast<Size>(t.tellg()); // report location (this is the length)
    QL_REQUIRE(length > 0, "File " << fileName << " is empty.");
    t.seekg(0, std::ios::beg);      // go back to the beginning
    _buffer = new char[length + 1]; // allocate memory for a buffer of appropriate dimension
    t.read(_buffer, length);        // read the whole file into the buffer
    _buffer[static_cast<int>(t.gcount())] = '\0';
    t.close(); // close file handle
    try {
        _doc->parse<0>(_buffer);
    } catch (rapidxml::parse_error& pe) {
        string where(pe.where<char>(), 30); // limit to first 30 chars.
        QL_FAIL("RapidXML Parse Error : " << pe.what() << ". where=" << where);
    }
}

XMLDocument::~XMLDocument() {
    if (_buffer != NULL)
        delete[] _buffer;
    if (_doc != NULL)
        delete _doc;
}

void XMLDocument::fromXMLString(const string& xmlString) {
    QL_REQUIRE(!_buffer, "XML Document is already loaded");
    Size length = xmlString.size();
    _buffer = new char[length + 1];
    strcpy(_buffer, xmlString.c_str());
    _buffer[length] = '\0';
    try {
        _doc->parse<0>(_buffer);
    } catch (rapidxml::parse_error& pe) {
        string where(pe.where<char>(), 30); // limit to first 30 chars.
        QL_FAIL("RapidXML Parse Error : " << pe.what() << ". where=" << where);
    }
}

XMLNode* XMLDocument::getFirstNode(const string& name) { return _doc->first_node(name == "" ? NULL : name.c_str()); }

void XMLDocument::appendNode(XMLNode* node) { _doc->append_node(node); }

void XMLDocument::toFile(const string& fileName) {
    std::ofstream ofs(fileName.c_str());
    ofs << *_doc;
    ofs.close();
}

string XMLDocument::toString() {
    ostringstream oss;
    oss << *_doc;
    return oss.str();
}

XMLNode* XMLDocument::allocNode(const string& nodeName) {
    XMLNode* n = _doc->allocate_node(node_element, allocString(nodeName));
    QL_REQUIRE(n, "Failed to allocate XMLNode for " << nodeName);
    return n;
}

XMLNode* XMLDocument::allocNode(const string& nodeName, const string& nodeValue) {
    XMLNode* n = _doc->allocate_node(node_element, allocString(nodeName), allocString(nodeValue));
    QL_REQUIRE(n, "Failed to allocate XMLNode for " << nodeName);
    return n;
}

char* XMLDocument::allocString(const string& str) {
    char* s = _doc->allocate_string(str.c_str());
    QL_REQUIRE(s, "Failed to allocate string for " << str);
    return s;
}

void XMLSerializable::fromFile(const string& filename) {
    XMLDocument doc(filename);
    fromXML(doc.getFirstNode(""));
}

void XMLSerializable::toFile(const string& filename) {
    XMLDocument doc;
    XMLNode* node = toXML(doc);
    doc.appendNode(node);
    doc.toFile(filename);
}

void XMLSerializable::fromXMLString(const string& xml) {
    ore::data::XMLDocument doc;
    doc.fromXMLString(xml);
    fromXML(doc.getFirstNode(""));
}

string XMLSerializable::toXMLString() {
    XMLDocument doc;
    XMLNode* node = toXML(doc);
    doc.appendNode(node);
    return doc.toString();
}

void XMLUtils::checkNode(XMLNode* node, const string& expectedName) {
    QL_REQUIRE(node, "XML Node is NULL (expected " << expectedName << ")");
    QL_REQUIRE(node->name() == expectedName,
               "XML Node name " << node->name() << " does not match expected name " << expectedName);
}

XMLNode* XMLUtils::addChild(XMLDocument& doc, XMLNode* parent, const string& name) {
    QL_REQUIRE(parent, "XML Parent Node is NULL (adding Child " << name << ")");
    XMLNode* node = doc.allocNode(name);
    parent->insert_node(0, node);
    return node;
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* n, const string& name, const char* value) {
    addChild(doc, n, name, string(value));
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* n, const string& name, const string& value) {
    if (value.size() == 0) {
        addChild(doc, n, name);
    } else {
        XMLNode* node = doc.allocNode(name, value);
        QL_REQUIRE(n, "XML Node is NULL (adding " << name << ")");
        n->insert_node(0, node);
    }
}

void XMLUtils::addChildAsCdata(XMLDocument& doc, XMLNode* n, const string& name, const string& value) {
    if (value.size() == 0) {
        addChild(doc, n, name);
    } else {
        QL_REQUIRE(n, "XML Node is NULL (adding " << name << ")");
        XMLNode* node = doc.allocNode(name);
        n->insert_node(0, node);
        XMLNode* cdata_node = doc.doc()->allocate_node(node_cdata);
        cdata_node->value(doc.allocString(value));
        QL_REQUIRE(cdata_node, "Failed to allocate cdata node for " << name);
        node->insert_node(0, cdata_node);
    }
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* n, const string& name, const string& value, const string& attrName,
                        const string& attr) {
    XMLNode* node;
    if (value.size() == 0) {
        node = addChild(doc, n, name);
    } else {
        node = doc.allocNode(name, value);
        QL_REQUIRE(n, "XML Node is NULL (adding " << name << ")");
        n->insert_node(0, node);
    }
    if (attrName != "" || attr != "") {
        XMLUtils::addAttribute(doc, node, attrName, attr);
    }
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* n, const string& name, Real value) {
    addChild(doc, n, name, convertToString(value));
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* n, const string& name, int value) {
    addChild(doc, n, name, convertToString(value));
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* n, const string& name, bool value) {
    string s = value ? "true" : "false";
    addChild(doc, n, name, s);
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* n, const string& name, const Period& value) {
    addChild(doc, n, name, convertToString(value));
}

void XMLUtils::addChild(XMLDocument& doc, XMLNode* parent, const string& name, const vector<Real>& values) {
    vector<string> strings(values.size());
    std::transform(values.begin(), values.end(), strings.begin(), [](Real x) { return convertToString2(x); });
    addChild(doc, parent, name, boost::algorithm::join(strings, ","));
}

void XMLUtils::addChildren(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
                           const string& firstName, const string& secondName, const map<string, string>& values) {
    QL_REQUIRE(parent, "XML Node is null (Adding " << names << ")");
    XMLNode* node = addChild(doc, parent, names);
    map<string, string>::const_iterator it;
    for (it = values.begin(); it != values.end(); ++it) {
        XMLNode* n = addChild(doc, node, name);
        QL_REQUIRE(n, "XML AllocNode failure (" << name << ")");
        addChild(doc, n, firstName, it->first);
        addChild(doc, n, secondName, it->second);
    }
}

string XMLUtils::getChildValue(XMLNode* node, const string& name, bool mandatory) {
    QL_REQUIRE(node, "XMLNode is NULL (was looking for child " << name << ")");
    xml_node<>* child = node->first_node(name.c_str());
    if (mandatory) {
        QL_REQUIRE(child, "Error: No XML Child Node " << name << " found.");
    }
    return child ? getNodeValue(child) : "";
}

Real XMLUtils::getChildValueAsDouble(XMLNode* node, const string& name, bool mandatory, double defaultValue) {
    string s = getChildValue(node, name, mandatory);
    return s == "" ? defaultValue : parseReal(s);
}

int XMLUtils::getChildValueAsInt(XMLNode* node, const string& name, bool mandatory, int defaultValue) {
    string s = getChildValue(node, name, mandatory);
    return s == "" ? defaultValue : parseInteger(s);
}

bool XMLUtils::getChildValueAsBool(XMLNode* node, const string& name, bool mandatory, bool defaultValue) {
    string s = getChildValue(node, name, mandatory);
    return s == "" ? defaultValue : parseBool(s);
}

Period XMLUtils::getChildValueAsPeriod(XMLNode* node, const string& name, bool mandatory, const Period& defaultValue) {
    string s = getChildValue(node, name, mandatory);
    return s == "" ? defaultValue : parsePeriod(s);
}

vector<string> XMLUtils::getChildrenValues(XMLNode* parent, const string& names, const string& name, bool mandatory) {
    vector<string> vec;
    xml_node<>* node = parent->first_node(names.c_str());
    if (mandatory) {
        QL_REQUIRE(node, "Error: No XML Node " << names << " found.");
    }
    if (node) {
        for (xml_node<>* child = node->first_node(name.c_str()); child; child = child->next_sibling(name.c_str()))
            vec.push_back(getNodeValue(child));
    }
    return vec;
}

vector<Real> XMLUtils::getChildrenValuesAsDoubles(XMLNode* node, const string& names, const string& name,
                                                  bool mandatory) {
    vector<string> vecS = getChildrenValues(node, names, name, mandatory);
    vector<Real> vecD(vecS.size());
    std::transform(vecS.begin(), vecS.end(), vecD.begin(), parseReal);
    return vecD;
}

vector<Real> XMLUtils::getChildrenValuesAsDoublesCompact(XMLNode* node, const string& name, bool mandatory) {
    string s = getChildValue(node, name, mandatory);
    return parseListOfValues(s, std::function<Real(string)>(&parseReal));
}

vector<Real> XMLUtils::getNodeValueAsDoublesCompact(XMLNode* node) {
    string s = getNodeValue(node);
    return parseListOfValues(s, std::function<Real(string)>(&parseReal));
}

vector<Period> XMLUtils::getChildrenValuesAsPeriods(XMLNode* node, const string& name, bool mandatory) {
    string s = getChildValue(node, name, mandatory);
    return parseListOfValues(s, std::function<Period(string)>(&parsePeriod));
}

vector<string> XMLUtils::getChildrenValuesAsStrings(XMLNode* node, const string& name, bool mandatory) {
    string s = getChildValue(node, name, mandatory);
    return parseListOfValues(s);
}

map<string, string> XMLUtils::getChildrenValues(XMLNode* parent, const string& names, const string& name,
                                                const string& firstName, const string& secondName, bool mandatory) {
    map<string, string> res;
    xml_node<>* node = parent->first_node(names.c_str());
    if (mandatory) {
        QL_REQUIRE(node, "Error: No XML Node " << names << " found.");
    }
    if (node) {
        for (xml_node<>* child = node->first_node(name.c_str()); child; child = child->next_sibling(name.c_str())) {
            string first = getChildValue(child, firstName, mandatory);
            string second = getChildValue(child, secondName, mandatory);
            // res[first] = second;
            res.insert(pair<string, string>(first, second));
        }
    }
    return res;
}

map<string, string> XMLUtils::getChildrenAttributesAndValues(XMLNode* parent, const string& names,
                                                             const string& attributeName, bool mandatory) {
    map<string, string> res;
    for (XMLNode* child = getChildNode(parent, names.c_str()); child;
         child = XMLUtils::getNextSibling(child, names.c_str())) {
        string first = getAttribute(child, attributeName);
        string second = getNodeValue(child);
        if (mandatory) {
            QL_REQUIRE(first != "", "empty attribute for " << names);
        }
        res.insert(pair<string, string>(first, second));
    }
    if (mandatory) {
        QL_REQUIRE(res.size() > 0, "Error: No XML Node " << names << " found.");
    }
    return res;
}

// returns first child node
XMLNode* XMLUtils::getChildNode(XMLNode* n, const string& name) {
    QL_REQUIRE(n, "XMLUtils::getChildNode(" << name << "): XML Node is NULL");
    return n->first_node(name == "" ? nullptr : name.c_str());
}

// return first node in the hierarchy of n that matches name, maybe n itself
XMLNode* XMLUtils::locateNode(XMLNode* n, const string& name) {
    QL_REQUIRE(n, "XMLUtils::locateNode(" << name << "): XML Node is NULL");
    if (n->name() == name)
        return n;
    else {
        const char* p = name.size() == 0 ? nullptr : name.c_str();
        XMLNode* node = n->first_node(p);
        QL_REQUIRE(node, "XML node with name " << name << " not found");
        return node;
    }
}

// append child to parent
void XMLUtils::appendNode(XMLNode* parent, XMLNode* child) {
    QL_REQUIRE(parent, "XMLUtils::appendNode() parent is NULL");
    QL_REQUIRE(child, "XMLUtils::appendNode() child is NULL");
    parent->append_node(child);
}

void XMLUtils::addAttribute(XMLDocument& doc, XMLNode* node, const string& attrName, const string& attrValue) {
    QL_REQUIRE(node, "XMLUtils::appendAttribute(" << attrName << "," << attrName << ") node is NULL");
    char* name = doc.allocString(attrName.c_str());
    char* value = doc.allocString(attrValue.c_str());
    node->append_attribute(doc.doc()->allocate_attribute(name, value));
}

string XMLUtils::getAttribute(XMLNode* node, const string& attrName) {
    QL_REQUIRE(node, "XMLUtils::getAttribute(" << attrName << ") node is NULL");
    xml_attribute<>* attr = node->first_attribute(attrName.c_str());
    if (attr && attr->value())
        return string(attr->value());
    else
        return "";
}

vector<XMLNode*> XMLUtils::getChildrenNodes(XMLNode* node, const string& name) {
    QL_REQUIRE(node, "XMLUtils::getChildredNodes(" << name << ") node is NULL");
    vector<XMLNode*> res;
    const char* p = name.size() == 0 ? nullptr : name.c_str();
    for (xml_node<>* c = node->first_node(p); c; c = c->next_sibling(p))
        res.push_back(c);
    return res;
}

string XMLUtils::getNodeName(XMLNode* node) {
    QL_REQUIRE(node, "XMLUtils::getNodeName(): XML Node is NULL");
    return node->name();
}

void XMLUtils::setNodeName(XMLDocument& doc, XMLNode* node, const string& name) {
    QL_REQUIRE(node, "XMLUtils::setNodeName(" << name << "): XML Node is NULL");
    char* nodeName = doc.allocString(name);
    node->name(nodeName);
}

XMLNode* XMLUtils::getNextSibling(XMLNode* node, const string& name) {
    QL_REQUIRE(node, "XMLUtils::getNextSibling(" << name << "): XML Node is NULL");
    return node->next_sibling(name == "" ? nullptr : name.c_str());
}

string XMLUtils::getNodeValue(XMLNode* node) {
    QL_REQUIRE(node, "XMLUtils::getNodeValue(): XML Node is NULL");
    // handle CDATA nodes
    XMLNode* n = node->first_node();
    if (n && n->type() == node_cdata)
        return n->value();
    // all other cases
    return node->value();
}

// template implementations

// FIXME only needed as long as we want convert2String2() here, can be removed once we have consolidated
// convert2String() and convert2String()
template <>
void XMLUtils::addChildren(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
                           const vector<Real>& values) {
    vector<string> strings(values.size());
    std::transform(values.begin(), values.end(), strings.begin(), [](Real x) { return convertToString2(x); });
    addChildren(doc, parent, names, name, strings);
}

template <class T>
void XMLUtils::addChildren(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
                           const vector<T>& values) {
    XMLNode* node = addChild(doc, parent, names);
    for (Size i = 0; i < values.size(); i++)
        addChild(doc, node, name, values[i]);
}

template <class T>
void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
                                         const vector<T>& values, const string& attrName, const vector<string>& attrs) {
    addChildrenWithAttributes(doc, parent, names, name, values, vector<string>{attrName},
                              vector<vector<string>>{attrs});
}

template <class T>
void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
                                         const vector<T>& values, const vector<string>& attrNames,
                                         const vector<vector<string>>& attrs) {
    QL_REQUIRE(attrNames.size() == attrs.size(),
               "attrNames size (" << attrNames.size() << ") must match attrs size (" << attrs.size() << ")");
    if (values.size() > 0) {
        for (auto const& attr : attrs) {
            QL_REQUIRE(values.size() == attr.size(), "Values / Attribute vector size mismatch");
        }
        QL_REQUIRE(parent, "XML Node is null (Adding " << names << ")");
        XMLNode* node = addChild(doc, parent, names);
        for (Size i = 0; i < values.size(); i++) {
            XMLNode* c = doc.allocNode(name, convertToString2(values[i]));
            QL_REQUIRE(c, "XML AllocNode failure (" << name << ")");
            QL_REQUIRE(node, "XML Node is NULL (" << name << ")");
            node->insert_node(0, c);
            for (Size j = 0; j < attrs.size(); ++j) {
                if (attrs[j][i] != "")
                    addAttribute(doc, c, attrNames[j], attrs[j][i]);
            }
        }
    }
}

template <class T>
void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names, const string& name,
                                                 const vector<T>& values, const string& attrName,
                                                 const vector<string>& attrs) {
    addChildrenWithOptionalAttributes(doc, n, names, name, values, vector<string>{attrName},
                                      vector<vector<string>>{attrs});
}

template <class T>
void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names, const string& name,
                                                 const vector<T>& values, const vector<string>& attrNames,
                                                 const vector<vector<string>>& attrs) {
    QL_REQUIRE(attrNames.size() == attrs.size(),
               "attrNames size (" << attrNames.size() << ") must match attrs size (" << attrs.size() << ")");
    for (auto const& attr : attrs)
        QL_REQUIRE(attr.empty() == attrs.front().empty(), "all attributes must be empty or non-empty at the same time");
    if (attrs.empty() || attrs.front().empty())
        XMLUtils::addChildren(doc, n, names, name, values);
    else
        XMLUtils::addChildrenWithAttributes(doc, n, names, name, values, attrNames, attrs);
}

vector<string> XMLUtils::getChildrenValuesWithAttributes(XMLNode* parent, const string& names, const string& name,
                                                         const string& attrName, vector<string>& attrs,
                                                         bool mandatory) {
    return getChildrenValuesWithAttributes<string>(
        parent, names, name, attrName, attrs, [](const string& x) { return x; }, mandatory);
}

template <class T>
vector<T> XMLUtils::getChildrenValuesWithAttributes(XMLNode* parent, const string& names, const string& name,
                                                    const string& attrName, vector<string>& attrs,
                                                    const std::function<T(string)> parser, bool mandatory) {
    std::vector<std::reference_wrapper<vector<string>>> attrs_v;
    attrs_v.push_back(attrs);
    return getChildrenValuesWithAttributes(parent, names, name, vector<string>{attrName}, attrs_v, parser, mandatory);
}

vector<string> XMLUtils::getChildrenValuesWithAttributes(XMLNode* parent, const string& names, const string& name,
                                                         const vector<string>& attrNames,
                                                         const vector<std::reference_wrapper<vector<string>>>& attrs,
                                                         bool mandatory) {
    return getChildrenValuesWithAttributes<string>(
        parent, names, name, attrNames, attrs, [](const string& x) { return x; }, mandatory);
}

template <class T>
vector<T> XMLUtils::getChildrenValuesWithAttributes(XMLNode* parent, const string& names, const string& name,
                                                    const vector<string>& attrNames,
                                                    const vector<std::reference_wrapper<vector<string>>>& attrs,
                                                    const std::function<T(string)> parser, bool mandatory) {
    QL_REQUIRE(parser, "XMLUtils::getChildrenValuesWithAttributes(): parser is null");
    QL_REQUIRE(attrNames.size() == attrs.size(),
               "attrNames size (" << attrNames.size() << ") must match attrs size (" << attrs.size() << ")");
    vector<T> vec;
    rapidxml::xml_node<>* node = parent->first_node(names.c_str());
    if (mandatory) {
        QL_REQUIRE(node, "Error: No XML Node " << names << " found.");
    }
    if (node) {
        for (rapidxml::xml_node<>* child = node->first_node(name.c_str()); child;
             child = child->next_sibling(name.c_str())) {

            std::string vstr = getNodeValue(child);
            vec.push_back(parser(vstr));

            for (Size i = 0; i < attrNames.size(); ++i) {
                xml_attribute<>* attr = child->first_attribute(attrNames[i].c_str());
                if (attr && attr->value())
                    ((vector<string>&)attrs[i]).push_back(attr->value());
                else
                    ((vector<string>&)attrs[i]).push_back("");
            }
        }
    }
    return vec;
}

/* Instantiate some template functions above for types T we want to support. Add instantiations for more types here,
   if required. This has two advantages over putting the templated version of the functions in the header file:
   - faster compile times
   - we do not need to include the rapid xml headers in xmlutils.hpp */

template void XMLUtils::addChildren(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
                                    const vector<string>& values);
// throws a warning currently, but that is ok, see the FIXME above in template <> void XMLUtils::addChildren(...)
// template void XMLUtils::addChildren(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
//                                     const vector<double>& values);
template void XMLUtils::addChildren(XMLDocument& doc, XMLNode* parent, const string& names, const string& name,
                                    const vector<bool>& values);

template void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names,
                                                  const string& name, const vector<string>& values,
                                                  const string& attrName, const vector<string>& attrs);
template void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names,
                                                  const string& name, const vector<double>& values,
                                                  const string& attrName, const vector<string>& attrs);
template void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names,
                                                  const string& name, const vector<bool>& values,
                                                  const string& attrName, const vector<string>& attrs);

template void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names,
                                                  const string& name, const vector<string>& values,
                                                  const vector<string>& attrNames, const vector<vector<string>>& attrs);
template void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names,
                                                  const string& name, const vector<double>& values,
                                                  const vector<string>& attrNames, const vector<vector<string>>& attrs);
template void XMLUtils::addChildrenWithAttributes(XMLDocument& doc, XMLNode* parent, const string& names,
                                                  const string& name, const vector<bool>& values,
                                                  const vector<string>& attrNames, const vector<vector<string>>& attrs);

template void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names,
                                                          const string& name, const vector<string>& values,
                                                          const string& attrName, const vector<string>& attrs);
template void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names,
                                                          const string& name, const vector<double>& values,
                                                          const string& attrName, const vector<string>& attrs);
template void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names,
                                                          const string& name, const vector<bool>& values,
                                                          const string& attrName, const vector<string>& attrs);

template void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names,
                                                          const string& name, const vector<string>& values,
                                                          const vector<string>& attrNames,
                                                          const vector<vector<string>>& attrs);
template void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names,
                                                          const string& name, const vector<double>& values,
                                                          const vector<string>& attrNames,
                                                          const vector<vector<string>>& attrs);
template void XMLUtils::addChildrenWithOptionalAttributes(XMLDocument& doc, XMLNode* n, const string& names,
                                                          const string& name, const vector<bool>& values,
                                                          const vector<string>& attrNames,
                                                          const vector<vector<string>>& attrs);

template vector<std::string> XMLUtils::getChildrenValuesWithAttributes(XMLNode* parent, const string& names,
                                                                       const string& name, const string& attrName,
                                                                       vector<string>& attrs,
                                                                       std::function<string(string)> parser,
                                                                       bool mandatory);

template vector<double> XMLUtils::getChildrenValuesWithAttributes(XMLNode* parent, const string& names,
                                                                  const string& name, const string& attrName,
                                                                  vector<string>& attrs,
                                                                  std::function<double(string)> parser, bool mandatory);
template vector<bool> XMLUtils::getChildrenValuesWithAttributes(XMLNode* parent, const string& names,
                                                                const string& name, const string& attrName,
                                                                vector<string>& attrs,
                                                                std::function<bool(string)> parser, bool mandatory);

} // namespace data
} // namespace ore

max_stars_repo_name: mrslezak/Engine

max_issues_repo_name: mrslezak/Engine

max_forks_repo_name: mrslezak/Engine


text:
/***************************************************************************
 *   Software License Agreement (BSD License)                              *
 *   Copyright (C) 2012 by Markus Bader <markus.bader@tuwien.ac.at>        *
 *                                                                         *
 *   Redistribution and use in source and binary forms, with or without    *
 *   modification, are permitted provided that the following conditions    *
 *   are met:                                                              *
 *                                                                         *
 *   1. Redistributions of source code must retain the above copyright     *
 *      notice, this list of conditions and the following disclaimer.      *
 *   2. Redistributions in binary form must reproduce the above copyright  *
 *      notice, this list of conditions and the following disclaimer in    *
 *      the documentation and/or other materials provided with the         *
 *      distribution.                                                      *
 *   3. Neither the name of the copyright holder nor the names of its      *
 *      contributors may be used to endorse or promote products derived    *
 *      from this software without specific prior written permission.      *
 *                                                                         *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS   *
 *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT     *
 *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS     *
 *   FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE        *
 *   COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,  *
 *   INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,  *
 *   BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;      *
 *   LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER      *
 *   CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT    *
 *   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY *
 *   WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE           *
 *   POSSIBILITY OF SUCH DAMAGE.                                           *
 ***************************************************************************/
#include <iostream>
#include <stdlib.h>
#include <signal.h>

#include "shmfw/variable.h"
#include <boost/program_options.hpp>

struct Prarmeters {
    std::string shm_memory_name;
    unsigned int shm_memory_size;
    std::string variable_name;
    bool loop;
    Prarmeters()
    : loop(true) {
    }
};
Prarmeters params;

Prarmeters readArgs ( int argc, char **argv ) {
    namespace po = boost::program_options;
    po::options_description desc ( "Allowed Parameters" );
    desc.add_options()
    ( "help", "get this help message" )
    ( "variable_name,n", po::value<std::string> ( &params.variable_name ), "shared variable name" )
    ( "shm_memory_name,m", po::value<std::string> ( &params.shm_memory_name )->default_value ( ShmFw::DEFAULT_SEGMENT_NAME() ), "shared memory segment name" )
    ( "shm_memory_size", po::value<unsigned int> ( &params.shm_memory_size )->default_value ( ShmFw::DEFAULT_SEGMENT_SIZE() ), "shared memory segment size" );

    po::variables_map vm;
    try {
        po::store ( po::parse_command_line ( argc, argv, desc ), vm );
    } catch ( const std::exception &ex ) {
        std::cout << desc << std::endl;;
        exit ( 1 );
    }
    po::notify ( vm );

    if ( vm.count ( "help" ) )  {
        std::cout << desc << std::endl;
        exit ( 1 );
    }
    return params;
}

void terminate ( int s ) {
    printf ( "Caught signal %d\n",s );
    fflush ( stdout );
    params.loop = false;
}

int main ( int argc, char *argv[] ) {

    readArgs ( argc, argv );
    signal ( SIGABRT,	terminate );
    signal ( SIGTERM,	terminate );

    ShmFw::HandlerPtr shmHdl = ShmFw::Handler::create ( params.shm_memory_name, params.shm_memory_size );
    if ( shmHdl->findName ( params.variable_name ) == NULL ) {
        std::cerr << "No variable: " << params.variable_name;
        std::cerr << " in share segment " << params.shm_memory_name;
        std::cerr << "!" << std::endl;
        return 0;
    }
    ShmFw::Var<ShmFw::Header> header( params.variable_name, shmHdl );
    boost::posix_time::ptime t0 = boost::posix_time::microsec_clock::local_time();
    unsigned int count = 0;
    while ( params.loop ) {           // Don't echo() while we do getch
	t0 = t0 + boost::posix_time::seconds(1);
	std::vector<boost::posix_time::time_duration> cycles;
	std::vector<boost::posix_time::time_duration> cycles_times;
	while(boost::posix_time::microsec_clock::local_time() < t0){
	  if(header.timed_wait(5000)){
	    cycles.push_back(boost::posix_time::microsec_clock::local_time() - header->timestampLocal());
	    header->dataProcessed();
	  } else {
	    std::cerr << "No signal within 5sec: " << std::endl;
	  }
	}
	if(cycles.size() > 0){
	  double sum = 0;
	  for(size_t i = 0; i < cycles.size(); i++){
	    sum += cycles[i].total_microseconds() / ( double ) boost::posix_time::seconds ( 1 ).total_microseconds();
	  }
	  double cycleAVR = sum / (double) cycles.size();

	  printf ( "%4i: %6.2fHz = %8.2fms :  %s", count, 1.0/cycleAVR, cycleAVR*1000, params.variable_name.c_str());
	  count++;
	}
    };
    exit ( 0 );

}

ext: cpp

max_stars_repo_name: ShmFw/shmfw

max_issues_repo_name: ShmFw/shmfw

max_forks_repo_name: ShmFw/shmfw


text:
/**
 * C++ implementation of Full Embeded Element beam model for MWA based on beam_full_EE.py script
 * and Sokolowski et al (2016) paper
 * Implemented by Marcin Sokolowski (May 2017) - marcin.sokolowski@curtin.edu.au
 * 20 May 2017 : Somewhat optimized by André Offringa.
 */

#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>

#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/filesystem.hpp>

#include <H5Cpp.h>

#include "beam2016implementation.h"
#include "system.h"

using namespace std;
using namespace H5;

#define nullptr NULL

// constants :
static const double deg2rad = M_PI/180.00;

// defines for debuging / logging - to be easily turned off on compilation level :
#define _PRINTF_LEVEL0 if(m_VerbLevel>0)printf
#define _PRINTF_LEVEL1 if(m_VerbLevel>1)printf
#define _PRINTF_LEVEL2 if(m_VerbLevel>2)printf

// timing function
static time_t get_dttm()
{
   long gm_time;
   time( &gm_time );
   return gm_time;
}

/*
  Print jones matrix in a format :

  j00 j01
  j10 j11
*/
void JonesMatrix::Print(const char* name,double az_deg,double za_deg)
{
/*   if( m_VerbLevel>3 ){
      printf("%s at (az,za)=(%.2f,%.2f) [deg] = \n",name,az_deg,za_deg);
      printf("\t%.8f + %.8fj     |     %.8f + %.8fj\n",j00.real(),j00.imag(),j01.real(),j01.imag());
      printf("\t%.8f + %.8fj     |     %.8f + %.8fj\n",j10.real(),j10.imag(),j11.real(),j11.imag());
   }*/
}

int Beam2016Implementation::m_VerbLevel=-1; // >=0 to enable talking
const double Beam2016Implementation::m_DefaultDelays[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; // default delays at zenith
const double Beam2016Implementation::m_DefaultAmps[]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};   // default amplitudes 1 for all the tile dipoles
const double Beam2016Implementation::_delayStep=435.0e-12; // beamformer step in pico-seconds

int Beam2016Implementation::has_freq(int freq_hz)
{
   for(size_t i=0;i<m_freq_list.size();i++){
      if( m_freq_list[i] == freq_hz ){
         return 1;
      }
   }

   return 0;
}

int Beam2016Implementation::find_closest_freq(int freq_hz)
{
   double min_diff=1e20;
   int best_idx=-1;

   for(size_t i=0;i<m_freq_list.size();i++){
      double diff = abs( m_freq_list[i] - freq_hz );
      if( diff < min_diff ){
         min_diff = diff;
         best_idx = i;
      }
   }

   if( best_idx >=0 ){
      return m_freq_list[best_idx];
   }

   return m_freq_list[0];
}

// local power calculations
// WARNING : please do not remove - the program works 2x slower when std::pow is used instead !!!
complex<double> Beam2016Implementation::power_complex( complex<double> val, int n )
{
   complex<double> ret=1;
   for(int i=0;i<n;i++){
      ret = ret*val;
   }

   return ret;
}

Beam2016Implementation::Beam2016Implementation( const double* delays, const double* amps ) :
  m_CalcModesLastFreqHz(-1),
  m_CalcModesLastDelays(nullptr),
  m_CalcModesLastAmps(nullptr),
  m_NormJones(1,1,1,1),
  m_NormFreqHz(-1),
  m_AntennaCount(N_ANT_COUNT),
  m_pH5File(nullptr)
{
	if( !delays ) {
		delays = m_DefaultDelays;
	}
	if( !amps ) {
		amps = m_DefaultAmps;
	}
	for(size_t i=0; i!=16; ++i)
	{
		_delays[i] = delays[i];
		_amps[i] = amps[i];
	}
	Read();
}

Beam2016Implementation::~Beam2016Implementation()
{
   delete m_pH5File;
   delete [] m_CalcModesLastDelays;
   delete [] m_CalcModesLastAmps;
}

//-------------------------------------------------------------------- Calculation of Jones matrix ------------------------------------------------------------------------------------------------------------------
void Beam2016Implementation::CalcJonesArray( vector< vector<double> >& azim_arr, vector< vector<double> >& za_arr, vector< vector<JonesMatrix> >& jones, int freq_hz_param, const double* delays, const double* amps, bool bZenithNorm )
{
  // convert AZIM -> FEKO PHI phi=90-azim or azim=90-phi :
  // python : phi_arr=math.pi/2-phi_arr #Convert to East through North (FEKO coords)
  for(size_t y=0;y<azim_arr.size();y++){
     vector<double>& image_row = azim_arr[y];

     for(size_t x=0;x<image_row.size();x++){
        image_row[x] = M_PI/2.00 - image_row[x];

        // phi_arr[phi_arr < 0] += 2*math.pi #360 wrap
        if( image_row[x] < 0 ){
           image_row[x] += 2*M_PI;
        }
     }
  }
  JonesMatrix::zeros( jones, (int)(azim_arr[0].size()), (int)(azim_arr.size()) );

  for(size_t y=0;y<azim_arr.size();y++){
     vector<double>& image_row = azim_arr[y];
      for(size_t x=0;x<image_row.size();x++){
          double azim_deg = image_row[x];
          double za_deg   = (za_arr[y])[x];

          (jones[y])[x] = CalcJones( azim_deg, za_deg, freq_hz_param, bZenithNorm );
      }
  }
}

JonesMatrix Beam2016Implementation::CalcJonesDirect( double az_rad, double za_rad )
{
  // convert AZIM -> FEKO PHI phi=90-azim or azim=90-phi :
  // python : phi_arr=math.pi/2-phi_arr #Convert to East through North (FEKO coords)
  JonesMatrix jones;
  double phi_rad = M_PI/2.00 - az_rad;

  CalcSigmas( phi_rad, za_rad, Q1_accum_X, Q2_accum_X, M_accum_X, N_accum_X, MabsM_X, Nmax_X, Cmn_X, 'X', jones );
  CalcSigmas( phi_rad, za_rad, Q1_accum_Y, Q2_accum_Y, M_accum_Y, N_accum_Y, MabsM_Y, Nmax_Y, Cmn_Y, 'Y', jones );

  // Swap test :
  JonesMatrix jones_tmp;
  jones_tmp = jones;

  jones.j00 = jones_tmp.j11;
  jones.j11 = jones_tmp.j00;
  jones.j01 = jones_tmp.j10;
  jones.j10 = jones_tmp.j01;

  return jones;
}

JonesMatrix Beam2016Implementation::CalcJones( double az_deg, double za_deg, int freq_hz_param, bool bZenithNorm )
{
	return CalcJones(az_deg, za_deg, freq_hz_param, _delays, _amps, bZenithNorm);
}

JonesMatrix Beam2016Implementation::CalcJones( double az_deg, double za_deg, int freq_hz, const double* delays, const double* amps, bool bZenithNorm )
{
   if( has_freq(freq_hz)<=0 ){
      freq_hz = find_closest_freq( freq_hz );
   }

   CalcModes( freq_hz , m_AntennaCount, delays, amps );
   JonesMatrix jones_ret = CalcJonesDirect( az_deg*deg2rad, za_deg*deg2rad );

   if( bZenithNorm ) {
      if( freq_hz != m_NormFreqHz ) {
         CalcZenithNormMatrix( freq_hz );
      }

      jones_ret.j00 = jones_ret.j00 / m_NormJones.j00;
      jones_ret.j01 = jones_ret.j01 / m_NormJones.j01;
      jones_ret.j10 = jones_ret.j10 / m_NormJones.j10;
      jones_ret.j11 = jones_ret.j11 / m_NormJones.j11;
   }

   jones_ret.Print("Jones",az_deg,za_deg);

   return jones_ret;
}

void Beam2016Implementation::CalcZenithNormMatrix( int freq_hz )
{
   if( freq_hz != m_NormFreqHz ){
      printf("INFO : calculating Jones matrix for frequency = %d Hz\n",freq_hz);

      // Azimuth angles at which Jones components are maximum (see beam_full_EE.py for comments):
      //  max_phis=[[math.pi/2,math.pi],[0,math.pi/2]] #phi where each Jones vector is max
      double j00_max_az= 90.00;
      double j01_max_az=180.00;
      double j10_max_az=  0.00;
      double j11_max_az= 90.00;

      JonesMatrix tmp_jones;

      // j00 :
      tmp_jones = CalcJones( j00_max_az, 0, freq_hz, m_DefaultDelays, m_DefaultAmps, false );
      m_NormJones.j00 = abs( tmp_jones.j00 );
      // m_NormJones.j00 =  tmp_jones.j00;

      // j01 :
      tmp_jones = CalcJones( j01_max_az, 0, freq_hz, m_DefaultDelays, m_DefaultAmps, false );
      m_NormJones.j01 = abs( tmp_jones.j01 );
      // m_NormJones.j01 = tmp_jones.j01;

      // j10 :
      tmp_jones = CalcJones( j10_max_az, 0, freq_hz, m_DefaultDelays, m_DefaultAmps, false );
      m_NormJones.j10 = abs( tmp_jones.j10 );
      // m_NormJones.j10 = tmp_jones.j10;

      // j11 :
      tmp_jones = CalcJones( j11_max_az, 0, freq_hz, m_DefaultDelays, m_DefaultAmps, false );
      m_NormJones.j11 = abs( tmp_jones.j11 );
      // m_NormJones.j11 = tmp_jones.j11;

      m_NormJones.Print("Zenith normalisation matrix",0,0);

      m_NormFreqHz = freq_hz;
   }else{
      printf("INFO : normalisation matrix already calculated for frequency = %d Hz\n",freq_hz);
   }
}

//-------------------------------------------------------------------- Calculation of Jones matrix components for given polarisation (eq. 4 and 5 in the Sokolowski et al (2017) paper --------------------------------
void Beam2016Implementation::CalcSigmas( double phi, double theta,
                                    vector< complex<double> >& Q1_accum, vector< complex<double> >& Q2_accum,
                                    vector<double>& M_accum, vector<double>& N_accum, vector<double>& MabsM, double Nmax, vector<double>& Cmn,
                                    char pol,
                                    JonesMatrix& jones_matrix )
{
   // int nmax = int( max(N_accum) );
   int nmax =  int( Nmax );

   //complex<double> complex_j(0,1);

   //double sin_theta = sin(theta);
   double cos_theta = cos(theta);
   double u = cos_theta;

   vector<double> P1sin_arr, P1_arr;

   P1sin( nmax, theta, P1sin_arr, P1_arr );
   _PRINTF_LEVEL0("DEBUG sizes P1_sin_arr.size = %d , P1_arr.size = %d, N.size = %d, M.size = %d, Q2.size = %d, Q1.size = %d\n",(int)(P1sin_arr.size()),(int)(P1_arr.size()),(int)(N_accum.size()),(int)(M_accum.size()),(int)(Q2_accum.size()),(int)(Q1_accum.size()));

   if( N_accum.size() != M_accum.size() ){
      printf("ERROR : size of N_accnum != M_accum ( %d != %d)\n",(int)(N_accum.size()),(int)(M_accum.size()));
      BOOST_ASSERT_MSG( (N_accum.size()==M_accum.size()),"Sizes of N_accum and M_accum must be equal");
   }
   _PRINTF_LEVEL0("M x N = %d x %d\n",(int)(M_accum.size()),(int)(N_accum.size()));

   complex<double> sigma_P(0,0),sigma_T(0,0);
   for(size_t i=0;i<N_accum.size();i++){
      double N = N_accum[i];
      int n    = int(N);

      double M = M_accum[i];
      //int    m = int(M);
      double m_abs_m = MabsM[i];

      double c_mn_sqr = (0.5*(2*N+1)*factorial_wrapper(N-abs(M))/factorial_wrapper(N+abs(M)));
      double c_mn = sqrt( c_mn_sqr );
      // Optimisation comment :
      // Possibly this might be a faster version - but does not seem so, so I live it as it was
      // MS tested it (2017-05-17), but does not seem to give a significant speed up, so for now the old version left
      // If tested again comment the 2 lines above and uncomment the line below ,
      // also verify that lines in CalcModes ( after comment "Intialisation of Cmn vector :") are un-commented - FOR NOW THEY ARE COMMENTED OUT NOT TO CALCULATE SOMETHING WHICH IS NOT USED
      // I've tested it on 2017-05-19 and the version with line below and lines at the end of CalcModes uncommented runs in ~9m30sec and the current one was ~9m40sec so I leave the current
      // version as I've tested it for longer time (MS), but it can be restored if some further optimisation is needed (but it will not be a breaktrough).
      // double c_mn = Cmn[i];

      complex<double> ejm_phi( cos(M*phi), sin(M*phi) );
      complex<double> phi_comp = ( ejm_phi*c_mn ) / ( sqrt(N*(N+1)) ) * m_abs_m;
      // printf("phi_comp[%d] = %e + %ej (c_mn=%.8f)\n",i,phi_comp.real(),phi_comp.imag(),c_mn);
      _PRINTF_LEVEL0("DEBUG : %.2f %.4f %.8f -> phi_comp = %e + %ej\n",M,phi,c_mn,phi_comp.real(),phi_comp.imag());
      _PRINTF_LEVEL0("DEBUG : Q1[%d] = %e + %ej , Q2[%d] = %e + %ej\n",int(i),Q1_accum[i].real(),Q1_accum[i].imag(),int(i),Q2_accum[i].real(),Q2_accum[i].imag());

      if( (n+1) >= int(j_power_cache.size()) ){
         char szError[1024];
         sprintf(szError,"ERROR : j_power_cache.size() = %d too few values as power %d requested\n",int(j_power_cache.size()),(n+1));
         BOOST_ASSERT_MSG( false , szError );
      }

//      complex<double> j_power_n = power_complex(complex_j,n); // WARNING : please do not remove - the program works 2x slower when std::pow(complex_j,n) is used instead !!!
      complex<double> j_power_n = j_power_cache[n];
      complex<double> E_theta_mn = j_power_n * ( P1sin_arr[i] * ( fabs(M) * Q2_accum[i]*u - M*Q1_accum[i] ) + Q2_accum[i]*P1_arr[i] );
      _PRINTF_LEVEL0("E_theta_mn[%d] = %e + %ej = (%e + %ej) * ( %e * ( %e * (%e + %ej) * %e - (%e)*(%e + %ej)) + (%e + %ej) * %e)\n",int(i),E_theta_mn.real(),E_theta_mn.imag(),j_power_n.real(),j_power_n.imag(),P1sin_arr[i],fabs(M),Q2_accum[i].real(),Q2_accum[i].imag(),u,M,Q1_accum[i].real(),Q1_accum[i].imag(),Q2_accum[i].real(),Q2_accum[i].imag(),P1_arr[i]);

//      complex<double> j_power_np1 = power_complex(complex_j,n+1); // WARNING : please do not remove - the program works 2x slower when std::pow(complex_j,n+1) is used instead !!!
      complex<double> j_power_np1 = j_power_cache[n+1];
      complex<double> E_phi_mn = j_power_np1 * ( P1sin_arr[i] * ( M*Q2_accum[i] - fabs(M)*Q1_accum[i]*u) - Q1_accum[i]*P1_arr[i] );
      _PRINTF_LEVEL0("E_phi_mn[%d] = %e + %ej = (%e + %ej) * ( %e * ( %e (%e + %ej) - %e (%e + %ej) %e ) + (%e + %ej) %e )\n",int(i),E_phi_mn.real(),E_phi_mn.imag(),j_power_np1.real(),j_power_np1.imag(),P1sin_arr[i],M,Q2_accum[i].real(),Q2_accum[i].imag(),fabs(M),Q2_accum[i].real(),Q2_accum[i].imag(),u,Q1_accum[i].real(),Q1_accum[i].imag(),P1_arr[i]);

       sigma_P = sigma_P + phi_comp * E_phi_mn;
       sigma_T = sigma_T + phi_comp * E_theta_mn;

       _PRINTF_LEVEL0("Sigma_P[%d] := (%.10f + %.10fj) \n",int(i),sigma_P.real(),sigma_P.imag());
       _PRINTF_LEVEL0("Sigma_T[%d] := (%.10f + %.10fj) \n",int(i),sigma_T.real(),sigma_T.imag());
   }

// NORMAL VERSION
   if( pol == 'X' ){
       jones_matrix.j00 = sigma_T;
       jones_matrix.j01 = -sigma_P; // as it is now in python code in sign_fix branch
   }else{
       jones_matrix.j10 = sigma_T;
       jones_matrix.j11 = -sigma_P; // as it is now in python code in sign_fix branch
   }

// swap test for wsclean :
/*   if( pol == 'X' ){
       jones_matrix.j10 = sigma_T;
       jones_matrix.j11 = -sigma_P; // as it is now in python code in sign_fix branch
   }else{
       jones_matrix.j00 = sigma_T;
       jones_matrix.j01 = -sigma_P; // as it is now in python code in sign_fix branch
   }   */

   if( m_VerbLevel > 0 ){
      jones_matrix.Print("Jones");
   }
}

//-------------------------------------------------------------------- Calculation of spherical harmonics coefficients (eq. 3-6 in the Sokolowski et al 2016 paper)  --------------------------------
// function comparing current parameters : frequency, delays and amplitudes with those previously used to calculate spherical waves coefficients (stored in the 3 variables :
// m_CalcModesLastFreqHz , , m_CalcModesLastAmps )
int Beam2016Implementation::IsCalcModesRequired( int freq_hz, int n_ant, const double* delays, const double* amps )
{
   if( freq_hz != m_CalcModesLastFreqHz ){
      return 1;
   }

   if( !m_CalcModesLastDelays || !m_CalcModesLastAmps ){
      return 1;
   }

   for(int i=0;i<n_ant;i++){
      if( delays[i] != m_CalcModesLastDelays[i] ){
         return 1;
      }
      if( amps[i] != m_CalcModesLastAmps[i] ){
         return 1;
      }
   }

   // printf(" Beam2016Implementation::IsCalcModesRequired = 0\n");
   return 0;
}

// function calculating coefficients for X and Y and storing parameters frequency, delays and amplitudes
void Beam2016Implementation::CalcModes( int freq_hz, size_t n_ant, const double* delays, const double* amps )
{
   if( IsCalcModesRequired(  freq_hz, n_ant, delays, amps )<=0  ){
      // if recalculation of modes is not required -> do not do it
      return;
   }

   Nmax_X = CalcModes( freq_hz , n_ant, delays, amps, 'X', Q1_accum_X, Q2_accum_X, M_accum_X, N_accum_X, MabsM_X, Cmn_X );
   // TODO : perhaps can be uncommented above and removed below :
   /*Nmax_X = max(N_accum_X);
   if( Nmax_X != Nmax_X_test ){
      printf("ERROR : Nmax_X_test=%.2f != Nmax_X=%.2f\n",Nmax_X_test,Nmax_X);
   }
   Nmax_X = Nmax_X_test;*/

   Nmax_Y = CalcModes( freq_hz , n_ant, delays, amps, 'Y', Q1_accum_Y, Q2_accum_Y, M_accum_Y, N_accum_Y, MabsM_Y, Cmn_Y );
   // TODO : perhaps can be uncommented above and removed below :
   /*Nmax_Y = max(N_accum_Y);
   if( Nmax_Y != Nmax_Y_test ){
      printf("ERROR : Nmax_X_test=%.2f != Nmax_X=%.2f\n",Nmax_Y_test,Nmax_Y);
   }
   Nmax_Y = Nmax_Y_test;*/

   // initialise powers of j
   double maxN = std::max(Nmax_X,Nmax_Y);
   complex<double> complex_j(0,1);
   j_power_cache.assign( int(maxN) + 20 , 1 );
   for(int i=1;i<int(j_power_cache.size());i++){
      j_power_cache[i] = j_power_cache[i-1] * complex_j;
   }

   m_CalcModesLastFreqHz = freq_hz;
   if( !m_CalcModesLastDelays ){
      m_CalcModesLastDelays = new double[n_ant];
   }
   memcpy(m_CalcModesLastDelays,delays,sizeof(double)*n_ant);

   if( !m_CalcModesLastAmps ){
     m_CalcModesLastAmps = new double[n_ant];
   }
   memcpy(m_CalcModesLastAmps,amps,sizeof(double)*n_ant);
}

// function calculating all coefficients Q1, Q2, N, M and derived MabsM, Nmax for a given polarisation (X or Y) - perhaps enum should be used here
double Beam2016Implementation::CalcModes( int freq_hz, size_t n_ant, const double* delays, const double* amp, char pol,
                                   vector< complex<double> >& Q1_accum, vector< complex<double> >& Q2_accum,
                                   vector<double>& M_accum, vector<double>& N_accum,
                                   vector<double>& MabsM, vector<double>& Cmn   )
{
   vector<double> phases(n_ant);
   double Nmax=0;
   M_accum.clear();
   N_accum.clear();
   MabsM.clear();
   Cmn.clear();

   int modes_size = m_Modes[0].size();
   _PRINTF_LEVEL0("Size(Modes) = %d\n",modes_size);
   Q1_accum.assign(modes_size, 0.0);
   Q2_accum.assign(modes_size, 0.0);

   for(size_t a=0;a<n_ant;a++){
      double phase = 2*M_PI*freq_hz*(-double(delays[a])*_delayStep);

      phases[a] = phase;

      // complex excitation voltage:
      // self.amps[pol]*np.exp(1.0j*phases)
      complex<double> phase_factor( cos(phase) , sin(phase) );

      complex<double> Vcplx = amp[a]*phase_factor;

      _PRINTF_LEVEL0("ANT_%d : Phase = %e , Vcplx = %e + %ej\n",int(a),phase,Vcplx.real(),Vcplx.imag());

      char Q_all_name[64];
      sprintf(Q_all_name,"%c%d_%d",pol,int(a)+1,freq_hz);

      // NO CACHE VERSION :
      vector< vector<double> > Q_all;
      ReadDataSet( Q_all_name, Q_all );

      _PRINTF_LEVEL0("Read dataaset %s has dimensions %d x %d\n",Q_all_name,(int)(Q_all.size()),(int)(Q_all[0].size()));

      size_t n_ant_coeff = Q_all[0].size();
      _PRINTF_LEVEL0("Number of coefficients for antenna = %d is %d\n",int(a),int(n_ant_coeff));

      vector<double> Ms1,Ns1,Ms2,Ns2;
      vector<double>&  m_Modes_Type = m_Modes[0];
      vector<double>&  m_Modes_M = m_Modes[1];
      vector<double>&  m_Modes_N = m_Modes[2];

      int bUpdateNAccum=0;
      vector<int> s1_list; // list of indexes where S=1 coefficients seat in array m_Modes ( and m_Modes_M and m_Modes_N )
      vector<int> s2_list; // list of indexes where S=2 coefficients seat in array m_Modes ( and m_Modes_M and m_Modes_N )
      for(size_t coeff=0;coeff<n_ant_coeff;coeff++){
         int mode_type = m_Modes_Type[coeff];

         if( mode_type <= 1 ){
            // s=1 modes :
            s1_list.push_back(coeff);

            Ms1.push_back( m_Modes_M[coeff] );
            Ns1.push_back( m_Modes_N[coeff] );
            if( m_Modes_N[coeff] > Nmax ){
               Nmax = m_Modes_N[coeff];
               bUpdateNAccum=1;
            }
         }else{
           // s=2 modes :
           s2_list.push_back(coeff);

           Ms2.push_back( m_Modes_M[coeff] );
           Ns2.push_back( m_Modes_N[coeff] );
         }
      }

      if( bUpdateNAccum > 0 ){
         N_accum = Ns1;
         M_accum = Ms1;
      }

      if( s1_list.size() == s2_list.size() && s2_list.size()==(n_ant_coeff/2) ){
         _PRINTF_LEVEL0("DEBUG : correct number of coefficients %d == %d == %d satisfied\n",(int)(s1_list.size()),(int)(s2_list.size()),int(n_ant_coeff/2));
      }else{
         char szError[1024];
         sprintf(szError,"ERROR : wrong number of coefficients for s1 and s2 condition %d = =%d == %d not satisfied!!!\n",(int)(s1_list.size()),(int)(s2_list.size()),int(n_ant_coeff/2));
         BOOST_ASSERT_MSG( s1_list.size() == s2_list.size() && s2_list.size()==(n_ant_coeff/2) , szError );
      }

      vector< std::complex<double> > Q1,Q2;
      vector<double>& Q_all_0 = Q_all[0];
      vector<double>& Q_all_1 = Q_all[1];
      int my_len_half=(n_ant_coeff/2);

      for(int i=0;i<my_len_half;i++){
         // calculate Q1:
         int s1_idx = s1_list[i];
         double s10_coeff = Q_all_0[s1_idx];
         double s11_coeff = Q_all_1[s1_idx];
         double arg = s11_coeff*deg2rad;
         complex<double> tmp( cos(arg) , sin(arg) );
         complex<double> q1_val = s10_coeff * tmp;
         Q1.push_back(q1_val);

         // calculate Q2:
         int s2_idx = s2_list[i];
         double s20_coeff = Q_all_0[s2_idx];
         double s21_coeff = Q_all_1[s2_idx];
         double arg2 = s21_coeff*deg2rad;
         complex<double> tmp2( cos(arg2) , sin(arg2) );
         complex<double> q2_val = s20_coeff * tmp2;
         Q2.push_back(q2_val);

         complex<double> debug_val = Q1_accum[i];

         Q1_accum[i] = Q1_accum[i] + q1_val*Vcplx;
         Q2_accum[i] = Q2_accum[i] + q2_val*Vcplx;

         if( m_VerbLevel>1 ){
            if( i==0 || i==1 ){
               printf("\ti=%d : Q1_accum := (%.14f + j%.14f) = (%.14f + j%.14f) + (%.14f + j%.14f)*(%.14f + j%.14f)\n",i,Q1_accum[i].real(),Q1_accum[i].imag(),debug_val.real(),debug_val.imag(),q1_val.real(),q1_val.imag(),Vcplx.real(),Vcplx.imag());
            }
         }
      }
      _PRINTF_LEVEL0("DEBUG Q2.size=%d vs. filled %d elements\n",(int)(Q2_accum.size()),my_len_half);
      print( Q1, Q1_accum, Q2, Q2_accum, my_len_half, N_accum, M_accum, Nmax );
   }

   _PRINTF_LEVEL0("DEBUG CalcModes(%c) : Q1[0] = %e + %ej , Q2[0] = %e + %ej\n",pol,Q1_accum[0].real(),Q1_accum[0].imag(),Q2_accum[0].real(),Q2_accum[0].imag());
   _PRINTF_LEVEL0("DEBUG TEST : Q1[0] = %e + %ej , Q2[0] = %e + %ej\n",Q1_accum[0].real(),Q1_accum[0].imag(),Q2_accum[0].real(),Q2_accum[0].imag());

   // Moved from CalcSigmas here :
   // Same as tragic python code:
   // MabsM=-M/np.abs(M)
   // MabsM[MabsM==np.NaN]=1 #for M=0, replace NaN with MabsM=1;
   // MabsM=(MabsM)**M
   for(size_t i=0;i<M_accum.size();i++){
      int m = int(M_accum[i]);

      double m_abs_m = 1;
      if( m>0 ){
         if( (m%2) != 0 ){
            m_abs_m = -1;
         }
      }

      MabsM.push_back( m_abs_m );
   }

   // Intialisation of Cmn vector :
   /*BOOST_ASSERT_MSG( (N_accum.size()==M_accum.size()),"Sizes of N_accum and M_accum must be equal when initialising Cmn vector");
   Cmn.assign( N_accum.size(), 0.00 );
   for(size_t i=0;i<N_accum.size();i++){
      double N = N_accum[i];
      double M = M_accum[i];

      if( fabs(M) <= fabs(N) ){  // only calculate really used coefficients, where m<=n - see eq. 3 in the Sokolowski et al (2016) paper :
         double c_mn_sqr = (0.5*(2*N+1)*factorial_wrapper(N-abs(M))/factorial_wrapper(N+abs(M)));
         double c_mn = sqrt( c_mn_sqr );

         Cmn[i] = c_mn;
      }else{
         char szError[1024];
         sprintf(szError,"Wrong value of m=%.2f or n=%.2f coefficient fabs(m) must be <= fabs(n)",M,N);
         printf("%s\n",szError);fflush(stdout);
         BOOST_ASSERT_MSG( false, szError);
      }
   }*/

   return Nmax;
}

//------------------------------------------------------------------------------------------------------ maths functions and wrappers ---------------------------------------------------------------------------------------
// factorial calculation :
double Beam2016Implementation::factorial_wrapper_base( unsigned n )
{
   return boost::math::factorial<double>(n);
}

// uses precalculated factorials stored in m_Factorial
double Beam2016Implementation::factorial_wrapper( unsigned n )
{
   if( m_Factorial.size() == 0 ){
      cache_factorial( 100 ); // pre-calculate maximum possible factorial used in the calculations, MS checked that 100 should be enough and pratical
   }
   if( n < 0 ){
      printf("WARNING : negative factorial required ??? n=%d\n",n);
   }

   if( n >= m_Factorial.size() ){
      printf("WARNING : required factorial %d! exceeds initially pre-calculated 100! -> re-calculating now (code optimisation required)\n",n);
      cache_factorial( n );
   }

   if( n < m_Factorial.size() ){
      return m_Factorial[n];
   }

   return factorial_wrapper_base( n );
}

void Beam2016Implementation::cache_factorial( unsigned max_n )
{
   time_t start_time = get_dttm();

   m_Factorial.clear();
   for(unsigned i=0;i<=max_n;i++){
      double fact = factorial_wrapper_base( i );
      m_Factorial.push_back( fact );
   }
   _PRINTF_LEVEL0("Initialisation of factorial %d! took %d seconds\n",max_n,(int)(get_dttm()-start_time));
}

// OUTPUT : returns list of Legendre polynomial values calculated up to order nmax :
int Beam2016Implementation::P1sin( int nmax, double theta, vector<double>& p1sin_out, vector<double>& p1_out )
{
	int size = nmax*nmax + 2*nmax;
	p1sin_out.resize(size);
	p1_out.resize(size);

	double sin_th, u;
	sincos(theta, &sin_th, &u);
	double delu=1e-6;

	vector<double> P, Pm1, Pm_sin, Pu_mdelu, Pm_sin_merged, Pm1_merged;
	P.reserve(nmax+1);
	Pm1.reserve(nmax+1);
	Pm_sin.reserve(nmax+1);
	Pu_mdelu.reserve(nmax+1);
	Pm_sin_merged.reserve(nmax*2+1);
	Pm1_merged.reserve(nmax*2+1);

	// Create a look-up table for the legendre polynomials
	// Such that legendre_table[ m * nmax + (n-1) ] = legendre(n, m, u)
	vector<double> legendre_table(nmax * (nmax+1));
	vector<double>::iterator legendre_iter = legendre_table.begin();
	for(int m=0; m!=nmax+1; ++m) {
		double leg0 = boost::math::legendre_p(0, m, u);
		double leg1 = boost::math::legendre_p(1, m, u);
		*legendre_iter = leg1;
			++legendre_iter;
		for(int n=2; n!=nmax+1; ++n) {
			if(m < n)
				*legendre_iter = boost::math::legendre_next(n-1, m, u, leg1, leg0);
			else if(m == n)
				*legendre_iter = boost::math::legendre_p(n, m, u);
			else
				*legendre_iter = 0.0;
			//std::cout << n << ',' << m << ' ' << leg1 << "~" << leg0 << "->" << *legendre_iter << '\n';
			//std::cout << boost::math::legendre_p(n, m, u) << '\n';
			leg0 = leg1;
			leg1 = *legendre_iter;
			++legendre_iter;
		}
	}

	for(int n=1; n<=nmax; n++) {
		P.resize(n+1);
		// Assign P[m] to legendre(n, m, u)
		// This is equal to:
		// lpmv(P, n , u);
		for(size_t m=0; m!=size_t(n)+1; ++m)
			P[m] = legendre_table[m*nmax + (n-1)];

		// skip first 1 and build table Pm1 (P minus 1 )
		Pm1.assign(P.begin()+1, P.end());
		Pm1.push_back(0);

		Pm_sin.assign(n+1, 0.0);
		if( u==1 || u==-1 ) {
			// In this case we take the easy approach and don't use
			// precalculated polynomials, since this path does not occur often.

			// TODO This path doesn't make sense.
			// I think Marcin should look at this; this path only occurs on polar positions so
			// is rare, but what is done here looks wrong: first calculate *all* polynomials,
			// then only use the pol for m=1. Pm_sin for indices 0 and >=2 are not calculated, this
			// seems not right.
			Pu_mdelu.resize(1);
			lpmv(Pu_mdelu, n, u-delu);

			// Pm_sin[1,0]=-(P[0]-Pu_mdelu[0])/delu #backward difference
			if( u == -1 )
				Pm_sin[1] = -(Pu_mdelu[0]-P[0])/delu; // #forward difference
			else
				Pm_sin[1] = -(P[0]-Pu_mdelu[0])/delu;
		}
		else {
			for(size_t i=0; i<P.size(); i++){
				Pm_sin[i] = P[i] / sin_th;
			}
		}

		Pm_sin_merged.assign(Pm_sin.rbegin(), Pm_sin.rend()-1);
		Pm_sin_merged.insert(Pm_sin_merged.end(), Pm_sin.begin(), Pm_sin.end());

		int ind_start=(n-1)*(n-1)+2*(n-1); // #start index to populate
		int ind_stop=n*n+2*n; //#stop index to populate

		// P_sin[np.arange(ind_start,ind_stop)]=np.append(np.flipud(Pm_sin[1::,0]),Pm_sin)
		int modified=0;
		for(int i=ind_start; i<ind_stop; i++){
			p1sin_out[i] = (Pm_sin_merged[modified]);
			modified++;
		}

		// P1[np.arange(ind_start,ind_stop)]=np.append(np.flipud(Pm1[1::,0]),Pm1)
		Pm1_merged.assign(Pm1.rbegin(), Pm1.rend()-1);
		Pm1_merged.insert(Pm1_merged.end(), Pm1.begin(), Pm1.end());
		modified=0;
		for(int i=ind_start;i<ind_stop;i++){
			p1_out[i] = Pm1_merged[modified];
			modified++;
		}
	}

	return nmax;
}

// Legendre polynomials :
void Beam2016Implementation::lpmv( vector<double>& output, int n, double x )
{
   for(size_t order=0; order<output.size(); order++) {
      double val = boost::math::legendre_p<double>(n, order, x);
      output[order] = val;
   }
}

//----------------------------------------------------------------------------------- HDF5 File interface and data structures for H5 data -----------------------------------------------------------------------------------
// This function goes thorugh all dataset names and records them info list of strings : Beam2016Implementation::m_obj_list
herr_t Beam2016Implementation::list_obj_iterate(hid_t loc_id, const char *name, const H5O_info_t *info, void *operator_data)
{
    string szTmp;
    Beam2016Implementation* pBeamModelPtr = ( Beam2016Implementation*)operator_data;
    BOOST_ASSERT_MSG( pBeamModelPtr!=NULL, "The pointer to  Beam2016Implementation class in Beam2016Implementation::list_obj_iterate must not be NULL");

    if (name[0] == '.'){         /* Root group, do not print '.' */
        _PRINTF_LEVEL2("  (Group)\n");
    }else
        switch (info->type) {
            case H5O_TYPE_GROUP:
                _PRINTF_LEVEL2("%s  (Group)\n", name);
                break;
            case H5O_TYPE_DATASET:
                _PRINTF_LEVEL2("%s  (Dataset)\n", name);
                szTmp = name;
                pBeamModelPtr->m_obj_list.push_back(szTmp);
                break;
            case H5O_TYPE_NAMED_DATATYPE:
                _PRINTF_LEVEL2("%s  (Datatype)\n", name);
                break;
            default:
                _PRINTF_LEVEL2("%s  (Unknown)\n", name);
        }

    return 0;
}

// Read dataset_name from H5 file
void Beam2016Implementation::ReadDataSet( const char* dataset_name, vector< vector<double> >& out_vector )
{
   _PRINTF_LEVEL0("\n\n\n\n\n\n---------------------------------------------------------------------------------------- Dataset = %s ----------------------------------------------------------------------------------------\n",dataset_name);

   DataSet modes = m_pH5File->openDataSet( dataset_name );
   hsize_t size=modes.getStorageSize();
   H5T_class_t type_class =modes.getTypeClass();
   size_t data_size=0;
   string type_class_str="unknown";
   H5std_string order_string;
   IntType intype;
   FloatType floattype;
   //H5T_order_t order;
   _PRINTF_LEVEL0("type_class = %d vs INT=%d, FLOAT=%d\n",type_class,H5T_INTEGER,H5T_FLOAT);
   out_vector.clear();

   switch( type_class ){
         case H5T_INTEGER :
            intype = modes.getIntType();

            /*
             * Get order of datatype and print message if it's a little endian.
            */
            //order = intype.getOrder( order_string );
            data_size = intype.getSize();
            type_class_str = "integer";
            break;

         case H5T_FLOAT :
            floattype = modes.getFloatType();
            data_size = floattype.getSize();
            type_class_str = "float";
            break;

         default :
            printf("ERROR : unknown type class = %d\n",type_class);
            break;
   }

   _PRINTF_LEVEL0("Storage size for modes = %d elements of class = %d (%s) of size %d bytes\n",(int)size,(int)type_class,type_class_str.c_str(),(int)data_size);

   DataSpace modes_dataspace = modes.getSpace();
   int rank = modes_dataspace.getSimpleExtentNdims();
   hsize_t dims_out[2];
   modes_dataspace.getSimpleExtentDims( dims_out, NULL);
   _PRINTF_LEVEL0("rank of modes = %d : %d x %d\n",(int)(rank),(int)(dims_out[0]),(int)(dims_out[1]));
   modes_dataspace.selectAll();

   float* data = new float[dims_out[0]*dims_out[1]];
   float** modes_data = new float*[dims_out[0]];
   for(size_t i=0;i<dims_out[0];i++){
      modes_data[i] = data + i*dims_out[1];
   }
   DataSpace memspace( rank, dims_out );
   modes.read( data, PredType::NATIVE_FLOAT, memspace, modes_dataspace );

   for(size_t i=0;i<dims_out[0];i++){
      vector<double> empty_vector;
      for(size_t j=0;j<dims_out[1];j++){
         empty_vector.push_back(modes_data[i][j]);

         _PRINTF_LEVEL0("%e ",modes_data[i][j]);
      }
      out_vector.push_back( empty_vector );

      _PRINTF_LEVEL0("\n-------------------------------\n");
   }

   if( m_VerbLevel > 0 ){
      printf("\n\n----- TEST -----\n\n");fflush(stdout);
      for(size_t i=0;i<out_vector.size();i++){
         vector<double>& line = out_vector[i];
         for(size_t j=0;j<line.size();j++){
            printf("%e ",modes_data[i][j]);
         }
         printf("\n");
         printf("-------------------------------\n");
      }
   }

   delete[] modes_data;
   delete[] data;
}

// Interface to HDF5 file format and structures to store H5 data
// Functions for reading H5 file and its datasets :
// Read data from H5 file, file name is specified in the object constructor
int Beam2016Implementation::Read()
{
   if( !m_pH5File ) {
      //if(!boost::filesystem::exists(_h5filename)) {
         string h5_test_path = DEFAULT_H5_FILE_PATH;
         h5_test_path += DEFAULT_H5_FILE;

         _PRINTF_LEVEL0("WARNING : file %s does not exist -> trying python path %s\n",_h5filename.c_str(),h5_test_path.c_str());
        std::string h5_path = System::FindPythonFilePath(h5_test_path);

         if(!boost::filesystem::exists(h5_path)){
            char szError[1024];

            sprintf(szError,"ERROR - cannot find the H5 file tested paths %s and %s -> cannot continue",_h5filename.c_str(),h5_path.c_str());
            printf("%s\n",szError);
            BOOST_ASSERT_MSG( false, szError );
         } else {
            // found H5 file and using it :
            printf("Found H5 file %s on %s -> using this path\n",_h5filename.c_str(),h5_path.c_str());
            _h5filename = h5_path.c_str();
         }
      //}

      m_pH5File = new H5File( _h5filename.c_str(), H5F_ACC_RDONLY );
   }else{
      _PRINTF_LEVEL0(" Beam2016Implementation::Read : file %s already read -> skipped\n",_h5filename.c_str());
      return 1;
   }

   if( m_pH5File ){
      hid_t group_id = m_pH5File->getId();
      _PRINTF_LEVEL0("Group ID = %d\n",(int)group_id);

      hid_t file_id = m_pH5File->getId();
      _PRINTF_LEVEL0("File ID = %d\n",(int)file_id);

      /* TODO :  not sure how to read attribute with the official HDF5 library ...
      if( H5Aexists( file_id, "VERSION" ) ){
         char szVersion[128];
         strcpy(szVersion,"TEST");
         //H5std_string szVersion;

         hid_t attr_id = H5Aopen(file_id, "VERSION", H5P_DEFAULT );
         hid_t attr_type = H5Aget_type(attr_id);
         herr_t err = H5Aread( attr_id, attr_type, (void*)szVersion );
         printf("Ids = %d -> %d -> type = %d\n",file_id,attr_id,attr_type);
         printf("Version of the %s file = %s (err = %d)\n",m_h5file.c_str(),szVersion,err);
      }else{
         printf("ERROR : attribute version does not exist\n");
      }*/

      m_obj_list.clear();
      m_freq_list.clear();
      herr_t status =  H5Ovisit (file_id, H5_INDEX_NAME, H5_ITER_NATIVE, list_obj_iterate, this);
      if( status < 0 ){
         BOOST_ASSERT_MSG( status>=0 , "H5Ovisit returned with negative value which indicates a critical error" );
      }

      int max_ant_idx=-1;
      for(size_t i=0;i<m_obj_list.size();i++){
          const char* key = m_obj_list[i].c_str();
          if(strstr(key,"X1_")){
              const char* szFreq = key+3;
              m_freq_list.push_back(atol(szFreq));
          }

          if( key[0] == 'X' ){
             int ant_idx=0,freq_hz=0;
             int scanf_ret = sscanf(key,"X%d_%d",&ant_idx,&freq_hz);
             if( scanf_ret == 2 ){
                if( ant_idx > max_ant_idx ){
                   max_ant_idx = ant_idx;
                }
             }
          }
       }
       m_AntennaCount = max_ant_idx; // number of antenna is read from the file
       if( m_AntennaCount != N_ANT_COUNT ){
          char szError[1024];
          sprintf(szError,"ERROR : number of simulated antennae = %d, the code is currently implemented for %d\n",m_AntennaCount,N_ANT_COUNT);
          BOOST_ASSERT_MSG( m_AntennaCount != N_ANT_COUNT, szError );
       }
       std::sort( m_freq_list.begin(), m_freq_list.end() );

       // show frequencies in the H5 file :
       if( m_VerbLevel >= 0 ){
          printf("Maximum antenna index in file = %d\n",max_ant_idx);
          printf("FREQUENCIES :\n");
          for(size_t i=0;i<m_freq_list.size();i++){
             printf("\t%d Hz\n",m_freq_list[i]);
          }
       }
       ReadDataSet("modes",m_Modes);
   }

   return 1;
}

//----------------------------------------------------------------------------------- debuging / logging functions ----------------------------------------------------------------------------------------------------------------
void Beam2016Implementation::print( vector< vector<double> >& arr, const char* name, int force )
{
   if( m_VerbLevel > 0 || force>0 ){
      printf("%s : \n",name);
      for(size_t i=0;i<arr.size();i++){
         vector<double>& row = arr[i];
         printf("\t");
         for(size_t j=0;j<row.size();j++){
            printf("%.2f ",row[j]);
         }
         printf("\n");
      }
      printf("\n");
   }
}

void Beam2016Implementation::print( vector<double>& arr, const char* name, int force )
{
   if( m_VerbLevel > 0 || force>0 ){
      printf("%s : ",name);
      for(size_t i=0;i<arr.size();i++){
         printf("%e ",arr[i]);
      }
      printf("\n");
   }
}

void Beam2016Implementation::print( vector<int>& arr, const char* name, int force )
{
   if( m_VerbLevel > 0 || force>0 ){
      printf("%s : ",name);
      for(size_t i=0;i<arr.size();i++){
         printf("%d ",arr[i]);
      }
      printf("\n");
   }
}

void Beam2016Implementation::print( vector< complex<double> >& Q1, vector< complex<double> >& Q1_accum,
                                    vector< complex<double> >& Q2, vector< complex<double> >& Q2_accum,
                                    int my_len_half,
                                    vector<double>& N_accum, vector<double>& M_accum,
                                    double Nmax )
{
      if( m_VerbLevel > 1 ){
         printf("Q1 = \n");
         for(int i=0;i<my_len_half;i++){
            complex<double>& q1_val = Q1[i];
            complex<double>& q1_accum_val = Q1_accum[i];

            printf("\tQ1[%d] = %e + j%e -> Q1_accum[%d] = %e + j%e\n",i,q1_val.real(),q1_val.imag(),i,q1_accum_val.real(),q1_accum_val.imag());
         }
         printf("\n\nQ2 = \n");
         for(int i=0;i<my_len_half;i++){
            complex<double>& q2_val = Q2[i];
            complex<double>& q2_accum_val = Q2_accum[i];

            printf("\tQ2[%d] = %e + j%e -> Q2_accum[%d] = %e + j%e\n",i,q2_val.real(),q2_val.imag(),i,q2_accum_val.real(),q2_accum_val.imag());
         }

         printf("\n\n");
         printf("N_max = %.2f\n",Nmax);
         printf("N_accum.size = %d , M_accum.size = %d\n",(int)(N_accum.size()),(int)(M_accum.size()));
         printf("N_accum =\n");
         for(size_t i=0;i<N_accum.size();i++){
            printf("\t");
            printf("%.0f ",N_accum[i]);
            if( ((i+1)%11)==0  ){
               printf("\n");
            }
         }
         printf("\n\n\n");

         printf("M_accum =\n");
         for(size_t i=0;i<M_accum.size();i++){
            printf("\t");
            printf("%.0f ",M_accum[i]);
            if( ((i+1)%11)==0  ){
               printf("\n");
            }
         }
         printf("\n");

      }
}

//----------------------------------------------------------------------------------- auxiliary functions for basic vector operations - TO BE REPLACED BY std calls ---------------------------------------------------------------

void JonesMatrix::zeros( vector< vector<JonesMatrix> >& jones, int x_size, int y_size )
{
   vector<JonesMatrix> zero_vector(x_size, JonesMatrix(0,0,0,0) );
   jones.assign(y_size, zero_vector);
}


text:
#include "adjusterwidget.hpp"

#include <stdexcept>
#include <string>

#include <boost/filesystem.hpp>

#include <QHBoxLayout>
#include <QLabel>
#include <QStyle>

CSVDoc::AdjusterWidget::AdjusterWidget (QWidget *parent)
    : QWidget (parent), mValid (false), mAction (ContentAction_Undefined)
{
    QHBoxLayout *layout = new QHBoxLayout (this);

    mIcon = new QLabel (this);

    layout->addWidget (mIcon, 0);

    mMessage = new QLabel (this);
    mMessage->setWordWrap (true);
    mMessage->setSizePolicy (QSizePolicy (QSizePolicy::Minimum, QSizePolicy::Minimum));

    layout->addWidget (mMessage, 1);

    setName ("", false);

    setLayout (layout);
}

void CSVDoc::AdjusterWidget::setAction (ContentAction action)
{
    mAction = action;
}

void CSVDoc::AdjusterWidget::setLocalData (const boost::filesystem::path& localData)
{
    mLocalData = localData;
}

boost::filesystem::path CSVDoc::AdjusterWidget::getPath() const
{
    if (!mValid)
        throw std::logic_error ("invalid content file path");

    return mResultPath;
}

bool CSVDoc::AdjusterWidget::isValid() const
{
    return mValid;
}

void CSVDoc::AdjusterWidget::setFilenameCheck (bool doCheck)
{
    mDoFilenameCheck = doCheck;
}

void CSVDoc::AdjusterWidget::setName (const QString& name, bool addon)
{
    QString message;

    mValid = (!name.isEmpty());

    if (!mValid)
    {
        message = "No name.";
    }
    else
    {
        boost::filesystem::path path (name.toUtf8().data());

        bool isLegacyPath = (path.extension() == ".esm" ||
                             path.extension() == ".esp");

        bool isFilePathChanged = (path.parent_path().string() != mLocalData.string());

        if (isLegacyPath)
            path.replace_extension (addon ? ".omwaddon" : ".omwgame");

        //if the file came from data-local and is not a legacy file to be converted,
        //don't worry about doing a file check.
        if (!isFilePathChanged && !isLegacyPath)
        {
            // path already points to the local data directory
            message = QString::fromUtf8 (("Will be saved as: " + path.string()).c_str());
            mResultPath = path;
        }
        //in all other cases, ensure the path points to data-local and do an existing file check
        else
        {
            // path points somewhere else or is a leaf name.
            if (isFilePathChanged)
                path = mLocalData / path.filename();

            message = QString::fromUtf8 (("Will be saved as: " + path.string()).c_str());
            mResultPath = path;

            if (boost::filesystem::exists (path))
            {
                /// \todo add an user setting to make this an error.
                message += "<p>A file with the same name already exists. If you continue, it will be overwritten.";
            }
        }
    }

    mMessage->setText (message);
    mIcon->setPixmap (style()->standardIcon (
        mValid ? QStyle::SP_MessageBoxInformation : QStyle::SP_MessageBoxWarning).
        pixmap (QSize (16, 16)));

    emit stateChanged (mValid);
}

max_stars_repo_name: Bodillium/openmw

max_issues_repo_name: Bodillium/openmw

max_forks_repo_name: Bodillium/openmw


text: #include <boost/asio/posix/stream_descriptor_service.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include "gl_core_4_3.h"
#include "SDL.h"
#include "SDL_opengl.h"
#include "SDL_timer.h"
#include <Eigen/Dense>

#include "loadshaders_4_3.h"

#include <cmath>
#include <iostream>
#include <stdint.h>

namespace {

constexpr GLubyte *bufferOffset(size_t bytes)
{
    return static_cast<GLubyte *>(0) + bytes;
}

const int kWindowWidth = 800;
const int kWindowHeight = 600;

enum VAO_IDs { Triangles, NumVAOs };
enum Buffer_IDs { ArrayBuffer, NumBuffers };
enum Attrib_IDs { vPosition = 0, mTransform = 1 };

GLuint VAOs[NumVAOs];
GLuint Buffers[NumBuffers];

float rotationAngle = 0;
uint32_t oldTime = 0;
uint32_t currentTime = 0;

const GLuint NumVertices = 6;

void init()
{
    glGenVertexArrays(NumVAOs, VAOs);
    glBindVertexArray(VAOs[Triangles]);

    Eigen::Vector2f vertices[NumVertices] = {
        { -0.90, -0.90 },
        {  0.85, -0.90 },
        { -0.90,  0.85 },
        {  0.90, -0.85 },
        {  0.90,  0.90 },
        { -0.85,  0.90 }
    };

    glGenBuffers(NumBuffers, Buffers);
    glBindBuffer(GL_ARRAY_BUFFER, Buffers[ArrayBuffer]);
    glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);

    ShaderInfo shaders[] = {
        { GL_VERTEX_SHADER, "../tests/hello/hello_4_3.vert" },
        { GL_FRAGMENT_SHADER, "../tests/hello/hello_4_3.frag" },
        { GL_NONE, NULL }
    };

    GLuint program = loadShaders(shaders);
    glUseProgram(program);

    glVertexAttribPointer(vPosition, 2, GL_FLOAT, GL_FALSE, 0, bufferOffset(0));
    glEnableVertexAttribArray(vPosition);
}

static const float PI = 2*3.1415926;
static const float TWO_PI = 2*PI;

// http://www.opengl.org/sdk/docs/man/xhtml/glBindFragDataLocation.xml
// http://www.opengl.org/sdk/docs/man/xhtml/glBindAttribLocation.xml

void display()
{
    glClear(GL_COLOR_BUFFER_BIT);

    oldTime = currentTime;
    currentTime = SDL_GetTicks();
    float deltaTime = (currentTime - oldTime) / 1000.0f;

    rotationAngle = fmod(rotationAngle + (deltaTime * (PI / 60.0f)), TWO_PI);

    Eigen::Matrix2f rotation;
    rotation << cos(rotationAngle), -sin(rotationAngle),
                sin(rotationAngle),  cos(rotationAngle);
    glUniformMatrix2fv(mTransform, 1, GL_FALSE, rotation.data());

    glBindVertexArray(VAOs[Triangles]);
    glDrawArrays(GL_TRIANGLES, 0, NumVertices);
}

} // namespace

int main(int argc, char *argv[])
{
    if (SDL_Init(SDL_INIT_VIDEO) != 0) {
        std::cerr << SDL_GetError() << std::endl;
        return -1;
    }

    SDL_GL_SetAttribute(SDL_GL_CONTEXT_MAJOR_VERSION, 4);
    SDL_GL_SetAttribute(SDL_GL_CONTEXT_MINOR_VERSION, 3);
    SDL_GL_SetAttribute(SDL_GL_CONTEXT_PROFILE_MASK, SDL_GL_CONTEXT_PROFILE_CORE);

    // Create an OpenGL capable window
    SDL_Window *window = SDL_CreateWindow(
        "Hello World",
        SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED,
        800, 600,
        SDL_WINDOW_OPENGL | SDL_WINDOW_RESIZABLE
    );
    if (!window) {
        std::cerr << SDL_GetError() << std::endl;
        SDL_Quit();
        return -1;
    }

    // Create an OpenGL context associated with the window.
    SDL_GLContext glcontext = SDL_GL_CreateContext(window);
    if (!glcontext) {
        std::cerr << SDL_GetError() << std::endl;
        SDL_Quit();
        return -1;
    }

    // Load extensions
    if (ogl_LoadFunctions() != ogl_LOAD_SUCCEEDED) {
        SDL_DestroyWindow(window);
        SDL_Quit();
        return -1;
    }

    // Initialize displayables
    init();

    SDL_Event event;
    bool keepRunning = true;
    while (keepRunning) {
        while (SDL_PollEvent(&event)) {
            if (event.type == SDL_QUIT) {
                keepRunning = false;
            }
        }

        // Render scene
        display();

        // Swap buffers
        SDL_GL_SwapWindow(window);
    }

    // Once finished with OpenGL functions, the SDL_GLContext can be deleted.
    SDL_GL_DeleteContext(glcontext);

    // Done! Close the window, clean-up and exit the program.
    SDL_DestroyWindow(window);
    SDL_Quit();

    return 0;
}

max_stars_repo_name: asgeir/old-school-projects

max_issues_repo_name: asgeir/old-school-projects

max_forks_repo_name: asgeir/old-school-projects


text:
/*
Copyright(c) 2014 Akihiro Nishimura

This software is released under the MIT License.
http://opensource.org/licenses/mit-license.php
*/

#define _SCL_SECURE_NO_WARNINGS

#include "mrlda.h"
#include "../helper/mapreduce_module.h"
#include "SigUtil/lib/calculation/ublas.hpp"
#include "SigUtil/lib/calculation/for_each.hpp"
#include "SigUtil/lib/calculation/assign_operation.hpp"
#include "SigUtil/lib/distance/norm.hpp"
#include "SigUtil/lib/tools/convergence.hpp"
#include "SigUtil/lib/file.hpp"
#include <boost/numeric/ublas/io.hpp>

namespace sigtm
{
const uint MrLDA:: max_local_iteration_ = 100;
const double MrLDA::global_convergence_threshold = 0.0001;
const double alpha_convergence_threshold = 0.001;
const uint max_alpha_update_iteration = 100;

const auto resume_info_fname = SIG_TO_FPSTR("mrlda_info");
const auto resume_alpha_fname = SIG_TO_FPSTR("mrlda_alpha");
const auto resume_gamma_fname = SIG_TO_FPSTR("mrlda_gamma");
const auto resume_phi_fname = SIG_TO_FPSTR("mrlda_phi");

void MrLDA::MapTask::process(value_type const& value, VectorK<double>& gamma, MatrixVK<double>& omega) const
{
	uint const vnum = value.vnum_;
	uint const knum = value.knum_;
	VectorK<double> const& alpha = *value.alpha_;
	MatrixKV<double> const& phi = *value.phi_;
	VectorV<uint> const& word_ct = *value.word_ct_;
	//VectorK<double> prev_gamma(value.knum_, 0);
	VectorK<double> sigma(value.knum_, 0);

	//while (!is_convergence(gamma, prev_gamma)){
	for (uint i = 0; i < max_local_iteration_; ++i){
		const auto exp_digamma = sig::map([](double g){ return std::exp(digamma(g)); }, gamma);

		for (uint v = 0; v<vnum; ++v){
			for (uint k = 0; k<knum; ++k){
				// update omega
				omega[v][k] = phi[k][v] * exp_digamma[k];
			}
			sig::normalize_dist(omega[v]);
			const uint wn = word_ct[v];
			sig::compound_assignment([wn](double& s, double p){ s += wn * p; }, sigma, omega[v]);
		}
		//prev_gamma = std::move(gamma);
		gamma = sig::plus(alpha, sigma);
	}
}

void MrLDA::init(bool resume)
{
	doc_word_ct_ = SIG_INIT_MATRIX(uint, D, V, 0);
	for (auto const& t : input_data_->tokens_){
		++doc_word_ct_[t.doc_id][t.word_id];
	}

	term_score_ = SIG_INIT_MATRIX(double, K, V, 0);

	auto base_pass = sig::modify_dirpass_tail(input_data_->getWorkingDirectory(), true);

	if (resume){
		auto load_info = sig::load_line(base_pass + resume_info_fname);
		if (sig::isJust(load_info)){
			auto info = sig::fromJust(load_info);
			total_iter_ct_ = std::stoul(info[0]);
		}

		auto tmp_alpha = std::move(alpha_);
		auto load_alpha = sig::load_num<double, VectorK<double>>(base_pass + resume_alpha_fname, " ");
		if (sig::isJust(load_alpha)){
			alpha_ = std::move(sig::fromJust(load_alpha));
			std::cout << "resume alpha" << std::endl;
		}
		else{
			alpha_ = std::move(tmp_alpha);
			std::cout << "resume alpha error : alpha is set as default" << std::endl;
		}
	}

	auto load_gamma = resume ? sig::load_num2d<double, MatrixDK<double>>(base_pass + resume_gamma_fname, " ") : nothing;
	if (sig::isJust(load_gamma)){
		gamma_ = std::move(sig::fromJust(load_gamma));
		std::cout << "resume gamma" << std::endl;
	}
	else{
		gamma_ = MatrixDK<double>(D_, VectorK<double>(K_));
		for (DocumentId d = 0; d<D_; ++d){
			for (TopicId k=0; k<K_; ++k){
				gamma_[d][k] = alpha_[k] + rand_d_();
			}
			bool ck = sig::normalize_dist(gamma_[d]);
		}
		if (resume) std::cout << "resume gamma error : gamma is set by random" << std::endl;
	}

	auto load_beta = resume ? sig::load_num2d<double, MatrixKV<double>>(base_pass + resume_phi_fname, " ") : nothing;
	if (sig::isJust(load_beta)){
		phi_ = std::move(sig::fromJust(load_beta));
		std::cout << "resume phi" << std::endl;
	}
	else{
		phi_ = MatrixKV<double>(K_, VectorV<double>(V_));
		for (TopicId k = 0; k<K_; ++k){
			for (WordId v = 0; v<V_; ++v){
				phi_[k][v] = eta_[k][v] + rand_d_();
			}
			bool ck = sig::normalize_dist(phi_[k]);
		}
		if (resume) std::cout << "resume phi error : phi is set by random" << std::endl;
	}
}

void MrLDA::saveResumeData() const
{
	std::cout << "save resume data... ";

	auto base_pass = input_data_->getWorkingDirectory();

	sig::save_num(alpha_, base_pass + resume_alpha_fname, " ");
	sig::save_num(gamma_, base_pass + resume_gamma_fname, " ");
	sig::save_num(phi_, base_pass + resume_phi_fname, " ");

	sig::clear_file(base_pass + resume_info_fname);
	sig::save_line(total_iter_ct_, base_pass + resume_info_fname, sig::WriteMode::append);

	std::cout << "completed" << std::endl;
}

double MrLDA::calcLiklihood(double term2, double term4) const
{
	return D_ * calcModule0(alpha_) + term2 + term3_ + term4;
}

void updateAlpha(VectorK<double>& alpha, VectorK<double> const& sufficient_statistics, const int D)
{
	const uint K = alpha.size();

	const auto calc_gradient = [D](double old_alpha, double suff_stat, double old_alpha_sum_digamma)
	{
		std::cout << "grad:" << D*old_alpha_sum_digamma << "," << -D*digamma(old_alpha) << "," << suff_stat << std::endl;
		return D * (old_alpha_sum_digamma - digamma(old_alpha)) + suff_stat;
	};

	const auto calc_hessian = [D](uint row, uint col, VectorK<double> const& old_alpha, double old_alpha_sum_trigamma)
	{
		return row == col ? D * (trigamma(old_alpha[col]) - old_alpha_sum_trigamma) : D * (-old_alpha_sum_trigamma);
	};

	auto regulate = [&](VectorK<double> const& delta, uint& decay) ->VectorK<double>
	{
		while (true){
			bool f = true;
			auto corr_delta = sig::mult(std::pow(0.8, decay), delta);

			std::cout << "delta:";
			for (uint i=0; i<K; ++i){
				std::cout << corr_delta[i] << ", ";
				if (alpha[i] < std::abs(corr_delta[i])){
					f = false;
					++decay;
					break;
				}
			}
			std::cout << std::endl;
			if (f) return corr_delta;
		}
	};

	uint iteration_count = 0;
	uint decay = 0;
	sig::ManageConvergence<VectorK<double>> convergence(alpha_convergence_threshold, sig::norm_L2);

	while (true) {
		const double alpha_sum = sig::sum(alpha);
		if (++iteration_count > max_alpha_update_iteration || convergence.update(alpha)) break;

		sig::vector_u<double> gradient(K);
		sig::matrix_u<double> hessian(K, K);
		auto alpha_sum_digamma = digamma(alpha_sum);
		auto alpha_sum_trigamma = trigamma(alpha_sum);

		for (uint i = 0; i < K; ++i){
			gradient(i) = calc_gradient(alpha[i], sufficient_statistics[i], alpha_sum_digamma);
			for (uint j = 0; j < K; ++j){
				hessian(i, j) = calc_hessian(i, j, alpha, alpha_sum_trigamma);
			}
		}

		auto delta = sig::from_vector_ublas(boost::numeric::ublas::prod(sig::fromJust(sig::invert_matrix(hessian)), gradient));
		auto corr_delta = regulate(delta, decay);

		sig::compound_assignment(
			[](double& a, double d){ a -= d; },
			alpha, corr_delta
		);

		std::cout << "alpha:";
		for (auto e : alpha) std::cout << e << ", ";
		std::cout << std::endl;
	}
	std::cout << "ct:" << iteration_count << std::endl;
}

/*
void updateAlpha(VectorK<double>& alpha, VectorK<double> const& sufficient_statistics, const int D)
{
	const uint K = alpha.size();
	uint iteration_count = 0;
	sig::ManageConvergence<VectorK<double>, sig::RelativeError, sig::Norm<1>> convergence(alpha_convergence_threshold, sig::norm_L1);

	VectorK<double> hessian(K);
	VectorK<double> gradient(K);
	VectorK<double> alpha_new(K);

	int decay = 0;
	while (true) {
		double sumG_H = 0;
		double sum1_H = 0;

		const double alpha_sum = sig::sum(alpha);
		if (++iteration_count > max_alpha_update_iteration || convergence.update(alpha)) break;

		for (int i = 0; i < K; i++) {
			gradient[i] =K * (digamma(alpha_sum) - digamma(alpha[i])) + sufficient_statistics[i];

			hessian[i] = -D * trigamma(alpha[i]);

			sumG_H += gradient[i] / hessian[i];
			sum1_H += 1 / hessian[i];
		}

		double z = D * trigamma(alpha_sum);
		double c = sumG_H / (1 / z + sum1_H);

		while (true) {
			bool singular = false;

			std::cout << "stepsize:";
			for (int i = 0; i < K; i++) {
				double step_size = std::pow(0.8, decay) * (gradient[i] - c) / hessian[i];
				std::cout << step_size << ", ";
				if (alpha[i] < std::abs(step_size)){
					// the current hessian matrix is singular
					singular = true;
					break;
				}
				alpha_new[i] = alpha[i] - step_size;
			}
			std::cout << std::endl;

			if (singular) {
				decay++;
				alpha_new = alpha;
			}
			else {
				break;
			}
		}
		alpha = alpha_new;
		std::cout << "alpha:" << alpha[0] << "," << alpha[1] << std::endl;
	}
	std::cout << "ct:" << iteration_count << std::endl;
}
*/

double MrLDA::iteration()
{
	auto updatePhi = [&](TopicId k, WordId v, double delta){
		phi_[k][v] = eta_[k][v] + delta;
	};

	double term2;
	VectorK<double> alpha_sufficient_statistics(K_, 0);

	mapreduce_->run<mapreduce::schedule_policy::cpu_parallel<mr_job>>(performance_result_);

	for (auto it = mapreduce_->begin_results(), end = mapreduce_->end_results(); it != end; ++it){
		if (std::get<0>(it->first) == ReduceKeyType::Lambda){
			updatePhi(std::get<1>(it->first), std::get<2>(it->first), it->second);
		}
		else if (std::get<0>(it->first) == ReduceKeyType::Alpha){
			alpha_sufficient_statistics[std::get<2>(it->first)] = it->second;
		}
		else{
			term2 = it->second;
		}
	}
	//updateAlpha(alpha_, alpha_sufficient_statistics, D_);

	for(auto& b : phi_) sig::normalize_dist(b);

	saveResumeData();

	return 0;//calcLiklihood(term2, sig::sum(lambda_, [&](VectorV<double> const& row){ return calcModule0(row); }));
}

void MrLDA::train(uint iteration_num, std::function<void(LDA const*)> callback)
{
	sig::ManageConvergenceSimple convergence(global_convergence_threshold);

	/*while (convergence.update(iteration())){
		std::cout << convergence.get_value() << std::endl;
	}*/

	for (uint i = 0; i < iteration_num; ++i, ++total_iter_ct_){
		iteration();
		calcTermScore(getPhi(), term_score_);
		save(Distribution::TOPIC, L"./test data");
		save(Distribution::TERM_SCORE, L"./test data");
		save(Distribution::DOCUMENT, L"./test data");
		callback(this);
		initMR();
	}
}

void MrLDA::save(Distribution target, FilepassString save_folder, bool detail) const
{
	save_folder = sig::modify_dirpass_tail(save_folder, true);

	switch (target){
	case Distribution::DOCUMENT:
		printTopic(getTheta(), input_data_->getInputFileNames(), save_folder + SIG_TO_FPSTR("document_mrlda"));
		break;
	case Distribution::TOPIC:
		printWord(getPhi(), std::vector<FilepassString>(), input_data_->words_, detail ? nothing : sig::Just<uint>(20), save_folder + SIG_TO_FPSTR("topic_mrlda"));
		break;
	case Distribution::TERM_SCORE:
		printWord(getTermScore(), std::vector<FilepassString>(), input_data_->words_, detail ? nothing : sig::Just<uint>(20), save_folder + SIG_TO_FPSTR("term-score_mrlda"));
		break;
	default:
		std::cout << "MrLDA::save error" << std::endl;
		getchar();
	}
}

auto MrLDA::getTheta() const->MatrixDK<double>
{
	MatrixDK<double> theta;

	for (DocumentId d = 0; d < D_; ++d){
		theta.push_back(getTheta(d));
	}
	return theta;
}

auto MrLDA::getTheta(DocumentId d_id) const->VectorK<double>
{
	double sum = sig::sum(gamma_[d_id]);
	// computed from the variational distribution
	return sig::map([sum](double g){ return g / sum; }, gamma_[d_id]);
}

auto MrLDA::getPhi() const->MatrixKV<double>
{
	MatrixKV<double> phi;

	for (TopicId k = 0; k < K_; ++k){
		phi.push_back(getPhi(k));
	}
	return phi;
}

auto MrLDA::getPhi(TopicId k_id) const->VectorV<double>
{
	/*double sum = sig::sum(lambda_[k_id]);
	// computed from the variational distribution
	return sig::map([sum](double l){ return l / sum; }, lambda_[k_id]);*/
	return phi_[k_id];
}

auto MrLDA::getWordOfTopic(Distribution target, uint return_word_num) const->VectorK< std::vector< std::tuple<std::wstring, double> > >
{
	VectorK< std::vector< std::tuple<std::wstring, double> > > result;

	for (TopicId k = 0; k < K_; ++k){
		result.push_back(getWordOfTopic(target, return_word_num, k));
	}
	return result;
}

auto MrLDA::getWordOfTopic(Distribution target, uint return_word_num, TopicId k_id) const->std::vector< std::tuple<std::wstring, double> >
{
	std::vector< std::tuple<std::wstring, double> > result;
	std::vector<double> df;

	if (target == Distribution::TOPIC) df = getPhi(k_id);
	else if (target == Distribution::TERM_SCORE) df = getTermScore(k_id);
	else{
		std::cout << "MrLDA::getWordOfTopic : Distributionが無効" << std::endl;
		return result;
	}
	return getTopWords(df, return_word_num, input_data_->words_);
}

auto MrLDA::getWordOfDocument(uint return_word_num) const->VectorD< std::vector< std::tuple<std::wstring, double> > >
{
	VectorD< std::vector< std::tuple<std::wstring, double> > > result;

	for (DocumentId d = 0; d<D_; ++d){
		result.push_back(getWordOfDocument(return_word_num, d));
	}
	return result;
}

auto MrLDA::getWordOfDocument(uint return_word_num, DocumentId d_id) const->std::vector< std::tuple<std::wstring, double> >
{
	std::vector< std::tuple<std::wstring, double> > result;
	auto top_wscore = getTermScoreOfDocument(getTheta(d_id), getTermScore());

	for (uint i = 0; i<return_word_num; ++i){
		result.push_back(std::make_tuple(*input_data_->words_.getWord(std::get<0>(top_wscore[i])), std::get<1>(top_wscore[i])));
	}
	return std::move(result);
}

}