text:
#include <chrono>
#include <ctime>
#include <iostream>
#include <fstream>
#include <vector>
#include <Eigen/Geometry>
using namespace std;

//default input values
string inputFile = "trajBig.pdb";
string outputFile = "output.txt";
int spheresAllocationFrame = 0;
double RMSDThreshPercent = 50;

//additional parameters
double sphereRadius = 8;

int SPHERES;
int ATOMS;
int FRAMES;
double highestSize;

//Atoms[<frame>][<atom>][<coordinate>]
vector<vector<vector<double>>> A;

//RMSDValue[<frame>][<sphere>]
vector<vector<double>> C;

//Maps sphere to CA; CAAtomNumber[<sphere>]
vector<int> sphereCA;

//Numer of atoms in [<sphere>]
vector<int> sphereSize;

struct sphere
{
    int frame;
    int index;
    double RMSD;
};

vector<vector<int>> sphereAtoms;
vector<vector<int>> atomToSphere;
vector<vector<double>> highestRMSD;
vector<sphere> highestRMSDbase;

bool operator <(const sphere& lhs, const sphere& rhs)
{
    return lhs.RMSD < rhs.RMSD;
}

vector<string> result;

//--------------------------------------------------

//reading data from input pdb file.
void readFile(string filename)
{
    string line;
    ifstream myfile(filename);
    if (myfile.is_open())
    {
        int frame = 0;
        int atom;
        SPHERES = 0;
        FRAMES = 0;
        ATOMS = 0;
        A = {};
        sphereCA = {};
        sphereSize = {};
        while (getline(myfile, line))
        {
            if (line[0] == 'M')
            {
                frame = stoi(line.substr(9, 5));
                frame--;
                A.push_back({});
                C.push_back({});
                FRAMES++;
            }
            else if (line[0] == 'A')
            {
                atom = stoi(line.substr(6, 5));
                atom--;
                A[frame].push_back({});
                A[frame][atom].push_back(stod(line.substr(30, 8)));
                A[frame][atom].push_back(stod(line.substr(38, 8)));
                A[frame][atom].push_back(stod(line.substr(46, 8)));
                if (frame == 0)
                {
                    ATOMS++;
                    if (line[14] == 'A' and line[13] == 'C')
                    {
                        sphereCA.push_back(atom);
                        sphereSize.push_back(0);
                        SPHERES++;
                    }
                }
            }
        }
        myfile.close();
    }
    else
    {
        cout << "Nie odnaleziono pliku!" << endl;
    }
}

//calculating distance between 2 atoms, used when allocating atoms to spheres.
double atomsDistanceCalc(int atom1, int atom2)
{
    double result = sqrt(pow(A[spheresAllocationFrame][atom1][0] - A[spheresAllocationFrame][atom2][0], 2) +
        pow(A[spheresAllocationFrame][atom1][1] - A[spheresAllocationFrame][atom2][1], 2) +
        pow(A[spheresAllocationFrame][atom1][2] - A[spheresAllocationFrame][atom2][2], 2));
    return result;
}

//allocating atoms into spheres, based on sphereRadius
void atomsAllocation()
{
    auto timestamp = std::chrono::system_clock::now();
    std::time_t timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    int checkpoint = 0;
    sphereAtoms = {};
    atomToSphere = {};
    for (int i = 0; i < SPHERES; i++)
    {
        sphereAtoms.push_back({});
    }
    for (int i = 0; i < ATOMS; i++)
    {
        atomToSphere.push_back({});
        for (int j = 0; j < SPHERES; j++)
        {
            if (atomsDistanceCalc(i, sphereCA[j]) <= sphereRadius)
            {
                atomToSphere[i].push_back(j);
                sphereAtoms[j].push_back(i);
                sphereSize[j]++;
            }
        }
        if ((i - checkpoint) >= 0.1 * ATOMS)
        {
            checkpoint = i;
            timestamp = std::chrono::system_clock::now();
            timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
            cout << "[RMSD] Atom number " << i + 1 << " of " << ATOMS << " allocated to spheres at " << std::ctime(&timestamp_print) << endl;
        }
    }
}

// Find3DAffineTransform is from oleg-alexandrov repository on github, available here https://github.com/oleg-alexandrov/projects/blob/master/eigen/Kabsch.cpp [as of 27.01.2022]
// Given two sets of 3D points, find the rotation + translation + scale
// which best maps the first set to the second.
// Source: http://en.wikipedia.org/wiki/Kabsch_algorithm

// The input 3D points are stored as columns.
Eigen::Affine3d Find3DAffineTransform(Eigen::Matrix3Xd in, Eigen::Matrix3Xd out) {

    // Default output
    Eigen::Affine3d A;
    A.linear() = Eigen::Matrix3d::Identity(3, 3);
    A.translation() = Eigen::Vector3d::Zero();

    if (in.cols() != out.cols())
        throw "Find3DAffineTransform(): input data mis-match";

    // First find the scale, by finding the ratio of sums of some distances,
    // then bring the datasets to the same scale.
    double dist_in = 0, dist_out = 0;
    for (int col = 0; col < in.cols() - 1; col++) {
        dist_in += (in.col(col + 1) - in.col(col)).norm();
        dist_out += (out.col(col + 1) - out.col(col)).norm();
    }
    if (dist_in <= 0 || dist_out <= 0)
        return A;
    double scale = dist_out / dist_in;
    out /= scale;

    // Find the centroids then shift to the origin
    Eigen::Vector3d in_ctr = Eigen::Vector3d::Zero();
    Eigen::Vector3d out_ctr = Eigen::Vector3d::Zero();
    for (int col = 0; col < in.cols(); col++) {
        in_ctr += in.col(col);
        out_ctr += out.col(col);
    }
    in_ctr /= in.cols();
    out_ctr /= out.cols();
    for (int col = 0; col < in.cols(); col++) {
        in.col(col) -= in_ctr;
        out.col(col) -= out_ctr;
    }

    // SVD
    Eigen::MatrixXd Cov = in * out.transpose();
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(Cov, Eigen::ComputeThinU | Eigen::ComputeThinV);

    // Find the rotation
    double d = (svd.matrixV() * svd.matrixU().transpose()).determinant();
    if (d > 0)
        d = 1.0;
    else
        d = -1.0;
    Eigen::Matrix3d I = Eigen::Matrix3d::Identity(3, 3);
    I(2, 2) = d;
    Eigen::Matrix3d R = svd.matrixV() * I * svd.matrixU().transpose();

    // The final transform
    A.linear() = scale * R;
    A.translation() = scale * (out_ctr - R * in_ctr);

    return A;
}

//superpose changes vector of atoms frame 2 to map atoms from frame 1 in the way to minimise RMSD between both frames
void superpose(const vector<vector<double>>& frame1, vector<vector<double>>& frame2)
{
    int atomsInSphere = frame1.size();
    Eigen::Matrix3Xd S1(3, atomsInSphere);
    Eigen::Matrix3Xd S2(3, atomsInSphere);
    for (int j = 0; j < atomsInSphere; j++)
    {
        for (int k = 0; k < 3; k++)
        {
            S1(k, j) = frame1[j][k];
            S2(k, j) = frame2[j][k];
        }
    }
    Eigen::Affine3d RT = Find3DAffineTransform(S2, S1);
    S2 = RT.linear() * S2;
    for (int j = 0; j < atomsInSphere; j++)
    {
        S2.block<3, 1>(0, j) += RT.translation();
    }
    for (int j = 0; j < atomsInSphere; j++)
    {
        for (int k = 0; k < 3; k++)
        {
            frame2[j][k] = S2(k, j);
        }
    }
}

//initializing vector for RMSD values
void initializeC()
{
    C = {};
    for (int i = 0; i < FRAMES; i++)
    {
        C.push_back({});
        for (int j = 0; j < SPHERES; j++)
            C[i].push_back(0);
    }
}

//calculating RMSD on spheres, on adjacent frames
void calculateRMSDSuperpose()
{
    initializeC();
    vector<vector<vector<double>>> sphereMatrix;
    vector<vector<double>> tempMatrix;
    double tempRMSD;
    auto timestamp = std::chrono::system_clock::now();
    std::time_t timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    int checkpoint = 0;
    for (int s = 0; s < SPHERES; s++)
    {
        int atomsInSphere = sphereAtoms[s].size();
        sphereMatrix = {};
        for (int i = 0; i < FRAMES; i++)
        {
            sphereMatrix.push_back({});
            for (int j = 0; j < atomsInSphere; j++)
            {
                sphereMatrix[i].push_back(A[i][sphereAtoms[s][j]]);
            }
            if (i > 0)
            {
                tempMatrix = sphereMatrix[i];
                superpose(sphereMatrix[i - 1], tempMatrix);
                for (int j = 0; j < atomsInSphere; j++)
                {
                    for (int k = 0; k < 3; k++)
                    {
                        tempRMSD = pow(tempMatrix[j][k] - sphereMatrix[i - 1][j][k], 2);
                        C[i - 1][s] += tempRMSD;
                    }
                }
                C[i - 1][s] /= ((long)atomsInSphere * (long)3);
                C[i - 1][s] = sqrt(C[i - 1][s]);
            }
        }
        if ((s - checkpoint) >= 0.1 * SPHERES)
        {
            checkpoint = s;
            timestamp = std::chrono::system_clock::now();
            timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
            cout << "[RMSD] RMSD calculated for " << s + 1 << " of " << SPHERES << " spheres at " << std::ctime(&timestamp_print) << endl;
        }
    }
}

//initializing vector for storing information about spheres with highest RMSD
void initializeHighestRMSD()
{
    highestRMSD = {};
    highestSize = (double)SPHERES / 100;
    highestSize *= (long)(FRAMES - 1);
    highestSize *= RMSDThreshPercent;
    for (int i = 0; i < highestSize; i++)
    {
        highestRMSD.push_back({});
        for (int j = 0; j < 3; j++)
            highestRMSD[i].push_back(-1);
    }
}

//choosing frames and spheres with highest RMSD, they will be colored later
void chooseHighestRMSD()
{
    initializeHighestRMSD();
    for (int i = 0; i < FRAMES - 1; i++)
    {
        for (int j = 0; j < SPHERES; j++)
        {
            for (unsigned int k = 0; k < highestRMSD.size(); k++)
            {
                if (C[i][j] > highestRMSD[k][2])
                {
                    for (int l = highestRMSD.size() - 1; l > k; l--)
                    {
                        for (int m = 0; m < 3; m++)
                        {
                            highestRMSD[l][m] = highestRMSD[l - 1][m];
                        }
                    }
                    highestRMSD[k][0] = i;
                    highestRMSD[k][1] = j;
                    highestRMSD[k][2] = C[i][j];
                    break;
                }
            }
        }
    }
}

//choosing frames and spheres with highest RMSD, they will be colored later
void chooseHighestRMSDquick()
{
    auto timestamp = std::chrono::system_clock::now();
    std::time_t timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    initializeHighestRMSD();
    int counter = 0;
    for (int i = 0; i < FRAMES - 1; i++)
    {
        for (int j = 0; j < SPHERES; j++)
        {
            highestRMSDbase.push_back(sphere());
            highestRMSDbase[counter].frame = i;
            highestRMSDbase[counter].index = j;
            highestRMSDbase[counter].RMSD = C[i][j];
            counter++;
        }
    }
    timestamp = std::chrono::system_clock::now();
    timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    cout << "[RMSD] RMSD started being sorted at " << std::ctime(&timestamp_print) << endl;
    sort(highestRMSDbase.begin(), highestRMSDbase.end());
    timestamp = std::chrono::system_clock::now();
    timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    cout << "[RMSD] RMSD sorted at " << std::ctime(&timestamp_print) << endl;
    for (unsigned int k = 0; k < highestRMSD.size(); k++)
    {
        highestRMSD[k][0] = highestRMSDbase[k].frame;
        highestRMSD[k][1] = highestRMSDbase[k].index;
        highestRMSD[k][2] = highestRMSDbase[k].RMSD;
    }
}

