text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>

#include <omp.h>

#include "util.h"

int main(int argc, char *argv[]) {
    char* input_fileName1 = argv[1];
    char* input_fileName2 = argv[2];
    int N_doc_bg = atoi(argv[3]);
    int N_kw = atoi(argv[4]);
    int N_obs = atoi(argv[5]);
    int N_iter = atoi(argv[6]);
    double sec_budget = atoi(argv[7]) / 100.0;
    char* output_fileName = argv[8];
    int N_doc = 480000 / 2;
    int fixed_padding = 2 * log(2) * (64 + log(N_doc) / log(2)) / sec_budget;
    printf("Padding: %d\n", fixed_padding);

    int* matrix = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* matrix_bg = (int*) malloc(sizeof(int) * N_kw * N_kw);
    int* matrix_padded = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* true_index = (int*) malloc(sizeof(int) * N_kw);
    int* permutation = (int*) malloc(sizeof(int) * N_obs);
    gsl_matrix* matrix_obs;

    for (int round = 0; round < 10; round++)
    {
        char input_fileName1_extend[40];
        char input_fileName2_extend[40];
        sprintf(input_fileName1_extend, "%s%d", input_fileName1, round);
        sprintf(input_fileName2_extend, "%s%d", input_fileName2, round);

        // Setup
        struct timeval tv1,tv2;
        gettimeofday(&tv1, NULL);
        read_matrix(&true_index, &matrix_bg, 1.0*N_doc/N_doc_bg, N_kw, input_fileName2_extend);
        read_matrix(&true_index, &matrix, 1.0, N_obs, input_fileName1_extend);
        gettimeofday(&tv2, NULL);
        printf("Reading done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
        fflush(stdout);

        for (int iter = 0; iter < 10; iter++)
        {
            printf("Run %d\n", iter);
            matrix_obs = gsl_matrix_alloc(N_obs, N_obs);

            gettimeofday(&tv1, NULL);
            pad_matrix(&matrix_padded, &matrix, N_obs, N_doc, sec_budget, fixed_padding);
            gettimeofday(&tv2, NULL);
            printf("Padding done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            gettimeofday(&tv1, NULL);
            observe_matrix(matrix_obs, &matrix_padded, N_obs);
            gettimeofday(&tv2, NULL);
            printf("Observed matrix generated: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            // Permute observed matrix randomly and attack
            gettimeofday(&tv1, NULL);
            attack(matrix_obs, &matrix_bg, &permutation, N_kw, N_obs, N_doc, sec_budget, fixed_padding, N_iter);
            gettimeofday(&tv2, NULL);
            printf("Main attack done: %f.\n",  (tv2.tv_usec - tv1.tv_usec) / 1000000.0 + (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            char output_fileName_full[40];
            sprintf(output_fileName_full, "%s%d-%d", output_fileName, round, iter);
            print_result(output_fileName_full, &permutation, &true_index, N_obs);
            //sprintf(output_fileName_full, "%s%d-%d-full", output_fileName, round, iter);
            //print_full_result(output_fileName_full, &permutation, &true_index, N_obs);
        }
    }

    free(matrix);
    free(matrix_padded);
    gsl_matrix_free(matrix_obs);
    return(0);
}

double log_score(int idx1, int idx2, gsl_matrix* matrix_obs, int** matrix, int** permutation, int N_kw, int N_doc, double sec_budget, int fixed_padding)
{
    int idx1_m = (*permutation)[idx1];
    int idx2_m = (*permutation)[idx2];

    if (idx1 == idx2)
    {
        int N_upper = 3 * sqrt((*matrix)[idx1_m*N_kw + idx2_m] * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]) / (double) N_doc);
        int N_lower = (*matrix)[idx1_m*N_kw + idx2_m] - N_upper;
        N_upper = (*matrix)[idx1_m*N_kw + idx2_m] + N_upper;
        int count_total = (int) gsl_matrix_get(matrix_obs, idx1, idx2);

        double score = 0;
        double prob = (*matrix)[idx1_m*N_kw + idx2_m] / (double) N_doc;
        for (int kk = N_lower; kk < N_upper+1; kk+=20)
            score += gsl_ran_binomial_pdf(kk, prob, N_doc) * gsl_ran_laplace_pdf(count_total - fixed_padding - 2*kk, 2.0/sec_budget);
        if (score == 0)
            return(-500.0);
    }

    int n1 = (int) gsl_matrix_get(matrix_obs, idx1, idx1) - (*matrix)[idx1_m*N_kw + idx1_m];
    int n2 = (int) gsl_matrix_get(matrix_obs, idx2, idx2) - (*matrix)[idx2_m*N_kw + idx2_m];

    double mean = 0.5 * gsl_matrix_get(matrix_obs, idx2, idx2) / N_doc * n1 + 0.5 * gsl_matrix_get(matrix_obs, idx1, idx1) / N_doc * n2;
    double var = 1.0 * (*matrix)[idx1_m*N_kw + idx2_m] / N_doc * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]);
    var += 0.25 * n1 / N_doc * gsl_matrix_get(matrix_obs, idx2, idx2) * (2 * N_doc - gsl_matrix_get(matrix_obs, idx2, idx2)) / N_doc;
    var += 0.25 * n2 / N_doc * gsl_matrix_get(matrix_obs, idx1, idx1) * (2 * N_doc - gsl_matrix_get(matrix_obs, idx1, idx1)) / N_doc;
    int count_total = (int) gsl_matrix_get(matrix_obs, idx1, idx2);

    int N_upper = 3 * sqrt((*matrix)[idx1_m*N_kw + idx2_m] * (N_doc - (*matrix)[idx1_m*N_kw + idx2_m]) / (double) N_doc);
    int N_lower = (*matrix)[idx1_m*N_kw + idx2_m] - N_upper;
    N_upper = (*matrix)[idx1_m*N_kw + idx2_m] + N_upper;

    double score = 0;
    double prob = (*matrix)[idx1_m*N_kw + idx2_m] / (double) N_doc;
    for (int kk = N_lower; kk < N_upper+1; kk+=20)
        score +=  gsl_ran_binomial_pdf(kk, prob, N_doc) * gsl_ran_gaussian_pdf(count_total - mean - kk, sqrt(var));

    if (score == 0)
        return(-500.0);
    return(log(20*score));
}

void initial_solution(int** permutation, gsl_matrix* matrix_obs, int** matrix, int fixed_padding, int N_kw, int N_obs)
{
    int* diff = (int*) malloc(sizeof(int) * N_kw);
    int* index = (int*) malloc(sizeof(int) * N_kw);
    int* occupied = (int*) malloc(sizeof(int) * N_kw);

    for (int ii = 0; ii < N_kw; ii++)
        occupied[ii] = -1;

    for (int ii = 0; ii < N_obs; ii++)
    {
        // reset index
        for (int jj = 0; jj < N_kw; jj++)
            index[jj] = jj;

        // compute differences
        for (int jj = 0; jj < N_kw; jj++)
            diff[jj] = abs(gsl_matrix_get(matrix_obs, ii, ii) - fixed_padding - 2*(*matrix)[jj*N_kw+jj]);

        // sorting
        for (int jj = 0; jj < N_kw-1; jj++)
        {
            for (int kk = 0; kk < N_kw-jj-1; kk++)
            {
                if (diff[kk] > diff[kk+1])
                {
                    int temp = diff[kk+1];
                    diff[kk+1] = diff[kk];
                    diff[kk] = temp;

                    temp = index[kk+1];
                    index[kk+1] = index[kk];
                    index[kk] = temp;
                }
            }
        }

        // use the smallest index available
        int idx = 0;
        while (occupied[idx] > 0)
            idx++;

        // update the data structures
        (*permutation)[ii] = index[idx];
        occupied[idx] = 1;
    }

}

void attack(gsl_matrix* matrix_obs, int** matrix, int** permutation, int N_kw, int N_obs, int N_doc, double sec_budget, int fixed_padding, int N_iter)
{
    // Initialise data structures
    double* score_matrix = (double*) malloc(sizeof(double) * N_obs * N_obs);
    double* score_row1 = (double*) malloc(sizeof(double) * N_obs);
    double* score_row2 = (double*) malloc(sizeof(double) * N_obs);

    int* permutation_tmp = (int*) malloc(sizeof(int) * N_obs);
    int* permutation_inv = (int*) malloc(sizeof(int) * N_kw);

    // Initialise permutations
    // Initial solution
    initial_solution(permutation, matrix_obs, matrix, fixed_padding, N_kw, N_obs);
    for (int ii = 0; ii < N_obs; ii++)
        permutation_tmp[ii] = (*permutation)[ii];
    for (int ii = 0; ii < N_kw; ii++)
        permutation_inv[ii] = -1;
    for (int ii = 0; ii < N_obs; ii++)
        permutation_inv[permutation_tmp[ii]] = ii;

    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    struct timeval tv1,tv2;
    gettimeofday(&tv1, NULL);

    // Compute initial score
    #pragma omp parallel for shared(score_matrix, matrix_obs, matrix)
    for (int ii = 0; ii < N_obs * N_obs; ii++)
        score_matrix[ii] = log_score((int) (ii / N_obs), ii % N_obs, matrix_obs, matrix, permutation, N_kw, N_doc, sec_budget, fixed_padding);
    gettimeofday(&tv2, NULL);
    printf("Initial score computed: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));

    // Iterations of simulated annealing
    double temp = (double) N_kw;
    int N_stuck = 0;
    for (int iter = 0; iter < N_iter; iter++)
    {
        /* Status code */
        if (iter % (N_iter / 10) == 0)
        {
            gettimeofday(&tv1, NULL);
            printf("Iteration: %d, %d, %d.\n", iter, N_stuck, (int) (tv1.tv_sec - tv2.tv_sec));
            fflush(stdout);
            gettimeofday(&tv2, NULL);
        }

	    // used to be 10k
        if (N_stuck >= N_kw * 20)
            iter = N_iter;

        /* Main code */
        int idx1, idx2;
        permutation_generation(&idx1, &idx2, &permutation_tmp, permutation, &permutation_inv, matrix_obs, matrix, fixed_padding, N_kw, N_obs, N_doc, sec_budget);
        if (idx1 == idx2)
        {
            N_stuck += 1;
                continue;
        }

        #pragma omp parallel for shared(score_row1)
        for (int ii = 0; ii < N_obs; ii++)
            score_row1[ii] = log_score(idx1, ii, matrix_obs, matrix, &permutation_tmp, N_kw, N_doc, sec_budget, fixed_padding);

        if (idx2 >= 0)
            #pragma omp parallel for shared(score_row2)
            for (int ii = 0; ii < N_obs; ii++)
                score_row2[ii] = log_score(idx2, ii, matrix_obs, matrix, &permutation_tmp, N_kw, N_doc, sec_budget, fixed_padding);

        double score_diff = 0;
        for (int ii = 0; ii < N_obs; ii++)
            score_diff += score_row1[ii];
        for (int ii = 0; ii < N_obs; ii++)
            score_diff -= score_matrix[idx1*N_obs + ii];
        if (idx2 >= 0)
        {
            for (int ii = 0; ii < N_obs; ii++)
                score_diff += score_row2[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_diff -= score_matrix[idx2*N_obs + ii];
        }

        // compute difference in score, with exponentiation
        score_diff = score_diff / temp;

        if (score_diff < -40)
            score_diff = 0;
        else if (score_diff > 0)
            score_diff = 1.01;
        else
            score_diff = exp(score_diff);

        if (gsl_ran_flat(r, 0, 1) < score_diff)
        {
            // Update the scores
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[idx1*N_obs + ii] = score_row1[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[ii*N_obs + idx1] = score_row1[ii];
            if (idx2 >= 0)
            {
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[idx2*N_obs + ii] = score_row2[ii];
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[ii*N_obs + idx2] = score_row2[ii];
            }

            // Update the permutation
            permutation_inv[(*permutation)[idx1]] = -1;
            (*permutation)[idx1] = permutation_tmp[idx1];
            permutation_inv[permutation_tmp[idx1]] = idx1;
            if (idx2 >= 0)
            {
                (*permutation)[idx2] = permutation_tmp[idx2];
                permutation_inv[permutation_tmp[idx2]] = idx2;
            }
            N_stuck = 0;
        }
        else
        {
            // Update the permutation
            permutation_tmp[idx1] = (*permutation)[idx1];
            if (idx2 >= 0)
                permutation_tmp[idx2] = (*permutation)[idx2];
            N_stuck += 1;
        }

        temp *= 0.995;
    }

    free(score_matrix);
    free(score_row1);
    free(score_row2);
    gsl_rng_free(r);
}

void print_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    for (int ii = 0; ii < N_obs; ii++)
    {
        //if (ii < 50)
        //    printf("%d, %d\n", (*permutation)[ii], (*true_index)[ii]);
        if ((*permutation)[ii] == (*true_index)[ii])
            count++;
    }

    fprintf(fp, "%d\n", count);
    fclose(fp);
    printf("Success: %d/%d.\n", count, N_obs);
}

void print_full_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    for (int ii = 0; ii < N_obs; ii++)
        if ((*permutation)[ii] == (*true_index)[ii])
            count++;

    fprintf(fp, "%d\n", count);
    for (int ii = 0; ii < N_obs; ii++)
        fprintf(fp, "%d,%d\n", (*permutation)[ii], (*true_index)[ii]);

    fclose(fp);
}

ext: c

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
///  @file rainbow-gen-usersk.c
///  @brief A command-line tool for generating the user secret keys.

#include <rainbow_keypair.h>
#include <stdlib.h>
#include <stdio.h>
#include <api.h>
#include <utils.h>
#include <string.h>
#include <blas.h>

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

    printf("%s\n", CRYPTO_ALGNAME);

    printf("msk size: %lu\n", CRYPTO_MASTER_SECRET_KEY_BYTES);
    printf("mpk size: %lu\n", CRYPTO_MASTER_PUBLIC_KEY_BYTES);
    printf("hash size: %d\n", _HASH_LEN);
    printf("signature size: %d\n\n", CRYPTO_BYTES);

    if (!(4 == argc || 5 == argc)) {
        printf("Usage:\n\n\trainbow-gen-usersk msk_file_name identity usk_file_name [upk_file_name]\n\n");
        return -1;
    }

    uint8_t *_msk = malloc(CRYPTO_MASTER_SECRET_KEY_BYTES);
    FILE *fp;
    unsigned r;

    fp = fopen(argv[1], "r");
    if (NULL == fp) {
        printf("fail to open master key file.\n");
        return -1;
    }
    r = byte_fget(fp, _msk, CRYPTO_MASTER_SECRET_KEY_BYTES);
    fclose(fp);
    if (CRYPTO_MASTER_SECRET_KEY_BYTES != r) {
        printf("fail to load key file.\n");
        return -1;
    }

    unsigned char *_identity = malloc((_ID + 1) / 2);
    r = 0;

    generate_identity_hash(_identity, (unsigned char *) argv[2], strlen(argv[2]));
    if (NULL == _identity) {
        printf("fail to create identity hash.\n");
        return -1;
    }
    printf("identity hash: ");
    for (unsigned i = 0; i < _ID; i++) {
        printf("%hhu ", gf16v_get_ele(_identity, i));
    }
    printf("\n\n");

    uint8_t *_usk = malloc(sizeof(usk_t));

    if (argc == 4) {
        //calculate usk with mpk and ID
        int re = calculate_usk((usk_t *) _usk, (msk_t *) _msk, _identity);
        if (0 != re) {
            printf("%s generate user-secret-key fails.\n", CRYPTO_ALGNAME);
            return -1;
        }
    } else {
        //calculate upk from usk (for debugging), if filename is given
        uint8_t *_upk = malloc(sizeof(upk_t));
        int re = calculate_usk_and_upk((usk_t *) _usk, (upk_t *) _upk, (msk_t *) _msk, _identity);
        if (0 != re) {
            printf("%s generate upk and usk fails.\n", CRYPTO_ALGNAME);
            return -1;
        }

        //write upk to disk
        fp = fopen(argv[4], "w+");
        if (NULL == fp) {
            printf("fail to write user public key file.\n");
            return -1;
        }
        byte_fdump(fp, CRYPTO_ALGNAME " user-public-key", _upk,
                   CRYPTO_USER_PUBLIC_KEY_BYTES); //upk speichern und beschriften
        fclose(fp);

        free(_upk);
    }

    //write usk to disk
    fp = fopen(argv[3], "w+");
    if (NULL == fp) {
        printf("fail to write user secret key file.\n");
        return -1;
    }
    byte_fdump(fp, CRYPTO_ALGNAME " user secret key", _usk, sizeof(usk_t)); //usk speichern und beschriften
    fclose(fp);

    //free usk and rest

    free(_msk);
    free(_identity);
    free(_usk);

    return 0;
}

max_stars_repo_name: mnmuser/ibs-rainbow

max_issues_repo_name: mnmuser/ibs-rainbow

max_forks_repo_name: mnmuser/ibs-rainbow


text:
#include <api/components.h>
#include <api/kvstore_itf.h>
#include <common/string_view.h>
#include <common/utils.h> /* MiB */
#include <gsl/pointers> /* not_null */

#include <cstddef> /* size_t */
#include <memory> /* unique_ptr */
#include <stdexcept> /* runtime_error */
#include <string>

namespace
{
	/* things which differ depending on the type of store used */
	struct custom_store
	{
		virtual ~custom_store() {}
		virtual std::size_t minimum_size(std::size_t s) const { return s; }
		virtual component::uuid_t factory() const = 0;
		virtual std::size_t presumed_allocation() const = 0;
		virtual bool uses_numa_nodes() const { return false; }
		virtual status_t rc_unknown_attribute() const { return E_NOT_SUPPORTED; }
		virtual status_t rc_percent_used() const { return S_OK; }
		virtual status_t rc_resize_locked() const { return E_LOCKED; }
		virtual status_t rc_attribute_key_null_ptr() const { return E_BAD_PARAM; }
		virtual status_t rc_attribute_key_not_found() const { return component::IKVStore::E_KEY_NOT_FOUND; }
		virtual status_t rc_attribute_hashtable_expansion() const { return S_OK; }
		virtual status_t rc_out_of_memory() const { return component::IKVStore::E_TOO_LARGE; }
		virtual status_t rc_atomic_update() const { return S_OK; }
		virtual status_t rc_allocate_pool_memory_size_0() const { return S_OK; }
		virtual bool swap_updates_timestamp() const { return false; }
	};

	struct custom_mapstore
		: public custom_store
	{
		virtual component::uuid_t factory() const override { return component::mapstore_factory; }
		std::size_t minimum_size(std::size_t s) const override { return std::max(std::size_t(8), s); }
		std::size_t presumed_allocation() const override { return 1ULL << DM_REGION_LOG_GRAIN_SIZE; }
		bool uses_numa_nodes() const override { return true; }
		status_t rc_unknown_attribute() const override { return E_INVALID_ARG; }
		status_t rc_percent_used() const override { return rc_unknown_attribute(); }
		status_t rc_resize_locked() const override { return E_INVAL; }
		status_t rc_attribute_key_null_ptr() const override { return E_INVALID_ARG; }
		status_t rc_attribute_key_not_found() const override { return rc_unknown_attribute(); }
		status_t rc_attribute_hashtable_expansion() const override { return rc_unknown_attribute(); }
		status_t rc_out_of_memory() const override { return E_INVAL; }
		status_t rc_atomic_update() const override { return E_NOT_SUPPORTED; }
		status_t rc_allocate_pool_memory_size_0() const override { return E_INVAL; }
	};

	struct custom_hstore
		: public custom_store
	{
		virtual component::uuid_t factory() const override { return component::hstore_factory; }
		std::size_t presumed_allocation() const override { return MiB(32); }
		bool swap_updates_timestamp() const override { return true; }
	};

	struct custom_hstore_cc
		: public custom_hstore
	{
		std::size_t presumed_allocation() const override { return MiB(32); }
	};

	custom_mapstore custom_mapstore_i{};
	custom_hstore custom_hstore_i{};
	custom_hstore_cc custom_hstore_cc_i{};

	const std::map<std::string, gsl::not_null<custom_store *>, std::less<>> custom_map =
	{
		{ "mapstore", &custom_mapstore_i },
		{ "hstore", &custom_hstore_i },
		{ "hstore-cc", &custom_hstore_cc_i },
		{ "hstore-mc", &custom_hstore_i },
		{ "hstore-mr", &custom_hstore_i },
		{ "hstore-mm", &custom_hstore_i },
	};
}

gsl::not_null<custom_store *> locate_custom_store(common::string_view store)
{
	const auto c_it = custom_map.find(store);
	if ( c_it == custom_map.end() )
    {
		throw std::runtime_error(common::format("store {} not recognized", store));
    }
	return c_it->second;
}

auto make_kvstore(
	common::string_view store
	, gsl::not_null<custom_store *> c
	, const component::IKVStore_factory::map_create & mc
) -> std::unique_ptr<component::IKVStore>
{
	using IKVStore_factory = component::IKVStore_factory;
	const std::string store_lib = "libcomponent-" + std::string(store) + ".so";
	auto *comp = component::load_component(store_lib.c_str(), c->factory());
	if ( ! comp )
    {
		throw std::runtime_error(common::format("failed to load component {}", store_lib));
    }

	const auto fact = make_itf_ref(static_cast<IKVStore_factory*>(comp->query_interface(IKVStore_factory::iid())));
	const auto kvstore = fact->create(0, mc);

	return std::unique_ptr<component::IKVStore>(kvstore);
}

ext: h

max_stars_repo_name: IBM/artemis

max_issues_repo_name: IBM/artemis

max_forks_repo_name: IBM/artemis


text:
#pragma once

#include <sstream>
#include <gsl/string_span>

/*
	String tokenizer that works like the ToEE vanilla one.
*/
class Tokenizer {
public:

	explicit Tokenizer(const std::string &input) : mIn(input) {
	}

	bool NextToken();

	bool IsQuotedString() const {
		return mTokenType == TokenType::QuotedString;
	}

	bool IsNumber() const {
		return mTokenType == TokenType::Number;
	}

	bool IsIdentifier() const {
		return mTokenType == TokenType::Identifier;
	}

	bool IsIdentifier(const char *identifier) const;

	const std::string &GetTokenText() const {
		return mTokenText;
	}

	const bool& GetEnableEscapes() const {
		return mEnableEscapes;
	}

	void SetEnableEscapes(bool enableEscapes) {
		mEnableEscapes = enableEscapes;
	}

	int GetTokenInt() const {
		return mTokenInt;
	}

	float GetTokenFloat() const {
		return (float) mTokenFloat;
	}

private:
	enum class TokenType {
		Number,
		QuotedString,
		Identifier,
		Unknown
	};

	std::istringstream mIn;
	std::string mTokenText;
	TokenType mTokenType = TokenType::Unknown;
	double mTokenFloat = 0;
	int mTokenInt = 0;

	std::string mLine; // Buffer for current line
	size_t mLinePos; // Current position within mLine
	int mLineNo = 0;
	bool mEnableEscapes = true;

	bool GetLine();
	bool LineHasMoreChars() const;

	bool ReadNumber();
	bool ReadQuotedString();
	bool ReadIdentifier();

	// Character based reading on the input source
	char PeekChar();
	void SkipChar();
	char TakeChar();
	void UngetChar();

	// Seeks past any control or space characters at current line pos
	void SkipSpaceAndControl();
	// Skips past comment at current line pos to end of line
	void SkipComment();
};

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Erhan Kenar $
// $Authors: Vipul Patel $
// --------------------------------------------------------------------------
#ifndef OPENMS_COMPARISON_SPECTRA_COMPAREFOURIERTRANSFORM_H
#define OPENMS_COMPARISON_SPECTRA_COMPAREFOURIERTRANSFORM_H

#include <OpenMS/COMPARISON/SPECTRA/PeakSpectrumCompareFunctor.h>
#include <OpenMS/KERNEL/StandardTypes.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>

namespace OpenMS
{
  /**
      @brief Compare Discrete Cosines value from a Fourier transformation, also known as Discrete Cosines Transformation

      The Direct Cosines Transformation based on the theory of the Fourier transformation.
      In this class the Fast Fourier Transformation(FFT) algorithm of the gsl library is used. FFT has a run-time complexity of
      n (log n). To get the Direct Cosines Transformation from a FFT there is preparation necessary. First the
      input data has to be mirrored. This is necessary, because FFT needs data which has a periodic nature.
      After the computation of FFT only the cosine values are important and stored in the Meta Data Array. So an inverse transformation of these values
      to get the original spectrum is not available.
      The comparison is done between two Meta Data Arrays, which contain the stored cosine values of their individual spectrum.
      The advantage of this method is how the comparison works. There is only one sum which has to be count, no multiplication is needed.

      Attention: only use the compare function, if the Spectrum was transformed earlier, else an error is going to appear.
      Only use this method of transformation, if you are sure there exists enough free memory. This is a fast estimation, but it only gives one or
      zero back.

      @htmlinclude OpenMS_CompareFouriertransform.parameters

      @ingroup SpectraComparison

  */

  class OPENMS_DLLAPI CompareFouriertransform :
    public PeakSpectrumCompareFunctor
  {
public:

    // @name Constructors and Destructors
    // @{
    /// default constructor
    CompareFouriertransform();

    /// copy constructor
    CompareFouriertransform(const CompareFouriertransform & source);

    /// destructor
    virtual ~CompareFouriertransform();
    // @}

    // @name Operators
    // @{
    /// assignment operator
    CompareFouriertransform & operator=(const CompareFouriertransform & source);

    /**
      @brief Dummy function

      This function only returns 0 for any given PeakSpectrum, please use the other compare operator function
    */
    double operator()(const PeakSpectrum &) const;
    /**
        @brief compare two PeakSpectrum by their Discrete Cosines Transformation.

        This function compares two given PeakSpectrum about their  Discrete Cosines Transformation.
        First, a transformation has to be calculated. Please use the function transform() in this class, before calling this
            function. The comparison works by summing the subtractions of each coefficient for all elements of both transformations. sum(_i=1)
            ^n x_i-y_i. If the sum is zero, both Spectra are identical in the real part and one is emitted, otherwise a zero.
        */
    double operator()(const PeakSpectrum & spec1, const PeakSpectrum & spec2) const;

    ///
    static PeakSpectrumCompareFunctor * create() { return new CompareFouriertransform(); }

    ///Returns the name used in the factory
    static const String getProductName()
    {
      return "CompareFouriertransform";
    }

    /**
    @brief calculate the Discrete Cosines Fourier Transformation.

    This function transforms a given PeakSpectrum to a Discrete Cosines Fourier
    Transformation. It stores only the part of the cosines of the FFT in the
    FloatDataArray which is a container from the PeakSpectrum.  Only call this
    function, if you are sure there is no other transformation done earlier
    over the same PeakSpectrum, because it isn't
    checked if there already exists a transformation.
    */
    void transform(PeakSpectrum & spec);
protected:
    /**
        @brief Search in the PeakSpectrum, if a Discrete Fourier transformation
        occurs, if not an error is going to be thrown, else the index
        of the occurrence is returned.

        This function gives back the position, where the transformation was
        saved in a FloatDataArray. If there is no entry, an error is thrown to
        indicate that a transformation has to be calculated before calling this
        comparison operator.
    */
    UInt searchTransformation_(const PeakSpectrum & spec) const;

  };

}
#endif /*OPENMS_COMPARISON_SPECTRA_COMPAREFOURIERTRANSFORM_H*/

max_stars_repo_name: open-ms/all-svn-branches

max_issues_repo_name: kreinert/OpenMS

max_forks_repo_name: kreinert/OpenMS


text:
#ifndef ESMCI_Mat_H
#define ESMCI_Mat_H

#include "ESMCI_Macros.h"
#include <stdexcept>
#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include <initializer_list>
#include <complex.h>
#include <cblas.h>
#include <lapacke.h>

namespace ESMCI{
  namespace MapperUtil{

    /* A Matrix class
     * The matrix class is used to create an n-dimensional matrix of values
     */
    template<typename T>
    class Matrix{
      public:
        Matrix(const std::vector<int> &dims, const std::vector<T> &data);
        Matrix(std::initializer_list<int> dims, std::initializer_list<T> data);
        Matrix(const std::vector<int> &dims, const T &val);

        /* Get the matrix dimensions */
        std::vector<int> get_dims(void ) const;
        /* Get a reference to the data strored in the matrix */
        const T *get_data_by_ref(void ) const;
        T *get_data_by_ref(void );
        /* Get the data stored in the matrix */
        std::vector<T> get_data(void ) const;

        /* Get the transpose of the matrix */
        Matrix<T> transpose(void ) const;
        /* Get the inverse of the matrix */
        Matrix<T> inv(void ) const;
        /* Get the pseudo inverse of the matrix */
        Matrix<T> pinv(void ) const;

        /* Returns true if the matrix, m, is equal to this matrix,
         * false otherwise. This function uses a tolerance, tol, to
         * compare the values stored in the matrix (unlike == that
         * does not provide a tolerance)
         */
        bool equals(const Matrix<T> &m, double tol) const;

        /* Compare if two matrices are equal */
        template<typename U>
        friend bool operator==(const Matrix<U> &lhs, const Matrix<U> &rhs);
        /* Compare if all values in the matrix is equal to val */
        template<typename U>
        friend bool operator==(const Matrix<U> &m, const U &val);
        template<typename U>
        friend bool operator==(const U &val, const Matrix<U> &m);
        /* Operator to add/subtract/multiply two matrices */
        template<typename U>
        friend Matrix<U> operator+(const Matrix<U> &lhs, const Matrix<U> &rhs);
        template<typename U>
        friend Matrix<U> operator-(const Matrix<U> &lhs, const Matrix<U> &rhs);
        template<typename U>
        friend Matrix<U> operator*(const Matrix<U> &lhs, const Matrix<U> &rhs);
        /* Operator to scale a matrix by a scalar */
        template<typename U, typename V>
        friend Matrix<U> operator*(const Matrix<U> &lhs, const V &rhs);
        template<typename U, typename V>
        friend Matrix<V> operator*(const U &lhs, const Matrix<V> &rhs);
        template<typename U>
        friend std::ostream& operator<<(std::ostream &ostr, const Matrix<U> &m);
      private:
        std::vector<int> dims_;
        std::vector<T> data_;
    };

    template<typename T>
    Matrix<T>::Matrix(const std::vector<int> &dims, const std::vector<T> &data):dims_(dims),data_(data)
    {}

    template<typename T>
    Matrix<T>::Matrix(std::initializer_list<int> dims, std::initializer_list<T> data):dims_(dims.begin(),dims.end()),data_(data.begin(),data.end())
    {}

    template<typename T>
    Matrix<T>::Matrix(const std::vector<int> &dims, const T& val):dims_(dims)
    {
      assert(dims.size() > 0);
      int data_sz = std::accumulate(dims.begin(), dims.end(), 1, std::multiplies<int>());
      assert(data_sz > 0);
      data_.resize(data_sz, val);
    }

    template<typename T>
    std::vector<int> Matrix<T>::get_dims(void ) const
    {
      return dims_;
    }

    template<typename T>
    const T *Matrix<T>::get_data_by_ref(void ) const
    {
      return ((data_.size() > 0 ) ? (&(data_[0])) : NULL);
    }

    template<typename T>
    T *Matrix<T>::get_data_by_ref(void )
    {
      return ((data_.size() > 0 ) ? (&(data_[0])) : NULL);
    }

    template<typename T>
    std::vector<T> Matrix<T>::get_data(void ) const
    {
      return data_;
    }

    template<typename T>
    Matrix<T> Matrix<T>::transpose(void ) const
    {
      assert(dims_.size() <= 2);
      int nelems = static_cast<int>(data_.size());
      if(dims_.size() == 1){
        return Matrix<T>(dims_, data_);
      }
      else{
        std::vector<int> dims(dims_.rbegin(), dims_.rend());
        if((dims_[0] == 1) || dims_[1] == 1){
          return Matrix<T>(dims, data_);
        }
        else{
          // FIXME: Use inplace transpose
          std::vector<T> data(data_.size(), static_cast<T>(0));
          int nrows = dims_[0];
          int ncols = dims_[1];
          int tnrows = ncols;
          int tncols = nrows;
          for(int idx = 0; idx < nelems; idx++){
            int ridx = idx/ncols;
            int cidx = idx - ridx * ncols;
            int tridx = cidx;
            int tcidx = ridx;
            int tidx = tridx * tncols + tcidx;
            data[tidx] = data_[idx];
          }
          return Matrix<T>(dims, data);
        }
      }
    }

    /* Calculate the inverse of the matrix using LAPACK */
    inline int LAPACK_Minv(int m, float *A)
    {
      lapack_int n = static_cast<lapack_int>(m);
      lapack_int *ipiv = (lapack_int *)calloc(n+1, sizeof(lapack_int));
      lapack_int info;
      /* Calculate LU decomposition of the matrix */
      info = LAPACKE_sgetrf(LAPACK_ROW_MAJOR, n, n, A, n, ipiv);
      /* info == 0 => success */
      if(info < 0){
        std::cout << "LAPACKE_sgetrf failed, the "
          << -info << "th arg in input array, A[" << -info << "] = "
          << A[-info] << " is invalid\n";
        return ESMF_FAILURE;
      }
      else if(info > 0){
        std::cout << "LAPACKE_sgetrf failed, the U["
          << info << "," << info << "] = 0, U is singular and div by zero can occur"
          << " if used to solve a system of equations\n";
        return ESMF_FAILURE;
      }
      /* Calculate the inverse of the matrix */
      info = LAPACKE_sgetri(LAPACK_ROW_MAJOR, n, A, n, ipiv);
      /* info == 0 => success */
      if(info < 0){
        std::cout << "LAPACKE_sgetri failed, the "
          << -info << "th arg in input array, A[" << -info << "] = "
          << A[-info] << " is invalid\n";
        return ESMF_FAILURE;
      }
      else if(info > 0){
        std::cout << "LAPACKE_sgetri failed, the U["
          << info << "," << info << "] = 0, the matrix is singular and its inverse "
          << "cannot be computed\n";
        return ESMF_FAILURE;
      }
      free(ipiv);

      return ESMF_SUCCESS;
    }

    /* Calculate the inverse of the matrix using LAPACK */
    inline int LAPACK_Minv(int m, double *A)
    {
      lapack_int n = static_cast<lapack_int>(m);
      lapack_int *ipiv = (lapack_int *)calloc(n+1, sizeof(lapack_int));
      lapack_int info;
      /* Calculate LU decomposition of the matrix */
      info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, A, n, ipiv);
      /* info == 0 => success */
      if(info < 0){
        std::cout << "LAPACKE_dgetrf failed, the "
          << -info << "th arg in input array, A[" << -info << "] = "
          << A[-info] << " is invalid\n";
        return ESMF_FAILURE;
      }
      else if(info > 0){
        std::cout << "LAPACKE_dgetrf failed, the U["
          << info << "," << info << "] = 0, U is singular and div by zero can occur"
          << " if used to solve a system of equations\n";
        return ESMF_FAILURE;
      }
      /* Calculate the inverse of the matrix */
      info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, A, n, ipiv);
      /* info == 0 => success */
      if(info < 0){
        std::cout << "LAPACKE_dgetri failed, the "
          << -info << "th arg in input array, A[" << -info << "] = "
          << A[-info] << " is invalid\n";
        return ESMF_FAILURE;
      }
      else if(info > 0){
        std::cout << "LAPACKE_dgetri failed, the U["
          << info << "," << info << "] = 0, the matrix is singular and its inverse "
          << "cannot be computed\n";
        return ESMF_FAILURE;
      }
      free(ipiv);

      return ESMF_SUCCESS;
    }

    /* Calculate the inverse of the matrix */
    template<typename T>
    Matrix<T> Matrix<T>::inv(void ) const
    {
      Matrix<T> res = *this;
      if(dims_.size() == 1){
        return res;
      }

      // We are only concerned about 2d matrices for now
      assert(dims_.size() == 2);
      assert(std::adjacent_find(dims_.cbegin(), dims_.cend(), std::not_equal_to<T>()) == dims_.cend());

      int ret = LAPACK_Minv(res.dims_[0], res.get_data_by_ref());
      //assert(ret == ESMF_SUCCESS);
      if(ret != ESMF_SUCCESS){
        std::cerr << "Finding matrix inv failed for : \n" << *this << "\n";
        std::string err_msg("LAPACK routine to find inv failed");
        throw std::runtime_error(err_msg);
      }

      return res;
    }

    /* Compute pseudo inverse = Inverse(A_transpose * A) * A_transpose */
    template<typename T>
    Matrix<T> Matrix<T>::pinv(void ) const
    {
      Matrix<T> trans = transpose();
      try{
        // Left inverse
        return ((trans * (*this)).inv() * trans);
      }
      catch(...){
        // Right inverse
        return (trans * ((*this) * trans).inv());
      }
    }

    template<typename T>
    bool Matrix<T>::equals(const Matrix<T> &m, double tol) const
    {
      std::vector<int> m_dims = m.get_dims();
      std::vector<T> m_data = m.get_data();

      if(m_dims != dims_){
        return false;
      }

      if(m_data.size() != data_.size()){
        return false;
      }

      for(typename std::vector<T>::const_iterator citer1 = m_data.cbegin(),
          citer2 = data_.cbegin();
          (citer1 != m_data.cend()) && (citer2 != data_.cend());
          ++citer1, ++citer2){
        if(fabs(*citer1 - *citer2) > tol){
          return false;
        }
      }

      return true;
    }

    template<typename T>
    bool operator==(const Matrix<T> &lhs, const Matrix<T> &rhs)
    {
      if(lhs.dims_ != rhs.dims_){
        return false;
      }
      if(lhs.data_ != rhs.data_){
        return false;
      }
      return true;
    }

    template<typename T>
    bool operator==(const Matrix<T> &m, const T& val)
    {
      for(typename std::vector<T>::const_iterator citer = m.data_.cbegin();
          citer != m.data_.cend(); ++citer){
        if(*citer != val){
          return false;
        }
      }
      return true;
    }

    template<typename T>
    bool operator==(const T &val, const Matrix<T> &m)
    {
      return (m == val);
    }

    template<typename T>
    Matrix<T> operator+(const Matrix<T> &lhs, const Matrix<T> &rhs)
    {
      assert(lhs.dims_.size() == rhs.dims_.size());
      assert(std::equal(lhs.dims_.begin(), lhs.dims_.end(), rhs.dims_.begin()));
      Matrix<T> res = lhs;
      std::transform(res.data_.begin(), res.data_.end(), rhs.data_.begin(),
                      res.data_.begin(), std::plus<T>());
      return res;
    }

    template<typename U, typename V>
    Matrix<U> operator+(const Matrix<U> &lhs, const V &rhs)
    {
      Matrix<U> res = lhs;
      std::transform(res.data_.begin(), res.data_.end(),
                      res.data_.begin(),
                      std::bind1st(std::plus<U>(), static_cast<U>(rhs)));
      return res;
    }

    template<typename U, typename V>
    Matrix<V> operator+(const U &lhs, const Matrix<V> &rhs)
    {
      return rhs + lhs;
    }

    template<typename T>
    Matrix<T> operator-(const Matrix<T> &lhs, const Matrix<T> &rhs)
    {
      assert(lhs.dims_.size() == rhs.dims_.size());
      assert(std::equal(lhs.dims_.begin(), lhs.dims_.end(), rhs.dims_.begin()));
      Matrix<T> res = lhs;
      std::transform(res.data_.begin(), res.data_.end(), rhs.data_.begin(),
                      res.data_.begin(), std::minus<T>());
      return res;
    }

    template<typename U, typename V>
    Matrix<U> operator-(const Matrix<U> &lhs, const V &rhs)
    {
      Matrix<U> res = lhs;
      std::transform(res.data_.begin(), res.data_.end(),
                      res.data_.begin(),
                      std::bind1st(std::minus<U>(), static_cast<U>(rhs)));
      return res;
    }

    /* Multiply two matrices using BLAS
     */
    inline int BLAS_Mmult(int m, int n, int k, int alpha,
                          const float *A, int lda, const float *B, int ldb,
                          int beta, float *C, int ldc)
    {
      /* cblas_sgemm has no return code ! */
      cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                  m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
      return ESMF_SUCCESS;
    }

    /* Multiply two matrices using BLAS
     */
    inline int BLAS_Mmult(int m, int n, int k, int alpha,
                          const double *A, int lda, const double *B, int ldb,
                          int beta, double *C, int ldc)
    {
      /* cblas_sgemm has no return code ! */
      cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                  m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
      return ESMF_SUCCESS;
    }

    /* Multiply two matrices, currently uses BLAS routines for matrix
     * multiplication
     */
    template<typename T>
    Matrix<T> operator*(const Matrix<T> &lhs, const Matrix<T> &rhs)
    {
      assert((lhs.dims_.size() > 0) && (rhs.dims_.size() > 0));
      // No tensor products for now
      assert((lhs.dims_.size() <= 2) && (rhs.dims_.size() <= 2));
      assert(lhs.dims_[lhs.dims_.size()-1] == rhs.dims_[0]);
      /* Determine inputs to the BLAS routine to perform
       * matrix multiplication
       */
      int m, n, k, alpha, lda, ldb, beta, ldc;
      alpha = 1;
      beta = 1;

      std::vector<int> res_dims;
      if(lhs.dims_.size() == 1){
        res_dims.push_back(lhs.dims_[0]);
        m = 1;
        k = 1;
        n = 1;
        lda = 1;
        ldb = 1;
        ldc = 1;
      }
      else{
        assert(lhs.dims_.size() == 2);
        res_dims.push_back(lhs.dims_[0]);
        res_dims.push_back(rhs.dims_[rhs.dims_.size()-1]);
        m = lhs.dims_[0];
        k = lhs.dims_[lhs.dims_.size() - 1];
        n = rhs.dims_[rhs.dims_.size() - 1];
        lda = k;
        ldb = n;
        ldc = n;
      }

      Matrix<T> res(res_dims, 0);
      const T *A = lhs.get_data_by_ref();
      const T *B = rhs.get_data_by_ref();
      T *C = res.get_data_by_ref();

      /* Use BLAS for multiplying the two matrices */
      int ret = BLAS_Mmult(m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
      assert(ret == ESMF_SUCCESS);

      return res;
    }

    template<typename U, typename V>
    Matrix<U> operator*(const Matrix<U> &lhs, const V &rhs)
    {
      Matrix<U> res = lhs;
      std::transform(res.data_.begin(), res.data_.end(),
                      res.data_.begin(),
                      std::bind1st(std::multiplies<U>(), static_cast<U>(rhs)));
      return res;
    }
    template<typename U, typename V>
    Matrix<V> operator*(const U &lhs, const Matrix<V> &rhs)
    {
      return rhs * lhs;
    }

    template<typename T>
    std::ostream& operator<<(std::ostream &ostr, const Matrix<T> &m)
    {
      std::vector<int> dim_wgts(m.dims_);
      //std::vector<int> dim_wgts(m.dims_.size(), 0);
      //std::partial_sum(m.dims_.rbegin(), m.dims_.rend(), dim_wgts.rbegin(), std::multiplies<int>());
      std::vector<int> dim_wgts_counter(dim_wgts.size(), 0);
      assert(dim_wgts.size() > 0);
      std::size_t dim_wgts_counter_idx = dim_wgts.size() - 1;
      for(typename std::vector<T>::const_iterator citer = m.data_.cbegin();
          citer != m.data_.cend(); ++citer){
        ostr << *citer << ", ";
        dim_wgts_counter[dim_wgts_counter_idx]++;
        if(dim_wgts_counter[dim_wgts_counter_idx] >=
            dim_wgts[dim_wgts_counter_idx]){
          while(dim_wgts_counter[dim_wgts_counter_idx]>=dim_wgts[dim_wgts_counter_idx]){
            dim_wgts_counter[dim_wgts_counter_idx] = 0;
            ostr << "\n";
            if(dim_wgts_counter_idx > 0){
              dim_wgts_counter[dim_wgts_counter_idx - 1]++;
              dim_wgts_counter_idx--;
            }
            else{
              break;
            }
          }

          dim_wgts_counter_idx = dim_wgts.size() - 1;
        }
      }
      return ostr;
    }

  } // namespace MapperUtil
} // namespace ESMCI

