text:
#include <time.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>

#include <uuid/uuid.h>

#include <openssl/sha.h>

#include <gsl/gsl_statistics_double.h>

#include "sig.h"
#include "pok.h"
#include "like.h"
#include "bn512.h"
#include "utils_like.h"

double print_trials_res(double * trials_res, int nb_trials)
{

    double mean = gsl_stats_mean(trials_res, 1, nb_trials);
    double stddev = gsl_stats_sd_m(trials_res, 1, nb_trials, mean);
    double confidence_lo = mean - ((1.96 * stddev) / sqrt(nb_trials));
    double confidence_hi = mean + ((1.96 * stddev) / sqrt(nb_trials));
    printf("    mean                   = %f\n", mean * 1000);
    printf("    stddev                 = %f\n", stddev * 1000);
    printf("    95%% confidence inteval = [%f, %f]\n\n", confidence_lo * 1000, confidence_hi * 1000);
    return mean * 1000;

}

void mesure_pairing(int nb_trials)
{

    printf("pairing : \n");

    mclBnG1 P;
    mclBnG2 Q;
    mclBnGT e;

    mclBn_init(MCL_BN462, MCLBN_COMPILED_TIME_VAR);

    mclBnG1_hashAndMapTo(&P, "abc", 3);
    mclBnG2_hashAndMapTo(&Q, "def", 3);

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        mclBn_pairing(&e, &P, &Q);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

}