//calculates result variable content, basen on highestRMSD
void calculateResult()
{
    result = {};
    int numberOfSpheresToColor = highestRMSD.size();
    vector<vector<int>> resultSpheres;
    for (int i = 0; i < FRAMES - 1; i++)
    {
        resultSpheres.push_back({});
    }
    for (int i = 0; i < numberOfSpheresToColor; i++)
    {
        resultSpheres[highestRMSD[i][0]].push_back(highestRMSD[i][1]);
    }
    string subResult;
    for (int i = 0; i < FRAMES - 1; i++)
    {
        subResult = "";
        for (int j = 0; j < numberOfSpheresToColor; j++)
        {
            if (highestRMSD[j][0] == i)
            {
                stringstream ss;
                string s;
                ss << sphereCA[highestRMSD[j][1]];
                ss >> s;
                subResult += s + " ";
            }
        }
        result.push_back(subResult);
    }
}

//saving results to an output file
void saveToFile(string filename)
{
    ofstream myfile(filename);
    if (myfile.is_open())
    {
        for (int i = 0; i < FRAMES - 1; i++)
        {
            if (!result[i].empty())
            {
                myfile << i << ' ' << result[i] << endl;
            }
        }
        myfile.close();
    }
}

//debug function for checking superpose results
void makePDB(string filename)
{
    vector<vector<vector<double>>> sphereMatrix;
    int atomsInSphere = sphereAtoms[0].size();
    sphereMatrix = {};
    for (int i = 0; i < 2; i++)
    {
        sphereMatrix.push_back({});
        for (int j = 0; j < atomsInSphere; j++)
        {
            sphereMatrix[i].push_back(A[i][sphereAtoms[0][j]]);
        }
        if (i > 0)
        {
            superpose(sphereMatrix[i - 1], sphereMatrix[i]);
        }
    }
    string line;
    ifstream myfilein(filename);
    filename += "out";
    ofstream myfileout("testBase.pdb");
    ofstream myfilesuper("testSuperpose.pdb");
    int counter = 0;
    if (myfilein.is_open() && myfileout.is_open())
    {
        int frame = 0;
        int atom;
        while (getline(myfilein, line))
        {
            if (line[0] == 'M')
            {
                myfileout << line << endl;
                myfilesuper << line << endl;
            }

            if (line[0] == 'E')
            {
                myfileout << line << endl;
                myfilesuper << line << endl;
                frame++;
                if (frame > 1)
                {
                    break;
                }
            }
            if (line[0] == 'A')
            {
                atom = stoi(line.substr(6, 5));
                atom--;
                for (unsigned int i = 0; i < sphereAtoms[0].size(); i++)
                {
                    if (sphereAtoms[0][i] == atom)
                    {
                        myfileout << line << endl;
                        if (frame == 1)
                        {
                            line.replace(32, 6, to_string(sphereMatrix[1][counter][0]), 0, 6);
                            line.replace(40, 6, to_string(sphereMatrix[1][counter][1]), 0, 6);
                            line.replace(48, 6, to_string(sphereMatrix[1][counter][2]), 0, 6);
                            myfilesuper << line << endl;
                            counter++;
                        }
                        else
                        {
                            myfilesuper << line << endl;
                        }
                    }
                }
            }
        }
        myfilein.close();
        myfileout.close();
        myfilesuper.close();
    }
}

int main(int argc, char* argv[])
{
    auto start = chrono::high_resolution_clock::now();
    if (argc == 5)
    {
        inputFile = string(argv[1]);
        outputFile = string(argv[2]);
        spheresAllocationFrame = stoi(argv[3]) - 1;
        RMSDThreshPercent = stod(argv[4]);
    }
    readFile(inputFile);
    if (spheresAllocationFrame < 0)
    {
        spheresAllocationFrame = 1;
    }
    else if (spheresAllocationFrame >= FRAMES)
    {
        spheresAllocationFrame = FRAMES - 1;
    }
    auto timestamp = std::chrono::system_clock::now();
    std::time_t timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    cout << "[RMSD] Input file: " << inputFile << " loaded at " << std::ctime(&timestamp_print) << endl;
    atomsAllocation();
    timestamp = std::chrono::system_clock::now();
    timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    cout << "[RMSD] Atoms allocated to spheres at " << std::ctime(&timestamp_print) << endl;
    calculateRMSDSuperpose();
    chooseHighestRMSDquick();
    calculateResult();
    saveToFile(outputFile);
    timestamp = std::chrono::system_clock::now();
    timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    cout << "[RMSD] Output file: " << outputFile << " saved at " << std::ctime(&timestamp_print) << endl;
    //makePDB(inputFile);
    auto stop = chrono::high_resolution_clock::now();
    chrono::duration<double> elapsed = stop - start;
    timestamp = std::chrono::system_clock::now();
    timestamp_print = std::chrono::system_clock::to_time_t(timestamp);
    cout << "[RMSD] Computation time: " << elapsed.count() << "s, ended at: " << std::ctime(&timestamp_print) << endl;
}

max_stars_repo_name: KrzysztofMularski/ProProtein

max_issues_repo_name: KrzysztofMularski/ProProtein

max_forks_repo_name: KrzysztofMularski/ProProtein


text:
/// Copyright 2012 by BBN Technologies Corp.
// All Rights Reserved.

// We could potentially do some caching here, if this is a bottleneck.

#include "Generic/common/leak_detection.h"

#include "Generic/actors/AWAKEActorInfo.h"
#include "Generic/actors/AWAKEDB.h"
#include "Generic/database/DatabaseConnection.h"
#include "Generic/common/SessionLogger.h"
#include "Generic/common/ParamReader.h"
#include "Generic/theories/EntityType.h"

#include <boost/make_shared.hpp>
#include <boost/foreach.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string.hpp>

namespace {
	Symbol PER = Symbol(L"PER");
	Symbol ORG = Symbol(L"ORG");
}

bool AWAKEActorInfo::_initialized = false;
AWAKEActorInfo_ptr AWAKEActorInfo::_actorInfo = AWAKEActorInfo_ptr();

struct AWAKEActorInfo::ActorRow {
	ActorRow() { };
	ActorRow(const std::wstring &canonicalName, Symbol entityType, Symbol entitySubtype, Symbol isoCode, double importance_score)
		: canonicalName(canonicalName), entityType(entityType), entitySubtype(entitySubtype), isoCode(isoCode), importance_score(importance_score) { }

	std::wstring canonicalName;
	Symbol entityType;
	Symbol entitySubtype;
	Symbol isoCode;
	double importance_score;
};

struct AWAKEActorInfo::AgentRow {
	AgentRow() {};
	AgentRow(const std::wstring &name) : name(name){}
	std::wstring name;
	std::vector<Symbol> associatedSectorCodes;
};

struct AWAKEActorInfo::Cache {
public:
	// These return refereneces to shared pointers.  The shared pointers will be null if the
	// value hasn't been added to the cache yet; in that case, simply assign a new value
	// to the shared pointer (which is returned by reference).
	boost::shared_ptr<std::vector<ActorId> > &getAssociatedActorIds(ActorId actorId, const char *date) {
		return lookup(associatedActorIds, actorId, date); }
	boost::shared_ptr<std::vector<ActorId> > &getAssociatedCountryActorIds(ActorId actorId, const char *date) {
		return lookup(associatedCountryActorIds, actorId, date); }
	boost::shared_ptr<std::vector<ActorId> > &getAssociatedLocationActorIds(ActorId actorId, const char *date) {
		return lookup(associatedLocationActorIds, actorId, date); }
	boost::shared_ptr<std::vector<SectorId> > &getAssociatedSectorIds(ActorId actorId, const char *date) {
		return lookup(associatedSectorIds, actorId, date); }
	boost::shared_ptr<std::vector<CountryId> > &getAssociatedCountryIds(ActorId actorId, const char *date) {
		return lookup(associatedCountryIds, actorId, date); }
	boost::shared_ptr<std::vector<Symbol> > &getAssociatedSectorCodes(ActorId actorId, const char *date) {
		return lookup(associatedSectorCodes, actorId, date); }

	Cache(): num_hits(0), num_misses(0) {}
	~Cache() {}

	const ActorRow& saveActorRow(DatabaseConnection::RowIterator &row, ActorRow& result) const {
		result.canonicalName = row.getCellAsWString(0);
		result.entityType = row.getCellAsSymbol(1);
		result.entitySubtype = row.getCellAsSymbol(2);
		result.isoCode = row.getCellAsSymbol(3);
		result.importance_score = row.getCellAsDouble(4, 0.0);
		return result;
	}

	// LT: don't just automatically search in the default database but look in the database
	// specified within the ActorId type target
	const ActorRow& getActorRow(ActorId target) {
		ActorId::HashMap<ActorRow>::iterator it = actorRows.find(target);
		if (it != actorRows.end()) {
			return (*it).second;
		} else {
			std::ostringstream query;
			query << "SELECT CanonicalName, EntityType, EntitySubtype, IsoCode, ImportanceScore"
				  << " FROM Actor a LEFT OUTER JOIN ActorIsoCode aic ON a.ActorId=aic.ActorId"
				  << " WHERE a.ActorId=" << target.getId();
			DatabaseConnection_ptr bbn_db(AWAKEDB::getDb(target));
			for (DatabaseConnection::RowIterator row = bbn_db->iter(query); row != bbn_db->end(); ++row) {
				if (actorRows.size() >= Cache::MAX_SIZE)
					actorRows.clear();
				return saveActorRow(row, actorRows[target]);
			}
			static ActorRow nullResult(L"UNKNOWN-ACTOR", Symbol(L"unknown"), Symbol(L"unknown"), Symbol(L"unknown"), 0.0);
			return nullResult;
		}
	}

	const AgentRow& getAgentRow(AgentId target) {
		if (agentRows.size()==0) readAgentTable();
		AgentId::HashMap<AgentRow>::iterator it = agentRows.find(target);
		if (it != agentRows.end()) {
			return (*it).second;
		} else {
			static AgentRow nullResult(L"UNKNOWN-AGENT");
			return nullResult;
		}
	}

	AgentId getAgentByName(Symbol name) {
		if (agentRows.size()==0) readAgentTable();
		Symbol::HashMap<AgentId>::iterator it = agentByName.find(name);
		if (it != agentByName.end()) {
			return (*it).second;
		} else {
			return AgentId();
		}
	}

	AgentId getDefaultPersonAgentId() {
		if (agentRows.size()==0) readAgentTable();
		return defaultPersonAgentId;
	}

	Symbol getDefaultPersonAgentCode() {
		if (agentRows.size()==0) readAgentTable();
		return defaultPersonAgentCode;
	}

private:
	size_t num_hits;
	size_t num_misses;

	ActorId::HashMap<ActorRow> actorRows;
	AgentId::HashMap<AgentRow> agentRows;
	Symbol::HashMap<AgentId> agentByName;

	typedef std::pair<ActorId, std::string> KeyType;
	struct HashOp {
		size_t operator()(const KeyType &v) const {
			return static_cast<size_t>(v.first.getId()); }
	};
	struct EqualOp {
		size_t operator()(const KeyType &v1, const KeyType &v2) const {
			return static_cast<size_t>((v1.first==v2.first) && (v1.second==v2.second));}
	};
	template<typename IdType>
	class IdVectorMap: public serif::hash_map<KeyType, boost::shared_ptr<std::vector<IdType> >, HashOp, EqualOp> {};

	IdVectorMap<ActorId> associatedActorIds;
	IdVectorMap<ActorId> associatedCountryActorIds;
	IdVectorMap<ActorId> associatedLocationActorIds;
	IdVectorMap<SectorId> associatedSectorIds;
	IdVectorMap<CountryId> associatedCountryIds;
	IdVectorMap<Symbol> associatedSectorCodes;

	template<typename IdType>
	boost::shared_ptr<std::vector<IdType> > &lookup(IdVectorMap<IdType> &idVectorMap, ActorId target, const char *date) {
		boost::shared_ptr<std::vector<IdType> > &result = idVectorMap[std::make_pair(target, date?date:"")];
		if ((!result) && (idVectorMap.size() > Cache::MAX_SIZE))
			idVectorMap.clear();
		return result;
	}

	static const size_t MAX_SIZE = 1000;

	AgentId defaultPersonAgentId;
	Symbol defaultPersonAgentCode;