#endif //ESMCI_Mat_H

max_stars_repo_name: joeylamcy/gchp

max_stars_repo_licenses: NCSA Apache-2.0 MIT

max_issues_repo_name: joeylamcy/gchp

max_issues_repo_licenses: NCSA Apache-2.0 MIT

max_forks_repo_name: joeylamcy/gchp

max_forks_repo_licenses: NCSA Apache-2.0 MIT


text:
/*
Copyright 2010-2011, D. E. Shaw Research.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

* Redistributions of source code must retain the above copyright
  notice, this list of conditions, and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions, and the following disclaimer in the
  documentation and/or other materials provided with the distribution.

* Neither the name of D. E. Shaw Research nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#ifndef __r123_compat_gslrng_dot_h__
#define __r123_compat_gslrng_dot_h__

#include <gsl/gsl_rng.h>
#include <string.h>

/**
   The macro:  GSL_CBRNG(NAME, CBRNGNAME)
   declares the necessary structs and  constants that define a
   gsl_rng_NAME type based on the counter-based RNG CBRNGNAME.  For example:

   Usage:

   @code
   #include <Random123/threefry.h>
   #include <Random123/conventional/gsl_cbrng.h>  // this file
   GSL_CBRNG(cbrng, threefry4x32); // creates gsl_rng_cbrng

   int main(int argc, char **argv){
       gsl_rng *r = gsl_rng_alloc(gsl_rng_cbrng);
       ... use r as you would use any other gsl_rng ...
    }
    @endcode

    It requires that NAME be the name of a CBRNG that follows the
    naming and stylistic conventions of the Random123 library.

    Note that wrapping a \ref CBRNG "counter-based PRNG" with a traditional API in
    this way obscures much of the power of the CBRNG API.
    Nevertheless, it may be of value to applications that are already
    coded to work with GSL random number generators, and that wish
    to use the RNGs in the Random123 library.

 */