void mesure_arithmetic_EC(int nb_trials)
{

    mclBnG1 P, xP;
    mclBnG2 Q, xQ;
    mclBnFr x;

    mclBn_init(MCL_BN462, MCLBN_COMPILED_TIME_VAR);

    mclBnFr_setByCSPRNG(&x);
    mclBnG1_hashAndMapTo(&P, "abc", 3);
    mclBnG2_hashAndMapTo(&Q, "def", 3);

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    printf("Scalar mult in G1 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        mclBnG1_mul(&xP, &P, &x);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }
    print_trials_res(trials_res, nb_trials);

    printf("Scalar mult in G2 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        mclBnG2_mul(&xQ, &Q, &x);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }
    print_trials_res(trials_res, nb_trials);

    printf("Add point in G1 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        mclBnG1_add(&xP, &P, &xP);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }
    print_trials_res(trials_res, nb_trials);

    printf("Add point in G2 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        mclBnG2_add(&xQ, &Q, &xQ);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }
    print_trials_res(trials_res, nb_trials);

}

void mesure_edd25519(int nb_trials)
{
    struct stat st = {0};
    size_t sig_len = ED25519_SIG_LENGTH;
    unsigned char * sig = (unsigned char *) malloc(sig_len * sizeof(unsigned char));;
    unsigned char msg[] = "Test sig";

    if (stat("./keys", &st) == -1) {
        mkdir("./keys", 0700);
    }

    sgen_ed25519("keys/pub.pem", "keys/priv.pem");

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    printf("sig ed25519 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ssig_ed25519("keys/priv.pem", msg, sizeof(msg), &sig, &sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

    printf("verif ed25519 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        sver_ed25519("keys/pub.pem", msg, sizeof(msg), sig, sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

    free(sig);

}

void mesure_sok(int nb_trials)
{

    const char * G2_basePoint_hexstr = "1 0257ccc85b58dda0dfb38e3a8cbdc5482e0337e7c1cd96ed61c913820408208f9ad2699bad92e0032ae1f0aa6a8b48807695468e3d934ae1e4df 1d2e4343e8599102af8edca849566ba3c98e2a354730cbed9176884058b18134dd86bae555b783718f50af8b59bf7e850e9b73108ba6aa8cd283 0a0650439da22c1979517427a20809eca035634706e23c3fa7a6bb42fe810f1399a1f41c9ddae32e03695a140e7b11d7c3376e5b68df0db7154e 073ef0cbd438cbe0172c8ae37306324d44d5e6b0c69ac57b393f1ab370fd725cc647692444a04ef87387aa68d53743493b9eba14cc552ca2a93a";

    char msg[] = "test";
    mclBnFr x, d;
    mclBnG2 Q, xQ, Rho;

    mclBnFr_setByCSPRNG(&x);
    mclBnG2_setStr(&Q, G2_basePoint_hexstr, strlen(G2_basePoint_hexstr), 16);
    mclBnG2_mul(&xQ, &Q, &x);

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    printf("sok_G2 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        sok_G2(&Q, &x, &xQ, (unsigned char *)msg, strlen(msg), &Rho, &d);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

    printf("sokver_G2 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        sokver_G2(&Q, &xQ, &Rho, &d, (unsigned char *)msg, strlen(msg));
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

}

void mesure_nipok(int nb_trials)
{

    const char * G1_basePoint_hexstr = "1 21a6d67ef250191fadba34a0a30160b9ac9264b6f95f63b3edbec3cf4b2e689db1bbb4e69a416a0b1e79239c0372e5cd70113c98d91f36b6980d 0118ea0460f7f7abb82b33676a7432a490eeda842cccfa7d788c659650426e6af77df11b8ae40eb80f475432c66600622ecaa8a5734d36fb03de";

    mclBnFr x, d;
    mclBnG1 P, xP, Rho;

    mclBnFr_setByCSPRNG(&x);
    mclBnG1_setStr(&P, G1_basePoint_hexstr, strlen(G1_basePoint_hexstr), 16);

    mclBnG1_mul(&xP, &P, &x);

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    printf("nipok_G1 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        nipok_G1(&P, &x, &xP, &Rho, &d);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

    printf("nipokver_G1 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        nipokver_G1(&P, &xP, &Rho, &d);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

}

void mesure_eqnipok(int nb_trials)
{
    const char * G1_basePoint_hexstr = "1 21a6d67ef250191fadba34a0a30160b9ac9264b6f95f63b3edbec3cf4b2e689db1bbb4e69a416a0b1e79239c0372e5cd70113c98d91f36b6980d 0118ea0460f7f7abb82b33676a7432a490eeda842cccfa7d788c659650426e6af77df11b8ae40eb80f475432c66600622ecaa8a5734d36fb03de";
    const char * G2_basePoint_hexstr = "1 0257ccc85b58dda0dfb38e3a8cbdc5482e0337e7c1cd96ed61c913820408208f9ad2699bad92e0032ae1f0aa6a8b48807695468e3d934ae1e4df 1d2e4343e8599102af8edca849566ba3c98e2a354730cbed9176884058b18134dd86bae555b783718f50af8b59bf7e850e9b73108ba6aa8cd283 0a0650439da22c1979517427a20809eca035634706e23c3fa7a6bb42fe810f1399a1f41c9ddae32e03695a140e7b11d7c3376e5b68df0db7154e 073ef0cbd438cbe0172c8ae37306324d44d5e6b0c69ac57b393f1ab370fd725cc647692444a04ef87387aa68d53743493b9eba14cc552ca2a93a";

    mclBnFr x, d;
    mclBnG1 P, xP, Rho;
    mclBnG2 Q, xQ, Sigma;

    mclBnFr_setByCSPRNG(&x);
    mclBnG1_setStr(&P, G1_basePoint_hexstr, strlen(G1_basePoint_hexstr), 16);
    mclBnG2_setStr(&Q, G2_basePoint_hexstr, strlen(G2_basePoint_hexstr), 16);

    mclBnG1_mul(&xP, &P, &x);

    mclBnG2_mul(&xQ, &Q, &x);

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    printf("eq_nipok_G1_G2 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        eq_nipok_G1_G2(&P, &xP, &Q, &xQ, &x, &Rho, &Sigma, &d);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

    printf("eq_nipokver_G1_G2 : \n");
    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        eq_nipokver_G1_G2(&P, &xP, &Q, &xQ, &Rho, &Sigma, &d);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    print_trials_res(trials_res, nb_trials);

}

double mesure_ake_a_get_mx(mclBnG1 * P, mclBnG2 * Q, mclBnFr * x, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_a_get_mx(P, Q, x, omega, omega_len, xP, xQ, x_ni);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_mx(mclBnG1 * P, mclBnG1 * xP, mclBnG2 * Q, mclBnG2 * xQ, unsigned char * omega, size_t omega_len, XY_ni * x_ni, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_mx(P, xP, Q, xQ, omega, omega_len, x_ni);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_ake_b_get_my(mclBnG2 * Q, mclBnFr * y, unsigned char * omega, size_t omega_len, mclBnG2 * yQ, XY_ni * y_ni, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_b_get_my(Q, y, omega, omega_len, yQ, y_ni);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_ake_b_get_sigma_Y_1(char * priv_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1, size_t sig_len, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_b_get_sigma_Y_1(priv_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_my(mclBnG2 * Q, mclBnG2 * yQ, unsigned char * omega, size_t omega_len, XY_ni * y_ni, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_my(Q, yQ, omega, omega_len, y_ni);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_sigma_Y_1(char * pub_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1, size_t sig_len, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_sigma_Y_1(pub_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_ake_a_get_sigma_X(char * priv_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1 , unsigned char * sigma_X, size_t sig_len, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_a_get_sigma_X(priv_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sigma_X, sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_sigma_X(char * pub_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1 , unsigned char * sigma_X, size_t sig_len, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_sigma_X(pub_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sigma_X, sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_ake_b_get_sigma_Y_2(char * priv_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1 , unsigned char * sigma_X, unsigned char * sigma_Y_2, size_t sig_len, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_b_get_sigma_Y_2(priv_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_sigma_Y_2(char * pub_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1 , unsigned char * sigma_X, unsigned char * sigma_Y_2, size_t sig_len, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_sigma_Y_2(pub_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_ake_a_get_shared_key(mclBnG1 * L_pk, mclBnG2 * yQ, mclBnFr * x, mclBnGT * ka, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_a_get_shared_key(L_pk, yQ, x, ka);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_ake_b_get_shared_key(mclBnG1 * L_pk, mclBnG2 * xQ, mclBnFr * y, mclBnGT * kb, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_b_get_shared_key(L_pk, xQ, y, kb);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_ake_O_get_sst(char * priv_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1 , unsigned char * sigma_X, unsigned char * sigma_Y_2, size_t sig_len, SST * sst, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        ake_O_get_sst(priv_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, sst);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
        free(sst->m);
        free(sst->sigma_O);
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_sst(char * pub_key_path, unsigned char * omega, size_t omega_len, mclBnG1 * xP, mclBnG2 * xQ, XY_ni * x_ni, mclBnG2 * yQ, XY_ni * y_ni, unsigned char * sigma_Y_1 , unsigned char * sigma_X, unsigned char * sigma_Y_2, size_t sig_len, SST * sst, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_sst(pub_key_path, omega, omega_len, xP, xQ, x_ni, yQ, y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, sst);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_L_ni_for_2L(mclBnG1 * P, mclBnG1 * L1_pk, Lambda_ni * l1_ni, mclBnG1 * L2_pk, Lambda_ni * l2_ni, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_L_ni(P, 4, L1_pk, l1_ni, L2_pk, l2_ni);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_tdgen_get_li_T1_for_2L(mclBnG1 * xP, mclBnG2 * yQ, mclBnFr * l1_sk, mclBnGT * l1_T1, mclBnFr * l2_sk, mclBnGT * l2_T1, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        tdgen_get_li_T1(xP, yQ, 4, l1_sk, l1_T1, l2_sk, l2_T1);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_tdgen_get_li_T2_for_2L(mclBnG1 * P, mclBnG1 * xP, mclBnG2 * yQ, mclBnG1 * L1_pk, mclBnFr * l1_sk, mclBnGT * l1_T1, Lambda_eq_ni * l1_T2, mclBnG1 * L2_pk, mclBnFr * l2_sk, mclBnGT * l2_T1, Lambda_eq_ni * l2_T2, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        tdgen_get_li_T2(P, xP, yQ, 8, L1_pk, l1_sk, l1_T1, l1_T2, L2_pk, l2_sk, l2_T1, l2_T2);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_verify_li_T2_for_2L(mclBnG1 * P, mclBnG1 * xP, mclBnG2 * yQ, mclBnG1 * L1_pk, mclBnGT * l1_T1, Lambda_eq_ni * l1_T2, mclBnG1 * L2_pk, mclBnGT * l2_T1, Lambda_eq_ni * l2_T2, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        verify_li_T2(P, xP, yQ, 6, L1_pk, l1_T1, l1_T2, L2_pk, l2_T1, l2_T2);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

double mesure_open_get_shared_key_for_2(mclBnGT * k, mclBnGT * l1_T1, mclBnGT * l2_T1, int nb_trials)
{

    double trials_res[nb_trials];
    clock_t begin;
    clock_t end;

    for(int i = 0; i < nb_trials; i++)
    {
        begin = clock();
        open_get_shared_key(k, 2, l1_T1, l2_T1);
        end = clock();
        trials_res[i] = (double)(end - begin) / CLOCKS_PER_SEC;
    }

    return print_trials_res(trials_res, nb_trials);

}

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

    printf("Run multiple times differents functions of the like potocol. All mesures are in ms.\n\n");

    // Mesure perf parameter
    int nb_trials = 100;
    if( argc == 2 )
        nb_trials = atoi(argv[1]);
    printf("Number of trials : %d\n\n", nb_trials);
    double total_time_A = 0.0;
    double total_time_B = 0.0;
    double total_time_O = 0.0;
    double total_time_Ver = 0.0;
    double total_time_TDGen = 0.0;
    double total_time_Open = 0.0;

    /*******************************************************************
     *                                                                 *
     *                      INDEPENDANT MESURES                        *
     *                                                                 *
     * *****************************************************************/

    mesure_pairing(nb_trials);
    mesure_arithmetic_EC(nb_trials);
    mesure_edd25519(nb_trials);
    mesure_sok(nb_trials);
    mesure_nipok(nb_trials);
    mesure_eqnipok(nb_trials);

    /*******************************************************************
     *                                                                 *
     *                          LIKE MESURES                           *
     *                                                                 *
     * *****************************************************************/

    // General Param
    mclBnG1 P, xP;
    mclBnG2 Q, xQ, yQ;
    mclBnFr x, y;
    XY_ni x_ni, y_ni;
    mclBnGT ka, kb, k, l1_T1, l2_T1;
    SST sst;
    Lambda_eq_ni l1_T2, l2_T2;

    // Id
    size_t id_len = sizeof(uuid_t);
    uuid_t id_A, id_B, id_L1, id_L2;
    uuid_generate_random(id_A);
    uuid_generate_random(id_B);
    uuid_generate_random(id_L1);
    uuid_generate_random(id_L2);
    size_t omega_len = 4 * id_len;
    unsigned char omega[omega_len];

    // Keys
    char pub_key_A[] = "keys/a_pub_key.pem";
    char priv_key_A[] = "keys/a_priv_key.pem";
    char pub_key_B[] = "keys/b_pub_key.pem";
    char priv_key_B[] = "keys/b_priv_key.pem";
    char pub_key_O[] = "keys/o_pub_key.pem";
    char priv_key_O[] = "keys/o_priv_key.pem";
    mclBnG1 L1_pk, L2_pk, L_pk;
    mclBnFr l1_sk, l2_sk;
    Lambda_ni l1_ni, l2_ni;

    // Sig buffer
    size_t sig_len = ED25519_SIG_LENGTH;
    unsigned char * sigma_Y_1 = (unsigned char *) malloc(sig_len * sizeof(unsigned char));
    unsigned char * sigma_X = (unsigned char *) malloc(sig_len * sizeof(unsigned char));
    unsigned char * sigma_Y_2 = (unsigned char *) malloc(sig_len * sizeof(unsigned char));

    // Setup
    setup(&P, &Q);

    //UKeyGen A
    u_o_key_gen(pub_key_A, priv_key_A);

    // UKeyGen B
    u_o_key_gen(pub_key_B, priv_key_B);

    // UKeyGen O
    u_o_key_gen(pub_key_O, priv_key_O);

    // LambdaKeyGen L1
    a_key_gen(&P, &l1_sk, &L1_pk, &l1_ni);

    // LambdaKeyGen L2
    a_key_gen(&P, &l2_sk, &L2_pk, &l2_ni);

    // AKE Precal
    verify_L_ni(&P, 4, &L1_pk, &l1_ni, &L2_pk, &l2_ni);
    ake_precalc_add_lipk(&L_pk, 2, &L1_pk, &L2_pk);
    ake_precalc_get_omega(omega, 8, id_A, id_len, id_B, id_len, id_L1, id_len, id_L2, id_len);

    // Mesure ake_a_get_mx
    printf("ake_a_get_mx : \n");
    double time_ake_a_get_mx = mesure_ake_a_get_mx(&P, &Q, &x, omega, omega_len, &xP, &xQ, &x_ni, nb_trials);

    // A : ake_a_get_mx
    ake_a_get_mx(&P, &Q, &x, omega, omega_len, &xP, &xQ, &x_ni);
    total_time_A += time_ake_a_get_mx;

    // Mesure verify_mx
    printf("verify_mx : \n");
    double time_verify_mx = mesure_verify_mx(&P, &xP, &Q, &xQ, omega, omega_len, &x_ni, nb_trials);

    // O : verify_mx
    verify_mx(&P, &xP, &Q, &xQ, omega, omega_len, &x_ni);
    total_time_O += time_verify_mx;

    // Mesure ake_b_get_my
    printf("ake_b_get_my : \n");
    double time_ake_b_get_my = mesure_ake_b_get_my(&Q, &y, omega, omega_len, &yQ, &y_ni, nb_trials);

    // Mesure ake_b_get_sigma_Y_1
    printf("ake_b_get_sigma_Y_1 : \n");
    double time_ake_b_get_sigma_Y_1 = mesure_ake_b_get_sigma_Y_1(priv_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sig_len, nb_trials);

    // B : ake_b_get_my, verify_mx, ake_b_get_sigma_Y_1
    verify_mx(&P, &xP, &Q, &xQ, omega, omega_len, &x_ni);
    ake_b_get_my(&Q, &y, omega, omega_len, &yQ, &y_ni);
    ake_b_get_sigma_Y_1(priv_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sig_len);
    total_time_B += time_verify_mx + time_ake_b_get_my + time_ake_b_get_sigma_Y_1;

    // Mesure verify_my
    printf("verify_my : \n");
    double time_verify_my = mesure_verify_my(&Q, &yQ, omega, omega_len, &y_ni, nb_trials);

    // Mesure verify_sigma_Y_1
    printf("verify_sigma_Y_1 : \n");
    double time_verify_sigma_Y_1 = mesure_verify_sigma_Y_1(pub_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sig_len, nb_trials);

    // O : verify_my, verify_sigma_Y_1
    verify_my(&Q, &yQ, omega, omega_len, &y_ni);
    verify_sigma_Y_1(pub_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sig_len);
    total_time_O += time_verify_my + time_verify_sigma_Y_1;

    // Mesure ake_a_get_sigma_X
    printf("ake_a_get_sigma_X : \n");
    double time_ake_a_get_sigma_X = mesure_ake_a_get_sigma_X(priv_key_A, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sig_len, nb_trials);

    // A : verify_my, verify_sigma_Y_1, ake_a_get_sigma_X
    verify_my(&Q, &yQ, omega, omega_len, &y_ni);
    verify_sigma_Y_1(pub_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sig_len);
    ake_a_get_sigma_X(priv_key_A, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sig_len);
    total_time_A += time_verify_my + time_verify_sigma_Y_1 + time_ake_a_get_sigma_X;

    // Mesure verify_sigma_X
    printf("verify_sigma_X : \n");
    double time_verify_sigma_X = mesure_verify_sigma_X(pub_key_A, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sig_len, nb_trials);

    // O : verify_sigma_X
    verify_sigma_X(pub_key_A, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sig_len);
    total_time_O += time_verify_sigma_X;

    // Mesure ake_b_get_sigma_Y_2
    printf("ake_b_get_sigma_Y_2 : \n");
    double time_ake_b_get_sigma_Y_2 = mesure_ake_b_get_sigma_Y_2(priv_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, nb_trials);

    // B : verify_sigma_X, ake_b_get_sigma_Y_2
    verify_sigma_X(pub_key_A, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sig_len);
    ake_b_get_sigma_Y_2(priv_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len);
    total_time_B += time_verify_sigma_X + time_ake_b_get_sigma_Y_2;

    // Mesure ake_a_get_shared_key
    printf("ake_a_get_shared_key : \n");
    double time_ake_a_get_shared_key = mesure_ake_a_get_shared_key(&L_pk, &yQ, &x, &ka, nb_trials);

    // A : ake_a_get_shared_key
    ake_a_get_shared_key(&L_pk, &yQ, &x, &ka);
    total_time_A += time_ake_a_get_shared_key;

    // Mesure ake_b_get_shared_key
    printf("ake_b_get_shared_key : \n");
    double time_ake_b_get_shared_key = mesure_ake_b_get_shared_key(&L_pk, &xQ, &y, &kb, nb_trials);

    // B : ake_b_get_shared_key
    ake_b_get_shared_key(&L_pk, &xQ, &y, &kb);
    total_time_B += time_ake_b_get_shared_key;

    // Mesure ake_b_get_sigma_Y_2
    printf("verify_sigma_Y_2 : \n");
    double time_verify_sigma_Y_2 = mesure_verify_sigma_Y_2(pub_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, nb_trials);

    // Mesure ake_O_get_sst
    printf("ake_O_get_sst : \n");
    double time_ake_O_get_sst = mesure_ake_O_get_sst(priv_key_O, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, &sst, nb_trials);

    // O : verify_sigma_Y_2, ake_O_get_sst
    verify_sigma_Y_2(pub_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len);
    ake_O_get_sst(priv_key_O, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, &sst);
    total_time_O += time_verify_sigma_Y_2 + time_ake_O_get_sst;

    // Mesure verify_sst
    printf("verify_sst : \n");
    double time_verify_sst = mesure_verify_sst(pub_key_O, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, &sst, nb_trials);

    // Mesure verify_L_ni
    printf("verify_L_ni (2 Authorities) : \n");
    double time_verify_L_ni = mesure_verify_L_ni_for_2L(&P, &L1_pk, &l1_ni, &L2_pk, &l2_ni, nb_trials);

    // Verification
    verify_L_ni(&P, 4, &L1_pk, &l1_ni, &L2_pk, &l2_ni);
    verify_mx(&P, &xP, &Q, &xQ, omega, omega_len, &x_ni);
    verify_my(&Q, &yQ, omega, omega_len, &y_ni);
    verify_sigma_Y_1(pub_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sig_len);
    verify_sigma_X(pub_key_A, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sig_len);
    verify_sigma_Y_2(pub_key_B, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len);
    verify_sst(pub_key_O, omega, omega_len, &xP, &xQ, &x_ni, &yQ, &y_ni, sigma_Y_1, sigma_X, sigma_Y_2, sig_len, &sst);
    total_time_Ver += time_verify_L_ni + time_verify_mx + time_verify_my + time_verify_sigma_Y_1 + time_verify_sigma_X + time_verify_sigma_Y_2 + time_verify_sst;

    // Mesure tdgen_get_li_T1
    printf("tdgen_get_li_T1 (2 Authorities) : \n");
    double time_tdgen_get_li_T1_for_2L = mesure_tdgen_get_li_T1_for_2L(&xP, &yQ, &l1_sk, &l1_T1, &l2_sk, &l2_T1, nb_trials);

    // Mesure tdgen_get_li_T2
    printf("tdgen_get_li_T2 (2 Authorities) : \n");
    double time_tdgen_get_li_T2_for_2L = mesure_tdgen_get_li_T2_for_2L(&P, &xP, &yQ, &L1_pk, &l1_sk, &l1_T1, &l1_T2, &L2_pk, &l2_sk, &l2_T1, &l2_T2, nb_trials);

    // TDGen
    tdgen_get_li_T1(&xP, &yQ, 4, &l1_sk, &l1_T1, &l2_sk, &l2_T1);
    tdgen_get_li_T2(&P, &xP, &yQ, 8, &L1_pk, &l1_sk, &l1_T1, &l1_T2, &L2_pk, &l2_sk, &l2_T1, &l2_T2);
    total_time_TDGen += time_tdgen_get_li_T1_for_2L + time_tdgen_get_li_T2_for_2L;

    // Mesure verify_li_T2
    printf("verify_li_T2 (2 Authorities) : \n");
    double time_verify_li_T2_for_2 = mesure_verify_li_T2_for_2L(&P, &xP, &yQ, &L1_pk, &l1_T1, &l1_T2, &L2_pk, &l2_T1, &l2_T2, nb_trials);

    // Mesure open_get_shared_key
    printf("open_get_shared_key (2 Authorities) : \n");
    double time_open_get_shared_key_for_2 = mesure_open_get_shared_key_for_2(&k, &l1_T1, &l2_T1, nb_trials);

    // Open
    verify_li_T2(&P, &xP, &yQ, 6, &L1_pk, &l1_T1, &l1_T2, &L2_pk, &l2_T1, &l2_T2);
    open_get_shared_key(&k, 2, &l1_T1, &l2_T1);
    total_time_Open += time_verify_li_T2_for_2 + time_open_get_shared_key_for_2;

    // Check key
    int res = mclBnGT_isEqual(&ka, &kb);
    if(res != 1)
    {
        printf("Keys not equal\n");
    }
    res = mclBnGT_isEqual(&ka, &k);
    if(res != 1)
    {
        printf("Keys not equal\n");
    }

    printf("Total CPU running time for A : %f\n", total_time_A);
    printf("Total CPU running time for B : %f\n", total_time_B);
    printf("Total CPU running time for O : %f\n", total_time_O);
    printf("Total CPU running time for Verification : %f\n", total_time_Ver);
    printf("Total CPU running time for TDGen (2 Authorities) : %f\n", total_time_TDGen);
    printf("Total CPU running time for Open (2 Authorities) : %f\n", total_time_Open);

    free(sigma_Y_1);
    free(sigma_X);
    free(sigma_Y_2);
    free(sst.m);
    free(sst.sigma_O);

    return 0;

}

max_stars_repo_name: AnonymLIKE/like

max_issues_repo_name: AnonymLIKE/like

max_forks_repo_name: AnonymLIKE/like


text:
#ifndef VKST_PLAT_FILE_HANDLE_H
#define VKST_PLAT_FILE_HANDLE_H

#include <plat/filesystem.h>
#include <turf/c/core.h>
#include <gsl.h>

namespace plat {

// A RAII wrapper around the C stdio file handle
class file_handle {
public:
  enum class open_modes : uint8_t {
    read = 0x1,
    write = 0x2,
    append = 0x4,
  }; // enum class open_modes

  static file_handle open(plat::filesystem::path const& path, open_modes mode,
                          std::error_code& ec) noexcept;

  file_handle() noexcept = default;

  std::FILE* get() noexcept { return _handle.get(); }
  operator std::FILE*() noexcept { return get(); }

  explicit operator bool() noexcept { return static_cast<bool>(_handle); }

  void reset() noexcept { _handle.reset(); }

private:
  file_handle(FILE* handle) noexcept : _handle{handle, &std::fclose} {}

  using file_ptr = gsl::unique_ptr<std::FILE, decltype(&std::fclose)>;
  file_ptr _handle{nullptr, &std::fclose};
}; // class file_handle

inline constexpr auto operator|(file_handle::open_modes a,
                                file_handle::open_modes b) noexcept {
  using U = std::underlying_type_t<file_handle::open_modes>;
  return static_cast<file_handle::open_modes>(static_cast<U>(a) |
                                              static_cast<U>(b));
}

inline constexpr auto operator&(file_handle::open_modes a,
                                file_handle::open_modes b) noexcept {
  using U = std::underlying_type_t<file_handle::open_modes>;
  return static_cast<file_handle::open_modes>(static_cast<U>(a) &
                                              static_cast<U>(b));
}

} // namespace plat

#endif // VKST_PLAT_FILE_HANDLE_H

ext: h

max_stars_repo_name: wesleygriffin/vkst

max_issues_repo_name: wesleygriffin/vkst

max_forks_repo_name: wesleygriffin/vkst


text:
// SPDX-License-Identifier: MIT
// The MIT License (MIT)
//
// Copyright (c) 2014-2018, Institute for Software & Systems Engineering
// Copyright (c) 2018-2019, Johannes Leupolz
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.

#ifndef PEMC_BASIC_RAW_MEMORY_H_
#define PEMC_BASIC_RAW_MEMORY_H_

#include <gsl/gsl_byte>
#include <memory>

namespace pemc {

using deleter_t = std::function<void(void*)>;

using unique_void_ptr = std::unique_ptr<void, deleter_t>;

// The following method can be used to reserve some bulk memory.
// The unique_void_ptr keeps track of the memory to avoid leaks.
// unique_ptr are not used, because they do not support void.
// Example of usage:
// > unique_void_ptr stateMemory = unique_void(::operator new(amount of memory )
// ); > gsl::byte* p_stateMemory = static_cast<gsl::byte*>(stateMemory.get());
unique_void_ptr unique_void(void* ptr);

/// <summary>
///   Compares the two buffers <paramref name="buffer1" /> and <paramref
///   name="buffer2" />, returning <c>true</c> when the buffers are equivalent.
/// </summary>
/// <param name="buffer1">The first buffer of memory to compare.</param>
/// <param name="buffer2">The second buffer of memory to compare.</param>
/// <param name="sizeInBytes">The size of the buffers in bytes.</param>
bool areBuffersEqual(gsl::byte* buffer1,
                     gsl::byte* buffer2,
                     size_t sizeInBytes);

/// <summary>
///   Copies <paramref name="sizeInBytes" />-many bytes from <paramref
///   name="source" /> to <see cref="destination" />.
/// </summary>
/// <param name="source">The first buffer of memory to compare.</param>
/// <param name="destination">The second buffer of memory to compare.</param>
/// <param name="sizeInBytes">The size of the buffers in bytes.</param>
void copyBuffers(gsl::byte* source, gsl::byte* destination, size_t sizeInBytes);

/// <summary>
///   Hashes the <paramref name="buffer" />.
/// </summary>
/// <param name="buffer">The buffer of memory that should be hashed.</param>
/// <param name="sizeInBytes">The size of the buffer in bytes.</param>
/// <param name="seed">The seed value for the hash.</param>
/// <remarks>See also https://en.wikipedia.org/wiki/MurmurHash (MurmurHash3
/// implementation)</remarks>
uint32_t hashBuffer(gsl::byte* buffer, size_t sizeInBytes, int32_t seed);

}  // namespace pemc

#endif  // PEMC_BASIC_RAW_MEMORY_H_

max_stars_repo_name: joleuger/pemc

max_issues_repo_name: joleuger/pemc

max_forks_repo_name: joleuger/pemc


text:
/* ################################################################################## */
/* ###                                                                            ### */
/* ###                                 sigvel.c                                   ### */
/* ###                                                                            ### */
/* ###   Original: hydra.c (public version of Gadget 2)                           ### */
/* ###   Author: Volker Springel                                                  ### */
/* ###                                                                            ### */
/* ###   Modified: February 2011                                                  ### */
/* ###   Author: Fabrice Durier                                                   ### */
/* ###                                                                            ### */
/* ################################################################################## */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef NOMPI
#include <mpi.h>
#endif
#include <gsl/gsl_math.h>
#include "allvars.h"
#include "proto.h"

/*! \file sigvel.c
 *  Re-Computation of SPH acceleration and maximum signal velocity
 *  for 'feedback' particles only.
 *
 *  This file is a modified version of hydra.c, where the SPH forces are
 *  computed, and where the rate of change of entropy due to the shock heating
 *  (via artificial viscosity) is computed.
 */

#ifdef TIMESTEP_UPDATE

static double hubble_a, atime, hubble_a2, fac_mu, fac_vsic_fix, a3inv, fac_egy;

#ifdef PERIODIC
static double boxSize, boxHalf;

#ifdef LONG_X
static double boxSize_X, boxHalf_X;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
static double boxSize_Y, boxHalf_Y;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
static double boxSize_Z, boxHalf_Z;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif

/*! This function is the driver routine for the update of the calculation
 *  of hydrodynamical force due to the feedback impact
 *  from either SNe or BHs on active gas particles.
 */
void get_sigvel(void)
{
  long long ntot, ntotleft;
  int i, j, k, n, ngrp, maxfill, source, ndone;
  int *nbuffer, *noffset, *nsend_local, *nsend, *numlist, *ndonelist;
  int level, sendTask, recvTask, nexport, place;
  double soundspeed_i;
  double tstart, tend, sumt, sumcomm;
  double timecomp = 0, timecommsumm = 0, timeimbalance = 0, sumimbalance;
#ifndef NOMPI
  MPI_Status status;
#endif

#ifdef PERIODIC
  boxSize = All.BoxSize;
  boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
  boxHalf_X = boxHalf * LONG_X;
  boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
  boxHalf_Y = boxHalf * LONG_Y;
  boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
  boxHalf_Z = boxHalf * LONG_Z;
  boxSize_Z = boxSize * LONG_Z;
#endif
#endif

  if(All.ComovingIntegrationOn)
    {
      /* Factors for comoving integration of hydro */
      hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
        + (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;

      hubble_a = All.Hubble * sqrt(hubble_a);
      hubble_a2 = All.Time * All.Time * hubble_a;

      fac_mu = pow(All.Time, 3 * (GAMMA - 1) / 2) / All.Time;

      fac_egy = pow(All.Time, 3 * (GAMMA - 1));

      fac_vsic_fix = hubble_a * pow(All.Time, 3 * GAMMA_MINUS1);

      a3inv = 1 / (All.Time * All.Time * All.Time);
      atime = All.Time;
    }
  else
    hubble_a = hubble_a2 = atime = fac_mu = fac_vsic_fix = a3inv = fac_egy = 1.0;

  /* `NumSphUpdate' gives the number of feedback particles on this processor that want a force update */
  for(n = 0, NumSphUpdate = 0; n < N_gas; n++)
    {
      if(P[n].Ti_endstep == All.Ti_Current && SphP[n].FeedbackFlag > 1)
	NumSphUpdate++;
    }

  numlist = malloc(NTask * sizeof(int) * NTask);

#ifndef NOMPI
  MPI_Allgather(&NumSphUpdate, 1, MPI_INT, numlist, 1, MPI_INT, GADGET_WORLD);
#else
    numlist[0] = NumSphUpdate;
#endif
  for(i = 0, ntot = 0; i < NTask; i++)
    ntot += numlist[i];
  free(numlist);

  if(ThisTask == 0 && ntot)
    printf("\t --->  Number of feedback particles = %ld \n\n", ntot);

  noffset = malloc(sizeof(int) * NTask);	/* offsets of bunches in common list */
  nbuffer = malloc(sizeof(int) * NTask);
  nsend_local = malloc(sizeof(int) * NTask);
  nsend = malloc(sizeof(int) * NTask * NTask);
  ndonelist = malloc(sizeof(int) * NTask);

  i = 0;			/* first particle for this task */
  ntotleft = ntot;		/* particles left for all tasks together */

  while(ntotleft > 0)
    {
      for(j = 0; j < NTask; j++)
	nsend_local[j] = 0;

      /* do local particles and prepare export list */
      tstart = second();
      for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeHydro - NTask; i++)
	if(P[i].Ti_endstep == All.Ti_Current && SphP[i].FeedbackFlag > 1)
	  {
	    ndone++;

	    for(j = 0; j < NTask; j++)
	      Exportflag[j] = 0;

	    get_sigvel_evaluate(i, 0);

	    for(j = 0; j < NTask; j++)
	      {
		if(Exportflag[j])
		  {
		    for(k = 0; k < 3; k++)
		      {
			HydroDataIn[nexport].Pos[k] = P[i].Pos[k];
			HydroDataIn[nexport].Vel[k] = SphP[i].VelPred[k];
		      }
                    HydroDataIn[nexport].Hsml = SphP[i].Hsml;
                    HydroDataIn[nexport].Mass = P[i].Mass;
                    HydroDataIn[nexport].DhsmlDensityFactor = SphP[i].DhsmlDensityFactor;
                    HydroDataIn[nexport].Density = SphP[i].Density;
                    HydroDataIn[nexport].Pressure = SphP[i].Pressure;
                    HydroDataIn[nexport].Timestep = P[i].Ti_endstep - P[i].Ti_begstep;

                    /* calculation of F1 */
                    soundspeed_i = sqrt(GAMMA * SphP[i].Pressure / SphP[i].Density);
                    HydroDataIn[nexport].F1 = fabs(SphP[i].DivVel) /
		      (fabs(SphP[i].DivVel) + SphP[i].CurlVel +
                       0.0001 * soundspeed_i / SphP[i].Hsml / fac_mu);

		    HydroDataIn[nexport].Index = i;
		    HydroDataIn[nexport].Task = j;
		    nexport++;
		    nsend_local[j]++;
		  }
	      }
	  }
      tend = second();
      timecomp += timediff(tstart, tend);

      qsort(HydroDataIn, nexport, sizeof(struct hydrodata_in), hydro_compare_key);

      for(j = 1, noffset[0] = 0; j < NTask; j++)
	noffset[j] = noffset[j - 1] + nsend_local[j - 1];

      tstart = second();
#ifndef NOMPI
      MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, GADGET_WORLD);
#else
    nsend[0] = nsend_local[0];
#endif
      tend = second();
      timeimbalance += timediff(tstart, tend);

      /* now do the particles that need to be exported */

      for(level = 1; level < (1 << PTask); level++)
	{
	  tstart = second();
	  for(j = 0; j < NTask; j++)
	    nbuffer[j] = 0;
	  for(ngrp = level; ngrp < (1 << PTask); ngrp++)
	    {
	      maxfill = 0;
	      for(j = 0; j < NTask; j++)
		{
		  if((j ^ ngrp) < NTask)
		    if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
		      maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
		}
	      if(maxfill >= All.BunchSizeHydro)
		break;

	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;
#ifndef NOMPI
	      if(recvTask < NTask)
		{
		  if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
		    {
		      /* get the particles */
		      MPI_Sendrecv(&HydroDataIn[noffset[recvTask]],
				   nsend_local[recvTask] * sizeof(struct hydrodata_in), MPI_BYTE,
				   recvTask, TAG_HYDRO_A,
				   &HydroDataGet[nbuffer[ThisTask]],
				   nsend[recvTask * NTask + ThisTask] * sizeof(struct hydrodata_in), MPI_BYTE,
				   recvTask, TAG_HYDRO_A, GADGET_WORLD, &status);
		    }
		}
#endif
	      for(j = 0; j < NTask; j++)
		if((j ^ ngrp) < NTask)
		  nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
	    }
	  tend = second();
	  timecommsumm += timediff(tstart, tend);

	  /* now do the imported particles */
	  tstart = second();
	  for(j = 0; j < nbuffer[ThisTask]; j++)
	    get_sigvel_evaluate(j, 1);
	  tend = second();
	  timecomp += timediff(tstart, tend);

	  /* do a block to measure imbalance */
	  tstart = second();
#ifndef NOMPI
	      MPI_Barrier(GADGET_WORLD);
#endif
	  tend = second();
	  timeimbalance += timediff(tstart, tend);

	  /* get the result */
	  tstart = second();
	  for(j = 0; j < NTask; j++)
	    nbuffer[j] = 0;
	  for(ngrp = level; ngrp < (1 << PTask); ngrp++)
	    {
	      maxfill = 0;
	      for(j = 0; j < NTask; j++)
		{
		  if((j ^ ngrp) < NTask)
		    if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
		      maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
		}
	      if(maxfill >= All.BunchSizeHydro)
		break;

	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;

#ifndef NOMPI
	      if(recvTask < NTask)
		{
		  if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
		    {
		      /* send the results */
		      MPI_Sendrecv(&HydroDataResult[nbuffer[ThisTask]],
				   nsend[recvTask * NTask + ThisTask] * sizeof(struct hydrodata_out),
				   MPI_BYTE, recvTask, TAG_HYDRO_B,
				   &HydroDataPartialResult[noffset[recvTask]],
				   nsend_local[recvTask] * sizeof(struct hydrodata_out),
				   MPI_BYTE, recvTask, TAG_HYDRO_B, GADGET_WORLD, &status);

		      /* add the result to the particles */
		      for(j = 0; j < nsend_local[recvTask]; j++)
			{
			  source = j + noffset[recvTask];
			  place = HydroDataIn[source].Index;

			  for(k = 0; k < 3; k++)
			    SphP[place].FeedAccel[k] += HydroDataPartialResult[source].Acc[k];

			  if(SphP[place].MaxSignalVel < HydroDataPartialResult[source].MaxSignalVel)
			    SphP[place].MaxSignalVel = HydroDataPartialResult[source].MaxSignalVel;
			}
		    }
		}
#endif
	      for(j = 0; j < NTask; j++)
		if((j ^ ngrp) < NTask)
		  nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
	    }
	  tend = second();
	  timecommsumm += timediff(tstart, tend);

	  level = ngrp - 1;
	}

#ifndef NOMPI
      MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, GADGET_WORLD);
#else
    ndonelist[0] = ndone;
#endif
      for(j = 0; j < NTask; j++)
	ntotleft -= ndonelist[j];
    }

  free(ndonelist);
  free(nsend);
  free(nsend_local);
  free(nbuffer);
  free(noffset);

  /* collect some timing information */
#ifndef NOMPI
  MPI_Reduce(&timecomp, &sumt, 1, MPI_DOUBLE, MPI_SUM, 0, GADGET_WORLD);
  MPI_Reduce(&timecommsumm, &sumcomm, 1, MPI_DOUBLE, MPI_SUM, 0, GADGET_WORLD);
  MPI_Reduce(&timeimbalance, &sumimbalance, 1, MPI_DOUBLE, MPI_SUM, 0, GADGET_WORLD);
#else
    sumt = timecomp;
    sumcomm = timecommsumm;
    sumimbalance = timeimbalance;
#endif
  if(ThisTask == 0)
    {
      All.CPU_HydCompWalk += sumt / NTask;
      All.CPU_HydCommSumm += sumcomm / NTask;
      All.CPU_HydImbalance += sumimbalance / NTask;
    }
}

/*! This function is the 'core' of the SPH force update. A target
 *  particle is specified which may either be local, or reside in the
 *  communication buffer.
 *  Note that only the acceleration due to feedback and the change
 *  of signal velocity are updated, not the actual variation of Entropy rate.
 */
void get_sigvel_evaluate(int target, int mode)
{
  int j, k, n, timestep, startnode, numngb;
  FLOAT *pos, *vel;
  FLOAT mass, h_i, dhsmlDensityFactor, rho, pressure, f1, f2;
  double acc[3], dtEntropy, maxSignalVel;
  double dx, dy, dz, dvx, dvy, dvz;
  double h_i2, hinv, hinv4;
  double p_over_rho2_i, p_over_rho2_j, soundspeed_i, soundspeed_j;
  double hfc, dwk_i, vdotr, vdotr2, visc, mu_ij, rho_ij, vsig;
  double h_j, dwk_j, r, r2, u, hfc_visc;

#ifndef NOVISCOSITYLIMITER
  double dt;
#endif

  if(mode == 0)
    {
      pos = P[target].Pos;
      vel = SphP[target].VelPred;
      h_i = SphP[target].Hsml;
      mass = P[target].Mass;
      dhsmlDensityFactor = SphP[target].DhsmlDensityFactor;
      rho = SphP[target].Density;
      pressure = SphP[target].Pressure;
      timestep = P[target].Ti_endstep - P[target].Ti_begstep;
      soundspeed_i = sqrt(GAMMA * pressure / rho);
      f1 = fabs(SphP[target].DivVel) /
        (fabs(SphP[target].DivVel) + SphP[target].CurlVel +
         0.0001 * soundspeed_i / SphP[target].Hsml / fac_mu);
    }
  else
    {
      pos = HydroDataGet[target].Pos;
      vel = HydroDataGet[target].Vel;
      h_i = HydroDataGet[target].Hsml;
      mass = HydroDataGet[target].Mass;
      dhsmlDensityFactor = HydroDataGet[target].DhsmlDensityFactor;
      rho = HydroDataGet[target].Density;
      pressure = HydroDataGet[target].Pressure;
      timestep = HydroDataGet[target].Timestep;
      soundspeed_i = sqrt(GAMMA * pressure / rho);
      f1 = HydroDataGet[target].F1;
    }

  /* initialize variables before SPH loop is started */
  acc[0] = acc[1] = acc[2] = dtEntropy = 0;
  maxSignalVel = 0;

  p_over_rho2_i = pressure / (rho * rho) * dhsmlDensityFactor;
  h_i2 = h_i * h_i;

  /* Now start the actual SPH computation for this particle */
  startnode = All.MaxPart;
  do
    {
      numngb = ngb_treefind_pairs(&pos[0], h_i, &startnode);

      for(n = 0; n < numngb; n++)
	{
	  j = Ngblist[n];

	  dx = pos[0] - P[j].Pos[0];
	  dy = pos[1] - P[j].Pos[1];
	  dz = pos[2] - P[j].Pos[2];

#ifdef PERIODIC			/*  find the closest image in the given box size  */
	  if(dx > boxHalf_X)
	    dx -= boxSize_X;
	  if(dx < -boxHalf_X)
	    dx += boxSize_X;
	  if(dy > boxHalf_Y)
	    dy -= boxSize_Y;
	  if(dy < -boxHalf_Y)
	    dy += boxSize_Y;
	  if(dz > boxHalf_Z)
	    dz -= boxSize_Z;
	  if(dz < -boxHalf_Z)
	    dz += boxSize_Z;
#endif
	  r2 = dx * dx + dy * dy + dz * dz;
	  h_j = SphP[j].Hsml;
	  if(r2 < h_i2 || r2 < h_j * h_j)
	    {
	      r = sqrt(r2);
	      if(r > 0)
		{
		  p_over_rho2_j = SphP[j].Pressure / (SphP[j].Density * SphP[j].Density);
		  soundspeed_j = sqrt(GAMMA * p_over_rho2_j * SphP[j].Density);
		  dvx = vel[0] - SphP[j].VelPred[0];
		  dvy = vel[1] - SphP[j].VelPred[1];
		  dvz = vel[2] - SphP[j].VelPred[2];
		  vdotr = dx * dvx + dy * dvy + dz * dvz;

		  if(All.ComovingIntegrationOn)
		    vdotr2 = vdotr + hubble_a2 * r2;
		  else
		    vdotr2 = vdotr;

                  if(r2 < h_i2)
                    {
                      hinv = 1.0 / h_i;
#ifndef  TWODIMS
                      hinv4 = hinv * hinv * hinv * hinv;
#else
                      hinv4 = hinv * hinv * hinv / boxSize_Z;
#endif
                      u = r * hinv;
                      if(u < 0.5)
                        dwk_i = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
                      else
                        dwk_i = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
                    }
                  else
                    {
                      dwk_i = 0;
                    }

                  if(r2 < h_j * h_j)
                    {
                      hinv = 1.0 / h_j;
#ifndef  TWODIMS
                      hinv4 = hinv * hinv * hinv * hinv;
#else
                      hinv4 = hinv * hinv * hinv / boxSize_Z;
#endif
                      u = r * hinv;
                      if(u < 0.5)
                        dwk_j = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
                      else
                        dwk_j = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
                    }
                  else
                    {
                      dwk_j = 0;
                    }

                  if(soundspeed_i + soundspeed_j > maxSignalVel)
                    maxSignalVel = soundspeed_i + soundspeed_j;

                  if(vdotr2 < 0)        /* ... artificial viscosity */
                    {
                      mu_ij = fac_mu * vdotr2 / r;      /* note: this is negative! */

                      vsig = soundspeed_i + soundspeed_j - 3 * mu_ij;

                      if(vsig > maxSignalVel)
                        maxSignalVel = vsig;

		      if(P[j].Ti_endstep == All.Ti_Current)
			if(vsig > SphP[j].MaxSignalVel)
			  SphP[j].MaxSignalVel = vsig;

                      rho_ij = 0.5 * (rho + SphP[j].Density);
                      f2 =
                        fabs(SphP[j].DivVel) / (fabs(SphP[j].DivVel) + SphP[j].CurlVel +
                                                0.0001 * soundspeed_j / fac_mu / SphP[j].Hsml);

                      visc = 0.25 * All.ArtBulkViscConst * vsig * (-mu_ij) / rho_ij * (f1 + f2);

                      /* .... end artificial viscosity evaluation */
#ifndef NOVISCOSITYLIMITER
                      /* make sure that viscous acceleration is not too large */
                      dt = imax(timestep, (P[j].Ti_endstep - P[j].Ti_begstep)) * All.Timebase_interval;
                      if(dt > 0 && (dwk_i + dwk_j) < 0)
                        {
                          visc = dmin(visc, 0.5 * fac_vsic_fix * vdotr2 /
                                      (0.5 * (mass + P[j].Mass) * (dwk_i + dwk_j) * r * dt));
                        }
#endif
                    }
                  else
                    visc = 0;

                  p_over_rho2_j *= SphP[j].DhsmlDensityFactor;

                  hfc_visc = 0.5 * P[j].Mass * visc * (dwk_i + dwk_j) / r;

                  hfc = hfc_visc + P[j].Mass * (p_over_rho2_i * dwk_i + p_over_rho2_j * dwk_j) / r;

                  acc[0] -= hfc * dx;
                  acc[1] -= hfc * dy;
                  acc[2] -= hfc * dz;

		  if(P[j].Ti_endstep == All.Ti_Current)
		    {
		      if(SphP[j].FeedbackFlag == 0)
			SphP[j].FeedbackFlag = 1;

		      if(maxSignalVel > SphP[j].MaxSignalVel)
			SphP[j].MaxSignalVel = maxSignalVel;

		      SphP[j].FeedAccel[0] -= hfc * dx;
		      SphP[j].FeedAccel[1] -= hfc * dy;
		      SphP[j].FeedAccel[2] -= hfc * dz;
		    }

		}
	    }
	}
      fflush(stdout);
    }
  while(startnode >= 0);

  /* Now collect the result at the right place */
  if(mode == 0)
    {
      SphP[target].MaxSignalVel = maxSignalVel;

      for(k = 0; k < 3; k++)
        SphP[target].FeedAccel[k] = acc[k];
    }
  else
    {
      HydroDataResult[target].MaxSignalVel = maxSignalVel;

      for(k = 0; k < 3; k++)
        HydroDataResult[target].Acc[k] = acc[k];
    }
}

#endif

ext: c

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
#ifndef buildin_40bf10d0_ba81_4f06_accf_8020d3698607_h
#define buildin_40bf10d0_ba81_4f06_accf_8020d3698607_h

#include <gslib\type.h>
#include <gslib\string.h>
#include <rathen\config.h>

__rathen_begin__

#define max_opr 4

struct oprnode
{
    int     prior;
    gchar   opr[max_opr];
    bool    unary;
};

struct oprinfo
{
    string      opr;
    bool        unary;

    oprinfo(const gchar* o, bool u): opr(o), unary(u) {}
    oprinfo(const string& s, bool u): opr(s), unary(u) {}
};

struct unary_oprinfo:
    public oprinfo
{
    unary_oprinfo(const gchar* o): oprinfo(o, true) {}
    unary_oprinfo(const string& s): oprinfo(s, true) {}
};

struct binary_oprinfo:
    public oprinfo
{
    binary_oprinfo(const gchar* o): oprinfo(o, false) {}
    binary_oprinfo(const string& s): oprinfo(s, false) {}
};

class opr_manager
{
public:
    /*
     * the operators:
       ( )
       ++ --
       - ! ~
       += -= *= /= %=
       * / %
       + -
       << >>
       < > <= >= == !=
       & | ^ && ||
       =
     */
    static const oprnode opr_list[];

private:
    opr_manager() {}

public:
    static opr_manager* get_singleton_ptr()
    {
        static opr_manager inst;
        return &inst;
    }
    const oprnode* get_unary_operator(const gchar* src) const;
    const oprnode* get_binary_operator(const gchar* src) const;
    const oprnode* get_unary_operator(const string& src) const { return get_unary_operator(src.c_str()); }
    const oprnode* get_binary_operator(const string& src) const { return get_binary_operator(src.c_str()); }
    const oprnode* get_operator(const oprinfo& info) const
    {
        return info.unary ? get_unary_operator(info.opr) :
            get_binary_operator(info.opr);
    }

    /*
     * Currently, only priority '1' was unary.
     * pass the 'opr' to this function, in case of any proper changes in the future.
     */
    bool is_unary(int prior, const gchar* opr) const { return prior == 1; }
    bool is_unary(int prior, const string& opr) const { return prior == 1; }
};

#define _opr_manager    opr_manager::get_singleton_ptr()

#define max_key 8

struct key_node
{
    int     pattern;
    gchar   key[max_key];
};

class key_manager
{
public:
    static const key_node key_list[];
    enum
    {
        /*
         * Respectively like:
         * stop;
         * exit(...);
         * else {...}
         * if(...) {...}
         */
        pat_err     = 0x00,
        pat_raw     = 0x01,
        pat_arg     = 0x02,
        pat_seg     = 0x04,
        pat_com     = 0x08,
    };

private:
    key_manager() {}

public:
    static key_manager* get_singleton_ptr()
    {
        static key_manager inst;
        return &inst;
    }
    const key_node* get_key(const gchar* src) const;
    const key_node* get_key(const string& src) const { return get_key(src.c_str()); }
    int get_key_pattern(const gchar* src) const;
    int get_key_pattern(const string& src) const { return get_key_pattern(src.c_str()); }
};

#define _key_manager    key_manager::get_singleton_ptr()

__rathen_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#ifndef _PARAMETERS_H_
#define _PARAMETERS_H_
#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include "../Headers/voxel.h"
#include <math.h>
#include <gsl/gsl_rng.h>

using namespace std;

class Parameters{
//Atributs
	public:
	int it; //Contador
	int nt;	// Nombre de passos total
   	int ntout; // Nombre de pasos per l'output
	int neq; //Number of equations
	float dt; // Pas de temps
	float t;

	string Dir_Output;
    string Dir_Input;
    string NumPar_Dat;

	string Data_Dat;
	string Init_Cond_Dat;
   	int nNodes;
   	int nNodes_inh;
    string Kex_Dat;
    string Kin_Dat;

	float ** Kex; //Coupling matrix for excitatory population
	float ** Kin; //Coupling matrix for inhibitory population
	voxel * V; //Voxel

	ofstream fout_data;
	ofstream fout_init_cond;

        public:
        Parameters();
        ~Parameters();

	//Mètodes

	void Read_Kex_Dat(void);
	void Read_Kin_Dat(void);
	void RefreshDelay(void);
	void Read_NumPar_Dat(void);

	void WriteData();
	void WriteInitialConditions(int);
	void openfiles();
	void closefiles();

        //private:

};

#endif


text:
//  SMCTC: rng.hh
//
//   Copyright Adam Johansen, 2008.
//
//   This file is part of SMCTC.
//
//   SMCTC 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 3 of the License, or
//   (at your option) any later version.
//
//   SMCTC 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 SMCTC.  If not, see <http://www.gnu.org/licenses/>.
//

//! \file
//! \brief Random number generation.
//!
//! This file contains the definitions for the smc::rng and smc::rnginfo class.
//! It wraps the random number generation facilities provided by the GSL and provides a convenient interfaces to access several of its more commonly-used features.

#ifndef __SMC_RNG_HH
#define __SMC_RNG_HH 1.0

extern "C" {
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
}

namespace smc {
  ///A gsl-rng information handling class (not templated)
  class gslrnginfo {
  private:
    ///This is a null terminated array of the available random number generators.
    const gsl_rng_type** typePtArray;
    ///The number of available random number generators.
    int nNumber;

  protected:
    gslrnginfo();

  public:
    ///Returns a reference to the sole static instance of this class.
    static gslrnginfo* GetInstance();

    ///Returns the number of available random number generators.
    int GetNumber();
    ///Returns the name of random number generator number nIndex.
    const char* GetNameByIndex(int nIndex);
    ///Returns a pointer to random number generator nIndex.
    const gsl_rng_type* GetPointerByIndex(int nIndex);
    ///Returns a pointer to the random number generator with name szName (or null if it doesn't exist).
    const gsl_rng_type* GetPointerByName(const char* szName);
  };

  ///The global application instance of the gslrnginfo class:
  extern gslrnginfo rngset;

  ///A random number generator class.

  ///    At present this serves as a wrapper for the gsl random number generation code.
  class rng {
  private:
    ///This is the type of random number generator underlying the class.
    const gsl_rng_type* type;
    ///This is a pointer to the internal workspace of the rng including its current state.
    gsl_rng* pWorkspace;

  public:
    ///Initialise the random number generator using default settings
    rng();
    ///Initialise the random number generator using the default seed for the type
    rng(const gsl_rng_type* Type);
    ///Initialise the random number generator using specified type and seed
    rng(const gsl_rng_type* Type,unsigned long int lSeed);

    ///Free the workspace allocated for random number generation
    ~rng();

    ///Provide access to the raw random number generator
    gsl_rng* GetRaw(void);

    ///Generate a multinomial random vector with parameters (n,w[1:k]) and store it in X
    void Multinomial(unsigned n, unsigned k, const double* w, unsigned* X);
    ///Returns a random integer generated uniformly between the minimum and maximum values specified
    long UniformDiscrete(long lMin, long lMax);

    ///Returns a random number generated from a Beta distribution with the specified parameters.
    double Beta(double da, double db);
    ///Returns a random number generated from a Cauchy distribution with the specified scale parameter.
    double Cauchy(double dScale);
    ///Returns a random number generated from an exponential distribution with the specified mean.
    double Exponential(double dMean);
    ///Return a random number generated from a gamma distribution with shape alpha and scale beta.
    double Gamma(double dAlpha, double dBeta);
    ///Returns a random number generated from a Laplacian distribution with the specified scale.
    double Laplacian(double dScale);
    ///Returns a random number generated from a Lognormal distribution of location mu and scale sigma
    double Lognormal(double dMu, double dSigma);
    ///Return a random number generated from a normal distribution with a specified mean and standard deviation
    double Normal(double dMean, double dStd);
    ///Return a random number generated from a standard normal distribution
    double NormalS(void);
    ///Returns a random number from a normal distribution, conditional upon it exceeding the specified threshold.
    double NormalTruncated(double dMean, double dStd, double dThreshold);
    ///Return a student-t random number generated with a specified number of degrees of freedom
    double StudentT(double dDF);
    ///Return a random number generated uniformly between dMin and dMax
    double Uniform(double dMin, double dMax);
    ///Returns a random number generated from the standard uniform[0,1) distribution
    double UniformS(void);
  };
}

#endif


text:
#include "lah.h"

#ifdef HAVE_LAPACK
#include <lapacke.h>

lah_Return lah_eigenValue(lah_mat *P, lah_value *S, lah_mat *Z)
{
    lapack_int res, nEig;
    lapack_int *isuppz;
    if (P == NULL || S == NULL || Z == NULL
        || P->nR != P->nC)
    {
        return lahReturnParameterError;
    }
    /* allocate space for the output parameters and workspace arrays */
    isuppz = malloc(2 * P->nR * sizeof(lapack_int));

    /* Data allocation */
    /* Sollte A-> data nur dreiecksmatrix enthalten ?*/
    res = SYEVR(LAH_LAPACK_LAYOUT, 'V', 'A', 'L', P->nR, P->data, P->nR,
                0, 0, 0, 0, LAMCH('S'), &nEig, S, Z->data, P->nR, isuppz);

    /* free support array */
    free(isuppz);
    if (res > 0)
        return lahReturnExternError;
    else if (res < 0)
        return lahReturnParameterError;
    else
        return lahReturnOk;
}
#endif


text:
/* modulo-bias - test for modulo bias in poor random number generation;
 * the bias does not appear where RAND_MAX is much larger than m. for a
 * three-sided dice roll, the last RAND_MAX to fail what I think is a
 * chi-squared 0.05 test is 257; given a system RAND_MAX of a much
 * larger value, this would likely only be of concern if the modulo m is
 * very, very large, or if an 8-bit system is being used. better results
 * are obtained where m evenly divides RAND_MAX; RAND_MAX of 100000 and
 * m = 4 only fails at 62 and below (assuming no errors on my part)
 *
 * actually, the exact numbers depend on the rand() implementation and
 * the seed; Mac OS X rand() produced the 257 for an arbitrary seed;
 * this version tries out multiple seeds (TRIALS) and emits what is
 * likely the maximum failure of the fitness test for that seed.
 * However! it still remains that if RAND_MAX is much larger than the
 * range, there is no (statistically significant) modulo bias */

#include <err.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

#include <gsl/gsl_cdf.h>

/* how many results to generate (3 == (0,1,2)) and how many tests to
 * run; TESTNUM will need to be high enough to distinguish signal from
 * the noise, but not too long that too much CPU is wasted */
#define RANGE 3
#define TESTNUM 99999

/* CHI^2 test foo, assuming uniform distribution of the RANGE values
 * over the TESTNUM runs, and 95% confidence */
#define EXPECTED TESTNUM / (double) RANGE
#define DEG_OF_FREE RANGE - 1
#define FITNESS 0.05

/* this may have to go up if larger RANGE are being tested, e.g. some
 * huge prime vs. some other huge prime, or otherwise scale in some
 * relation to RANGE and randmax */
#define NOBIAS_RUN 100

/* how many (hopefully) different random seeds to do the test under */
#define TRIALS 500

#define randof(x) ((rand() % x) % RANGE)

unsigned int results[RANGE];

int main(void)
{
    unsigned int i, maxfail, nobias, randmax, testnum, trialnum;
    double sq;

#ifdef __OpenBSD__
    if (pledge("stdio", NULL) == -1)
        err(1, "pledge failed");
#endif

    //setvbuf(stdout, (char *)NULL, _IOLBF, (size_t) 0);

    for (trialnum = 0; trialnum < TRIALS; trialnum++) {
        /* pseudorandom within pseudorandom, per Dune */
        srand(arc4random());

        nobias = 0;

        for (randmax = RANGE; randmax < UINT_MAX; randmax += 1) {
            for (testnum = 0; testnum < TESTNUM; testnum++) {
                results[randof(randmax)]++;
            }
            sq = 0.0;
            for (i = 0; i < RANGE; i++) {
                sq += pow((double) results[i] - EXPECTED, 2) / EXPECTED;
                results[i] = 0; // reset for the next test
            }
            if (gsl_cdf_chisq_Q(sq, DEG_OF_FREE) <= FITNESS) {
                maxfail = randmax;
                nobias = 0;
            } else {
                if (nobias++ >= NOBIAS_RUN) {
                    break;
                }
            }
        }
        printf("%d\n", maxfail);
    }

    exit(EXIT_SUCCESS);
}

max_stars_repo_name: thrig/scripts

max_issues_repo_name: thrig/scripts

max_forks_repo_name: thrig/scripts


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira

 This program 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, 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., 51 Franklin Street, Fifth Floor, Boston, MA
 02110-1301, USA.

 $Id: minimizer_low.c 11822 2014-02-25 21:18:45Z dstrubbe $
*/

#include <config.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_multimin.h>

#include "string_f.h"

/* A. First, the interface to the gsl function that calculates the minimum
   of an N-dimensional function, with the knowledge of the function itself
   and its gradient. */

/* This is a type used to communicate with Fortran; func_d is the type of the
   interface to a Fortran subroutine that calculates the function and its
   gradient. */
typedef void (*func_d)(const int*, const double*, double*, const int*, double*);
typedef struct{
  func_d func;
} param_fdf_t;

/* Compute both f and df together. */
static void
my_fdf (const gsl_vector *v, void *params, double *f, gsl_vector *df)
{
  double *x, *gradient, ff2;
  int i, dim, getgrad;
  param_fdf_t * p;

  p = (param_fdf_t *) params;

  dim = v->size;
  x = (double *)malloc(dim*sizeof(double));
  gradient = (double *)malloc(dim*sizeof(double));

  for(i=0; i<dim; i++) x[i] = gsl_vector_get(v, i);
  getgrad = (df == NULL) ? 0 : 1;

  p->func(&dim, x, &ff2, &getgrad, gradient);

  if(f != NULL)
    *f = ff2;

  if(df != NULL){
    for(i=0; i<dim; i++)
      gsl_vector_set(df, i, gradient[i]);
  }

  free(x); free(gradient);
}

static double
my_f (const gsl_vector *v, void *params)
{
  double val;

  my_fdf(v, params, &val, NULL);
  return val;
}

/* The gradient of f, df = (df/dx, df/dy). */
static void
my_df (const gsl_vector *v, void *params, gsl_vector *df)
{
  my_fdf(v, params, NULL, df);
}

typedef void (*print_f_ptr)(const int*, const int*, const double*, const double*, const double*, const double*);

int FC_FUNC_(oct_minimize, OCT_MINIMIZE)
     (const int *method, const int *dim, double *point, const double *step, const double *line_tol,
      const double *tolgrad, const double *toldr, const int *maxiter, func_d f,
      const print_f_ptr write_info, double *minimum)
{
  int iter = 0;
  int status;
  double maxgrad, maxdr;
  int i;
  double * oldpoint;
  double * grad;

  const gsl_multimin_fdfminimizer_type *T = NULL;
  gsl_multimin_fdfminimizer *s;
  gsl_vector *x;
  gsl_vector *absgrad, *absdr;
  gsl_multimin_function_fdf my_func;

  param_fdf_t p;

  p.func = f;

  oldpoint = (double *) malloc(*dim * sizeof(double));
  grad     = (double *) malloc(*dim * sizeof(double));

  my_func.f = &my_f;
  my_func.df = &my_df;
  my_func.fdf = &my_fdf;
  my_func.n = *dim;
  my_func.params = (void *) &p;

  /* Starting point */
  x = gsl_vector_alloc (*dim);
  for(i=0; i<*dim; i++) gsl_vector_set (x, i, point[i]);

  /* Allocate space for the gradient */
  absgrad = gsl_vector_alloc (*dim);
  absdr = gsl_vector_alloc (*dim);

  //GSL recommends line_tol = 0.1;
  switch(*method){
  case 1:
    T = gsl_multimin_fdfminimizer_steepest_descent;
    break;
  case 2:
    T = gsl_multimin_fdfminimizer_conjugate_fr;
    break;
  case 3:
    T = gsl_multimin_fdfminimizer_conjugate_pr;
    break;
  case 4:
    T = gsl_multimin_fdfminimizer_vector_bfgs;
    break;
  case 5:
    T = gsl_multimin_fdfminimizer_vector_bfgs2;
    break;
  }

  s = gsl_multimin_fdfminimizer_alloc (T, *dim);

  gsl_multimin_fdfminimizer_set (s, &my_func, x, *step, *line_tol);
  do
    {
      iter++;
      for(i=0; i<*dim; i++) oldpoint[i] = point[i];

      /* Iterate */
      status = gsl_multimin_fdfminimizer_iterate (s);

      /* Get current minimum, point and gradient */
      *minimum = gsl_multimin_fdfminimizer_minimum(s);
      for(i=0; i<*dim; i++) point[i] = gsl_vector_get(gsl_multimin_fdfminimizer_x(s), i);
      for(i=0; i<*dim; i++) grad[i] = gsl_vector_get(gsl_multimin_fdfminimizer_gradient(s), i);

      /* Compute convergence criteria */
      for(i=0; i<*dim; i++) gsl_vector_set(absdr, i, fabs(point[i]-oldpoint[i]));
      maxdr = gsl_vector_max(absdr);
      for(i=0; i<*dim; i++) gsl_vector_set(absgrad, i, fabs(grad[i]));
      maxgrad = gsl_vector_max(absgrad);

      /* Print information */
      write_info(&iter, dim, minimum, &maxdr, &maxgrad, point);

      /* Store infomation for next iteration */
      for(i=0; i<*dim; i++) oldpoint[i] = point[i];

      if (status)
        break;

      if ( (maxgrad <= *tolgrad) || (maxdr <= *toldr) ) status = GSL_SUCCESS;
      else status = GSL_CONTINUE;
    }
  while (status == GSL_CONTINUE && iter <= *maxiter);

  if(status == GSL_CONTINUE) status = 1025;

  gsl_multimin_fdfminimizer_free (s);
  gsl_vector_free (x); gsl_vector_free(absgrad); gsl_vector_free(absdr);

  free(oldpoint);
  free(grad);

  return status;
}

/* B. Second, the interface to the gsl function that calculates the minimum
   of a one-dimensional function, with the knowledge of the function itself,
   but not its gradient. */

/* This is a type used to communicate with Fortran; func1 is the type of the
   interface to a Fortran subroutine that calculates the function and its
   gradient. */
typedef void (*func1)(const double*, double*);
typedef struct{
  func1 func;
} param_f1_t;

double fn1(double x, void * params)
{
  param_f1_t * p = (param_f1_t* ) params;
  double fx;
  p->func(&x, &fx);
  return fx;
}

void FC_FUNC_(oct_1dminimize, OCT_1DMINIMIZE)(double *a, double *b, double *m, func1 f, int *status)
{
  int iter = 0;
  int max_iter = 100;
  const gsl_min_fminimizer_type *T;
  gsl_min_fminimizer *s;
  gsl_function F;
  param_f1_t p;

  p.func = f;

  F.function = &fn1;
  F.params = (void *) &p;

  T = gsl_min_fminimizer_brent;
  s = gsl_min_fminimizer_alloc (T);

  *status = gsl_min_fminimizer_set (s, &F, *m, *a, *b);

  gsl_set_error_handler_off();

  do
    {
      iter++;
      *status = gsl_min_fminimizer_iterate (s);

      *m = gsl_min_fminimizer_x_minimum (s);
      *a = gsl_min_fminimizer_x_lower (s);
      *b = gsl_min_fminimizer_x_upper (s);

      *status = gsl_min_test_interval (*a, *b, 0.00001, 0.0);

      /*if (*status == GSL_SUCCESS) printf ("Converged:\n");*/
      /*printf ("%5d [%.7f, %.7f] %.7f \n", iter, *a, *b,*m);*/
    }
  while (*status == GSL_CONTINUE && iter < max_iter);
  gsl_min_fminimizer_free(s);

}

/* C. Third, the interface to the gsl function that calculates the minimum
   of an N-dimensional function, with the knowledge of the function itself,
   but not its gradient. */

/* This is a type used to communicate with Fortran; funcn is the type of the
   interface to a Fortran subroutine that calculates the function and its
   gradient. */
typedef void (*funcn)(int*, double*, double*);
typedef struct{
  funcn func;
} param_fn_t;

double fn(const gsl_vector *v, void * params)
{
  double val;
  double *x;
  int i, dim;
  param_fn_t * p;

  p = (param_fn_t *) params;
  dim = v->size;
  x = (double *)malloc(dim*sizeof(double));

  for(i=0; i<dim; i++) x[i] = gsl_vector_get(v, i);
  p->func(&dim, x, &val);

  free(x);
  return val;
}

typedef void (*print_f_fn_ptr)(const int*, const int*, const double*, const double*, const double*);

int FC_FUNC_(oct_minimize_direct, OCT_MINIMIZE_DIRECT)
     (const int *method, const int *dim, double *point, const double *step,
      const double *toldr, const int *maxiter, funcn f,
      const print_f_fn_ptr write_info, double *minimum)
{
  int iter = 0, status, i;
  double size;

  const gsl_multimin_fminimizer_type *T = NULL;
  gsl_multimin_fminimizer *s = NULL;
  gsl_vector *x, *ss;
  gsl_multimin_function my_func;

  param_fn_t p;
  p.func = f;

  my_func.f = &fn;
  my_func.n = *dim;
  my_func.params = (void *) &p;

  /* Set the initial vertex size vector */
  ss = gsl_vector_alloc (*dim);
  gsl_vector_set_all (ss, *step);

  /* Starting point */
  x = gsl_vector_alloc (*dim);
  for(i=0; i<*dim; i++) gsl_vector_set (x, i, point[i]);

  switch(*method){
  case 6:
    T = gsl_multimin_fminimizer_nmsimplex;
    break;
  }

  s = gsl_multimin_fminimizer_alloc (T, *dim);
  gsl_multimin_fminimizer_set (s, &my_func, x, ss);

  do
    {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);

      if(status) break;

      *minimum = gsl_multimin_fminimizer_minimum(s);
      for(i=0; i<*dim; i++) point[i] = gsl_vector_get(gsl_multimin_fminimizer_x(s), i);

      size = gsl_multimin_fminimizer_size (s);
      status = gsl_multimin_test_size (size, *toldr);

      write_info(&iter, dim, minimum, &size, point);

    }
  while (status == GSL_CONTINUE && iter < *maxiter);

  if(status == GSL_CONTINUE) status = 1025;

  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free(s);
  return status;
}

max_stars_repo_name: gimunu/octopus-metric

max_issues_repo_name: gimunu/octopus-metric

max_forks_repo_name: gimunu/octopus-metric


text:
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define STRASSEN_LOWER_BOUND 128
#define NAN -1

typedef float REAL_NUMBER;
typedef struct
{
    int nrows;
    int ncols;
    REAL_NUMBER* data;
} matrix;

matrix read_matrix(const char* file_name) {
    FILE* fp = fopen(file_name, "r");
    matrix m;
    m.nrows = 0;
    m.ncols = 0;

    // get ncols and nrows
    REAL_NUMBER tempf = 0;
    char tempc = 1;
    while (tempc != '\n') {
        fscanf(fp, "%f", &tempf);
        fscanf(fp, "%c", &tempc);
        m.ncols++;
    }
    m.nrows++;
    while (fscanf(fp, "%c", &tempc) != EOF) {
        if (tempc == '\n') {
            m.nrows++;
        }
    }

    // read data
    rewind(fp);
    m.data = (REAL_NUMBER*)malloc(m.ncols * m.nrows * sizeof(REAL_NUMBER));
    REAL_NUMBER* pread = m.data;
    for (int i = 0; i < m.ncols * m.nrows; i++) {
        // fscanf(fp, "%f", &m.data[i]);
        // fscanf(fp, "%f", m.data + i);
        fscanf(fp, "%f", pread);
        pread++;
    }

    fclose(fp);
    return m;
}

void print_matrix(matrix* m, const char* file_name) {
    FILE* fp = fopen(file_name, "w");

    for (int i = 0; i < m->ncols * m->nrows; i++) {
        if (i > 0 && i % m->ncols == 0) {
            fputc('\n', fp);
        }
        fprintf(fp, "%g ", m->data[i]);
    }
    fputc('\n', fp);  // last blank line

    fclose(fp);
}

matrix create_matrix(int nrows, int ncols, REAL_NUMBER fill) {
    matrix m = {
        .nrows = nrows,
        .ncols = ncols,
        .data = (REAL_NUMBER*)malloc(nrows * ncols * sizeof(REAL_NUMBER))};

    if (fill != NAN) {
        for (int i = 0; i < nrows * ncols; i++) {
            m.data[i] = fill;
        }
    }

    return m;
}

int main(int argc, char const* argv[]) {
    if (argc != 4) {
        printf("Wrong number of arguments.");
        exit(EXIT_FAILURE);
    }

    const enum CBLAS_ORDER Order = CblasRowMajor;
    const enum CBLAS_TRANSPOSE TransA = CblasNoTrans;
    const enum CBLAS_TRANSPOSE TransB = CblasNoTrans;
    const float alpha=1;
    const float beta=0;

    double start, end;

    start = clock();

    matrix m1 = read_matrix(argv[1]);
    matrix m2 = read_matrix(argv[2]);

    end = clock();
    printf("Read file time: %fs\n", (double)(end - start) / CLOCKS_PER_SEC);

    start = clock();

    matrix res = create_matrix(m1.nrows, m2.ncols, NAN);
    cblas_sgemm(Order, TransA, TransB, m1.nrows, m2.ncols, m1.ncols, alpha, m1.data, m1.ncols, m2.data, m2.ncols, beta, res.data, res.ncols);

    end = clock();
    printf("Multiplication time: %fs\n", (double)(end - start) / CLOCKS_PER_SEC);

    print_matrix(&res, argv[3]);
    return 0;
}

// gcc -I D:\OpenBLAS-0.3.18-x64\include -L D:\OpenBLAS-0.3.18-x64\lib .\matmul_OpenBLAS.c -lopenblas -o matmul_OpenBLAS


text:
#ifndef CIMPLE_POLYTOPE_LIBRARY_CIMPLE_H
#define CIMPLE_POLYTOPE_LIBRARY_CIMPLE_H

#include <math.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_matrix.h>
#include "cimple_gsl_library_extension.h"
#include <gurobi_c.h>
#include "setoper.h"
#include <cdd.h>
#include "cimple_auxiliary_functions.h"
#include "cimple_minksum_wrapper.h"

/**
 * H left side of polytope (sometimes noted A or L)
 * G right side of polytope (sometimes noted b or M)
 *
 *      H.x <= G
 *
 * Chebyshev center is a possible definition  of the "center" of the polytope.
 */
typedef struct polytope{

    gsl_matrix * H;
    gsl_vector * G;
    double *chebyshev_center;

}polytope;

/**
 * @brief "Constructor" Dynamically allocates the space a polytope needs
 * @param k H.size1 == G.size
 * @param n H.size2
 * @return
 */
struct polytope *polytope_alloc(size_t k, size_t n);

/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the polytope
 * @param polytope
 */
void polytope_free(polytope *polytope);

/**
 * Subdivision of abstract state containing additionally to the polytope also safe mode instructions
 * The array of polytopes "polytope **safe_mode" contains N polytopes the system has to go through to reach the invariant set
 */
typedef struct cell{

    polytope **safe_mode;
    polytope *polytope_description;

}cell;

/**
 * @brief "Constructor" Dynamically allocates the space a cell needs
 * @param k cell.polytope.H.size1 == G.size
 * @param n cell.polytope.HH.size2
 * @return
 */
struct cell *cell_alloc(size_t k,
                        size_t n,
                        int time_horizon);

/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the cell
 * @param cell
 */
void cell_free(cell *cell);

/**
 * Convex region of several (cells_count) polytopes (array of polytopes**)
 * hull_over_polytopes is the convex hull of polytopes in that abstract state
 * hull.A = NULL if only one polytope exists in that region
 */
typedef struct abstract_state{

    polytope *convex_hull;

    int cells_count;
    cell **cells;

    int transitions_in_count;
    struct abstract_state** transitions_in;
    int transitions_out_count;
    struct abstract_state** transitions_out;

    int distance_invariant_set;
    struct abstract_state * next_state;
    cell *invariant_set;

}abstract_state;

/**
 * @brief "Constructor" Dynamically allocates the space a region of polytope needs
 *
 * Allocates memory space according to the cells_count and their respective sizes
 *
 * @param k Array with the number of rows of each polytope dim([cells_count])
 * @param k_hull number of rows of convex hull polytope of the region
 * @param n system_dynamics size (e.g. s_dyn.A.size1)
 * @param cells_count
 * @return
 */
struct abstract_state *abstract_state_alloc(size_t *k,
                                            size_t k_hull,
                                            size_t n,
                                            int transitions_in_count,
                                            int transitions_out_count,
                                            int cells_count,
                                            int time_horizon);

/**
 * @brief "Destructor" Deallocates the dynamically allocated memory of the region of polytopes
 * @param abstract_state
 */
void abstract_state_free(abstract_state * abstract_state);

/**
 * @brief Converts two C arrays to a polytope consistent of a left side matrix (i.e. H) and right side vector (i.e. G)
 * @param polytope empty polytope with allocated memory
 * @param k number of rows of polytope (H.size1 or G.size)
 * @param n dimension of polytope = system_dynamics size (e.g. s_dyn.A.size1)
 * @param left_side C array to be converted to left side matrix
 * @param right_side C array to be converted to right side vector
 * @param name name of polytope, will be displayed to user terminal to inform about successfull initialization
 */
void polytope_from_arrays(polytope *polytope,
                          double *left_side,
                          double *right_side,
                          double *cheby,
                          char*name);

/**
 * @brief Converts a polytope in gsl form to cdd constraint form
 * @param original
 * @param err
 * @return
 */
dd_PolyhedraPtr polytope_to_cdd(polytope *original,
                                dd_ErrorType *err);

/**
 * @brief Converts a polytope in cdd constraint form to gsl form
 * @param original
 * @return
 */
polytope* cdd_to_polytope(dd_PolyhedraPtr *original);

/**
 * @brief Generate a polytope representing a scaled unit cube
 * @param scale
 * @param dimensions
 * @return gsl polytope
 */
polytope * polytope_scaled_unit_cube(double scale,
                                     int dimensions);

/**
 * @brief Checks whether a state is in a certain polytope
 * @param polytope
 * @param x state to be checked
 * @return 0 if state is not in polytope or 1 if it is
 */
bool polytope_check_state(polytope *polytope, gsl_vector *x);

/**
 * @brief Check whether P1 \ subset P2
 * @param P1
 * @param P2
 * @return
 */
bool polytope_is_subset(polytope *P1,
                        polytope *P2);

/**
 * @brief Unite inequalities of P1 and P2 in new polytope and remove redundancies
 * @param P1
 * @param P2
 * @return
 */
polytope * polytope_unite_inequalities(polytope *P1,
                                       polytope *P2);

/**
 * @brief Project gsl polytope of n+ dimensions to the first n dimensions
 * @param original
 * @param n
 * @return gsl polytope
 */
polytope* polytope_projection(polytope * original,
                              size_t n);
/**
 * @brief Multiplication of a polytope with a matrix
 * @param original
 * @param scale
 * @return gsl polytope
 */
polytope * polytope_linear_transform(polytope *original,
                                     gsl_matrix *scale);

/**
 * @brief Remove redundancies from gsl polytope inequalities
 * @param original
 * @return
 */
polytope * polytope_minimize(polytope *original);

/**
 * @brief Compute Minkowski sum of two polytopes
 * @param P1
 * @param P2
 * @return
 */
polytope * polytope_minkowski(polytope *P1,
                              polytope *P2);

/**
 * @brief Compute Pontryagin difference of two polytopes C=A-B s.t.:
 * A-B = {c \in A-B| c+b \in A, \forall b \in B}
 * @param P1
 * @param P2
 * @return
 */
polytope * polytope_pontryagin(polytope* P1,
                               polytope* P2);

/**
 * @brief Set up constraints in quadratic problem for GUROBI
 * @param constraints
 * @param model
 * @param N
 * @return
 */
int polytope_to_constraints_gurobi(polytope *constraints,
                                   GRBmodel *model,
                                   size_t N);

/**
 * @brief Generate a polytope representing a scaled unit cube
 * @param scale
 * @param dimensions
 * @return cdd polytope
 */
dd_PolyhedraPtr cdd_scaled_unit_cube(double scale,
                                     int dimensions);

/**
 * @brief Project cdd polytope of n+ dimensions to the first n dimensions
 * @param original
 * @param n
 * @return cdd polytope
 */
void cdd_projection(dd_PolyhedraPtr *original,
                    dd_PolyhedraPtr *new,
                    size_t n,
                    dd_ErrorType *err);

/**
 * @brief Remove redundancies from cdd polytope inequalities
 * @param original
 * @return cdd polytope
 */
dd_PolyhedraPtr cdd_minimize(dd_PolyhedraPtr *original,
                             dd_ErrorType *err);

#endif //CIMPLE_POLYTOPE_LIBRARY_CIMPLE_H

max_stars_repo_name: shaesaert/TuLiPXML

max_issues_repo_name: shaesaert/TuLiPXML

max_forks_repo_name: shaesaert/TuLiPXML


text:
/***********************************
V1 IMPLEMENTATION USING RANDOM POINTS
************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <cblas.h>
#include <string.h>
#include <mpi.h>

#include "../inc/knnring.h"

int main(int argc, char *argv[]){
  int p_rank;
  int p_size;
  int n;
  int d;
  int k;
  double *X;

  //INITIALIZE MPI COMMUNICATIONS
  MPI_Init(&argc,&argv);
  if(argc<4){
    printf("ERROR: GIVE CORRECT NUMBERS OF ARGUMENTS n,d,k\n");
    exit(1);
  }else{
    n = atoi(argv[1]);
    d = atoi(argv[2]);
    k = atoi(argv[3]);
  }
  MPI_Comm_size(MPI_COMM_WORLD, &p_size);
  MPI_Comm_rank(MPI_COMM_WORLD, &p_rank);
  MPI_Status status;                      // STATUS COMMUNICATION
  int tag = 1;                            //COMMUNICATION TAG

  struct timespec ts_start;
	struct timespec ts_end;
  double time;

  if(p_rank == 0){

    X = (double *)malloc(n*d*sizeof(double));

    for(int i=1;i<p_size;i++){
      create_points(X,n*d);
      MPI_Send(X,n*d,MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
    }

    //compute points for myself
    create_points(X,n*d);

    clock_gettime(CLOCK_MONOTONIC, &ts_start);
    knnresult results = distrAllkNN(X,n,d,k);
    clock_gettime(CLOCK_MONOTONIC, &ts_end);
    time = 1000000*(double)(ts_end.tv_sec-ts_start.tv_sec)+(double)(ts_end.tv_nsec-ts_start.tv_nsec)/1000;
    printf("TOTAL V1 time: %lf\n",time);

    //WRITE RESULTS TO FILE
    char str[200];
    snprintf(str,sizeof(str),"results/dist_knn_%d_%d_%d.txt",n,d,k);
    //write_to_file(str,results.nidx,results.ndist,n*k,results.k);

    for(int p=1;p<p_size;p++){
      MPI_Recv(results.nidx, n*k,    MPI_INT,p,tag,MPI_COMM_WORLD, &status);
      MPI_Recv(results.ndist,n*k, MPI_DOUBLE,p,tag,MPI_COMM_WORLD, &status);
      //write_to_file(str,results.nidx,results.ndist,n*k,results.k);
    }

  }else{
    //receive points
    X = (double *)malloc(n*d*sizeof(double));
    MPI_Recv(X,n*d,MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);

    clock_gettime(CLOCK_MONOTONIC, &ts_start);
    knnresult results = distrAllkNN(X,n,d,k);
    clock_gettime(CLOCK_MONOTONIC, &ts_end);
    time = 1000000*(double)(ts_end.tv_sec-ts_start.tv_sec)+(double)(ts_end.tv_nsec-ts_start.tv_nsec)/1000;
    //printf("TOTAL V1 time: %lf from process %d\n",time,p_rank);

    MPI_Send(results.nidx, n*k, MPI_INT   ,0,tag,MPI_COMM_WORLD);
    MPI_Send(results.ndist,n*k, MPI_DOUBLE,0,tag,MPI_COMM_WORLD);
  }

  //execution is over
  MPI_Finalize();
  printf("DONE ....=>EXITING\n");
  return 0;
}


text:
#include <math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include "gauss.h"

#define SQRT_2_PI     2.506628274631      /* sqrt(2*PI) */
#define LOG_SQRT_2_PI 0.918938533204673   /* log(sqrt(2*PI)) */

gauss_t gauss_init_prec(double prec_mean, double prec){
	gauss_t r;
	r.prec_mean = prec_mean;
	r.prec      = prec;
	r.var       = 1.0 / prec;
	r.mean      = prec_mean / prec;
	r.stddev    = sqrt(r.var);
	return r;
}

gauss_t gauss_init_var(double mean, double var){
	gauss_t r;
	r.var       = var;
	r.mean      = mean;
	r.stddev    = sqrt(var);
	r.prec      = 1.0 / var;
	r.prec_mean = r.prec * r.mean;
	return r;
}

gauss_t gauss_init_std(double mean, double stddev){
	gauss_t r;
	r.var       = stddev * stddev;
	r.mean      = mean;
	r.stddev    = stddev;
	r.prec      = 1.0 / r.var;
	r.prec_mean = r.prec * r.mean;
	return r;
}

double gauss_log_product_normal(gauss_t a, gauss_t b){
	if(a.prec == 0 || b.prec == 0) return 0.0;

	double var_sum   = a.var + b.var;
	double mean_diff = a.mean - b.mean;

	return -LOG_SQRT_2_PI - (log(var_sum)/2.0) - ((mean_diff * mean_diff) / (2.0 * var_sum));
}

double gauss_log_ratio_normal(gauss_t a, gauss_t b){
	if(a.prec == 0 || b.prec == 0) return 0.0;

	double var_diff  = a.var - b.var;
	double mean_diff = a.mean - b.mean;

	return log(b.var) + LOG_SQRT_2_PI - (log(var_diff)/2.0) + ((mean_diff * mean_diff) / (2.0 * var_diff));
}

double gauss_cdf_P(gauss_t a, double x){
	return gsl_cdf_gaussian_P(x - a.mean, a.stddev);
}

double gauss_cdf_Q(gauss_t a, double x){
	return gsl_cdf_gaussian_Q(x - a.mean, a.stddev);
}

gauss_t gauss_add(gauss_t a, gauss_t b){
	gauss_t r;
	r.mean      = a.mean + b.mean;
	r.var       = a.var  + b.var;
	r.prec      = 1.0 / r.var;
	r.prec_mean = r.prec * r.mean;
	r.stddev    = sqrt(r.var);
	return r;
}

gauss_t gauss_sub(gauss_t a, gauss_t b){
	gauss_t temp = b;
	temp.mean = -temp.mean;
	return gauss_add(a, temp);
}

double gauss_abs_diff(gauss_t a, gauss_t b){
	double prec_mean_diff = fabs(a.prec_mean - b.prec_mean);
	double prec_diff = sqrt(fabs(a.prec - b.prec));
	return (prec_mean_diff > prec_diff) ? prec_mean_diff : prec_diff;
}

gauss_t gauss_div(gauss_t a, gauss_t b){
	double prec      = a.prec      - b.prec;
	double prec_mean = a.prec_mean - b.prec_mean;
	return gauss_init_prec(prec_mean, prec);
}

gauss_t gauss_mul(gauss_t a, gauss_t b){
	double prec      = a.prec      + b.prec;
	double prec_mean = a.prec_mean + b.prec_mean;
	return gauss_init_prec(prec_mean, prec);
}

max_stars_repo_name: cosmicturtle/gamesetmatch

max_issues_repo_name: cosmicturtle/gamesetmatch

max_forks_repo_name: cosmicturtle/gamesetmatch


text:
#include "sparse-matrix.h"
#include "sparse-solver.h"
#include "sparse-solver-private.h"

#include <stdlib.h>
#include <string.h>
#if HAVE_MKL_LAPACK
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif
#include <assert.h>

void
sparse_lapack_init (SparseSolver *solver)
{

}

void
sparse_lapack_clear (SparseSolver *solver)
{

}

int
sparse_lapack_solve (SparseSolver *solver,
                     SparseMatrix *A,
                     void         *x,
                     const void   *b)
{
    lapack_int ret;

    assert (A->storage == SPARSE_MATRIX_STORAGE_DENSE);
    assert (A->shape[0] == A->shape[1]);

    lapack_int *pivot = malloc (A->shape[0] * sizeof (lapack_int));

    if (A->type == SPARSE_MATRIX_TYPE_DOUBLE)
    {
        memcpy (x, b, A->shape[1] * sizeof (double));
        ret = LAPACKE_dgesv (LAPACK_ROW_MAJOR,
                             A->shape[0],
                             1,
                             A->data,
                             A->shape[1],
                             pivot,
                             (double*) x,
                             1);
    }
    else if (A->type == SPARSE_MATRIX_TYPE_FLOAT)
    {
        memcpy (x, b, A->shape[1] * sizeof (float));
        ret = LAPACKE_sgesv (LAPACK_ROW_MAJOR,
                             A->shape[0],
                             1,
                             A->data,
                             A->shape[1],
                             pivot,
                             (float*) x,
                             1);
    }
    else
    {
        assert (0);
    }

    free (pivot);

    return ret == 0;
}

max_stars_repo_name: major-lab/mirbooking

max_issues_repo_name: major-lab/mirbooking

max_forks_repo_name: major-lab/mirbooking


text:
#ifndef OPENMC_TALLIES_FILTER_H
#define OPENMC_TALLIES_FILTER_H

#include <cstdint>
#include <string>
#include <unordered_map>

#include "pugixml.hpp"
#include <gsl/gsl>

#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/memory.h"
#include "openmc/particle.h"
#include "openmc/tallies/filter_match.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Modifies tally score events.
//==============================================================================

class Filter
{
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors, factory functions

  Filter();
  virtual ~Filter();

  //! Create a new tally filter
  //
  //! \tparam T Type of the filter
  //! \param[in] id  Unique ID for the filter. If none is passed, an ID is
  //!    automatically assigned
  //! \return Pointer to the new filter object
  template<typename T>
  static T* create(int32_t id = -1);

  //! Create a new tally filter
  //
  //! \param[in] type  Type of the filter
  //! \param[in] id  Unique ID for the filter. If none is passed, an ID is
  //!    automatically assigned
  //! \return Pointer to the new filter object
  static Filter* create(const std::string& type, int32_t id = -1);

  //! Create a new tally filter from an XML node
  //
  //! \param[in] node XML node
  //! \return Pointer to the new filter object
  static Filter* create(pugi::xml_node node);

  //! Uses an XML input to fill the filter's data fields.
  virtual void from_xml(pugi::xml_node node) = 0;

  //----------------------------------------------------------------------------
  // Methods

  virtual std::string type() const = 0;

  //! Matches a tally event to a set of filter bins and weights.
  //!
  //! \param[in] p Particle being tracked
  //! \param[in] estimator Tally estimator being used
  //! \param[out] match will contain the matching bins and corresponding
  //!   weights; note that there may be zero matching bins
  virtual void
  get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const = 0;

  //! Writes data describing this filter to an HDF5 statepoint group.
  virtual void
  to_statepoint(hid_t filter_group) const
  {
    write_dataset(filter_group, "type", type());
    write_dataset(filter_group, "n_bins", n_bins_);
  }

  //! Return a string describing a filter bin for the tallies.out file.
  //
  //! For example, an `EnergyFilter` might return the string
  //! "Incoming Energy [0.625E-6, 20.0)".
  virtual std::string text_label(int bin) const = 0;

  //----------------------------------------------------------------------------
  // Accessors

  //! Get unique ID of filter
  //! \return Unique ID
  int32_t id() const { return id_; }

  //! Assign a unique ID to the filter
  //! \param[in]  Unique ID to assign. A value of -1 indicates that an ID should
  //!   be automatically assigned
  void set_id(int32_t id);

  //! Get number of bins
  //! \return Number of bins
  int n_bins() const { return n_bins_; }

  gsl::index index() const { return index_; }

  //----------------------------------------------------------------------------
  // Data members

protected:
  int n_bins_;
private:
  int32_t id_ {C_NONE};
  gsl::index index_;
};

//==============================================================================
// Global variables
//==============================================================================

namespace model {
  extern "C" int32_t n_filters;
  extern std::unordered_map<int, int> filter_map;
  extern vector<unique_ptr<Filter>> tally_filters;
}

//==============================================================================
// Non-member functions
//==============================================================================

//! Make sure index corresponds to a valid filter
int verify_filter(int32_t index);

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_H

max_stars_repo_name: cjwyett/openmc

max_issues_repo_name: cjwyett/openmc

max_forks_repo_name: cjwyett/openmc