	// LT: iterate through all databases to find the actorDBName
	// remove actorDBName--instead create AgentID using db_name from for loop
	void readAgentTable() {

		defaultPersonAgentId = AgentId();
		defaultPersonAgentCode = Symbol(L"UNKNOWN-AGENT");

		std::ostringstream query;
		query << "SELECT AgentId, AgentName FROM Agent";
		std::ostringstream sector_query;
		sector_query << "SELECT AgentId, SectorName FROM AgentSectorLink ASL, Sector S WHERE ASL.SectorId = S.SectorId";
		DatabaseConnectionMap dbs = AWAKEDB::getNamedDbs();
		for (DatabaseConnectionMap::iterator i = dbs.begin(); i != dbs.end(); ++i) {
			for (DatabaseConnection::RowIterator row = i->second->iter(query); row!=i->second->end(); ++row) {
				AgentId agentId(row.getCellAsInt32(0), i->first);
				AgentRow& agentRow = agentRows[agentId];
				agentRow.name = row.getCellAsWString(1);
				agentByName[Symbol(agentRow.name)] = agentId;
				if (agentRow.name == L"Civilian" || agentRow.name == L"Citizen") {
					defaultPersonAgentId = agentId;
					defaultPersonAgentCode = Symbol(agentRow.name);
				}
			}
			// Read in the agent/sector mappings as well.
			// It's OK to use SectorName as the unique "code" here, since the DB constrains this value to be unique
			for (DatabaseConnection::RowIterator row = i->second->iter(sector_query); row!=i->second->end(); ++row) {
				AgentId agentId(row.getCellAsInt32(0), i->first);
				agentRows[agentId].associatedSectorCodes.push_back(row.getCellAsSymbol(1));
			}
		}
		// Add an entry for the NULL agent id.
		agentRows[AgentId()] = AgentRow(L"UNKNOWN-AGENT");

		if (defaultPersonAgentId.isNull()) {
			throw UnexpectedInputException("AWAKEActorInfo::readAgentTable", "No default person agent ('Civilian' or 'Citizen') in dictionary");
		}
	}
};

// LT: this is fine to leave as default for now because we probably won't need link table from more than one database
AWAKEActorInfo::AWAKEActorInfo(): ActorInfo(), _affiliation_link_type_id(-1), _country_link_type_id(-1), _location_link_type_id(-1) {
	DatabaseConnection_ptr bbn_db(AWAKEDB::getDefaultDb());

	std::ostringstream query1;
	query1 << "SELECT LinkTypeId FROM LinkType"
		   << " WHERE LinkType='Affiliation'";
	for (DatabaseConnection::RowIterator row = bbn_db->iter(query1); row != bbn_db->end(); ++row)
		_affiliation_link_type_id = row.getCellAsInt32(0);

	std::ostringstream query2;
	query2 << "SELECT LinkTypeId FROM LinkType"
		   << " WHERE LinkType='Country'";
	for (DatabaseConnection::RowIterator row = bbn_db->iter(query2); row != bbn_db->end(); ++row)
		_country_link_type_id = row.getCellAsInt32(0);

	std::ostringstream query3;
	query3 << "SELECT LinkTypeId FROM LinkType"
		   << " WHERE LinkType='Location'";
	for (DatabaseConnection::RowIterator row = bbn_db->iter(query3); row != bbn_db->end(); ++row)
		_location_link_type_id = row.getCellAsInt32(0);

	//LT: removing references to _actorDBName
	//_actorDBName = ParamReader::getParam(Symbol(L"bbn_actor_db_name"));
	//if (_actorDBName.is_null())
	//	_actorDBName = Symbol(L"default");

	_cache = _new AWAKEActorInfo::Cache();

	loadActorPatterns();
}

ActorInfo_ptr AWAKEActorInfo::getAWAKEActorInfo() {
	if (!_initialized) {
		_actorInfo = boost::make_shared<AWAKEActorInfo>();
		_initialized = true;
	}
	return _actorInfo;
}

AWAKEActorInfo::~AWAKEActorInfo() {
	delete _cache;
};

CountryId AWAKEActorInfo::getCountryId(ActorId target) {
	if (!isACountry(target))
		return CountryId();

	return CountryId(target.getId(), target.getDbName());
}

bool AWAKEActorInfo::mightBeALocation(ActorId target) {
	if (target.isNull()) return false;
	return _cache->getActorRow(target).entityType == EntityType::getGPEType().getName() ||
		   _cache->getActorRow(target).entityType == EntityType::getLOCType().getName();
}

Symbol AWAKEActorInfo::getActorCode(ActorId target) {
	std::wstringstream str;
	str << target.getId();
	return Symbol(str.str());
}

Symbol AWAKEActorInfo::getAgentCode(AgentId target) {
	std::wstringstream str;
	str << target.getId();
	return Symbol(str.str());
}

boost::shared_ptr<std::vector<ActorId> > AWAKEActorInfo::getAssociatedActorIds(ActorId target, std::vector<int> link_types, const char *date) {
	boost::shared_ptr<std::vector<ActorId> > result = boost::make_shared<std::vector<ActorId> >();

	// LT: replaced getDefaultDb with getDb(target)
	DatabaseConnection_ptr bbn_db(AWAKEDB::getDb(target));
	std::ostringstream query;

	query << "SELECT DISTINCT RightActorId FROM ActorLink"
		<< " WHERE LeftActorId=" << target.getId();
	if (date)
		query << " AND ((StartDate is null) or (StartDate < " << bbn_db->toDate(date) << "))"
		<< " AND ((EndDate is null) or (EndDate > " << bbn_db->toDate(date) << "))";
	if (link_types.size() > 0) {
		query << " AND LinkTypeId in (";
		for (std::vector<int>::iterator iter = link_types.begin(); iter != link_types.end(); iter++) {
			if (iter != link_types.begin())
				query << ", ";
			query << (*iter);
		}
		query << ")";
	}

	for (DatabaseConnection::RowIterator row = bbn_db->iter(query); row!=bbn_db->end(); ++row) {
		result->push_back(ActorId(row.getCellAsInt32(0), target.getDbName()));
	}

	return result;
}

std::vector<ActorId> AWAKEActorInfo::getAssociatedActorIds(ActorId target, const char *date) {
	boost::shared_ptr<std::vector<ActorId> >& result = _cache->getAssociatedActorIds(target, date);

	if (!result) {
		std::vector<int> link_types;
		result = getAssociatedActorIds(target, link_types, date); // This adds to the cache
	}

	return (*result);
}

std::vector<ActorId> AWAKEActorInfo::getAssociatedCountryActorIds(ActorId target, const char *date) {
	boost::shared_ptr<std::vector<ActorId> >& result = _cache->getAssociatedCountryActorIds(target, date);

	if (!result) {
		std::vector<int> link_types;
		link_types.push_back(_country_link_type_id);
	    result = getAssociatedActorIds(target, link_types, date); // This adds to the cache
	}

	return (*result);
}

std::vector<ActorId> AWAKEActorInfo::getAssociatedLocationActorIds(ActorId target, const char *date) {
	boost::shared_ptr<std::vector<ActorId> >& result = _cache->getAssociatedLocationActorIds(target, date);

	if (!result) {
		std::vector<int> link_types;
		link_types.push_back(_location_link_type_id);
		result = getAssociatedActorIds(target, link_types, date); // This adds to the cache
	}

	return (*result);
}

std::vector<CountryId> AWAKEActorInfo::getAssociatedCountryIds(ActorId target, const char *date) {
	boost::shared_ptr<std::vector<CountryId> >& result = _cache->getAssociatedCountryIds(target, date);

	if (!result) {
		result = boost::make_shared<std::vector<CountryId> >();

	    std::vector<ActorId> actorIds = getAssociatedCountryActorIds(target, date);
		BOOST_FOREACH(ActorId aid, actorIds) {
			result->push_back(CountryId(aid.getId(), aid.getDbName()));
		}
	}

	return *result;
}

std::vector<SectorId> AWAKEActorInfo::getAssociatedSectorIds(ActorId target, const char *date) {
	boost::shared_ptr<std::vector<SectorId> >& result = _cache->getAssociatedSectorIds(target, date);

	if (!result) {
		// LT: replaced getDefaultDb with getDb(target)
		DatabaseConnection_ptr bbn_db(AWAKEDB::getDb(target));
		result = boost::make_shared<std::vector<SectorId> >();
		std::ostringstream query;

		query << "SELECT DISTINCT SectorId FROM ActorSectorLink"
			  << " WHERE ActorId=" << target.getId();
		if (date)
			query << " AND ((StartDate is null) or (StartDate < " << bbn_db->toDate(date) << "))"
			      << " AND ((EndDate is null) or (EndDate > " << bbn_db->toDate(date) << "))";

		for (DatabaseConnection::RowIterator row = bbn_db->iter(query); row != bbn_db->end(); ++row) {
			result->push_back(SectorId(row.getCellAsInt32(0), target.getDbName()));
		}
	}
	return *result;
}

std::vector<Symbol> AWAKEActorInfo::getAssociatedSectorCodes(ActorId target, const char *date, int min_frequency) {
	boost::shared_ptr<std::vector<Symbol> >& result = _cache->getAssociatedSectorCodes(target, date);

	if (!result) {
		// LT: replaced getDefaultDb with getDb(target)
		DatabaseConnection_ptr bbn_db(AWAKEDB::getDb(target));
		result = boost::make_shared<std::vector<Symbol> >();
		std::vector<int> idResults;
		std::ostringstream query;

		query << "SELECT DISTINCT s.SectorName, s.SectorId FROM ActorSectorLink asl, Sector s"
			  << " WHERE ActorId=" << target.getId()
			  << " AND asl.SectorId=s.SectorId";
		if (date)
			query << " AND ((asl.StartDate is null) or (asl.StartDate < " << bbn_db->toDate(date) << "))"
			<< " AND ((asl.EndDate is null) or (asl.EndDate > " << bbn_db->toDate(date) << "))";

		for (DatabaseConnection::RowIterator row = bbn_db->iter(query); row != bbn_db->end(); ++row) {
			result->push_back(row.getCellAsSymbol(0));
			idResults.push_back(row.getCellAsInt32(1));
		}

		// Get parent sectors as well, (and parents of parents, etc.)
		BOOST_FOREACH(int id, idResults) {
			while (id != -1) {
				std::ostringstream parentSectorQuery;
				parentSectorQuery
					<< "SELECT s2.SectorName, s2.SectorId FROM Sector s1, Sector s2"
					<< " WHERE s1.SectorId=" << id
					<< " AND s1.ParentSectorId=s2.SectorId";
				id = -1;
				for (DatabaseConnection::RowIterator row = bbn_db->iter(parentSectorQuery); row != bbn_db->end(); ++row) {
					Symbol parentSectorName = row.getCellAsSymbol(0);
					if (std::find(result->begin(), result->end(), parentSectorName) != result->end())
						result->push_back(parentSectorName);
					id = row.getCellAsInt32(1);
				}
			}
		}
	}
	return *result;
}

std::wstring AWAKEActorInfo::getSectorName(SectorId target) {
	std::ostringstream query;
	query << "SELECT SectorName from Sector where SectorId=" << target;
	// LT: replaced getDefaultDb with getDb(target)
	DatabaseConnection_ptr bbn_db(AWAKEDB::getDb(target));
	for (DatabaseConnection::RowIterator row = bbn_db->iter(query); row != bbn_db->end(); ++row) {
		return row.getCellAsWString(0);
	}
	return L"UNKNOWN-SECTOR";
}

std::wstring AWAKEActorInfo::getActorName(ActorId target) {
	if (target.isNull()) return L"UNKNOWN-ACTOR";
	return _cache->getActorRow(target).canonicalName;
}

Symbol AWAKEActorInfo::getIsoCodeForActor(ActorId target) {
	if (target.isNull()) return Symbol(L"unknown");
	return _cache->getActorRow(target).isoCode;
}

Symbol AWAKEActorInfo::getEntityTypeForActor(ActorId target) {
	if (target.isNull()) return Symbol();
	return _cache->getActorRow(target).entityType;
}

bool AWAKEActorInfo::isACountry(ActorId target) {
	if (target.isNull()) return false;
	return _cache->getActorRow(target).entityType == EntityType::getGPEType().getName() &&
		   _cache->getActorRow(target).entitySubtype == Symbol(L"Nation");
}

bool AWAKEActorInfo::isAnIndividual(ActorId target) {
	if (target.isNull()) return false;
	return _cache->getActorRow(target).entityType ==  EntityType::getPERType().getName();
}

bool AWAKEActorInfo::isAGPE(ActorId target) {
	if (target.isNull()) return false;
	return _cache->getActorRow(target).entityType == EntityType::getGPEType().getName();
}

bool AWAKEActorInfo::isAPerson(ActorId target) {
	if (target.isNull()) return false;
	return _cache->getActorRow(target).entityType == EntityType::getPERType().getName();
}