#define GSL_CBRNG(NAME, CBRNGNAME)                                      \
const gsl_rng_type *gsl_rng_##NAME;                                     \
                                                                        \
typedef struct{                                                         \
    CBRNGNAME##_ctr_t ctr;                                                   \
    CBRNGNAME##_ctr_t r;                                                     \
    CBRNGNAME##_key_t key;                                                   \
    int elem;                                                           \
} NAME##_state;                                                         \
                                                                        \
static unsigned long int NAME##_get(void *vstate){                      \
    NAME##_state *st = (NAME##_state *)vstate;                          \
    const int N=sizeof(st->ctr.v)/sizeof(st->ctr.v[0]);                 \
    if( st->elem == 0 ){                                                \
        ++st->ctr.v[0];                                                 \
        if( N>1 && st->ctr.v[0] == 0 ) ++st->ctr.v[1];                  \
        if( N>2 && st->ctr.v[1] == 0 ) ++st->ctr.v[2];                  \
        if( N>3 && st->ctr.v[2] == 0 ) ++st->ctr.v[3];                  \
        st->r = CBRNGNAME(st->ctr, st->key);                                 \
        st->elem = N;                                                   \
    }                                                                   \
    return st->r.v[--st->elem];                                         \
}                                                                       \
                                                                        \