bool AWAKEActorInfo::isAnOrganization(ActorId target) {
	if (target.isNull()) return false;
	return _cache->getActorRow(target).entityType == EntityType::getORGType().getName();
}

bool AWAKEActorInfo::isAFacility(ActorId target) {
	if (target.isNull()) return false;
	return _cache->getActorRow(target).entityType == EntityType::getFACType().getName();
}

void AWAKEActorInfo::loadActorPatterns() {
	double min_confidence = ParamReader::getOptionalFloatParamWithDefaultValue("minimum_actor_string_confidence", 0);
	//LT: changed to be looping over all dbs instead of just getting default
	//DatabaseConnection_ptr bbn_db = AWAKEDB::getDefaultDb();
	DatabaseConnectionMap dbs = AWAKEDB::getNamedDbs();
	for (DatabaseConnectionMap::iterator i = dbs.begin(); i != dbs.end(); ++i) {
		double min_per_actor_importance_score =
			ParamReader::getOptionalFloatParamWithDefaultValue("min_per_actor_importance_score", 0.0);
		double min_org_actor_importance_score =
			ParamReader::getOptionalFloatParamWithDefaultValue("min_org_actor_importance_score", 0.0);

		std::vector<std::string> allowedSources;
		std::vector<std::string> disallowedSources;
		if (ParamReader::hasParam("actor_string_allowed_sources") && ParamReader::hasParam("actor_string_disallowed_sources")) {
			throw UnexpectedInputException("AWAKEActorInfo::loadActorPatterns", "Only one of actor_string_allowed_sources or actor_string_disallowed_sources may be specified");
		}
		if (ParamReader::hasParam("actor_string_allowed_sources")) {
			allowedSources = ParamReader::getStringVectorParam("actor_string_allowed_sources");
		} else if (ParamReader::hasParam("actor_string_disallowed_sources")) {
			disallowedSources = ParamReader::getStringVectorParam("actor_string_disallowed_sources");
		}

		std::ostringstream query;
		query << "SELECT a.ActorId, a.EntityType, s.ActorStringId, s.String, s.Acronym, s.Confidence"
			  << ", s.RequiresContext, a.ImportanceScore";

		// Maintain backwards compatibility with databases that don't have an ActorStringSource table
		if (allowedSources.size() != 0 || disallowedSources.size() != 0)
			query << ", asr.OriginalSourceElement";

		query << " FROM Actor a, ActorString s";

		if (allowedSources.size() != 0 || disallowedSources.size() != 0)
			query << ", ActorStringSource asr";

		query << " WHERE a.ActorId = s.ActorId";

		if (allowedSources.size() != 0 || disallowedSources.size() != 0)
			query << " AND s.ActorStringId = asr.ActorStringId";

		if (min_confidence > 0)
			query << " AND s.Confidence > " << min_confidence;

		bool printed = false;
		for (DatabaseConnection::RowIterator row = i->second->iter(query); row != i->second->end(); ++row) {
			ActorPattern *ap = _new ActorPattern();
			Symbol entityType = Symbol(row.getCellAsWString(1));
			double importance_score = row.getCellAsDouble(7);

			if (allowedSources.size() != 0 || disallowedSources.size() != 0) {
				std::string source = row.getCellAsString(8);
				std::vector<std::string> stringSources;
				boost::split(stringSources, source, boost::is_any_of(" ,"));
				if (allowedSources.size() != 0) {
					bool found_match = false;
					BOOST_FOREACH(std::string as, allowedSources) {
						if (std::find(stringSources.begin(), stringSources.end(), as) != stringSources.end()) {
							found_match = true;
						}
					}
					if (!found_match)
						continue;
				}
				if (disallowedSources.size() != 0) {
					bool found_match = false;
					BOOST_FOREACH(std::string as, disallowedSources) {
						if (std::find(stringSources.begin(), stringSources.end(), as) != stringSources.end()) {
							found_match = true;
						}
					}
					if (found_match)
						continue;
				}
			}
			if (min_per_actor_importance_score > 0.0 && entityType == PER && importance_score < min_per_actor_importance_score)
				continue;
			if (min_org_actor_importance_score > 0.0 && entityType == ORG && importance_score < min_org_actor_importance_score)
				continue;
			// LT: replaced _actorDBName with i->first
			ap->actor_id = ActorId(row.getCellAsInt32(0), i->first);
			ap->entityTypeSymbol = entityType;
			ap->pattern_id = ActorPatternId(row.getCellAsInt32(2), i->first);

			std::wstring str = row.getCellAsWString(3);
			std::vector<std::wstring> words;
			boost::split(words, str, boost::is_any_of(" "));
			ap->pattern = std::vector<Symbol>();
			BOOST_FOREACH(std::wstring word, words) {
				ap->pattern.push_back(Symbol(word));
				std::transform(word.begin(), word.end(), word.begin(), towlower);
				ap->lcPattern.push_back(Symbol(word));
			}

			int acronym_int = row.getCellAsBool(4);
			if (acronym_int)
				ap->acronym = true;
			else
				ap->acronym = false;

			std::string confidenceStr = row.getCellAsString(5);
			ap->confidence = (float)::atof(confidenceStr.c_str());

			ap->requires_context = false;
			int context_int = row.getCellAsBool(6);
			if (context_int)
				ap->requires_context = true;
			else
				ap->requires_context = false;

			ap->lcString = ActorPattern::getNameFromSymbolList(ap->lcPattern);

			_patterns.push_back(ap);
		}
	}
}

ActorId AWAKEActorInfo::getActorIdForGeonameId(std::wstring &geonameid) {
	std::map<std::wstring, ActorId>::iterator iter = _geonameIdToActorIdCache.find(geonameid);
	if (iter != _geonameIdToActorIdCache.end())
		return (*iter).second;

	std::ostringstream query;
	query << "SELECT ActorId"
		  << " FROM Actor WHERE GeonameId=" << DatabaseConnection::quote(geonameid);
	//LT: changed to be looping over all dbs instead of just getting default
	//DatabaseConnection_ptr bbn_db(AWAKEDB::getDefaultDb());
	DatabaseConnectionMap dbs = AWAKEDB::getNamedDbs();
	for (DatabaseConnectionMap::iterator i = dbs.begin(); i != dbs.end(); ++i) {
		for (DatabaseConnection::RowIterator row = i->second->iter(query); row != i->second->end(); ++row) {
			// LT: replaced _actorDBName with i->first
			ActorId actor_id(row.getCellAsInt32(0), i->first);
			_geonameIdToActorIdCache[geonameid] = actor_id;
			return actor_id;
		}
	}
	_geonameIdToActorIdCache[geonameid] = ActorId();
	return ActorId();
}

double AWAKEActorInfo::getImportanceScoreForActor(ActorId target) {
	if (target.isNull()) return 0.0;
	return _cache->getActorRow(target).importance_score;
}

ActorId AWAKEActorInfo::getActorByName(std::wstring name) {
	// Inefficient, but is currently only used to get the United States actor
	std::ostringstream query;
	query << "SELECT ActorId"
		  << " FROM Actor WHERE CanonicalName=" << DatabaseConnection::quote(name);
	//LT: changed to be looping over all dbs instead of just getting default
	//DatabaseConnection_ptr bbn_db(AWAKEDB::getDefaultDb());
	DatabaseConnectionMap dbs = AWAKEDB::getNamedDbs();
	for (DatabaseConnectionMap::iterator i = dbs.begin(); i != dbs.end(); ++i) {
		for (DatabaseConnection::RowIterator row = i->second->iter(query); row != i->second->end(); ++row) {
			// LT: replaced _actorDBName with i->first
			ActorId actor_id(row.getCellAsInt32(0), i->first);
			return actor_id;
		}
	}
	return ActorId();
}

AgentId AWAKEActorInfo::getAgentByName(Symbol name) {
	return _cache->getAgentByName(name);
}

std::wstring AWAKEActorInfo::getAgentName(AgentId target) {
	return _cache->getAgentRow(target).name;
}

bool AWAKEActorInfo::isRestrictedToCountryActors(AgentId target) {
	// Always false for AWAKE agents
	return false;
}

std::vector<Symbol> AWAKEActorInfo::getAssociatedSectorCodes(AgentId target) {
	return _cache->getAgentRow(target).associatedSectorCodes;
}

AgentId AWAKEActorInfo::getDefaultPersonAgentId() {
	return _cache->getDefaultPersonAgentId();
}

Symbol AWAKEActorInfo::getDefaultPersonAgentCode() {
	return _cache->getDefaultPersonAgentCode();
}

bool AWAKEActorInfo::isCountryActorName(std::wstring name) {
	std::transform(name.begin(), name.end(), name.begin(), towlower);

	if (_countryActorNameCache.find(name) != _countryActorNameCache.end())
		return _countryActorNameCache[name];

	std::ostringstream query;
	query << "SELECT ActorId"
		  << " FROM Actor WHERE EntitySubtype='Nation'"
		  << " AND lower(CanonicalName)=" << DatabaseConnection::quote(name);

	//LT: changed to be looping over all dbs instead of just getting default
	//DatabaseConnection_ptr bbn_db(AWAKEDB::getDefaultDb());
	DatabaseConnectionMap dbs = AWAKEDB::getNamedDbs();
	for (DatabaseConnectionMap::iterator i = dbs.begin(); i != dbs.end(); ++i) {
		for (DatabaseConnection::RowIterator row = i->second->iter(query); row != i->second->end(); ++row) {
			_countryActorNameCache[name] = true;
			return true;
		}
	}

	_countryActorNameCache[name] = false;
	return false;
}

std::vector<ActorPattern *>& AWAKEActorInfo::getPatterns() {
	return _patterns;
}

ext: cpp

max_stars_repo_name: BBN-E/serif

max_issues_repo_name: BBN-E/serif

max_forks_repo_name: BBN-E/serif


text:
// Copyright (c) 2019 The Unit-e developers
// Distributed under the MIT software license, see the accompanying
// file COPYING or http://www.opensource.org/licenses/mit-license.php.
#include <coins.h>

#include <test/test_unite.h>
#include <boost/test/unit_test.hpp>

BOOST_FIXTURE_TEST_SUITE(coins_view_db_tests, BasicTestingSetup)

BOOST_AUTO_TEST_CASE(ccoins_view_cache_clear_coins)
{
  auto view_db = MakeUnique<CCoinsViewDB>(0, true, true);
  auto view_backend = MakeUnique<CCoinsViewCache>(view_db.get());  // extra layer
  auto view_cache = MakeUnique<CCoinsViewCache>(view_backend.get());

  CCoinsMap coins;
  for (size_t i = 0; i < 5; ++i) {
    COutPoint point(uint256S("aa"), static_cast<uint32_t>(i));
    Coin coin(CTxOut(1, CScript()), 1, TxType::REGULAR);
    CCoinsCacheEntry entry(std::move(coin));
    entry.flags |= CCoinsCacheEntry::DIRTY;
    coins.emplace(point, std::move(entry));
  }
  view_db->BatchWrite(coins, uint256S("aa"), snapshot::SnapshotHash());

  size_t total = 0;
  std::unique_ptr<CCoinsViewCursor> cursor(view_db->Cursor());
  while (cursor->Valid()) {
    ++total;
    cursor->Next();
  }
  BOOST_CHECK_EQUAL(total, 5); // sanity check

  view_cache->ClearCoins();
  std::unique_ptr<CCoinsViewCursor> empty_cursor(view_db->Cursor());
  BOOST_CHECK(!empty_cursor->Valid());
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: frolosofsky/unit-e

max_issues_repo_name: frolosofsky/unit-e

max_forks_repo_name: frolosofsky/unit-e


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

// Standard includes
#include <fstream>
#include <iomanip>
#include <iostream>

// Third party includes
#include <boost/algorithm/string/trim.hpp>
#include <boost/lexical_cast.hpp>

// VOTCA includes
#include <votca/tools/getline.h>
#include <votca/tools/rangeparser.h>
#include <votca/tools/table.h>
#include <votca/tools/tokenizer.h>

// Local VOTCA includes
#include "votca/csg/imcio.h"

namespace votca {
namespace csg {

using group_matrix = Eigen::MatrixXd;
using namespace std;

void imcio_write_dS(const string &file, const tools::Table &dS,
                    const std::list<Index> *list) {
  // write the dS
  ofstream out_dS;
  out_dS.open(file);
  out_dS << setprecision(8);
  if (!out_dS) {
    throw runtime_error(string("error, cannot open file ") + file);
  }

  if (list == nullptr) {
    for (Index i = 0; i < dS.size(); ++i) {
      out_dS << dS.x(i) << " " << dS.y(i) << endl;
    }
  } else {
    for (Index i : *list) {
      out_dS << dS.x(i) << " " << dS.y(i) << endl;
    }
  }

  out_dS.close();
  cout << "written " << file << endl;
}

void imcio_write_matrix(const string &file, const Eigen::MatrixXd &gmc,
                        const std::list<Index> *list) {
  ofstream out_A;
  out_A.open(file);
  out_A << setprecision(8);

  if (!out_A) {
    throw runtime_error(string("error, cannot open file ") + file);
  }

  if (list == nullptr) {
    for (Index i = 0; i < gmc.rows(); ++i) {
      for (Index j = 0; j < gmc.cols(); ++j) {
        out_A << gmc(i, j) << " ";
      }
      out_A << endl;
    }
  } else {
    for (Index i : *list) {
      for (Index j : *list) {
        out_A << gmc(i, j) << " ";
      }
      out_A << endl;
    }
  }
  out_A.close();
  cout << "written " << file << endl;
}

void imcio_write_index(
    const string &file,
    const std::vector<std::pair<std::string, tools::RangeParser> > &ranges) {
  // write the index
  ofstream out_idx;
  out_idx.open(file);

  if (!out_idx) {
    throw runtime_error(string("error, cannot open file ") + file);
  }

  for (const auto &range : ranges) {
    out_idx << range.first << " " << range.second << endl;
  }
  out_idx.close();
  cout << "written " << file << endl;
}

Eigen::MatrixXd imcio_read_matrix(const std::string &filename) {
  std::ifstream intt;
  intt.open(filename);
  if (!intt) {
    throw std::runtime_error(std::string("error, cannot open file ") +
                             filename);
  }

  std::string line;
  std::vector<double> result;
  Index numrows = 0;
  size_t numcols = 0;
  while (tools::getline(intt, line)) {
    if (line[0] == '#') {
      continue;
    }
    votca::tools::Tokenizer tok(line, " \t");
    std::vector<std::string> tokens = tok.ToVector();
    if (numrows == 0) {
      numcols = tokens.size();
    } else if (numcols != tokens.size()) {
      throw std::runtime_error(
          "Matrix has not the same number of entries in each row.");
    }
    for (const std::string &s : tokens) {
      result.push_back(std::stod(s));
    }
    numrows++;
  }
  intt.close();

  return Eigen::Map<Eigen::MatrixXd>(result.data(), numrows, numcols);
}

std::vector<std::pair<std::string, tools::RangeParser> > imcio_read_index(
    const string &filename) {
  ifstream in;
  in.open(filename);
  if (!in) {
    throw runtime_error(string("error, cannot open file ") + filename);
  }

  std::vector<std::pair<std::string, tools::RangeParser> > indeces;
  string line;
  // read till the first data line
  while (tools::getline(in, line)) {
    // remove comments and xmgrace stuff
    line = line.substr(0, line.find("#"));
    line = line.substr(0, line.find("@"));

    boost::trim(line);
    size_t found = line.find(" ");
    if (found == string::npos) {
      throw runtime_error(string("wrong format in ") + filename);
    }

    string name = line.substr(0, found);
    string range = line.substr(found);
    tools::RangeParser rp;
    rp.Parse(range);
    indeces.push_back(std::pair<std::string, tools::RangeParser>(name, rp));
  }
  in.close();
  return indeces;
}

}  // namespace csg
}  // namespace votca

ext: cc

max_stars_repo_name: BuildJet/csg

max_issues_repo_name: BuildJet/csg

max_forks_repo_name: BuildJet/csg


text:
////////////////////////////////////////////////////////////////////////////////
//  Copyright (c) 2010 Bryce Lelbach
//
//  Distributed under the Boost Software License, Version 1.0. (See accompanying
//  file BOOST_LICENSE_1_0.rst or copy at http://www.boost.org/LICENSE_1_0.txt)
////////////////////////////////////////////////////////////////////////////////

#include <phxpr/config.hpp>

#include <iterator>

#include <boost/assert.hpp>

#include <prana/utree/predicates.hpp>

#include <phxpr/evaluator.hpp>
#include <phxpr/predicate.hpp>
#include <phxpr/exception.hpp>

using boost::fusion::at_c;

namespace phxpr {

evaluator::result_type
evaluator::execute (evaluator::range_type const& range) {
  // {{{
  typedef evaluator::range_type::const_iterator iterator;

  iterator it = range.begin(), end = range.end();

  if (it == end)
    BOOST_THROW_EXCEPTION(invalid_operator_expression(utree(range)));

  utree f;

  // Evaluate the operator expresion.
  if (prana::is_utree_container(*it)) {
    f = evaluate(*it, *this);

    if (prana::recursive_which(f) != utree_type::function_type)
      BOOST_THROW_EXCEPTION(procedure_call_or_macro_use_expected(*it));
  }

  // Evaluate the named operator.
  else if (prana::recursive_which(*it) == utree_type::symbol_type) {
    // TODO: syntax checks
    if (*it == utree(spirit::utf8_symbol_type("lambda"))) {
      iterator formals = ++it;
      iterator body = ++it;
      f = make_lambda_expression(*formals, range_type(body, end));
      phxpr::shared_ptr<runtime_environment> new_env
        = phxpr::make_shared<runtime_environment>();
      return new_env->invoke(f);
    }

    else if (*it == utree(spirit::utf8_symbol_type("variable"))) {
      iterator identifier = ++it;
      iterator value = ++it;

      if (++it != end)
        BOOST_THROW_EXCEPTION(invalid_variable_definition
          (utree(range_type(range.begin(), value))));

      // Module-level variable definition.
      BOOST_ASSERT(frame == 0);
      return make_module_level_variable(*identifier, *value);
    }

    // TODO: syntax checks
    else if (*it == utree(spirit::utf8_symbol_type("if"))) {
      iterator test = ++it;
      iterator then = ++it;
      iterator else_ = ++it;

      bool r = predicate(evaluate(*test, *this));

      if (r)
        return evaluate(*then, *this);

      // (if <false> <then>) returns unspecified (aka invalid)
      if (else_ == end)
        return utree();

      return evaluate(*else_, *this);
    }

    else {
      phxpr::shared_ptr<utree> pf = variables->lookup(*it);

      if (!pf)
        BOOST_THROW_EXCEPTION(identifier_not_found(*it));

      f = *pf;
    }
  }

  else
    BOOST_THROW_EXCEPTION(procedure_call_or_macro_use_expected(*it));

  BOOST_ASSERT(unsigned(f.tag()) <= global_procedure_table->size());

  // Load the lazy argument's signature from the gpt.
  signature const& sig = (*global_procedure_table)[f.tag()];

  if ((at_c<3>(sig) == function_type::thunk) ||
      (at_c<3>(sig) == function_type::lambda))
  {
    phxpr::shared_ptr<runtime_environment> new_env
      = phxpr::make_shared<runtime_environment>();
    f = new_env->invoke(f);
  }

  // Invoke nullary procedures.
  if (++it == end) {
    phxpr::shared_ptr<runtime_environment> new_env
      = phxpr::make_shared<runtime_environment>();
    return new_env->invoke(f);
  }

  // Invoke non-nullary procedures.
  const displacement env_size = at_c<4>(sig) + std::distance(it, end);
  phxpr::shared_ptr<runtime_environment::upvalue_type> storage
    = phxpr::make_shared<runtime_environment::upvalue_type>
      (env_size, sheol::fill);

  // Evaluate the arguments.
  for (std::size_t i = 0; it != end; ++it)
    (*storage)[i++] = evaluate(*it, *this);

  phxpr::shared_ptr<runtime_environment> new_env
    = phxpr::make_shared<runtime_environment>(storage);
  return new_env->invoke(f);
} // }}}

} // phxpr

max_stars_repo_name: brycelelbach/prana

max_issues_repo_name: brycelelbach/prana

max_forks_repo_name: brycelelbach/prana


text:
/// @file
///
/// XXX Notes on program XXX
///
/// @copyright (c) 2011 CSIRO
/// Australia Telescope National Facility (ATNF)
/// Commonwealth Scientific and Industrial Research Organisation (CSIRO)
/// PO Box 76, Epping NSW 1710, Australia
/// atnf-enquiries@csiro.au
///
/// This file is part of the ASKAP software distribution.
///
/// The ASKAP software distribution is free software: you can redistribute it
/// and/or modify it under the terms of the GNU General Public License as
/// published by the Free Software Foundation; either version 2 of the License,
/// or (at your option) any later version.
///
/// This program is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with this program; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
///
/// @author XXX XXX <XXX.XXX@csiro.au>
///
#include <askap_analysis.h>
#include <patternmatching/CatalogueMatcher.h>

#include <askap/AskapLogging.h>
#include <askap/AskapError.h>

#include <patternmatching/Triangle.h>
#include <patternmatching/Point.h>
#include <patternmatching/PointCatalogue.h>
#include <patternmatching/MatchingUtilities.h>
#include <casainterface/CasaInterface.h>

#include <Common/ParameterSet.h>
#include <boost/shared_ptr.hpp>
#include <casacore/images/Images/ImageInterface.h>
#include <casacore/images/Images/ImageOpener.h>
#include <casacore/images/Images/FITSImage.h>
#include <casacore/images/Images/MIRIADImage.h>
#include <casacore/casa/Arrays/Vector.h>
#include <casacore/coordinates/Coordinates/CoordinateSystem.h>
#include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
#include <casacore/casa/Quanta.h>

#include <vector>

ASKAP_LOGGER(logger, ".cataloguematching");

namespace askap {

namespace analysis {

namespace matching {

CatalogueMatcher::CatalogueMatcher()
{
}

CatalogueMatcher::CatalogueMatcher(const LOFAR::ParameterSet& parset)
{
    itsReferenceImage = parset.getString("referenceImage", "");

    LOFAR::ParameterSet subset = parset.makeSubset("source.");
    subset.add("referenceImage", itsReferenceImage);
    itsSrcCatalogue = PointCatalogue(subset);

    subset = parset.makeSubset("reference.");
    subset.add("referenceImage", itsReferenceImage);
    itsRefCatalogue = PointCatalogue(subset);

    itsMatchFile = parset.getString("matchFile", "matches.txt");
    itsMissFile = parset.getString("missFile", "misses.txt");
    itsPositionUnits = casa::Unit(parset.getString("positionUnits", "deg"));
    if (!parset.isDefined("epsilon")) {
        ASKAPTHROW(AskapError, "The epsilon parameter must be provided.");
    }
    std::string epsilonString = parset.getString("epsilon");
    casa::Quantity q;
    casa::Quantity::read(q, epsilonString);
    itsEpsilon = q.getValue(itsPositionUnits);
    itsEpsilonUnits = q.getUnit();
    this->convertEpsilon();
    ASKAPLOG_DEBUG_STR(logger, "Requested epsilon value was " << epsilonString <<
                       ", which is " << itsEpsilon << " " << itsPositionUnits.getName());
    // itsEpsilon = parset.getDouble("epsilon");
    if (itsEpsilon < 0) {
        ASKAPTHROW(AskapError, "The epsilon parameter must be positive.");
    }
    itsMeanDx = itsMeanDy = 0.;
    itsSourceSummaryFile = parset.getString("srcSummaryFile", "match-summary-sources.txt");
    itsReferenceSummaryFile = parset.getString("refSummaryFile", "match-summary-reference.txt");
}

void CatalogueMatcher::convertEpsilon()
{
    if (itsReferenceImage != "") {
        const boost::shared_ptr<ImageInterface<Float> > imagePtr =
            analysisutilities::openImage(itsReferenceImage);

        int dirCooNum = imagePtr->coordinates().findCoordinate(casa::Coordinate::DIRECTION);
        casa::DirectionCoordinate dirCoo = imagePtr->coordinates().directionCoordinate(dirCooNum);

        ASKAPLOG_DEBUG_STR(logger, "Converting epsilon from " << itsEpsilon <<
                           " " << itsPositionUnits.getName());
        ASKAPASSERT(dirCoo.worldAxisUnits()[0] == dirCoo.worldAxisUnits()[1]);

        casa::Quantity eps(itsEpsilon, itsPositionUnits);
        double epsInWorldUnits = eps.getValue(dirCoo.worldAxisUnits()[0]);
        double pixScale = sqrt(fabs(dirCoo.increment()[0] *  dirCoo.increment()[1]));
        itsEpsilon = epsInWorldUnits / pixScale;
        itsEpsilonUnits.setValue(itsEpsilonUnits.getValue() / pixScale);
        ASKAPLOG_DEBUG_STR(logger, "Now have epsilon = " << itsEpsilon);
    }

}

bool CatalogueMatcher::read()
{

    bool filesOK = true;

    if (!itsSrcCatalogue.read()) {
        ASKAPLOG_FATAL_STR(logger, "Could not read Source Catalogue");
        return false;
    }
    if (!itsRefCatalogue.read()) {
        ASKAPLOG_FATAL_STR(logger, "Could not read Reference Catalogue");
        return false;
    }

    filesOK = itsSrcCatalogue.pointList().size() > 0 && itsRefCatalogue.pointList().size() > 0;
    if (itsSrcCatalogue.pointList().size() == 0) {
        ASKAPLOG_ERROR_STR(logger, "Could not read source catalogue from " <<
                           itsSrcCatalogue.filename());
    }
    if (itsRefCatalogue.pointList().size() == 0) {
        ASKAPLOG_ERROR_STR(logger, "Could not read source catalogue from " <<
                           itsRefCatalogue.filename());
    }
    if (filesOK) {
        ASKAPLOG_INFO_STR(logger, "Size of source pixel list = " <<
                          itsSrcCatalogue.pointList().size() <<
                          " and triangle list = " <<
                          itsSrcCatalogue.triangleList().size());
        ASKAPLOG_INFO_STR(logger, "Size of reference pixel list = " <<
                          itsRefCatalogue.pointList().size() <<
                          " and triangle list = " <<
                          itsRefCatalogue.triangleList().size());
    }

    // if(!itsRefCatalogue.crudeMatch(itsSrcCatalogue.fullPointList(), itsEpsilon))
    //      ASKAPLOG_WARN_STR(logger, "Crude matching failed! Using full reference point list");

    return filesOK;
}

//**************************************************************//

void CatalogueMatcher::findMatches()
{

    itsMatchingTriList = matchLists(itsSrcCatalogue.triangleList(),
                                    itsRefCatalogue.triangleList(),
                                    itsEpsilon);

    trimTriList(itsMatchingTriList);
    ASKAPLOG_INFO_STR(logger, "Found " << itsMatchingTriList.size() << " matches");
    itsNumInitialMatches = 0;

    if (itsMatchingTriList.size() > 0) {
        itsMatchingPixList = vote(itsMatchingTriList);
        itsNumInitialMatches = itsMatchingPixList.size();
        ASKAPLOG_INFO_STR(logger, "After voting, have found " <<
                          itsMatchingPixList.size() << " matching points");

        itsSenseMatch = (itsMatchingTriList[0].first.isClockwise() ==
                         itsMatchingTriList[0].second.isClockwise());

        if (itsSenseMatch) {
            ASKAPLOG_INFO_STR(logger, "The two lists have the same sense.");
        } else {
            ASKAPLOG_INFO_STR(logger, "The two lists have the opposite sense.");
        }
    }
}

//**************************************************************//

void CatalogueMatcher::zeroOffsetMatch()
{
    size_t srcSize = itsSrcCatalogue.pointList().size();
    size_t refSize = itsRefCatalogue.pointList().size();
    ASKAPLOG_DEBUG_STR(logger, "Performing zero-offset match of lists of size " <<
                       srcSize << " and " << refSize);

    std::vector<bool> srcMatched(srcSize, false);
    std::vector<bool> refMatched(refSize, false);
    int nmatch = 0;

    std::sort(itsSrcCatalogue.pointList().begin(), itsSrcCatalogue.pointList().end());
    std::sort(itsRefCatalogue.pointList().begin(), itsRefCatalogue.pointList().end());

    for (size_t s = 0; s < itsSrcCatalogue.pointList().size(); s++) {

        for (size_t r = 0; r < itsRefCatalogue.pointList().size() && !srcMatched[s]; r++) {

            if (!refMatched[r]) {
                if (itsSrcCatalogue.pointList()[s].sep(itsRefCatalogue.pointList()[r]) <
                        itsEpsilon) {

                    itsMatchingPixList.push_back(
                        std::pair<Point, Point>(itsSrcCatalogue.pointList()[s],
                                                itsRefCatalogue.pointList()[r]));

                    refMatched[r] = true;
                    srcMatched[s] = true;
                    nmatch++;
                }
            }
        }
    }

    itsSenseMatch = true;

    ASKAPLOG_DEBUG_STR(logger, "Matched " << nmatch << " pairs of points");

}

//**************************************************************//

void CatalogueMatcher::findOffsets()
{
    std::vector<double> dx(itsMatchingPixList.size(), 0.);
    std::vector<double> dy(itsMatchingPixList.size(), 0.);

    double sensescale = itsSenseMatch ? -1. : 1.;
    for (size_t i = 0; i < itsMatchingPixList.size(); i++) {
        dx[i] = itsMatchingPixList[i].first.x() - itsMatchingPixList[i].second.x();
        dy[i] = itsMatchingPixList[i].first.y() + sensescale * itsMatchingPixList[i].second.y();
    }

    itsMeanDx = itsMeanDy = 0.;
    for (size_t i = 0; i < itsMatchingPixList.size(); i++) {
        itsMeanDx += dx[i];
        itsMeanDy += dy[i];
    }
    itsMeanDx /= double(itsMatchingPixList.size());
    itsMeanDy /= double(itsMatchingPixList.size());

    // Don't need these elsewhere, so make local only
    double rmsDx = 0., rmsDy = 0.;
    for (size_t i = 0; i < itsMatchingPixList.size(); i++) {
        rmsDx += (dx[i] - itsMeanDx) * (dx[i] - itsMeanDx);
        rmsDy += (dy[i] - itsMeanDy) * (dy[i] - itsMeanDy);
    }
    rmsDx = sqrt(rmsDx / (double(itsMatchingPixList.size() - 1)));
    rmsDy = sqrt(rmsDy / (double(itsMatchingPixList.size() - 1)));

    std::stringstream ss;
    ss << "Offsets between the two are dx = "
       << casa::Quantity(itsMeanDx, itsPositionUnits).getValue(itsEpsilonUnits)
       << " +- " << casa::Quantity(rmsDx, itsPositionUnits).getValue(itsEpsilonUnits)
       << " dy = "
       << casa::Quantity(itsMeanDy, itsPositionUnits).getValue(itsEpsilonUnits)
       << " +- " << casa::Quantity(rmsDy, itsPositionUnits).getValue(itsEpsilonUnits);
    ASKAPLOG_INFO_STR(logger, ss.str());
}

//**************************************************************//

void CatalogueMatcher::addNewMatches()
{

    if (itsNumInitialMatches > 0) {

        this->rejectMultipleMatches();
        const float matchRadius = 3.;
        std::vector<Point>::iterator src, ref;
        std::vector<std::pair<Point, Point> >::iterator match;

        for (src = itsSrcCatalogue.fullPointList().begin();
                src < itsSrcCatalogue.fullPointList().end();
                src++) {

            bool isMatch = false;
            match = itsMatchingPixList.begin();

            for (; match < itsMatchingPixList.end() && !isMatch; match++) {
                isMatch = (src->ID() == match->first.ID());
            }

            if (!isMatch) {
                float minOffset = 0.;
                int minRef = -1;

                for (ref = itsRefCatalogue.fullPointList().begin();
                        ref < itsRefCatalogue.fullPointList().end();
                        ref++) {

                    float offset = hypot(src->x() - ref->x() - itsMeanDx,
                                         src->y() - ref->y() - itsMeanDy);

                    if (offset < matchRadius * itsEpsilon) {
                        if ((minRef == -1) || (offset < minOffset)) {
                            minOffset = offset;
                            minRef = int(ref - itsRefCatalogue.fullPointList().begin());
                        }
                    }
                }

                if (minRef >= 0) { // there was a match within errors
                    ref = itsRefCatalogue.fullPointList().begin() + minRef;
                    std::pair<Point, Point> newMatch(*src, *ref);
                    itsMatchingPixList.push_back(newMatch);
                }
            }
        }

        this->rejectMultipleMatches();

        ASKAPLOG_INFO_STR(logger, "Total number of matches = " << itsMatchingPixList.size());
    }
}

//**************************************************************//

void CatalogueMatcher::rejectMultipleMatches()
{

    if (itsMatchingPixList.size() < 2) return;

    std::vector<std::pair<Point, Point> >::iterator alice, bob;
    alice = itsMatchingPixList.begin();

    while (alice < itsMatchingPixList.end() - 1) {
        bool aliceGone = false;
        bob = alice + 1;

        while (bob < itsMatchingPixList.end() && !aliceGone) {
            bool bobGone = false;
            if (alice->second.ID() == bob->second.ID()) {
                // alice & bob have the same reference source
                double df_alice = alice->first.flux() - alice->second.flux();
                double df_bob   = bob->first.flux() - bob->second.flux();

                if (fabs(df_alice) < fabs(df_bob)) {
                    itsMatchingPixList.erase(bob);
                    bobGone = true;
                } else {
                    itsMatchingPixList.erase(alice);
                    aliceGone = true;
                }
            }
            if (!bobGone) bob++;
        }

        if (!aliceGone) alice++;
    }
}

//**************************************************************//

void CatalogueMatcher::outputMatches()
{
    std::ofstream fout(itsMatchFile.c_str());
    if (fout.is_open()) {
        std::vector<std::pair<Point, Point> >::iterator match;
        int prec = 3;
        size_t width = 0;
        for (match = itsMatchingPixList.begin();
                match < itsMatchingPixList.end();
                match++) {
            prec = std::max(prec, int(ceil(log10(1. / match->first.flux()))) + 1);
            width = std::max(width, match->first.ID().size());
            width = std::max(width, match->second.ID().size());
        }

        fout.setf(std::ios::fixed);
        unsigned int ct = 0;
        char matchType;

        for (match = itsMatchingPixList.begin();
                match < itsMatchingPixList.end();
                match++) {

            if (ct++ < itsNumInitialMatches) matchType = '1';
            else matchType = '2';

            casa::Quantity sep(match->first.sep(match->second), itsPositionUnits);

            fout << std::setw(3) << matchType << " "
                 << std::setw(width) << match->first.ID() << " "
                 << std::setw(width) << match->second.ID() << " "
                 << std::setw(8)  << std::setprecision(6)
                 << sep.getValue(itsEpsilonUnits) << "\n";

        }

        fout.close();
    } else {
        ASKAPLOG_ERROR_STR(logger, "Could not open match file " << itsMatchFile);
    }

}

//**************************************************************//

void CatalogueMatcher::outputMisses()
{
    std::ofstream fout(itsMissFile.c_str());
    if (fout.is_open()) {
        fout.setf(std::ios::fixed);
        std::vector<Point>::iterator pt;
        std::vector<std::pair<Point, Point> >::iterator match;

        size_t width = 0;
        for (pt = itsRefCatalogue.fullPointList().begin();
                pt < itsRefCatalogue.fullPointList().end();
                pt++) {
            width = std::max(width, pt->ID().size());
        }

        for (pt = itsRefCatalogue.fullPointList().begin();
                pt < itsRefCatalogue.fullPointList().end();
                pt++) {

            bool isMatch = false;

            for (match = itsMatchingPixList.begin();
                    match < itsMatchingPixList.end() && !isMatch;
                    match++) {
                isMatch = (pt->ID() == match->second.ID());
            }

            if (!isMatch) {
                fout << "R "
                     << std::setw(width) << pt->ID() << " "
                     << std::setw(10) << std::setprecision(3) << pt->x()  << " "
                     << std::setw(10) << std::setprecision(3) << pt->y() << " "
                     << std::setw(10) << std::setprecision(8) << pt->flux()  << "\n";
            }
        }

        width = 0;
        for (pt = itsSrcCatalogue.fullPointList().begin();
                pt < itsSrcCatalogue.fullPointList().end();
                pt++) {
            width = std::max(width, pt->ID().size());
        }

        for (pt = itsSrcCatalogue.fullPointList().begin();
                pt < itsSrcCatalogue.fullPointList().end();
                pt++) {

            bool isMatch = false;

            for (match = itsMatchingPixList.begin();
                    match < itsMatchingPixList.end() && !isMatch;
                    match++) {
                isMatch = (pt->ID() == match->first.ID());
            }

            if (!isMatch) {
                fout << "S "
                     << std::setw(width) << pt->ID() << " "
                     << std::setw(10) << std::setprecision(3) << pt->x()  << " "
                     << std::setw(10) << std::setprecision(3) << pt->y()  << " "
                     << std::setw(10) << std::setprecision(8) << pt->flux() << "\n";
            }
        }
    } else {
        ASKAPLOG_ERROR_STR(logger, "Could not open miss file " << itsMissFile);
    }

}

//**************************************************************//

void CatalogueMatcher::outputSummary()
{

    if (itsSourceSummaryFile != "") {
        this->outputSummary(itsSrcCatalogue, "src", itsSourceSummaryFile);
    }

    if (itsReferenceSummaryFile != "") {
        this->outputSummary(itsRefCatalogue, "ref", itsReferenceSummaryFile);
    }

}

void CatalogueMatcher::outputSummary(PointCatalogue &cat, std::string whichOne, std::string filename)
{
    size_t width = 0;
    std::vector<Point>::iterator pt;
    std::vector<std::pair<Point, Point> >::iterator mpair;
    std::vector<std::pair<Point, Point> >::iterator match;

    if ( (whichOne!="src") && (whichOne!="ref") ){
        ASKAPTHROW(AskapError, "Incorrect value for 'whichOne' in CatalogueMatcher::outputSummary");
    }

    for (match = itsMatchingPixList.begin();
            match < itsMatchingPixList.end();
            match++) {
        width = std::max(width, match->first.ID().size());
        width = std::max(width, match->second.ID().size());
    }

    std::ofstream fout(filename.c_str());
    if (fout.is_open()) {
        for (pt = cat.fullPointList().begin();
                pt < cat.fullPointList().end();
                pt++) {

            std::string noMatch = "---";
            std::string matchID = noMatch;
            for (mpair = itsMatchingPixList.begin();
                 (mpair < itsMatchingPixList.end()) && (matchID==noMatch);
                 mpair++) {
                if ( (whichOne == "src") && (pt->ID() == mpair->first.ID()) ) {
                    matchID = mpair->second.ID();
                }
                else if ( (whichOne == "ref") && (pt->ID() == mpair->second.ID()) ){
                    matchID = mpair->first.ID();
                }
            }
            fout << std::setw(width) << pt->ID() << " "
                 << std::setw(width) << matchID << " "
                 << std::setw(10) << std::setprecision(7) << pt->x()  << " "
                 << std::setw(10) << std::setprecision(7) << pt->y() << " "
                 << std::setw(10) << std::setprecision(8) << pt->flux()  << "\n";
        }
        fout.close();
    } else {
        ASKAPLOG_ERROR_STR(logger, "Could not open summary file " << filename);
    }

}

}
}
}

max_stars_repo_name: rtobar/askapsoft

max_issues_repo_name: ATNF/askapsoft

max_forks_repo_name: ATNF/askapsoft


text:
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Unit Test

// Copyright (c) 2019 Tinko Bartels, Berlin, Germany.

// Contributed and/or modified by Tinko Bartels,
//   as part of Google Summer of Code 2019 program.

// This file was modified by Oracle on 2021.
// Modifications copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#include <geometry_test_common.hpp>

#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/strategy/cartesian/in_circle_robust.hpp>

template <typename P>
void test_all()
{
    typedef bg::strategy::in_circle::in_circle_robust<double, 2> inc2;
    typedef bg::strategy::in_circle::in_circle_robust<double, 1> inc1;
    typedef bg::strategy::in_circle::in_circle_robust<double, 0> inc0;

    P col1(0.0, 0.0), col2(1.0, 0.0), col3(0.0, 1.0);
    P in(0.5,0.5) , on(1.0, 0.0), out(-0.5, -0.5);
    int in2 = inc2::apply(col1, col2, col3, in);
    BOOST_CHECK_GT(in2, 0);
    int in1 = inc1::apply(col1, col2, col3, in);
    BOOST_CHECK_GT(in1, 0);
    int in0 = inc0::apply(col1, col2, col3, in);
    BOOST_CHECK_GT(in0, 0);

    int on2 = inc2::apply(col1, col2, col3, on);
    BOOST_CHECK_EQUAL(on2, 0);
    int on1 = inc1::apply(col1, col2, col3, on);
    BOOST_CHECK_EQUAL(on1, 0);
    int on0 = inc0::apply(col1, col2, col3, on);
    BOOST_CHECK_EQUAL(on0, 0);

    int out2 = inc2::apply(col1, col2, col3, out);
    BOOST_CHECK_GT(0, out2);
    int out1 = inc1::apply(col1, col2, col3, out);
    BOOST_CHECK_GT(0, out1);
    int out0 = inc0::apply(col1, col2, col3, out);
    BOOST_CHECK_GT(0, out0);

    P hard1(0, 0), hard2(1e20, 0), hard3(0, 1e20);
    P inhard(0.5, 0.5);
    int hardr  = inc2::apply(hard1, hard2, hard3, inhard);
    BOOST_CHECK_GT(hardr, 0);
}

int test_main(int, char* [])
{
    test_all<bg::model::d2::point_xy<double> >();
    return 0;
}

max_stars_repo_name: jhypolite/geometry

max_issues_repo_name: jhypolite/geometry

max_forks_repo_name: jhypolite/geometry


text: #include <boost/preprocessor/library.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2009, 2014 Jose Aparicio

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#ifndef quantlib_recursive_loss_model_hpp
#define quantlib_recursive_loss_model_hpp

#include <ql/experimental/credit/constantlosslatentmodel.hpp>
#include <ql/experimental/credit/defaultlossmodel.hpp>

#if defined(__GNUC__) && (((__GNUC__ == 4) && (__GNUC_MINOR__ >= 8)) || (__GNUC__ > 4))
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
#endif
#include <boost/bind.hpp>
#if defined(__GNUC__) && (((__GNUC__ == 4) && (__GNUC_MINOR__ >= 8)) || (__GNUC__ > 4))
#pragma GCC diagnostic pop
#endif
#include <map>
#include <algorithm>

namespace QuantLib {

    /*! Recursive STCDO default loss model for a heterogeneous pool of names.
    The pool names are heterogeneous in their default probabilities, notionals
    and recovery rates. Correlations are given by the latent model.
    The recursive pricing algorithm used here is described in Andersen, Sidenius
    and Basu; "All your hedges in one basket", Risk, November 2003, pages 67-72

        Notice that using copulas other than Gaussian it is only an
        approximation (see remark on p.68).

        \todo Make the loss unit equal to some small fraction depending on the
        portfolio loss weights (notionals and recoveries). As it is now this
        is ok for pricing but not for risk metrics. See the discussion in O'Kane
        18.3.2
        \todo Intengrands should all use the inverted probabilities for
        performance instead of calling the copula inversion with the same vals.
    */
    template<class copulaPolicy>
    class RecursiveLossModel : public DefaultLossModel {
    public:
        RecursiveLossModel(
            const ext::shared_ptr<ConstantLossLatentmodel<copulaPolicy> >& m,
// nope! use max common divisor. See O'Kane. Or give both options at least.
            Size nbuckets  = 1)
        : copula_(m), nBuckets_(nbuckets), wk_() { }
      private:
          /*!
          @param pDefDate Vector of unconditional default probabilities for each
          live name (at the current evaluation date). This is passed instead of
          the date for performance reasons (if in the future other magnitudes
          -e.g. lgd- are contingent on the date they shouldd be passed too).
          */
        Disposable<std::map<Real, Probability> > conditionalLossDistrib(
            const std::vector<Probability>& pDefDate,
            const std::vector<Real>& mktFactor) const;
        Real expectedConditionalLoss(const std::vector<Probability>& pDefDate, //<< never used!!
            const std::vector<Real>& mktFactor) const;
        Disposable<std::vector<Real> > conditionalLossProb(
            const std::vector<Probability>& pDefDate,
            //const Date& date,
            const std::vector<Real>& mktFactor) const;
        //versions using the P-inverse, deprecate the former
        Disposable<std::map<Real, Probability> > conditionalLossDistribInvP(
            const std::vector<Real>& pDefDate,
            //const Date& date,
            const std::vector<Real>& mktFactor) const;
        Real expectedConditionalLossInvP(const std::vector<Real>& pDefDate,
            //const Date& date,
            const std::vector<Real>& mktFactor) const;
    protected:
        void resetModel();
    public:
        /*  Expected tranche Loss calculation.
            This is computed from the first equation on page 70 (not numbered)
            Notice that while we want to compute:
            \f[
            EL(t) = \sum_{l_k}l_k P(l;t) =
              \sum_{l_k}l_k \int P(l_k;t|\omega) d\omega q(\omega)
            \f]
            One can invert the sumation and the integral order to:
            \f[
            EL(t) = \int\,q(\omega)\,d\omega\,\sum_{l_k}\,l_k\,P(l_k;t|\omega) =
              \int\,q(\omega)\,d\omega\,EL(t|\omega)
            \f]
            and this is the way it is integrated here. The recursion formula
            makes it easier this way.
        */
       Real expectedTrancheLoss(const Date& date) const;
       Disposable<std::vector<Real> > lossProbability(const Date& date) const;
       // REMEBER THIS HAS TO BE MOVED TO A DISTRIBUTION OBJECT.............
       Disposable<std::map<Real, Probability> > lossDistribution(
           const Date& d) const;
       // INTEGRATE THEN SEARCH RATHER THAN SEARCH AND THEN INTEGRATE:
       // Here I am not using a search because the point might not be attainable
       //  (loss distrib is not continuous)
       Real percentile(const Date& d, Real percentile) const;
       Real expectedShortfall(const Date& d, Real perctl) const;
    protected:
        const ext::shared_ptr<ConstantLossLatentmodel<copulaPolicy> > copula_;
    private:
        // loss model descriptor members
        const Size nBuckets_;
        mutable std::vector<Real> wk_;
        mutable Real lossUnit_;
        //! name to name factor. In the single factor copula:
        //    correl = beta * beta
        // When constructing through a single correlation number the factor is
        //   taken to be the positive swuare root of this number in the copula.
        ////////in the latent model now: mutable std::vector<Real> oneFactorCorrels_;
        // cached remaining basket magnitudes:
        mutable Real attachAmount_,
            detachAmount_,
            notional_;
        mutable Size remainingBsktSize_;
        mutable std::vector<Real> notionals_;
    };

    typedef RecursiveLossModel<GaussianCopulaPolicy> RecursiveGaussLossModel;

    // Inlines ------------------------------------------------

    template<class CP>
    inline Real RecursiveLossModel<CP>::expectedTrancheLoss(
        const Date& date) const
    {
/*
        std::map<Real, Probability> dist = lossDistribution(date);

        Real expLoss = 0.;
        std::map<Real, Probability>::iterator distIt = dist.begin();

        while(distIt != dist.end()) {
            Real loss = distIt->first * lossUnit_;
            loss = std::max(std::min(loss, detachAmount_)-attachAmount_, 0.);
            // MIN MAX BUGS ....??
            expLoss += loss * distIt->second;
            distIt++;
        }
        return expLoss ;

    ///////////////////////////////////////////////////////////////////////

        // calculate inverted unconditional Ps first so we save the inversion:
        // TO DO : turn to STL algorithm code
        std::vector<Probability> uncDefProb =
            basket_->remainingProbabilities(date);

        return copula_->integratedExpectedValue(
            boost::function<Real (const std::vector<Real>& v1)>(
                boost::bind(
                    &RecursiveLossModel::expectedConditionalLoss,
                    this,
                    boost::cref(uncDefProb),
                    _1)
                )
            );
            */
/**/
        std::vector<Probability> uncDefProb =
            basket_->remainingProbabilities(date);
        std::vector<Real> invProb;
        for(Size i=0; i<uncDefProb.size(); ++i)
           invProb.push_back(copula_->inverseCumulativeY(uncDefProb[i], i));
           ///  invProb.push_back(CP::inverseCumulativeY(uncDefProb[i], i));//<-static call
        return copula_->integratedExpectedValue(
            boost::function<Real (const std::vector<Real>& v1)>(
                boost::bind(
                    &RecursiveLossModel::expectedConditionalLossInvP,
                    this,
                    boost::cref(invProb),
                    _1)
                )
            );

    }

    template<class CP>
    inline Disposable<std::vector<Real> >
        RecursiveLossModel<CP>::lossProbability(const Date& date) const {

        std::vector<Probability> uncDefProb =
            basket_->remainingProbabilities(date);
        return copula_->integratedExpectedValue(
            boost::function<Disposable<std::vector<Real> > (const std::vector<Real>& v1)>(
                boost::bind(
                    &RecursiveLossModel::conditionalLossProb,
                    this,
                    boost::cref(uncDefProb),
                    _1)
                )
            );
    }

    // -------------------------------------------------------------------

    template<class CP>
    void RecursiveLossModel<CP>::resetModel() {
        // basket update:
        notionals_ = basket_->remainingNotionals();
        notional_  = basket_->remainingNotional();
        attachAmount_ = basket_->remainingAttachmentAmount();
        detachAmount_ = basket_->remainingDetachmentAmount();
        // model parameters:
        remainingBsktSize_ = notionals_.size();

        copula_->resetBasket(basket_.currentLink());

        std::vector<Real> lgdsTmp, lgds;
        for(Size i=0; i<remainingBsktSize_; ++i)
            lgds.push_back(notionals_[i]*(1.-copula_->recoveries()[i]));
        lgdsTmp = lgds;
        ///////////////std::remove(lgds.begin(), lgds.end(), 0.);
        lgds.erase(std::remove(lgds.begin(), lgds.end(), 0.), lgds.end());
        lossUnit_ = *(std::min_element(lgds.begin(), lgds.end()))
            / nBuckets_;
        for(Size i=0; i<remainingBsktSize_; ++i)
            wk_.push_back(std::floor(lgdsTmp[i]/lossUnit_ + .5));
    }

    // make it return a distribution object?
    template<class CP>
    Disposable<std::map<Real, Probability> >
        RecursiveLossModel<CP>::lossDistribution(const Date& d) const
    {
        std::map<Real, Probability> distrib;
        std::vector<Real> values  = lossProbability(d);
        Real sum = 0.;
        for(Size i=0; i<values.size(); ++i) {
            distrib.insert(std::make_pair<Real, Probability>(i * lossUnit_,
                sum + values[i]));
            sum += values[i];
        }
        return distrib;
    }

    // Integrate then search rather than search and then integrate?
    // Here I am not using a search because the point might be not attainable
    //   (loss distrib is not continuous)
    template<class CP>
    Real RecursiveLossModel<CP>::percentile(const Date& d,
        Real percentile) const
    {
        std::map<Real, Probability> dist = lossDistribution(d);

        if(dist.begin()->second >=1.) return dist.begin()->first;

        // deterministic case (e.g. date requested is todays date)
        if(dist.size() == 1) return dist.begin()->first;

        if(percentile == 1.) return dist.rbegin()->second;
        if(percentile == 0.) return dist.begin()->second;
        std::map<Real, Probability>::const_iterator itdist = dist.begin();
        while (itdist->second <= percentile) ++itdist;
        Real valPlus = itdist->second;
        Real xPlus   = itdist->first;
        --itdist;  //we're never 1st or last, because of tests above
        Real valMin  = itdist->second;
        Real xMin    = itdist->first;

        // return xPlus-(xPlus-xMin)*(valPlus-percentile)/(valPlus-valMin);
        Real portfLoss =  xPlus-(xPlus-xMin)*(valPlus-percentile)
            /(valPlus-valMin);
        return //remainingNotional_ *
            std::min(std::max(portfLoss - attachAmount_, 0.),
                detachAmount_ - attachAmount_);/////(detach_ - attach_);
    }

    template<class CP>
    Real RecursiveLossModel<CP>::expectedShortfall(const Date& d,
        Real perctl) const
    {
        if(d == Settings::instance().evaluationDate()) return 0.;
        std::map<Real, Probability> distrib = lossDistribution(d);

        std::map<Real, Probability>::iterator itNxt, itDist =
            distrib.begin();
        for(; itDist != distrib.end(); ++itDist)
            if(itDist->second >= perctl) break;
        itNxt = itDist;
        --itDist; // what if we are on the first one?!!!

        // One could linearly triangulate the exact point and get extra
        // precission on the first(broken) period.
        if(itNxt != distrib.end()) {
            Real lossNxt = std::min(std::max(itNxt->first - attachAmount_,
                0.), detachAmount_ - attachAmount_);
            Real lossHere = std::min(std::max(itDist->first - attachAmount_,
                0.), detachAmount_ - attachAmount_);

            Real val =  lossNxt - (itNxt->second - perctl) *
                (lossNxt - lossHere) / (itNxt->second - itDist->second);
            Real suma = (itNxt->second - perctl) * (lossNxt + val) * .5;
            ++itDist; ++itNxt;
            do{
                lossNxt = std::min(std::max(itNxt->first - attachAmount_,
                    0.), detachAmount_ - attachAmount_);
                lossHere = std::min(std::max(itDist->first - attachAmount_,
                    0.), detachAmount_ - attachAmount_);
                suma += .5 * (lossHere + lossNxt) * (itNxt->second -
                    itDist->second);
                ++itDist; ++itNxt;
            }while(itNxt != distrib.end());
            return suma / (1.-perctl);
        }
        return 0.;// well, we are in error....  fix: FAIL
    }

    template<class CP>
    Disposable<std::map<Real, Probability> >
        RecursiveLossModel<CP>::conditionalLossDistrib(
            const std::vector<Probability>& pDefDate,
            //const Date& date,
            const std::vector<Real>& mktFactor) const
    {
        //eq. 10 p.68
        //attainable losses distribution, recursive algorithm
        const std::vector<Probability>& uncDefProb = pDefDate;// alias, remove

        std::map<Real, Probability> pIndepDistrib;
        ////////  K=0
        pIndepDistrib.insert(std::make_pair(0., 1.));
        for(Size iName=0; iName<remainingBsktSize_; ++iName) {
            Probability pDef =
                copula_->conditionalDefaultProbability(uncDefProb[iName], iName,
                                                mktFactor);
            ////// iterate on all possible losses in the distribution:
            std::map<Real, Probability> pDistTemp;
            std::map<Real, Probability>::iterator distIt =
                pIndepDistrib.begin();
            while(distIt != pIndepDistrib.end()) {
              ///   update prob if this name does not default
                std::map<Real, Probability>::iterator matchIt
                    = pDistTemp.find(distIt->first);
                if(matchIt != pDistTemp.end()) {
                    matchIt->second += distIt->second * (1.-pDef);
                }else{
                    pDistTemp.insert(std::make_pair(distIt->first,
                        distIt->second * (1.-pDef)));
                }
              ////   and if it does
                matchIt = pDistTemp.find(distIt->first + wk_[iName]);
                if(matchIt != pDistTemp.end()) {
                    matchIt->second += distIt->second * pDef;
                }else{
                    pDistTemp.insert(std::make_pair(
                        distIt->first+wk_[iName], distIt->second * pDef));
                }
                ++distIt;
            }
           /////  copy back
            pIndepDistrib = pDistTemp;
        }
        /* Apply tranche limits now .... mind you this could be done outside*/
        ////  to be done....
        return pIndepDistrib;
    }

    template<class CP>
    Disposable<std::map<Real, Probability> >
        // twice?! rewrite one in terms of the other, this is a duplicate!
        RecursiveLossModel<CP>::conditionalLossDistribInvP(
            const std::vector<Real>& invpDefDate,
            //const Date& date,
            const std::vector<Real>& mktFactor) const
    {
        // eq. 10 p.68
        // attainable losses distribution, recursive algorithm

        std::map<Real, Probability> pIndepDistrib;
        // K=0
        pIndepDistrib.insert(std::make_pair(0., 1.));
        for(Size iName=0; iName<remainingBsktSize_; ++iName) {
            Probability pDef =
                copula_->conditionalDefaultProbabilityInvP(invpDefDate[iName],
                    iName, mktFactor);

            // iterate on all possible losses in the distribution:
            std::map<Real, Probability> pDistTemp;
            std::map<Real, Probability>::iterator distIt =
                pIndepDistrib.begin();
            while(distIt != pIndepDistrib.end()) {
                // update prob if this name does not default
                std::map<Real, Probability>::iterator matchIt
                    = pDistTemp.find(distIt->first);
                if(matchIt != pDistTemp.end()) {
                    matchIt->second += distIt->second * (1.-pDef);
                }else{
                    pDistTemp.insert(std::make_pair(distIt->first,
                        distIt->second * (1.-pDef)));
                }
                // and if it does
                matchIt = pDistTemp.find(distIt->first + wk_[iName]);
                if(matchIt != pDistTemp.end()) {
                    matchIt->second += distIt->second * pDef;
                }else{
                    pDistTemp.insert(std::make_pair(
                        distIt->first+wk_[iName], distIt->second * pDef));
                }
                ++distIt;
            }
            // copy back
            pIndepDistrib = pDistTemp;
        }
        /* Apply tranche limits now .... mind you this could be done outside*/
        return pIndepDistrib;
    }

    /*
    Bugs here???. The max min on the tranche looks
    wrong. It is better to have a tranche function since that way we can avoid
    adding up losses over all the posible losses rather than just over the
    tranche limits.
    */
    //! Portfolio loss conditional to the market factor value
    template<class CP>
    Real RecursiveLossModel<CP>::expectedConditionalLoss(
        const std::vector<Probability>& pDefDate,
        //const Date& date,
        const std::vector<Real>& mktFactor) const
    {
        std::map<Real, Probability> pIndepDistrib =
            conditionalLossDistrib(pDefDate, mktFactor);

        // get the expected value subject to the value of the market
        //   factor.
        Real expLoss = 0.;
        //---------------------------------------------------------------
        /* This is the original (easy to read) loop which I have partially
             unroll below to take profit of the fact that once we go over
             the tranche top the loss amount is fixed:
        */
        std::map<Real, Probability>::iterator distIt =
            pIndepDistrib.begin();

        while(distIt != pIndepDistrib.end()) {
            Real loss = distIt->first * lossUnit_;
     //       loss = std::max(std::min(loss, detachAmount_)-attachAmount_, 0.);
            loss = std::min(std::max(loss - attachAmount_, 0.),
                detachAmount_ - attachAmount_);
            // MIN MAX BUGS ....??
            expLoss += loss * distIt->second;
            ++distIt;
        }
        return expLoss ;
    }

    template<class CP>
    // again, I am duplicating code.
    Real RecursiveLossModel<CP>::expectedConditionalLossInvP(
                                 const std::vector<Real>& invPDefDate,
                                 //const Date& date,
                                 const std::vector<Real>& mktFactor) const
    {
        std::map<Real, Probability> pIndepDistrib =
            conditionalLossDistribInvP(invPDefDate, mktFactor);

        // get the expected value subject to the value of the market
        //   factor.
        Real expLoss = 0.;
        //---------------------------------------------------------------
        /* This is the original (easy to read) loop which I have partially
             unroll below to take profit of the fact that once we go over
             the tranche top the loss amount is fixed:
        */
        std::map<Real, Probability>::iterator distIt =
            pIndepDistrib.begin();

        while(distIt != pIndepDistrib.end()) {
            Real loss = distIt->first * lossUnit_;
   //         loss = std::max(std::min(loss, detachAmount_)-attachAmount_, 0.);
            loss = std::min(std::max(loss - attachAmount_, 0.),
                detachAmount_ - attachAmount_);
            // MIN MAX BUGS ....???
            expLoss += loss * distIt->second;
            ++distIt;
        }
        return expLoss ;
    }

    template<class CP>
    Disposable<std::vector<Real> > RecursiveLossModel<CP>::conditionalLossProb(
        const std::vector<Probability>& pDefDate,
        //const Date& date,
        const std::vector<Real>& mktFactor) const
    {
        std::map<Real, Probability> pIndepDistrib =
            conditionalLossDistrib(pDefDate, mktFactor);

        std::vector<Real> results;
        std::map<Real, Probability>::iterator distIt = pIndepDistrib.begin();
        while(distIt != pIndepDistrib.end()) {
            //Real loss = distIt->first * loss_unit_
            //                    ;
            //loss = std::max(std::min(loss,
            //    results_.xMax)-results_.xMin, 0.);
            //expLoss += loss * distIt->second;

            results.push_back(distIt->second);
             ++distIt;
        }
        return results;
    }

}

#endif

max_stars_repo_name: tlapfai/My-Quantlib

max_issues_repo_name: tlapfai/My-Quantlib

max_forks_repo_name: tlapfai/My-Quantlib