static double NAME##_get_double (void * vstate);                        \
                                                                        \
static void NAME##_set(void *vstate, unsigned long int s){              \
    NAME##_state *st = (NAME##_state *)vstate;                          \
    st->elem = 0;                                                       \
    /* Assume that key and ctr have an array member, v,                 \
       as if they are r123arrayNxW.  If not, this will fail             \
       to compile.  In particular, this macro fails to compile          \
       when the underlying CBRNG requires use of keyinit */             \
    memset(&st->ctr.v[0], 0, sizeof(st->ctr.v));                        \
    memset(&st->key.v[0], 0, sizeof(st->key.v));                        \
    /* GSL 1.15 documentation says this about gsl_rng_set:              \
         Note that the most generators only accept 32-bit seeds, with higher \
         values being reduced modulo 2^32.  For generators with smaller \
         ranges the maximum seed value will typically be lower.         \
     so we won't jump through any hoops here to deal with               \
     high bits if sizeof(unsigned long) > sizeof(uint32_t). */          \
    st->key.v[0] = s;                                                   \
}                                                                       \
                                                                        \
static const gsl_rng_type NAME##_type = {                               \
    #NAME,                                                              \
    ~0UL>>((R123_W(CBRNGNAME##_ctr_t)>=8*sizeof(unsigned long))? 0 : (8*sizeof(unsigned long) - R123_W(CBRNGNAME##_ctr_t))),     \
    0,                                                                  \
    sizeof(NAME##_state),                                               \
    &NAME##_set,                                                        \
    &NAME##_get,                                                        \
    &NAME##_get_double                                                  \
};                                                                      \
                                                                        \
static double                                                           \
NAME##_get_double (void * vstate)                                       \
{                                                                       \
    return NAME##_get (vstate)/(double)NAME##_type.max;                 \
}                                                                       \
                                                                        \
const gsl_rng_type *gsl_rng_##NAME = &NAME##_type

#endif

max_stars_repo_name: maumueller/rehashing

max_issues_repo_name: maumueller/rehashing

max_forks_repo_name: maumueller/rehashing


text:
/*! \file allvars.h
 *  \brief declares global variables and structures.
 *
 *  This file declares all global variables and structures. Further variables should be added here, and declared as
 *  \e \b extern. The actual existence of these variables is provided by the file \ref allvars.cxx. To produce
 *  \ref allvars.cxx from \ref allvars.h, do the following:
 *
 *     \arg Erase all \#define's, typedef's, and enum's
 *     \arg add \#include "allvars.h", delete the \#ifndef ALLVARS_H conditional
 *     \arg delete all keywords 'extern'
 *     \arg delete all struct definitions enclosed in {...}, e.g.
 *        "extern struct global_data_all_processes {....} All;"
 *        becomes "struct global_data_all_processes All;"
 */

#ifndef ALLVARS_H
#define ALLVARS_H
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <cmath>
#include <string>
#include <getopt.h>
#include <sys/stat.h>
#include <sys/timeb.h>

//gsl stuff
#include <gsl/gsl_heapsort.h>

//nbody code
#include <NBody.h>
#include <NBodyMath.h>
#include <KDTree.h>

#ifdef USEOPENMP
#include <omp.h>
#endif

#ifdef USEMPI
#include <mpi.h>
///Includes external global variables used for MPI version of code
#include "mpivar.h"
#endif

#include "../../src/gadgetitems.h"
#include "../../src/tipsy_structs.h"

//for morphological classifications (where tidal tails, or completely disrupted)
#define NUMMORPHCLASS 3;
#define QTIDAL 0.60
#define STIDAL 0.40
#define QLIMTIDAL 0.30
#define SLIMTIDAL 0.25
#define VFTIDAL 0.85
#define VFLIMTIDAL 0.50
#define V2TIDAL 0.6
#define V3TIDAL 0.4
#define V2LIMTIDAL 0.35
#define V3LIMTIDAL 0.30

#define ETIDAL 0.90
#define ELIMTIDAL 0.50

//number of bands in which to calculate magnitude
#define NBANDS 9

///define the particle types
//@{
#define NPARTTYPES 5
#define ALLTYPE 4
#define GASTYPE 0
#define DMTYPE 1
#define STARTYPE 2
#define BHTYPE 3
//@}

///\todo need to alter how stats are stored if one can update the number of quantiles stored
///the number of values that is used to characterize a bulk distribution
///where idea is 0=mean,1=sd,2-4=quantiles of 0.16,0.5,0.84
#define NDISTRIB 5
///max number of quantiles
#define MAXNQUANTS 11

///max number of radial bins
#define NRAD 1000

/// Structures and external variables

///for stars,gas, define lower limits for Z and sfr
#define MINGASZ 1e-7
#define MINSTARZ 1e-7
#define MINSFR 1e-3

///define sphi io block points for gadget
//@{
///self energy
#define sphioblock_u 0
///density
#define sphioblock_d 1
#ifdef SPHCOOLING
///mean molecular weight
#define sphioblock_mmw 2
///neutral hydrogen abundance
#define sphioblock_nha 3
#endif
///smoothing lengths
#define sphioblock_h
//@}

///define different reference frames
//@{
///cm ref
#define ICMFRAME 1
///particle with largest potential
#define IPOTFRAME 2
///most bound particle based in input data
#define IMBFRAME 3
//@}

///define omp number of elements threshold
#define ompunbindnum 1000
#define omppropnum 50000

///minimum number of particles for properties to be calculated
#define MINNPROP 10
///minimum number of particles for a profile to be calculated
#define MINNPROFILE 100
///minimum mass fraction of candidate structure for enclosed quantities to be useful
#define MINMASSFRAC 0.05

///some useful constants for gas
//@{
///mass of helium relative to hydrogen
#define M_HetoM_H 4.0026
//@}

using namespace std;
using namespace Math;
using namespace NBody;

///for potential calculation
struct PotInfo
{
    ///bucket size for tree calculation
    int BucketSize;
    ///for tree potential calculation;
    Double_t TreeThetaOpen;
    ///gravitational softening (here simple plummer)
    Double_t eps;
    ///\todo When calculating the thermal internal energy of gas particles
    ///I should probably include some temperature floors and ceilings
    ///and perhaps some sound speed constrains. But I should produce
    ///properties of the baryonic material for all, gravitationally self-bound and grav+therm bound

    PotInfo(){
        BucketSize=8;
        TreeThetaOpen=0.7;
        eps=0.01;
    }
};

///Options structure
struct Options
{
    ///filenames
    //@{
    ///input particle data and output name
    char *fname,*outname;
    ///halo catalog base name
    char *halocatname;
    ///optional config file name
    char *configname;
    ///optional stellar synthesis model file name
    char *magname;
    //@}
    ///number of files per particle data snapshot, if zero then tipsy file
    int nfiles;
    ///number of files per halo catalog
    int nfileshalocat;
    ///number of sph blocks in gadget file
    int sphnio;
    ///defining units and cosmology
    //@{
    ///length,m,v,grav units pressure unit, temperature
    Double_t L, M, V, G, PUnit, TUnit;
    ///scalbodydata/utils/haloanalcode/e factor, Hubunit, h and cosmology
    Double_t a,H,h, Omega_m, Omega_Lambda,rhoc;
    ///softening length and period (comove)
    Double_t eps, p;
    ///background density and virial level
    Double_t rhobg, virlevel;
    int comove;
    ///for gas, mass fraction of hydrogen
    Double_t Xh;
    ///for gas, boltzmann constant
    Double_t Boltzmann;
    ///for gas, number of hydrogen atoms per sph particle
    Double_t N_H;
    ///for gas, correction for helium to mass of sph particle /mh ie: Number of particles = M_p / (M_H*Xh+M_He*(1-Xh)) or even M_H*Xh+M_He*(1-Xh)*Y_Li+M_Li*(1-Xh)*(1-Y_Li) ...
    Double_t N_He_corr;
    ///for gas, relation of specific internal energy to Temperature (u=3/2*k_B*(N_H*(Xh+(1-Xh)*N_He)*T
    Double_t utoT;
    ///for gas, temperature threshold when calculating X-ray based spectroscopic-like temperature)
    Double_t Tempthreshold;
    ///\todo may want some specific stuff for stars too
    //@}

    ///store number of particles
    Int_t npart[NPARTTYPES];

    ///structure that contains variables for potentialcalc
    PotInfo pinfo;
    Double_t error;

    ///quantile points for characterization of distribution
    int nquants;
    Double_t quants[MAXNQUANTS];

    ///minimum number a halo has to have to be analyzed
    Int_t minnum;
    ///Kinetic/Potential Energy ratio
    Double_t TVratio;
    ///number of mpi files that were written in parallel that must be read for group_catalog files
    int mpinum;

    ///for calculaion of CM
    Double_t cmfrac, cmadjustfac;
    PotInfo uinfo;

    ///for calculating kinetic energy to a reference frame
    int reframe;

    ///if no mass is stored
    Double_t MassValue;

    ///if zoom simulations, store effective resolution
    long int Neff;

    ///for verbose output
    int verbose;

    ///store whether potential has been already calculated as object has been moved to particle with deepest potential well.
    int ipotcalc;

    ///used to calculate kinetic reference frame about particle at bottom of potential well
    Int_t Nbref;
    Double_t fracbref;

    ///input STF/VELOCIraptor format flags
    //@{
    int iseparatefiles,ibinaryin;
    //@}

    ///output format
    int ibinaryout;

    Options()
    {
        fname=outname=halocatname=NULL;
        nfiles=1;
        sphnio=2;
        nfileshalocat=1;

        L = 1.0;
        M = 1.0;
        V = 1.0;
        PUnit=5.588e-14;//(assumes T in kelvin) makes pressure into ergs/cm^3 assuming standard gadget units (this is boltzmann's constant / proton mass *Massunit_g/(Lengthunit_cm)^3
        TUnit=1.0;
        MassValue=1.0;

        p = 0.0;
        a = 1.0;
        H = 0.1;
        h = 1.0;
        Omega_m = 1.0;
        Omega_Lambda = 0.0;
        rhobg = 1.0;
        virlevel = 50.;
        comove=0;
        G = 43.01;//for G in units of (km/s)^2 Mpc/1e10 Msun
        rhoc=27.753;//for same units as above, of course missing factors of h^2
        G = 43021.1349;//for G in units of (km/s)^2 kpc/1e10 Msun
        rhoc=2.7753e-8;

        Xh=0.76;
        Boltzmann=log(6.94107285)-70.0*log(10.0); //ln of boltzmann in units of K^(-1) (km/s)^2 (1e10 Msun)
        N_H=log(1.1892)+67.0*log(10.0); //number of hydrogen atoms per 10^10 Msun
        N_He_corr=-log(Xh+M_HetoM_H*(1.0-Xh));
        utoT=log(1.5)+Boltzmann+N_H+N_He_corr;
        Tempthreshold=7.0*log(10.0);

        error=1e-3;

        nquants=3;
        quants[0]=0.16;quants[1]=0.5;quants[2]=0.84;
        for (int i=0;i<NPARTTYPES;i++) npart[i]=0;

        minnum=20;
        TVratio=-1.0;
        mpinum=0;

        cmfrac=0.1;
        cmadjustfac=0.7;
        Neff=-1;

        reframe=IPOTFRAME;
        verbose=0;

        ipotcalc=1;
        Nbref=10;
        fracbref=0.05;

        iseparatefiles=0;
        ibinaryin=0;
        ibinaryout=0;

    }
};

///Radial profile data in spherical shells
struct RadPropData
{
    ///To store radial profiles
    //@{
    Int_t nbins;
    int ptype;
    Double_t *numinbin;
    Double_t *radval;
    Double_t *Menc;
    Coordinate *velenc;
    Matrix *sigmavenc;
    Coordinate *Jenc;
    Double_t *qenc;
    Double_t *senc;
    Matrix *eigvecenc;
    ///primarily for sph particles
    Double_t *Den;
    ///gas/star particle info
    Double_t *Tempenc;
    Double_t *Tdispenc;
    Double_t *Temenc;
    Double_t *Tslenc;
#ifdef RADIATIVE
    Double_t *Zmetenc;
    Double_t *Tageenc;
#endif
//@}
    RadPropData(){
        //for (int i=0;i<NRAD;i++) {numinbin[i]=radval[i]=Menc[i]=0.;Jenc[i]=Coordinate(0,0,0);velenc[i]=Coordinate(0,0,0);qenc[i]=senc[i]=1.0;}
        ptype=DMTYPE;
        numinbin=NULL;radval=NULL;Menc=NULL;velenc=NULL;sigmavenc=NULL;Jenc=NULL;qenc=NULL;senc=NULL;eigvecenc=NULL;
        Den=NULL;Tempenc=NULL;Tdispenc=NULL;Temenc=NULL;Tslenc=NULL;
#ifdef RADIATIVE
        Zmetenc=NULL;
        Tageenc=NULL;
#endif
    }
    void SetBins(int N) {
        nbins=N;
        if(numinbin!=NULL) delete[] numinbin;numinbin=new Double_t[nbins];
        if(radval!=NULL) delete[] radval;radval=new Double_t[nbins];
        if(Menc!=NULL) delete[] Menc;Menc=new Double_t[nbins];
        if(velenc!=NULL) delete[] velenc;velenc=new Coordinate[nbins];
        if(sigmavenc!=NULL) delete[] sigmavenc;sigmavenc=new Matrix[nbins];
        if(Jenc!=NULL) delete[] Jenc;Jenc=new Coordinate[nbins];
        if(qenc!=NULL) delete[] qenc;qenc=new Double_t[nbins];
        if(senc!=NULL) delete[] senc;senc=new Double_t[nbins];
        if(eigvecenc!=NULL) delete[] eigvecenc;eigvecenc=new Matrix[nbins];
        if(Tempenc!=NULL) delete[] Tempenc;Tempenc=new Double_t[nbins];
        if(Tdispenc!=NULL) delete[] Tdispenc;Tdispenc=new Double_t[nbins];
        if(Temenc!=NULL) delete[] Temenc;Temenc=new Double_t[nbins];
        if(Tslenc!=NULL) delete[] Tslenc;Tslenc=new Double_t[nbins];
        if(Den!=NULL) delete[] Den;Den=new Double_t[nbins];
#ifdef RADIATIVE
        if(Zmetenc!=NULL) delete[] Zmetenc;Zmetenc=new Double_t[nbins];
        if(Tageenc!=NULL) delete[] Tageenc;Tageenc=new Double_t[nbins];
#endif
        for (int i=0;i<nbins;i++) {numinbin[i]=radval[i]=Menc[i]=0;Jenc[i]=Coordinate(0,0,0);velenc[i]=Coordinate(0,0,0);qenc[i]=senc[i]=1.0;}
        for (int i=0;i<nbins;i++) {Tempenc[i]=Tdispenc[i]=Den[i]=Temenc[i]=Tslenc[i]=0.;}
    }
};

///Cylindrical profile data
struct CylPropData
{
    ///reference axis of cylinder
    Coordinate zref;
    ///To store cylindrical profiles
    //@{
    Int_t nbins;
    Double_t *numinbin;
    Double_t *radval;
    Double_t *Menc;
    Coordinate *Jenc;
    Coordinate *velenc;
    Matrix *sigmavenc;
    Double_t *zmean;
    Double_t *Rmean;
    Double_t *zdisp;
    Double_t *Rdisp;
    Double_t *qenc;
    Double_t *senc;
    Matrix *eigvecenc;
    ///primarily for sph particles
    Double_t *Den;
    ///gas/star particle info
    Double_t *Tempenc;
    Double_t *Tdispenc;
    Double_t *Temenc;
    Double_t *Tslenc;
#ifdef RADIATIVE
    Double_t *Zmetenc;
    Double_t *Tageenc;
#endif
    //@}
    CylPropData(){
        numinbin=NULL;radval=NULL;Menc=NULL;velenc=NULL;sigmavenc=NULL;Jenc=NULL;zmean=NULL;Rmean=NULL;zdisp=NULL;Rdisp=NULL;qenc=NULL;senc=NULL;eigvecenc=NULL;
        Den=NULL;Tempenc=NULL;Tdispenc=NULL;Temenc=NULL;Tslenc=NULL;
#ifdef RADIATIVE
        Zmetenc=NULL;
        Tageenc=NULL;
#endif
    }
    void SetBins(int N) {
        nbins=N;
        if(numinbin!=NULL) delete[] numinbin;numinbin=new Double_t[nbins];
        if(radval!=NULL) delete[] radval;radval=new Double_t[nbins];
        if(Menc!=NULL) delete[] Menc;Menc=new Double_t[nbins];
        if(velenc!=NULL) delete[] velenc;velenc=new Coordinate[nbins];
        if(sigmavenc!=NULL) delete[] sigmavenc;sigmavenc=new Matrix[nbins];
        if(Jenc!=NULL) delete[] Jenc;Jenc=new Coordinate[nbins];
        if(zmean!=NULL) delete[] zmean;zmean=new Double_t[nbins];
        if(Rmean!=NULL) delete[] Rmean;Rmean=new Double_t[nbins];
        if(zdisp!=NULL) delete[] zdisp;zdisp=new Double_t[nbins];
        if(Rdisp!=NULL) delete[] Rdisp;Rdisp=new Double_t[nbins];
        if(qenc!=NULL) delete[] qenc;qenc=new Double_t[nbins];
        if(senc!=NULL) delete[] senc;senc=new Double_t[nbins];
        if(eigvecenc!=NULL) delete[] eigvecenc;eigvecenc=new Matrix[nbins];
        if(Tempenc!=NULL) delete[] Tempenc;Tempenc=new Double_t[nbins];
        if(Tdispenc!=NULL) delete[] Tdispenc;Tdispenc=new Double_t[nbins];
        if(Temenc!=NULL) delete[] Temenc;Temenc=new Double_t[nbins];
        if(Tslenc!=NULL) delete[] Tslenc;Tslenc=new Double_t[nbins];
        if(Den!=NULL) delete[] Den;Den=new Double_t[nbins];
#ifdef RADIATIVE
        if(Zmetenc!=NULL) delete[] Zmetenc;Zmetenc=new Double_t[nbins];
        if(Tageenc!=NULL) delete[] Tageenc;Tageenc=new Double_t[nbins];
#endif
        for (int i=0;i<nbins;i++) {numinbin[i]=radval[i]=Menc[i]=0.;Jenc[i]=velenc[i]=Coordinate(0,0,0);qenc[i]=senc[i]=1.0;zmean[i]=Rmean[i]=zdisp[i]=Rdisp[i]=0.;}
        for (int i=0;i<nbins;i++) {Tempenc[i]=Tdispenc[i]=Den[i]=Temenc[i]=Tslenc[i]=0.;}
    }
};

///Based on PropData in \ref allvars.h of STF/VELOCIraptor
///Here the data also stores the type of particles it corresponds to
struct PropData
{
    ///type of structure, dm, gas, star, etc
    Int_t ptype;
    ///number of particles
    Int_t num;
    ///centre of mass
    Coordinate gcm, gcmvel;
    ///Position of most bound particle
    Coordinate gpos, gvel;
    ///physical properties regarding mass, size
    Double_t gmass,gsize,gMvir,gRvir,gRmbp,gmaxvel,gRmaxvel,gMmaxvel;
    ///physical properties for shape/mass distribution
    Double_t gq,gs;
    Matrix geigvec;
    ///physical properties for velocity
    Double_t gsigma_v;
    Matrix gveldisp;
    ///physical properties for dynamical state
    Double_t Efrac,Pot,T;
    ///physical properties for dynamical state for specific particle types
    Double_t Efractyped[NPARTTYPES],Pottyped[NPARTTYPES],Ttyped[NPARTTYPES];
    ///physical properties for angular momentum
    Coordinate gJ;

    ///metallicity, temperature, pressure,etc
    Double_t Temp[NDISTRIB];
#ifdef RADIATIVE
    Double_t Zmet[NDISTRIB];
    ///stellar age
    Double_t Tage[NDISTRIB];
    ///luminosity, which for gas which would probably be Infrared and radio bands
    Double_t LUM[NBANDS];
#endif

    ///To store cylindrical profiles
    RadPropData radprofile;
    ///To store cylindrical profiles
    CylPropData cylprofile;

    PropData(){
        num=0;
        gmass=gsize=gRmaxvel=Efrac=Pot=T=gcm[0]=gcm[1]=gcm[2]=gcmvel[0]=gcmvel[1]=gcmvel[2]=0.;
        gveldisp=Matrix(0.);
        gq=gs=1.0;
        gRmbp=0.0;

        for (int i=0;i<NDISTRIB;i++)Temp[i]=0;
#ifdef RADIATIVE
        for (int i=0;i<NDISTRIB;i++)Zmet[i]=Tage[i]=0.;
        for (int i=0;i<NBANDS;i++)LUM[i]=0;
#endif
    }
};

///used to store halo particle id data
struct HaloParticleData{
    ///store the halo id
    long unsigned haloID;
    ///number of particles, eventually used to store bound number
    long unsigned NumberofParticles;
    ///number of particles of different types, eventually used to store bound number
    Int_t NumofType[NPARTTYPES];
    ///number of all particles
    long unsigned AllNumberofParticles;
    ///number of all particles of different types
    Int_t AllNumofType[NPARTTYPES];

    ///offset of halo's particle in list of all particles in structures
    Int_t noffset;
    ///array of particle ids
    long unsigned *ParticleID;
    HaloParticleData(long unsigned numinhalo=0){
        AllNumberofParticles=NumberofParticles=numinhalo;
        if (NumberofParticles>0) {
            ParticleID=new long unsigned[NumberofParticles];
        }
        for (int i=0;i<NPARTTYPES;i++) AllNumofType[i]=NumofType[i]=0;
    }
    void Alloc(long unsigned ninhalos=0){
        if (AllNumberofParticles>0){
            delete[] ParticleID;
        }
        AllNumberofParticles=NumberofParticles=ninhalos;
        if (AllNumberofParticles>0) {
            ParticleID=new long unsigned[NumberofParticles];
        }
    }
    ~HaloParticleData(){
        if (NumberofParticles>0){
            delete[] ParticleID;
        }
    }
};

///External types for tipsy
//@{
#define TGASTYPE 0;
#define TDARKTYPE 1;
#define TSTARTYPE 2;
//@}

///note that for grouped particles type is 10+TYPE.
#define SUBSTRUCTTYPE 10;

#endif

max_stars_repo_name: broukema/VELOCIraptor-STF

max_issues_repo_name: broukema/VELOCIraptor-STF

max_forks_repo_name: broukema/VELOCIraptor-STF


text:
// Copyright (c) 2015-2016, Massachusetts Institute of Technology
// Copyright (c) 2016-2017 Sandia Corporation
// Copyright (c) 2017 NTESS, LLC.

// This file is part of the Compressed Continuous Computation (C3) Library
// Author: Alex A. Gorodetsky
// Contact: alex@alexgorodetsky.com

// All rights reserved.

// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:

// 1. Redistributions of source code must retain the above copyright notice,
//    this list of conditions and the following disclaimer.

// 2. Redistributions in binary form must reproduce the above copyright notice,
//    this list of conditions and the following disclaimer in the documentation
//    and/or other materials provided with the distribution.

// 3. Neither the name of the copyright holder nor the names of its contributors
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.

// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

//Code

#ifndef LINALG_H
#define LINALG_H

#ifdef __APPLE__
    #include <Accelerate/Accelerate.h>
    /* #include "/System/Library/Frameworks/Accelerate.framework/Versions/Current/Frameworks/vecLib.framework/Headers/clapack.h" */

    #define dgetri_(X, Y, Z, A , B, C, D ) \
            ( dgetri_( (__CLPK_integer *) X, Y, (__CLPK_integer *) Z, (__CLPK_integer *) A, \
                        B, (__CLPK_integer *)C , (__CLPK_integer *) D) )
    #define dgetrf_(X, Y, Z, A ,B, C ) \
            ( dgetrf_( (__CLPK_integer *) X,(__CLPK_integer *) Y, Z, (__CLPK_integer *) A, \
                       (__CLPK_integer *) B, (__CLPK_integer *)C ))
    #define dorgqr_(X,Y,Z,A,B,C,D,E,F) \
            ( dorgqr_( (__CLPK_integer *) X, (__CLPK_integer *) Y, (__CLPK_integer *) Z, A, \
                       (__CLPK_integer *) B, C , D, (__CLPK_integer *) E, (__CLPK_integer *)F) )
    #define dgeqrf_(X, Y, Z, A , B, C, D, E ) \
            ( dgeqrf_( (__CLPK_integer *) X,  (__CLPK_integer *)Y, Z, (__CLPK_integer *) A, \
                        B, C, (__CLPK_integer *) D, (__CLPK_integer *) E) )

    #define dorgrq_(X,Y,Z,A,B,C,D,E,F) \
            ( dorgrq_( (__CLPK_integer *) X, (__CLPK_integer *) Y, (__CLPK_integer *) Z, A, \
                       (__CLPK_integer *) B, C , D, (__CLPK_integer *) E, (__CLPK_integer *)F) )
    #define dgerqf_(X, Y, Z, A , B, C, D, E ) \
            ( dgerqf_( (__CLPK_integer *) X,  (__CLPK_integer *)Y, Z, (__CLPK_integer *) A, \
                        B, C, (__CLPK_integer *) D, (__CLPK_integer *) E) )

    #define dgesdd_(X,Y,Z,A,B,C,D,E,F,G,H,I,J,K) \
            ( dgesdd_(X, (__CLPK_integer *)Y, (__CLPK_integer *)Z, A, (__CLPK_integer *) B,\
                C,D,(__CLPK_integer *) E, F, (__CLPK_integer *) G, H, (__CLPK_integer *) I,\
                (__CLPK_integer *) J, (__CLPK_integer *) K ) )

    #define dgebal_(X, Y, Z, A , B, C, D, E ) \
            ( dgebal_( X,  (__CLPK_integer *)Y, Z, (__CLPK_integer *) A, \
                   (__CLPK_integer *) B, (__CLPK_integer *)C, D, (__CLPK_integer *) E) )

    #define dpotrf_(X,Y,Z,A,B) \
            ( dpotrf_(X, (__CLPK_integer *)Y, Z, (__CLPK_integer *) A, (__CLPK_integer *) B ))

    #define dpotri_(X,Y,Z,A,B) \
            ( dpotri_(X, (__CLPK_integer *)Y, Z, (__CLPK_integer *) A, (__CLPK_integer *) B ))

    #define dtrtri_(X,Y,Z,A,B,C) \
            ( dtrtri_(X,Y, (__CLPK_integer *)Z, A, (__CLPK_integer *) B, (__CLPK_integer *) C ))

    #define dgesv_(X,Y,Z,A,B,C,D,E) \
            ( dgesv_((__CLPK_integer *)X, (__CLPK_integer *)Y, Z, (__CLPK_integer *) A, (__CLPK_integer *) B, \
                C, (__CLPK_integer *) D, (__CLPK_integer *) E))

    #define dhseqr_(X, Y, Z, A , B, C, D, E,F,G,H,I,J,K ) \
            ( dhseqr_( X, Y, (__CLPK_integer *)Z, (__CLPK_integer *)A, \
                   (__CLPK_integer *) B, C, (__CLPK_integer *)D, E, F, G, (__CLPK_integer *) H, \
                   I, (__CLPK_integer *) J, (__CLPK_integer *)K ) )

    #define dsyev_(A,B,C,D,E,F,G,H,J) \
        ( dsyev_(A,B,(__CLPK_integer *) C,D,(__CLPK_integer *) E,F, \
            G, (__CLPK_integer *) H, (__CLPK_integer *) J) )

    #define dgeev_(X, Y, Z, A , B, C, D, E,F,G,H,I,J,K ) \
            ( dgeev_( X, Y, (__CLPK_integer *)Z, A, \
                   (__CLPK_integer *) B, C, D, E, (__CLPK_integer *) F, G, (__CLPK_integer *) H, \
                   I, (__CLPK_integer *) J, (__CLPK_integer *)K ) )

    #define dstegr_(X,Y,Z,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q) \
            ( dstegr_(X,Y,(__CLPK_integer *)Z,A,B,C,D,\
                   (__CLPK_integer *)E,(__CLPK_integer *)F,G,(__CLPK_integer *)H, \
                   I,J,(__CLPK_integer *)K,(__CLPK_integer *)L,M,(__CLPK_integer *)N, \
                   (__CLPK_integer *)O,(__CLPK_integer *)P,(__CLPK_integer *)Q))

    #define dgelsd_(X,Y,Z,A,B,C,D,E,F,G,H,I,J,K)                         \
            ( dgelsd_( (__CLPK_integer *)X, (__CLPK_integer *)Y, (__CLPK_integer *)Z, \
                       (__CLPK_doublereal *)A, (__CLPK_integer *)B, (__CLPK_doublereal *)C,  \
                       (__CLPK_integer *)D, (__CLPK_doublereal *) E, \
	                   (__CLPK_doublereal *)F, (__CLPK_integer *)G, (__CLPK_doublereal *)H, \
                       (__CLPK_integer *)I, (__CLPK_integer *)J, (__CLPK_integer *)K))
#else
    /* #include <gsl/gsl_cblas.h> */
    #include <cblas.h>

void dgetri_(int * X, double *Y, int * Z, int * A, double *B, int *C , int * D);
void dgetrf_(int * X,int * Y, double*Z, int * A, int * B, int *C );
void dorgqr_(int * X, int * Y, int * Z, double *A,int * B, double *C , double *D, int * E, int *F);
void dgeqrf_(int *X, int *Y, double *Z, int * A, double * B, double * C, int * D, int * E);
void dgeev_(char * X, char *Y, int * Z, double *A, int * B, double * C, double *D, double *E,
            int * F, double *G, int * H, double *, int *, int *K);
void dgebal_(char *X,  int *Y, double *Z, int * A, int * B, int *C, double *D, int * E);
void dhseqr_(char *X, char *Y, int *Z, int *A,int * B, double *C, int *D,
             double *E, double *F, double *G, int * H, double *, int *, int *K );
void dorgrq_(int * X, int * Y, int * Z, double *A,int * B, double *C , double *D, int * E, int *F);
void dgerqf_(int * X,  int *Y, double *Z, int * A, double *B, double *C, int * D, int * E);
void dgesdd_(char *X, int * Y, int *Z, double *A, int * B,double * C,double * D,
             int * E,double * F, int * G, double *H, int *, int *, int * K );
void dstev_(char *X,  int *Y, double *Z, double *A, double *B, int *C, double *D, int * E);
void dsyev_(char *A,char *B,int * C,double *D,int * E,double *F, double *G, int* H, int * J);
void dpotrf_(char *X, int*Y, double *Z, int * A, int * B );
void dpotri_(char *X, int*Y, double *Z, int * A, int* B);
void dtrtri_(char *X,char*Y, int *Z, double *A, int * B, int * C);
void dgesv_(int *X, int *Y, double*Z, int * A, int * B,double*C, int * D, int * E);
void dgelsd_(int *X, int *Y, int *Z, double *A, int *B, double *C,
             int *D, double * E, double *F, int *G, double *H, int *,int *, int *K);

#endif

#include "matrix_util.h"

void c3linalg_multiple_vec_mat(size_t, size_t, size_t, const double *, size_t,
                               const double *, size_t, double *,size_t);
void c3linalg_multiple_mat_vec(size_t, size_t, size_t, const double *, size_t,
                               const double *, size_t, double *,size_t);
int qr(size_t, size_t, double *, size_t);
void rq_with_rmult(size_t, size_t, double *, size_t, size_t, size_t, double *, size_t);
void svd(size_t, size_t, size_t, double *, double *, double *, double *);
size_t truncated_svd(size_t, size_t, size_t, double *, double **, double **, double **, double);
size_t pinv(size_t, size_t, size_t, double *, double *, double);

double norm2(double *, int);
double norm2diff(double *, double *, int);
double mean(double *, size_t);
double mean_size_t(size_t *, size_t);

struct mat * kron(const struct mat *, const struct mat *);
void kron_col(int, int, double *, int, int, int, double *, int, double *, int);
void vec_kron(size_t, size_t, double *, size_t, size_t, size_t,
        double *, size_t, double *, double, double *);
void vec_kronl(size_t, size_t, double *, size_t, size_t, size_t,
        double *, size_t, long double *, double, long double *);

// decompositions
struct fiber_list{
    size_t index;
    double * vals;
    struct fiber_list * next;
};
struct fiber_info{
    size_t nfibers;
    struct fiber_list * head;
};
void AddFiber(struct fiber_list **, size_t, double *, size_t);
int IndexExists(struct fiber_list *, size_t);
double * getIndex(struct fiber_list *, size_t);
void DeleteFiberList(struct fiber_list **);

struct sk_decomp {
    size_t n;
    size_t m;
    size_t rank;
    size_t * rows_kept;
    size_t * cols_kept;
    size_t cross_rank;
    double * cross_inv;
    struct fiber_info * row_vals;
    struct fiber_info * col_vals;
    int success;
};

void init_skf(struct sk_decomp **, size_t, size_t, size_t);
void sk_decomp_to_full(struct sk_decomp *, double *);
void free_skf(struct sk_decomp **);

/* int comp_pivots(const double *, int, int, int *); */
int maxvol_rhs(const double *, size_t, size_t, size_t *, double *); //
int skeleton(double *, size_t, size_t, size_t, size_t *, size_t *, double);
int skeleton_func(double (*A)(int,int, int, void*), void *, size_t,
                size_t, size_t, size_t *, size_t *, double);
int
skeleton_func2(int (*Ap)(double *, double, size_t, size_t, double *,
                void *),
                   void *, struct sk_decomp **,  double *, double *,
                   double);
void linear_ls(size_t, size_t, double *, double *, double *);
#endif


text:
#ifndef VKST_PLAT_FS_NOTIFY_WIN32_H
#define VKST_PLAT_FS_NOTIFY_WIN32_H

#ifndef VKST_PLAT_FS_NOTIFY_H
#error "Include fs_notify.h only"
#endif

#include "fs_notify.h"
#include <gsl.h>
#include <vector>

namespace plat {

class fs_notify final : public impl::fs_notify<fs_notify> {
public:
  class watch {
  public:
    OVERLAPPED overlapped{};
    plat::filesystem::path path{};
    notify_delegate delegate{};
    bool recursive{false};
    HANDLE handle{INVALID_HANDLE_VALUE};
    std::array<BYTE, 32 * 1024> buffer{};
    watch_id id{UINT32_MAX};
    bool stop{false};

    bool refresh(bool clear = false) noexcept;
    void notify(plat::filesystem::path changed_path, DWORD action) noexcept;

    watch(plat::filesystem::path p, notify_delegate d, bool r) noexcept
    : path{std::move(p)}, delegate{std::move(d)}, recursive{r} {}

    watch() = default;
    watch(watch const&) = delete;
    watch& operator=(watch const&) = delete;
    ~watch() noexcept;
  }; // struct watch

private:
  std::vector<gsl::unique_ptr<watch>> _watches;

  watch_id do_add(plat::filesystem::path path,
                  impl::fs_notify<fs_notify>::notify_delegate delegate,
                  bool recursive, std::error_code& ec) noexcept;

  void do_remove(watch_id id) noexcept;

  void do_tick() noexcept;

  friend class impl::fs_notify<fs_notify>;
}; // class fs_notify

} // namespace plat

#endif // VKST_PLAT_FS_NOTIFY_WIN32_H

max_stars_repo_name: wesleygriffin/vkst

max_issues_repo_name: wesleygriffin/vkst

max_forks_repo_name: wesleygriffin/vkst