text:
/* linalg/test_cod.c
 *
 * Copyright (C) 2017 Patrick Alken
 *
 * 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

static int test_COD_decomp_eps(const gsl_matrix * m, const double eps, const char *desc);
static int test_COD_decomp(gsl_rng *r);
static int test_COD_lssolve_eps(const gsl_matrix * m, const double * actual, const double eps, const char *desc);
static int test_COD_lssolve(void);
static int test_COD_lssolve2_eps(const double lambda, const gsl_matrix * A, const gsl_vector * b, const double eps, const char *desc);
static int test_COD_lssolve2(gsl_rng * r);

/* create a matrix of a given rank */
static int
create_rank_matrix(const size_t rank, gsl_matrix * m, gsl_rng * r)
{
  const size_t M = m->size1;
  const size_t N = m->size2;
  size_t i;
  gsl_vector *u = gsl_vector_alloc(M);
  gsl_vector *v = gsl_vector_alloc(N);

  gsl_matrix_set_zero(m);

  /* add several rank-1 matrices together */
  for (i = 0; i < rank; ++i)
    {
      create_random_vector(u, r);
      create_random_vector(v, r);
      gsl_blas_dger(1.0, u, v, m);
    }

  gsl_vector_free(u);
  gsl_vector_free(v);

  return GSL_SUCCESS;
}

static int
test_COD_decomp_eps(const gsl_matrix * m, const double eps, const char *desc)
{
  int s = 0;
  size_t i, j, M = m->size1, N = m->size2;
  size_t rank;

  gsl_matrix * QRZT = gsl_matrix_alloc(M, N);
  gsl_matrix * Q = gsl_matrix_alloc(M, M);
  gsl_matrix * R = gsl_matrix_alloc(M, N);
  gsl_matrix * QR = gsl_matrix_alloc(M, N);
  gsl_matrix * Z = gsl_matrix_alloc(N, N);
  gsl_vector * tau_Q = gsl_vector_alloc(GSL_MIN(M, N));
  gsl_vector * tau_Z = gsl_vector_alloc(GSL_MIN(M, N));
  gsl_vector * work = gsl_vector_alloc(N);
  gsl_matrix * lhs = gsl_matrix_alloc(M, N);
  gsl_matrix * rhs = gsl_matrix_alloc(M, N);

  gsl_permutation * perm = gsl_permutation_alloc(N);

  gsl_matrix_memcpy(QRZT, m);

  s += gsl_linalg_COD_decomp(QRZT, tau_Q, tau_Z, perm, &rank, work);
  s += gsl_linalg_COD_unpack(QRZT, tau_Q, tau_Z, rank, Q, R, Z);

  /* compute lhs = m P */
  gsl_matrix_memcpy(lhs, m);
  gsl_permute_matrix(perm, lhs);

  /* compute rhs = Q R Z^T */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, Q, R, 0.0, QR);
  gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, QR, Z, 0.0, rhs);

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          double aij = gsl_matrix_get(rhs, i, j);
          double bij = gsl_matrix_get(lhs, i, j);

          gsl_test_rel(aij, bij, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i, j, aij, bij);
        }
    }

  gsl_permutation_free (perm);
  gsl_vector_free(work);
  gsl_vector_free(tau_Q);
  gsl_vector_free(tau_Z);
  gsl_matrix_free(QRZT);
  gsl_matrix_free(lhs);
  gsl_matrix_free(rhs);
  gsl_matrix_free(QR);
  gsl_matrix_free(Q);
  gsl_matrix_free(R);
  gsl_matrix_free(Z);

  return s;
}

static int
test_COD_decomp(gsl_rng *r)
{
  int s = 0;
  size_t N;
  gsl_matrix *m;

  /* test COD decomposition on Hilbert matrices */
  for (N = 1; N <= 12; ++N)
    {
      m = gsl_matrix_alloc(N, N);
      create_hilbert_matrix2(m);

      test_COD_decomp_eps(m, 256.0 * N * GSL_DBL_EPSILON, "COD_decomp hilbert");

      gsl_matrix_free(m);
    }

  /* build some matrices of a given rank and test */

  m = gsl_matrix_alloc(100, 50);
  create_rank_matrix(26, m, r);
  test_COD_decomp_eps(m, 1.0e2 * GSL_DBL_EPSILON, "COD_decomp rank 26");
  gsl_matrix_free(m);

  m = gsl_matrix_alloc(550, 200);
  create_rank_matrix(176, m, r);
  test_COD_decomp_eps(m, 1.0e3 * GSL_DBL_EPSILON, "COD_decomp rank 176");
  gsl_matrix_free(m);

  test_COD_decomp_eps(m35, 1.0e1 * GSL_DBL_EPSILON, "COD_decomp m(3,5)");
  test_COD_decomp_eps(m53, 1.0e1 * GSL_DBL_EPSILON, "COD_decomp m(5,3)");

  test_COD_decomp_eps(s35, 1.0e1 * GSL_DBL_EPSILON, "COD_decomp s(3,5)");
  test_COD_decomp_eps(s53, 1.0e1 * GSL_DBL_EPSILON, "COD_decomp s(5,3)");

  test_COD_decomp_eps(vander2, 1.0e1 * GSL_DBL_EPSILON, "COD_decomp vander(2)");
  test_COD_decomp_eps(vander3, 1.0e1 * GSL_DBL_EPSILON, "COD_decomp vander(3)");
  test_COD_decomp_eps(vander4, 1.0e1 * GSL_DBL_EPSILON, "COD_decomp vander(4)");
  test_COD_decomp_eps(vander12, 1e-3, "COD_decomp vander(12)"); /* FIXME: large tolerance needed */

  return s;
}

static int
test_COD_lssolve_eps(const gsl_matrix * m, const double * actual, const double eps, const char *desc)
{
  int s = 0;
  size_t i, M = m->size1, N = m->size2;

  gsl_vector * lhs = gsl_vector_alloc(M);
  gsl_vector * rhs = gsl_vector_alloc(M);
  gsl_matrix * QRZT  = gsl_matrix_alloc(M, N);
  gsl_vector * tau_Q = gsl_vector_alloc(GSL_MIN(M, N));
  gsl_vector * tau_Z = gsl_vector_alloc(GSL_MIN(M, N));
  gsl_vector * work = gsl_vector_alloc(N);
  gsl_vector * x = gsl_vector_alloc(N);
  gsl_vector * r = gsl_vector_alloc(M);
  gsl_vector * res = gsl_vector_alloc(M);
  gsl_permutation * perm = gsl_permutation_alloc(N);
  size_t rank;

  gsl_matrix_memcpy(QRZT, m);

  for (i = 0; i < M; i++)
    gsl_vector_set(rhs, i, i + 1.0);

  s += gsl_linalg_COD_decomp(QRZT, tau_Q, tau_Z, perm, &rank, work);
  s += gsl_linalg_COD_lssolve(QRZT, tau_Q, tau_Z, perm, rank, rhs, x, res);

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      gsl_test_rel(xi, actual[i], eps,
                   "%s (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                   desc, M, N, i, xi, actual[i]);
    }

  /* compute residual r = b - m x */
  if (M == N)
    {
      gsl_vector_set_zero(r);
    }
  else
    {
      gsl_vector_memcpy(r, rhs);
      gsl_blas_dgemv(CblasNoTrans, -1.0, m, x, 1.0, r);
    }

  for (i = 0; i < N; i++)
    {
      double r1 = gsl_vector_get(res, i);
      double r2 = gsl_vector_get(r, i);

      if (fabs(r2) < 1.0e3 * GSL_DBL_EPSILON)
        {
          gsl_test_abs(r1, r2, 10.0 * eps,
                       "%s res (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i, r1, r2);
        }
      else
        {
          gsl_test_rel(r1, r2, eps,
                       "%s res (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i, r1, r2);
        }
    }

  gsl_vector_free(r);
  gsl_vector_free(res);
  gsl_vector_free(x);
  gsl_vector_free(tau_Q);
  gsl_vector_free(tau_Z);
  gsl_matrix_free(QRZT);
  gsl_vector_free(rhs);
  gsl_vector_free(lhs);
  gsl_vector_free(work);
  gsl_permutation_free(perm);

  return s;
}

static int
test_COD_lssolve(void)
{
  int s = 0;

  test_COD_lssolve_eps(m53, m53_lssolution, 1.0e4 * GSL_DBL_EPSILON, "COD_lssolve m(5,3)");

  test_COD_lssolve_eps(hilb2, hilb2_solution, 1.0e2 * GSL_DBL_EPSILON, "COD_lssolve hilbert(2)");
  test_COD_lssolve_eps(hilb3, hilb3_solution, 1.0e2 * GSL_DBL_EPSILON, "COD_lssolve hilbert(3)");
  test_COD_lssolve_eps(hilb4, hilb4_solution, 1.0e4 * GSL_DBL_EPSILON, "COD_lssolve hilbert(4)");

  test_COD_lssolve_eps(vander2, vander2_solution, 1.0e1 * GSL_DBL_EPSILON, "COD_lssolve vander(2)");
  test_COD_lssolve_eps(vander3, vander3_solution, 1.0e1 * GSL_DBL_EPSILON, "COD_lssolve vander(3)");
  test_COD_lssolve_eps(vander4, vander4_solution, 1.0e2 * GSL_DBL_EPSILON, "COD_lssolve vander(4)");

  /* rank-1 least squares problem from the 'lin2' test dataset for multifit_nlinear */
  {
    const size_t M = 20;
    const size_t N = 5;

    /* unique minimum norm solution from "Factorize" matlab package */
    const double x_ls[] = { 1.818181818181817e-02, 3.636363636363636e-02, 5.454545454545454e-02,
                            7.272727272727272e-02, 9.090909090909088e-02 };

    gsl_matrix *m = gsl_matrix_alloc(M, N);
    size_t i, j;

    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            gsl_matrix_set(m, i, j, (i + 1.0) * (j + 1.0));
          }
      }

    test_COD_lssolve_eps(m, x_ls, 1.0e2 * GSL_DBL_EPSILON, "COD_lssolve lin2");

    gsl_matrix_free(m);
  }

  return s;
}

/* solve: min ||b - A x||^2 + lambda^2 ||x||^2 */
static int
test_COD_lssolve2_eps(const double lambda, const gsl_matrix * A, const gsl_vector * b, const double eps, const char *desc)
{
  int s = 0;
  size_t i, M = A->size1, N = A->size2;

  gsl_vector * lhs = gsl_vector_alloc(M);
  gsl_matrix * QRZT  = gsl_matrix_alloc(M, N);
  gsl_vector * tau_Q = gsl_vector_alloc(GSL_MIN(M, N));
  gsl_vector * tau_Z = gsl_vector_alloc(GSL_MIN(M, N));
  gsl_vector * work = gsl_vector_alloc(N);
  gsl_vector * x = gsl_vector_alloc(N);
  gsl_vector * x_aug = gsl_vector_alloc(N);
  gsl_vector * r = gsl_vector_alloc(M);
  gsl_vector * res = gsl_vector_alloc(M);
  gsl_permutation * perm = gsl_permutation_alloc(N);
  size_t rank;

  /* form full rank augmented system B = [ A ; lambda*I_N ], f = [ rhs ; 0 ] and solve with QRPT */
  {
    gsl_vector_view v;
    gsl_matrix_view m;
    gsl_permutation *p = gsl_permutation_alloc(N);
    gsl_matrix * B = gsl_matrix_calloc(M + N, N);
    gsl_vector * f = gsl_vector_calloc(M + N);
    gsl_vector * tau = gsl_vector_alloc(N);
    gsl_vector * residual = gsl_vector_alloc(M + N);
    int signum;

    m = gsl_matrix_submatrix(B, 0, 0, M, N);
    gsl_matrix_memcpy(&m.matrix, A);

    m = gsl_matrix_submatrix(B, M, 0, N, N);
    v = gsl_matrix_diagonal(&m.matrix);
    gsl_vector_set_all(&v.vector, lambda);

    v = gsl_vector_subvector(f, 0, M);
    gsl_vector_memcpy(&v.vector, b);

    /* solve: [ A ; lambda*I ] x_aug = [ b ; 0 ] */
    gsl_linalg_QRPT_decomp(B, tau, p, &signum, work);
    gsl_linalg_QRPT_lssolve(B, tau, p, f, x_aug, residual);

    gsl_permutation_free(p);
    gsl_matrix_free(B);
    gsl_vector_free(f);
    gsl_vector_free(tau);
    gsl_vector_free(residual);
  }

  gsl_matrix_memcpy(QRZT, A);

  s += gsl_linalg_COD_decomp(QRZT, tau_Q, tau_Z, perm, &rank, work);

  {
    gsl_matrix *S = gsl_matrix_alloc(rank, rank);
    gsl_vector *workr = gsl_vector_alloc(rank);

    s += gsl_linalg_COD_lssolve2(lambda, QRZT, tau_Q, tau_Z, perm, rank, b, x, res, S, workr);

    gsl_matrix_free(S);
    gsl_vector_free(workr);
  }

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(x_aug, i);
      gsl_test_rel(xi, yi, eps,
                   "%s (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                   desc, M, N, i, xi, yi);
    }

  /* compute residual r = b - A x */
  if (M == N)
    {
      gsl_vector_set_zero(r);
    }
  else
    {
      gsl_vector_memcpy(r, b);
      gsl_blas_dgemv(CblasNoTrans, -1.0, A, x, 1.0, r);
    }

  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(res, i);
      double yi = gsl_vector_get(r, i);

      gsl_test_rel(xi, yi, sqrt(eps),
                   "%s res (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                   desc, M, N, i, xi, yi);
    }

  gsl_vector_free(r);
  gsl_vector_free(res);
  gsl_vector_free(x);
  gsl_vector_free(x_aug);
  gsl_vector_free(tau_Q);
  gsl_vector_free(tau_Z);
  gsl_matrix_free(QRZT);
  gsl_vector_free(lhs);
  gsl_vector_free(work);
  gsl_permutation_free(perm);

  return s;
}

static int
test_COD_lssolve2(gsl_rng * r)
{
  int s = 0;
  gsl_matrix *m;
  gsl_vector *b;
  size_t M, N;
  const double lambda = 2.3;

  M = 100;
  N = 50;
  m = gsl_matrix_alloc(M, N);
  b = gsl_vector_alloc(M);
  create_rank_matrix(26, m, r);
  create_random_vector(b, r);
  test_COD_lssolve2_eps(lambda, m, b, 1.0e4 * M * GSL_DBL_EPSILON, "COD_lssolve2 rank 26");
  gsl_matrix_free(m);
  gsl_vector_free(b);

  M = 500;
  N = 450;
  m = gsl_matrix_alloc(M, N);
  b = gsl_vector_alloc(M);
  create_rank_matrix(278, m, r);
  create_random_vector(b, r);
  test_COD_lssolve2_eps(lambda, m, b, 1.0e6 * M * GSL_DBL_EPSILON, "COD_lssolve2 rank 278");
  gsl_matrix_free(m);
  gsl_vector_free(b);

  return s;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#pragma once

#include <string>

#include <gsl/gsl_complex.h>

#include "Expression.h"

struct NumericalExpression: public Expression {

    NumericalExpression(Scanner::Type kind);

    expression derivative(const std::string& var) override;
    expression integrate(const std::string& var) override;
};

class NumExpression: public NumericalExpression {
    double real, imag;

    NumExpression(double real, double imag);
    NumExpression(const std::string&);

    public:
        static expression construct(double real, double imag = 0);
        static expression construct(const gsl_complex& z);
        static expression construct(const std::string&);

        bool isNumber() const override { return true; }

        gsl_complex complex(const Variables& vars = emptyVars) const override;

        EXPRESSION_OVERRIDES

};

class HexExpression: public NumericalExpression {
    unsigned long long num;

    HexExpression(unsigned long long num);
    HexExpression(const std::string&);

    public:
        static expression construct(unsigned long long num);
        static expression construct(const std::string&);

        bool isNumber() const override { return true; }

        EXPRESSION_OVERRIDES

};

class BinExpression: public NumericalExpression {
    unsigned long long num;

    BinExpression(unsigned long long num);
    BinExpression(const std::string&);

    public:
        static expression construct(unsigned long long num);
        static expression construct(const std::string&);

        bool isNumber() const override { return true; }

        EXPRESSION_OVERRIDES

};

ext: h

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
#include "AsynchronousUpdate.h"
#include <stdio.h>
#include <stdlib.h>
#include "constants.h"
#include <gsl/gsl_rng.h>
#include <omp.h>

int getBool(gsl_rng* r)
{
	return gsl_rng_uniform_int(r,2);
}

void singleStateUpdate(int n, int* state, int** topology, int seed, int* fixed_nodes)
{
	//Perform Asynchronous Updates
	gsl_rng* asyncer = gsl_rng_alloc(gsl_rng_ranlxs2);
	gsl_rng_set(asyncer,seed);

	for (int i = 0; i < ITER_ASYNCHRONOUS; ++i)
	{
		int node_to_update = gsl_rng_uniform_int(asyncer,n);
		if(fixed_nodes[node_to_update] == NORMAL_FLAG)
		{
			int val_temp = 0;

			for (int j = 0; j < n; ++j)
			{
				val_temp+= state[j]*topology[node_to_update][j];
			}

			if(val_temp > 0)
			{
				state[node_to_update] = 1;
			}
			else if(val_temp<0)
			{
				state[node_to_update] = LOWER_EXPRESSION;
			}
		}
	}
	//Free the rng
	gsl_rng_free(asyncer);
}

void getFinalState(int n, int state[n], gsl_rng* rangen, int** topology, int* fixed_nodes, int i)
{
	#pragma omp parallel for
	for (int j = 0; j < n; ++j)
	{
		if(fixed_nodes[j]==NORMAL_FLAG)
		{
			#pragma omp critical(rand)
			{
				state[j] = getBool(rangen);
			}
			if(state[j]==0)
			{
				state[j] = LOWER_EXPRESSION;
			}
		}
		else{
			state[j] = fixed_nodes[j];
		}
	}

	singleStateUpdate(n,state,topology,i+1,fixed_nodes);
}

asynclist* updateAsyncList(int n, int** topology, int* fixed_nodes, int height, int seed_init)
{
	if(height)
	{
		if(height%2)
		{
			gsl_rng* rangen = gsl_rng_alloc(gsl_rng_ranlxs2);
			gsl_rng_set(rangen,seed_init);
			asynclist* to_merge[2];
			int seed[2] = {gsl_rng_get(rangen),gsl_rng_get(rangen)};
			int i;
			#pragma omp parallel for shared(n,topology,fixed_nodes,seed,height) private(i)
			for (i = 0; i < 2; ++i)
			{
				to_merge[i] = updateAsyncList(n,topology,fixed_nodes,height-1,seed[i]);
			}
			return merge_asynclist(to_merge[0],to_merge[1],n);
		}
		else
		{
			gsl_rng* rangen = gsl_rng_alloc(gsl_rng_taus2);
			gsl_rng_set(rangen,seed_init);
			asynclist* to_merge[2];
			int seed[2] = {gsl_rng_get(rangen),gsl_rng_get(rangen)};
			int i;
			#pragma omp parallel for shared(n,topology,fixed_nodes,seed,height) private(i)
			for (i = 0; i < 2; ++i)
			{
				to_merge[i] = updateAsyncList(n,topology,fixed_nodes,height-1,seed[i]);
			}
			return merge_asynclist(to_merge[0],to_merge[1],n);
		}
	}
	else{
		asynclist* list = (asynclist*)malloc(sizeof(asynclist));
		init_asynclist(list);
		int state[n];

		gsl_rng* rangen = gsl_rng_alloc(gsl_rng_ranlxs2);
		gsl_rng_set(rangen,seed_init);

		for (int i = 0; i < n; ++i)
		{
			if(fixed_nodes[i] == NORMAL_FLAG)
			{
				state[i] = gsl_rng_uniform_int(rangen,2);
				if(state[i]==0)
				{
					state[i] = LOWER_EXPRESSION;
				}
			}
			else
			{
				state[i] = fixed_nodes[i];
			}
		}

		singleStateUpdate(n,state,topology,gsl_rng_get(rangen),fixed_nodes);
		list->n_elements = 1;
		asynclistnode* node = (asynclistnode*)malloc(sizeof(asynclistnode));
		if(!node)
		{
			exit(10);
		}
		node->n_occurances = 1;
		node->state = (int*)malloc(n*sizeof(int));
		if(!(node->state))
		{
			exit(11);
		}
		for (int i = 0; i < n; ++i)
		{
			node->state[i] = state[i];
		}
		node->next = NULL;
		node->prev = NULL;
		list->next = node;
		return list;
	}
}

void updateAsyncTree(int n, int** topology, int* fixed_nodes, base* stable)
{
	int iter = (1<<N_SAMPLES);

	gsl_rng* rangen = gsl_rng_alloc(gsl_rng_ranlxs2);
	gsl_rng_set(rangen,0);

	#pragma omp parallel for
	for (int i = 0; i < iter; ++i)
	{
		/*
		gsl_rng* rangen = gsl_rng_alloc(gsl_rng_ranlxs2);
		gsl_rng_set(rangen,i);
		*/
		//int* state=(int*)malloc(n*sizeof(int));
		int state[n];
		getFinalState(n,state,rangen,topology,fixed_nodes,i);

		#pragma omp critical(addToTree)
		{
			add_node(stable,n,state);
		}

		//free(state);
	}
	gsl_rng_free(rangen);
}

ext: c

max_stars_repo_name: SamarthH/CBool

max_issues_repo_name: SamarthH/CBool

max_forks_repo_name: SamarthH/CBool


text:
/* matrix/gsl_matrix_uchar.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * 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 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_UCHAR_H__
#define __GSL_MATRIX_UCHAR_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_uchar.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned char * data;
  gsl_block_uchar * block;
  int owner;
} gsl_matrix_uchar;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_view;

typedef _gsl_matrix_uchar_view gsl_matrix_uchar_view;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_const_view;

typedef const _gsl_matrix_uchar_const_view gsl_matrix_uchar_const_view;

/* Allocation */

gsl_matrix_uchar *
gsl_matrix_uchar_alloc (const size_t n1, const size_t n2);

gsl_matrix_uchar *
gsl_matrix_uchar_calloc (const size_t n1, const size_t n2);

gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_block (gsl_block_uchar * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_matrix (gsl_matrix_uchar * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

gsl_vector_uchar *
gsl_vector_uchar_alloc_row_from_matrix (gsl_matrix_uchar * m,
                                        const size_t i);

gsl_vector_uchar *
gsl_vector_uchar_alloc_col_from_matrix (gsl_matrix_uchar * m,
                                        const size_t j);

void gsl_matrix_uchar_free (gsl_matrix_uchar * m);

/* Views */

_gsl_matrix_uchar_view
gsl_matrix_uchar_submatrix (gsl_matrix_uchar * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_uchar_view
gsl_matrix_uchar_row (gsl_matrix_uchar * m, const size_t i);

_gsl_vector_uchar_view
gsl_matrix_uchar_column (gsl_matrix_uchar * m, const size_t j);

_gsl_vector_uchar_view
gsl_matrix_uchar_diagonal (gsl_matrix_uchar * m);

_gsl_vector_uchar_view
gsl_matrix_uchar_subdiagonal (gsl_matrix_uchar * m, const size_t k);

_gsl_vector_uchar_view
gsl_matrix_uchar_superdiagonal (gsl_matrix_uchar * m, const size_t k);

_gsl_matrix_uchar_view
gsl_matrix_uchar_view_array (unsigned char * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_uchar_view
gsl_matrix_uchar_view_array_with_tda (unsigned char * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector (gsl_vector_uchar * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector_with_tda (gsl_vector_uchar * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_submatrix (const gsl_matrix_uchar * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_row (const gsl_matrix_uchar * m,
                            const size_t i);

_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_column (const gsl_matrix_uchar * m,
                               const size_t j);

_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_diagonal (const gsl_matrix_uchar * m);

_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subdiagonal (const gsl_matrix_uchar * m,
                                    const size_t k);

_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_superdiagonal (const gsl_matrix_uchar * m,
                                      const size_t k);

_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array (const unsigned char * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array_with_tda (const unsigned char * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector (const gsl_vector_uchar * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector_with_tda (const gsl_vector_uchar * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

unsigned char   gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j);
void    gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x);

unsigned char * gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j);
const unsigned char * gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j);

void gsl_matrix_uchar_set_zero (gsl_matrix_uchar * m);
void gsl_matrix_uchar_set_identity (gsl_matrix_uchar * m);
void gsl_matrix_uchar_set_all (gsl_matrix_uchar * m, unsigned char x);

int gsl_matrix_uchar_fread (FILE * stream, gsl_matrix_uchar * m) ;
int gsl_matrix_uchar_fwrite (FILE * stream, const gsl_matrix_uchar * m) ;
int gsl_matrix_uchar_fscanf (FILE * stream, gsl_matrix_uchar * m);
int gsl_matrix_uchar_fprintf (FILE * stream, const gsl_matrix_uchar * m, const char * format);

int gsl_matrix_uchar_memcpy(gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);
int gsl_matrix_uchar_swap(gsl_matrix_uchar * m1, gsl_matrix_uchar * m2);

int gsl_matrix_uchar_swap_rows(gsl_matrix_uchar * m, const size_t i, const size_t j);
int gsl_matrix_uchar_swap_columns(gsl_matrix_uchar * m, const size_t i, const size_t j);
int gsl_matrix_uchar_swap_rowcol(gsl_matrix_uchar * m, const size_t i, const size_t j);
int gsl_matrix_uchar_transpose (gsl_matrix_uchar * m);
int gsl_matrix_uchar_transpose_memcpy (gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);

unsigned char gsl_matrix_uchar_max (const gsl_matrix_uchar * m);
unsigned char gsl_matrix_uchar_min (const gsl_matrix_uchar * m);
void gsl_matrix_uchar_minmax (const gsl_matrix_uchar * m, unsigned char * min_out, unsigned char * max_out);

void gsl_matrix_uchar_max_index (const gsl_matrix_uchar * m, size_t * imax, size_t *jmax);
void gsl_matrix_uchar_min_index (const gsl_matrix_uchar * m, size_t * imin, size_t *jmin);
void gsl_matrix_uchar_minmax_index (const gsl_matrix_uchar * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

int gsl_matrix_uchar_isnull (const gsl_matrix_uchar * m);
int gsl_matrix_uchar_ispos (const gsl_matrix_uchar * m);
int gsl_matrix_uchar_isneg (const gsl_matrix_uchar * m);

int gsl_matrix_uchar_add (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
int gsl_matrix_uchar_sub (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
int gsl_matrix_uchar_mul_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
int gsl_matrix_uchar_div_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
int gsl_matrix_uchar_scale (gsl_matrix_uchar * a, const double x);
int gsl_matrix_uchar_add_constant (gsl_matrix_uchar * a, const double x);
int gsl_matrix_uchar_add_diagonal (gsl_matrix_uchar * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_uchar_get_row(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t i);
int gsl_matrix_uchar_get_col(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t j);
int gsl_matrix_uchar_set_row(gsl_matrix_uchar * m, const size_t i, const gsl_vector_uchar * v);
int gsl_matrix_uchar_set_col(gsl_matrix_uchar * m, const size_t j, const gsl_vector_uchar * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
unsigned char
gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
unsigned char *
gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (unsigned char *) (m->data + (i * m->tda + j)) ;
}

extern inline
const unsigned char *
gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const unsigned char *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_UCHAR_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* spmatrix/gsl_spmatrix_float.h
 *
 * Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2018 Patrick Alken
 *
 * 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_SPMATRIX_COMPLEX_FLOAT_H__
#define __GSL_SPMATRIX_COMPLEX_FLOAT_H__

#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_bst.h>
#include <gsl/gsl_vector_complex_float.h>
#include <gsl/gsl_matrix_complex_float.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/*
 * COO format:
 *
 * If data[n] = A_{ij}, then:
 *   i = A->i[n]
 *   j = A->p[n]
 *
 * Compressed column format (CSC):
 *
 * If data[n] = A_{ij}, then:
 *   i = A->i[n]
 *   A->p[j] <= n < A->p[j+1]
 * so that column j is stored in
 * [ data[p[j]], data[p[j] + 1], ..., data[p[j+1] - 1] ]
 *
 * Compressed row format (CSR):
 *
 * If data[n] = A_{ij}, then:
 *   j = A->i[n]
 *   A->p[i] <= n < A->p[i+1]
 * so that row i is stored in
 * [ data[p[i]], data[p[i] + 1], ..., data[p[i+1] - 1] ]
 */

typedef struct
{
  size_t size1;              /* number of rows */
  size_t size2;              /* number of columns */

  /* i (size nzmax) contains:
   *
   * COO/CSC: row indices
   * CSR: column indices
   */
  int *i;

  float *data;               /* matrix elements of size nzmax */

  /*
   * COO: p[n] = column number of element data[n]
   * CSC: p[j] = index in data of first non-zero element in column j
   * CSR: p[i] = index in data of first non-zero element in row i
   */
  int *p;

  size_t nzmax;              /* maximum number of matrix elements */
  size_t nz;                 /* number of non-zero values in matrix */

  gsl_bst_workspace *tree;   /* binary tree structure */
  gsl_spmatrix_pool *pool;   /* memory pool for binary tree nodes */
  size_t node_size;          /* size of individual tree node in bytes */

  /*
   * workspace of size 2*MAX(size1,size2)*MAX(sizeof(float),sizeof(int))
   * used in various routines
   */
  union
    {
      void *work_void;
      int *work_int;
      float *work_atomic;
    } work;

  int sptype;                /* sparse storage type */
  size_t spflags;            /* GSL_SPMATRIX_FLG_xxx */
} gsl_spmatrix_complex_float;

/*
 * Prototypes
 */

/* allocation / initialization */

gsl_spmatrix_complex_float * gsl_spmatrix_complex_float_alloc (const size_t n1, const size_t n2);
gsl_spmatrix_complex_float * gsl_spmatrix_complex_float_alloc_nzmax (const size_t n1, const size_t n2,
                                                                     const size_t nzmax, const int sptype);
void gsl_spmatrix_complex_float_free (gsl_spmatrix_complex_float * m);
int gsl_spmatrix_complex_float_realloc (const size_t nzmax, gsl_spmatrix_complex_float * m);
size_t gsl_spmatrix_complex_float_nnz (const gsl_spmatrix_complex_float * m);
const char * gsl_spmatrix_complex_float_type (const gsl_spmatrix_complex_float * m);
int gsl_spmatrix_complex_float_set_zero (gsl_spmatrix_complex_float * m);
int gsl_spmatrix_complex_float_tree_rebuild (gsl_spmatrix_complex_float * m);

/* compress */

int gsl_spmatrix_complex_float_csc (gsl_spmatrix_complex_float * dest, const gsl_spmatrix_complex_float * src);
int gsl_spmatrix_complex_float_csr (gsl_spmatrix_complex_float * dest, const gsl_spmatrix_complex_float * src);
gsl_spmatrix_complex_float * gsl_spmatrix_complex_float_compress (const gsl_spmatrix_complex_float * src, const int sptype);
gsl_spmatrix_complex_float * gsl_spmatrix_complex_float_compcol (const gsl_spmatrix_complex_float * src);
gsl_spmatrix_complex_float * gsl_spmatrix_complex_float_ccs (const gsl_spmatrix_complex_float * src);
gsl_spmatrix_complex_float * gsl_spmatrix_complex_float_crs (const gsl_spmatrix_complex_float * src);

/* copy */

int gsl_spmatrix_complex_float_memcpy (gsl_spmatrix_complex_float * dest, const gsl_spmatrix_complex_float * src);

/* file I/O */

int gsl_spmatrix_complex_float_fprintf (FILE * stream, const gsl_spmatrix_complex_float * m, const char * format);
gsl_spmatrix_complex_float * gsl_spmatrix_complex_float_fscanf (FILE * stream);
int gsl_spmatrix_complex_float_fwrite (FILE * stream, const gsl_spmatrix_complex_float * m);
int gsl_spmatrix_complex_float_fread (FILE * stream, gsl_spmatrix_complex_float * m);

/* get/set */

gsl_complex_float gsl_spmatrix_complex_float_get (const gsl_spmatrix_complex_float * m, const size_t i, const size_t j);
int gsl_spmatrix_complex_float_set (gsl_spmatrix_complex_float * m, const size_t i, const size_t j, const gsl_complex_float x);
gsl_complex_float * gsl_spmatrix_complex_float_ptr (const gsl_spmatrix_complex_float * m, const size_t i, const size_t j);

/* operations */

int gsl_spmatrix_complex_float_scale (gsl_spmatrix_complex_float * m, const gsl_complex_float x);
int gsl_spmatrix_complex_float_scale_columns (gsl_spmatrix_complex_float * m, const gsl_vector_complex_float * x);
int gsl_spmatrix_complex_float_scale_rows (gsl_spmatrix_complex_float * m, const gsl_vector_complex_float * x);
int gsl_spmatrix_complex_float_add (gsl_spmatrix_complex_float * c, const gsl_spmatrix_complex_float * a, const gsl_spmatrix_complex_float * b);
int gsl_spmatrix_complex_float_d2sp (gsl_spmatrix_complex_float * T, const gsl_matrix_complex_float * A);
int gsl_spmatrix_complex_float_sp2d (gsl_matrix_complex_float * A, const gsl_spmatrix_complex_float * S);

/* properties */

int gsl_spmatrix_complex_float_equal (const gsl_spmatrix_complex_float * a, const gsl_spmatrix_complex_float * b);

/* swap */

int gsl_spmatrix_complex_float_transpose (gsl_spmatrix_complex_float * m);
int gsl_spmatrix_complex_float_transpose2 (gsl_spmatrix_complex_float * m);
int gsl_spmatrix_complex_float_transpose_memcpy (gsl_spmatrix_complex_float * dest, const gsl_spmatrix_complex_float * src);

__END_DECLS

#endif /* __GSL_SPMATRIX_COMPLEX_FLOAT_H__ */

max_stars_repo_name: vinej/sml

max_issues_repo_name: vinej/sml

max_forks_repo_name: vinej/sml


text:
//
// Created by ZL on 2019-06-13.
//

#ifndef ZQCNNDEMO_ZQCNN_MTCNN_NCHWC4_H
#define ZQCNNDEMO_ZQCNN_MTCNN_NCHWC4_H

#include "ZQ_CNN_Net_NCHWC.h"
#include "ZQ_CNN_MTCNN_NCHWC.h"
#include <vector>
#include <iostream>
#include "opencv2/opencv.hpp"
#include "ZQ_CNN_CompileConfig.h"
#include <cblas.h>

class MTCNNNCHWC {
public:
    MTCNNNCHWC(const std::string &model_path);

    std::vector<ZQ::ZQ_CNN_BBox> detect(const std::string &img_path);

    std::vector<ZQ::ZQ_CNN_BBox> detectMat(cv::Mat &image0);

private:
    ZQ::ZQ_CNN_MTCNN_NCHWC *mtcnn;
    int thread_num = 1;
};

#endif //ZQCNNDEMO_ZQCNN_MTCNN_NCHWC4_H


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <map>
#include "DrawableGameComponent.h"
#include "FullScreenRenderTarget.h"
#include "FullScreenQuad.h"
#include "MatrixHelper.h"

namespace Library
{
	class Texture2D;
	class TexturedModelMaterial;
}

namespace Rendering
{
	class DistortionMaskingDemo final : public Library::DrawableGameComponent
	{
	public:
		DistortionMaskingDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		DistortionMaskingDemo(const DistortionMaskingDemo&) = delete;
		DistortionMaskingDemo(DistortionMaskingDemo&&) = default;
		DistortionMaskingDemo& operator=(const DistortionMaskingDemo&) = default;
		DistortionMaskingDemo& operator=(DistortionMaskingDemo&&) = default;
		~DistortionMaskingDemo();

		bool DrawCutoutModeEnabled() const;
		void SetDrawCutoutModeEnabled(bool enabled);
		void ToggledDrawCutoutModeEnabled();

		float DisplacementScale() const;
		void SetDisplacementScale(float displacementScale);

		virtual void Initialize() override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		struct PixelCBufferPerObject
		{
			float DisplacementScale{ 1.0f };
			DirectX::XMFLOAT3 Padding;
		};

		enum class DistortionShaderClass
		{
			Cutout = 0,
			Composite,
			NoDistortion
		};

		inline static const DirectX::XMVECTORF32 CutoutZeroColor{ DirectX::Colors::Black };

		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		std::shared_ptr<Library::TexturedModelMaterial> mTexturedModelMaterial;
		Library::FullScreenRenderTarget mSceneRenderTarget;
		Library::FullScreenRenderTarget mCutoutRenderTarget;
		Library::FullScreenQuad mFullScreenQuad;
		winrt::com_ptr<ID3D11Buffer> mPixelCBufferPerObject;
		PixelCBufferPerObject mPixelCBufferPerObjectData;
		std::shared_ptr<Library::Texture2D> mDistortionMap;
		std::map<DistortionShaderClass, winrt::com_ptr<ID3D11ClassInstance>> mShaderClassInstances;
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		bool mDrawCutoutModeEnabled{ false };
		bool mUpdateMaterial{ true };
	};
}


text:
#pragma once
/*
 * (C) Copyright 2020-2021 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/*! \defgroup ioda_cxx_types Type System
 * \brief The data type system
 * \ingroup ioda_cxx_api
 *
 * @{
 * \file Type.h
 * \brief Interfaces for ioda::Type and related classes. Implements the type system.
 */
#include <array>
#include <cstring>
#include <functional>
#include <gsl/gsl-lite.hpp>
#include <memory>
#include <string>
#include <typeindex>
#include <typeinfo>
#include <vector>

#include "ioda/Types/Type_Provider.h"
#include "ioda/Exception.h"
#include "ioda/defs.h"

namespace ioda {
class Type;

/// Basic pre-defined types (Python convenience wrappers)
/// \see py_ioda.cpp
/// \note Names here do not match the python equivalents. The
///   Python names match numpy's definitions.
enum class BasicTypes {
  undefined_,  ///< Internal use only
  float_,
  double_,
  ldouble_,
  char_,
  short_,
  ushort_,
  int_,
  uint_,
  lint_,
  ulint_,
  llint_,
  ullint_,
  int32_,
  uint32_,
  int16_,
  uint16_,
  int64_,
  uint64_,
  bool_,
  str_
};

namespace detail {
IODA_DL size_t COMPAT_strncpy_s(char* dest, size_t destSz, const char* src, size_t srcSz);

class Type_Backend;

template <class Type_Implementation = Type>
class Type_Base {
  friend class ::ioda::Type;
  std::shared_ptr<Type_Backend> backend_;

protected:
  ::ioda::detail::Type_Provider* provider_;

  /// @name General Functions
  /// @{

  Type_Base(std::shared_ptr<Type_Backend> b, ::ioda::detail::Type_Provider* p)
      : backend_(b), provider_(p) {}

  /// Get the type provider.
  inline detail::Type_Provider* getTypeProvider() const { return provider_; }

  /*
  /// \brief Convenience function to check a type.
  /// \param DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \returns True if the type matches
  /// \returns False (0) if the type does not match
  /// \throws if an error occurred.
  template <class DataType>
  bool isA() const {
    Type templateType = Types::GetType_Wrapper<DataType>::GetType(getTypeProvider());

    return isA(templateType);
  }
  /// Hand-off to the backend to check equivalence
  virtual bool isA(Type lhs) const;

  /// Python compatability function
  inline bool isA(BasicTypes dataType) const { return isA(Type(dataType, getTypeProvider())); }
  */
public:
  virtual ~Type_Base() {}
  std::shared_ptr<Type_Backend> getBackend() const { return backend_; }
  bool isValid() const { return (backend_.use_count() > 0); }

  /// \brief Get the size of a single element of a type, in bytes.
  /// \details This function is paired with the read and write functions to allow you to
  /// read and write data in a type-agnostic manner.
  /// This size report is a bit complicated when variable-length strings are encountered.
  /// In these cases, the size of the string pointer is returned.
  virtual size_t getSize() const;

  /// @}
};
}  // namespace detail

/// \brief Represents the "type" (i.e. integer, string, float) of a piece of data.
/// \ingroup ioda_cxx_types
///
/// Generally, you do not have to use this class directly. Attributes and Variables have
/// templated functions that convert your type into the type used internally by ioda.
/// \see Types::GetType and Types::GetType_Wrapper for functions that produce these types.
class IODA_DL Type : public detail::Type_Base<> {
public:
  Type();
  Type(std::shared_ptr<detail::Type_Backend> b, std::type_index t);
  Type(BasicTypes, gsl::not_null<::ioda::detail::Type_Provider*> t);

  virtual ~Type();

  /// @name Type-querying functions
  /// @{

  /// @deprecated This function is problematic since we cannot query a type properly
  /// when loading from a file.
  std::type_index getType() const { return as_type_index_; }
  inline std::type_index operator()() const { return getType(); }
  inline std::type_index get() const { return getType(); }

  /// @}
private:
  std::type_index as_type_index_;
};

namespace detail {

/// Backends inherit from this and they provide their own functions.
/// Lots of std::dynamic_cast, unfortunately.
class IODA_DL Type_Backend : public Type_Base<> {
public:
  virtual ~Type_Backend();
  size_t getSize() const override;

protected:
  Type_Backend();
};

}  // namespace detail

/// \brief Defines the type system used for manipulating IODA objects.
namespace Types {

// using namespace ioda::Handles;

/// \brief Convenience struct to determine if a type can represent a string.
/// \ingroup ioda_cxx_types
/// \todo extend to UTF-8 strings, as HDF5 supports these. No support for UTF-16, but conversion
/// functions may be applied.
/// \todo Fix for "const std::string".
template <typename T>
struct is_string : public std::integral_constant<
                     bool, std::is_same<char*, typename std::decay<T>::type>::value
                             || std::is_same<const char*, typename std::decay<T>::type>::value> {};
/// \brief Convenience struct to determine if a type can represent a string.
/// \ingroup ioda_cxx_types
template <>
struct is_string<std::string> : std::true_type {};

/// Useful compile-time definitions.
namespace constants {
/// \note Different than ObsSpace variable-length dimension. This is for a Type.
constexpr size_t _Variable_Length = 0;
// constexpr int _Not_An_Array_type = -3;
}  // namespace constants

/// \brief For fundamental, non-string types.
/// \ingroup ioda_cxx_types
template <class DataType, int Array_Type_Dimensionality = 0>
Type GetType(gsl::not_null<const ::ioda::detail::Type_Provider*> t,
             std::initializer_list<Dimensions_t> Adims                   = {},
             typename std::enable_if<!is_string<DataType>::value>::type* = 0) {
  if (Array_Type_Dimensionality <= 0)
    throw Exception(
      "Bad assertion / unsupported fundamental type at the frontend side "
      "of the ioda type system.", ioda_Here());
  else
    return t->makeArrayType(Adims, typeid(DataType[]), typeid(DataType));
}
/// \brief For fundamental string types. These are either constant or variable length arrays.
/// Separate handling elsewhere.
/// \ingroup ioda_cxx_types
template <class DataType, int String_Type_Length = constants::_Variable_Length>
Type GetType(gsl::not_null<const ::ioda::detail::Type_Provider*> t,
             std::initializer_list<Dimensions_t>                        = {},
             typename std::enable_if<is_string<DataType>::value>::type* = 0) {
  return t->makeStringType(String_Type_Length, typeid(DataType));
}

// This macro just repeats a long definition

/// @def IODA_ADD_FUNDAMENTAL_TYPE
/// Macro that defines a "fundamental type" that needs to be supported
/// by the backend. These match C++11.
/// \ingroup ioda_cxx_types
/// \see https://en.cppreference.com/w/cpp/language/types
/// \since C++11: we use bool, short int, unsigned short int,
///   int, unsigned int, long int, unsigned long int,
///   long long int, unsigned long long int,
///   signed char, unsigned char, char,
///   wchar_t, char16_t, char32_t,
///   float, double, long double.
/// \since C++20: we also add char8_t.
#define IODA_ADD_FUNDAMENTAL_TYPE(x)                                                               \
  template <>                                                                                \
  inline Type GetType<x, 0>(gsl::not_null<const ::ioda::detail::Type_Provider*> t,                 \
                            std::initializer_list<Dimensions_t>, void*) {                          \
    return t->makeFundamentalType(typeid(x));                                                      \
  }

IODA_ADD_FUNDAMENTAL_TYPE(bool);
IODA_ADD_FUNDAMENTAL_TYPE(short int);
IODA_ADD_FUNDAMENTAL_TYPE(unsigned short int);
IODA_ADD_FUNDAMENTAL_TYPE(int);
IODA_ADD_FUNDAMENTAL_TYPE(unsigned int);
IODA_ADD_FUNDAMENTAL_TYPE(long int);
IODA_ADD_FUNDAMENTAL_TYPE(unsigned long int);
IODA_ADD_FUNDAMENTAL_TYPE(long long int);
IODA_ADD_FUNDAMENTAL_TYPE(unsigned long long int);
IODA_ADD_FUNDAMENTAL_TYPE(signed char);
IODA_ADD_FUNDAMENTAL_TYPE(unsigned char);
IODA_ADD_FUNDAMENTAL_TYPE(char);
IODA_ADD_FUNDAMENTAL_TYPE(wchar_t);
IODA_ADD_FUNDAMENTAL_TYPE(char16_t);
IODA_ADD_FUNDAMENTAL_TYPE(char32_t);
// IODA_ADD_FUNDAMENTAL_TYPE(char8_t); // C++20
IODA_ADD_FUNDAMENTAL_TYPE(float);
IODA_ADD_FUNDAMENTAL_TYPE(double);
IODA_ADD_FUNDAMENTAL_TYPE(long double);

#undef IODA_ADD_FUNDAMENTAL_TYPE

/*
/// Used in an example. Incomplete.
/// \todo Pop off std::array as a 1-D object
template<> inline Type GetType<std::array<int,2>, 0>
        (gsl::not_null<const ::ioda::detail::Type_Provider*> t,
        std::initializer_list<Dimensions_t>, void*) {
        return t->makeArrayType({2}, typeid(std::array<int,2>), typeid(int)); }

*/

/*
template <class DataType, int Array_Type_Dimensionality = 0>
Type GetType(
        gsl::not_null<const ::ioda::detail::Type_Provider*> t,
        std::initializer_list<Dimensions_t> Adims = {},
        typename std::enable_if<!is_string<DataType>::value>::type* = 0);
template <class DataType, int String_Type_Length = constants::_Variable_Length>
Type GetType(
        gsl::not_null<const ::ioda::detail::Type_Provider*> t,
        typename std::enable_if<is_string<DataType>::value>::type* = 0);
        */

/// \brief Wrapper struct to call GetType. Needed because of C++ template rules.
/// \ingroup ioda_cxx_types
/// \see ioda::Attribute, ioda::Has_Attributes, ioda::Variable, ioda::Has_Variables
template <class DataType,
          int Length = 0>  //, typename = std::enable_if_t<!is_string<DataType>::value>>
struct GetType_Wrapper {
  static Type GetType(gsl::not_null<const ::ioda::detail::Type_Provider*> t) {
    /// \note Currently breaks array types, but these are not yet used.
    return ::ioda::Types::GetType<DataType, Length>(t, {Length});
  }
};
/// \ingroup ioda_cxx_types
typedef std::function<Type(gsl::not_null<const ::ioda::detail::Type_Provider*>)>
  TypeWrapper_function;
/*
template <class DataType, int Length = 0, typename = std::enable_if_t<is_string<DataType>::value>>
struct GetType_Wrapper {
        Type GetType(gsl::not_null<const ::ioda::detail::Type_Provider*> t) const {
                // string split
                return ::ioda::Types::GetType<DataType, Length>(t);
        }
};
*/

// inline Encapsulated_Handle GetTypeFixedString(Dimensions_t sz);
}  // namespace Types
}  // namespace ioda

/// @}

max_stars_repo_name: Gibies/ioda

max_issues_repo_name: Gibies/ioda

max_forks_repo_name: Gibies/ioda


text:
#ifndef THREADEDNESS_INCLUDE
#define THREADEDNESS_INCLUDE

#include <string>

#include <gsl/string_span>

#include "config/commandLineOptions.h"
#include "config/fleetingOptionsInterface.h"
#include "config/variablesMap.h"
#include "core/task.h"

namespace execHelper::plugins {
using Jobs = config::Jobs_t;

const gsl::czstring<> JOBS_KEY = "jobs";

/**
 * \brief Extends the functionality to include the _jobs_ config parameter and processes this parameter, using the --jobs flag
 */
struct JobsLong {
    /**
     * Adds the variables for this functionality to the given variables map
     *
     * @param[in] options       The fleeting options to take into account
     * @param[out] variables    The variables map to add the variables to
     */
    static void
    getVariables(config::VariablesMap& variables,
                 const config::FleetingOptionsInterface& options) noexcept;

    /**
     * Applies the given variables to the task
     *
     * @param[in] variables The variables map to use
     * @param[out] task The task with the given variables map applied to it
     */
    inline static void apply(core::Task& task,
                             const config::VariablesMap& variables) noexcept {
        task.append(
            {"--jobs", std::to_string(*(variables.get<Jobs>(JOBS_KEY)))});
    }
};

/**
 * \brief Extends the functionality to include the _jobs_ config parameter and processes this parameter, using the -j flag
 */
struct JobsShort {
    /*! @copydoc JobsLong::getVariables(config::VariablesMap&, const config::FleetingOptionsInterface&)
     */
    inline static void
    getVariables(config::VariablesMap& variables,
                 const config::FleetingOptionsInterface& options) noexcept {
        JobsLong::getVariables(variables, options);
    }

    /*! @copydoc JobsLong::apply(core::Task&, const config::VariablesMap&)
     */
    inline static void apply(core::Task& task,
                             const config::VariablesMap& variables) noexcept {
        task.append({"-j", std::to_string(*(variables.get<Jobs>(JOBS_KEY)))});
    }
};

} // namespace execHelper::plugins

#endif /* THREADEDNESS_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
#pragma once

#include <gsl.h>
#include <boost/optional.hpp>

#include <mitkPythonService.h>

#include <QObject>
#include <QSet>
#include <QString>
#include <QVariant>
#include <QVariantList>
#include <QVector>

#include <SolidData.h>
#include <FaceIdentifier.h>
#include <ISolverStudyData.h>
#include <MeshData.h>
#include <VesselForestData.h>
#include <PCMRIData.h>

#include "PythonSolverSetupServiceActivator.h"
#include <PythonSolverSetupManager.h>
#include "PythonSolverSetupServiceExports.h"

namespace crimson
{

class   FaceTypeQtWrapper : public QObject
{
    Q_OBJECT
public:
    Q_ENUMS(FaceType);
    enum FaceType {
        ftCapInflow = FaceIdentifier::ftCapInflow,
        ftCapOutflow = FaceIdentifier::ftCapOutflow,
        ftWall = FaceIdentifier::ftWall,

        ftUndefined = FaceIdentifier::ftUndefined
    };
};

class   ArrayDataTypeQtWrapper : public QObject
{
    Q_OBJECT
public:
    Q_ENUMS(ArrayDataType);
    enum class ArrayDataType {
        Int = 0,
        Double = 1,
    };
};

class   PyFaceIdentifier
{
public:
    static PythonQtObjectPtr cppToPy(const FaceIdentifier& cppFaceId)
    {
        auto pythonService = PythonSolverSetupServiceActivator::getPythonService();

        pythonService->Execute("import CRIMSONCore");

        QVariantList parentSolidIndices;
        for (const auto& parentSolidIndex : cppFaceId.parentSolidIndices) {
            parentSolidIndices << QString::fromStdString(parentSolidIndex);
        }

        auto result = PythonQtObjectPtr{pythonService->GetPythonManager()->mainContext().call(
            "CRIMSONCore.FaceIdentifier.FaceIdentifier",
            QVariantList{static_cast<FaceTypeQtWrapper::FaceType>(cppFaceId.faceType), parentSolidIndices})};
        Ensures(!result.isNull());
        return result;
    }

    static boost::optional<FaceIdentifier> pyToCpp(PythonQtObjectPtr pyFaceId)
    {
        FaceIdentifier cppFaceId;

        QVariant faceIdVariant = pyFaceId.getVariable("faceType");
        if (!faceIdVariant.canConvert<int>()) {
            MITK_ERROR << "Face identifier not recognized - wrong faceType (expected int)";
            return boost::none;
        }
        cppFaceId.faceType = static_cast<FaceIdentifier::FaceType>(faceIdVariant.toInt());

        QVariant parentSolidIndicesVariant = pyFaceId.getVariable("parentSolidIndices");
        if (!parentSolidIndicesVariant.canConvert<QVariantList>()) {
            MITK_ERROR << "Face identifier not recognized - wrong parentSolidIndices (expected list of strings)";
            return boost::none;
        }

        for (const auto& parentSolidIndexVariant : parentSolidIndicesVariant.toList()) {
            if (!parentSolidIndexVariant.canConvert<QString>()) {
                MITK_ERROR << "Face identifier not recognized - wrong item in parentSolidIndices (expected string)";
                return boost::none;
            }

            cppFaceId.parentSolidIndices.insert(parentSolidIndexVariant.toString().toStdString());
        }

        return cppFaceId;
    }
};

class   Utils : public QObject
{
    Q_OBJECT
public slots:
    void static_Utils_reloadAll() { crimson::PythonSolverSetupServiceActivator::reloadPythonSolverSetups(true); }

    void static_Utils_logError(QString message) { MITK_ERROR << message.toStdString(); }

    void static_Utils_logWarning(QString message) { MITK_WARN << message.toStdString(); }

    void static_Utils_logInformation(QString message) { MITK_INFO << message.toStdString(); }
};

class   SolidDataQtWrapper : public QObject
{
    Q_OBJECT
public:
	SolidDataQtWrapper()
        : _brep(nullptr)
    {
    }
	SolidDataQtWrapper(gsl::not_null<const SolidData*> brep)
        : _brep(brep)
    {
    }

	SolidDataQtWrapper(const SolidDataQtWrapper& other)
        : _brep(other._brep)
    {
    }

public slots:

    QVector<double> getFaceNormal(PythonQtObjectPtr pyFaceIdentifier) const
    {
        Expects(_brep != nullptr);
        auto cppFaceIdOptional = PyFaceIdentifier::pyToCpp(pyFaceIdentifier);

        if (!cppFaceIdOptional) {
            return {};
        }

        auto normal = _brep->getFaceNormal(cppFaceIdOptional.get());
        return QVector<double>{normal[0], normal[1], normal[2]};
    }

    double getDistanceToFaceEdge(PythonQtObjectPtr pyFaceIdentifier, double x, double y, double z) const
    {
        Expects(_brep != nullptr);
        auto cppFaceIdOptional = PyFaceIdentifier::pyToCpp(pyFaceIdentifier);

        if (!cppFaceIdOptional) {
            return 0;
        }

        auto pos = mitk::Point3D{};
        pos[0] = x;
        pos[1] = y;
        pos[2] = z;
        return _brep->getDistanceToFaceEdge(cppFaceIdOptional.get(), pos);
    }

    PythonQtObjectPtr getFaceIdentifierForModelFace(int faceIndex) const
    {
        Expects(_brep != nullptr);
        auto identifierOpt = _brep->getFaceIdentifierMap().getFaceIdentifierForModelFace(faceIndex);
        return identifierOpt ? PyFaceIdentifier::cppToPy(identifierOpt.get()) : nullptr;
    }

    int getNumberOfModelFaces() const
    {
        Expects(_brep != nullptr);
        return _brep->getFaceIdentifierMap().getNumberOfModelFaces();
    }

    int getNumberOfFaceIdentifiers() const
    {
        Expects(_brep != nullptr);
        return _brep->getFaceIdentifierMap().getNumberOfFaceIdentifiers();
    }

    int faceIdentifierIndex(PythonQtObjectPtr pyFaceIdentifier) const
    {
        Expects(_brep != nullptr);

        auto cppFaceIdOptional = PyFaceIdentifier::pyToCpp(pyFaceIdentifier);

        if (!cppFaceIdOptional) {
            return 0;
        }

        return _brep->getFaceIdentifierMap().faceIdentifierIndex(cppFaceIdOptional.get());
    }

    PythonQtObjectPtr getFaceIdentifier(int uniqueFaceId) const
    {
        Expects(_brep != nullptr);
        return PyFaceIdentifier::cppToPy(_brep->getFaceIdentifierMap().getFaceIdentifier(uniqueFaceId));
    }

    QVector<int> getModelFacesForFaceIdentifier(PythonQtObjectPtr pyFaceIdentifier) const
    {
        Expects(_brep != nullptr);
        auto cppFaceIdOptional = PyFaceIdentifier::pyToCpp(pyFaceIdentifier);

        if (!cppFaceIdOptional) {
            return {};
        }

        return QVector<int>::fromStdVector(_brep->getFaceIdentifierMap().getModelFacesForFaceIdentifier(cppFaceIdOptional.get()));
    }

private:
	const SolidData* _brep;
};

class   MeshDataQtWrapper : public QObject
{
    Q_OBJECT
public:
    MeshDataQtWrapper()
        : _mesh(nullptr)
    {
    }
    MeshDataQtWrapper(gsl::not_null<const MeshData*> mesh)
        : _mesh(mesh)
    {
    }

    MeshDataQtWrapper(const MeshDataQtWrapper& other)
        : _mesh(other._mesh)
    {
    }

public slots:

    int getNNodes() const
    {
        Expects(_mesh != nullptr);
        return _mesh->getNNodes();
    }

    int getNEdges() const
    {
        Expects(_mesh != nullptr);
        return _mesh->getNEdges();
    }

    int getNFaces() const
    {
        Expects(_mesh != nullptr);
        return _mesh->getNFaces();
    }

    int getNElements() const
    {
        Expects(_mesh != nullptr);
        return _mesh->getNElements();
    }

    QVector<double> getNodeCoordinates(int nodeIndex) const
    {
        Expects(_mesh != nullptr);
        auto coords = _mesh->getNodeCoordinates(nodeIndex);
        return QVector<double>{coords[0], coords[1], coords[2]};
    }

    QVector<int> getElementNodeIds(int elementIndex) const
    {
        Expects(_mesh != nullptr);
        return QVector<int>::fromStdVector(_mesh->getElementNodeIds(elementIndex));
    }

    QVector<int> getAdjacentElements(int elementIndex) const
    {
        Expects(_mesh != nullptr);
        return QVector<int>::fromStdVector(_mesh->getAdjacentElements(elementIndex));
    }

    QVector<int> getNodeIdsForFace(PythonQtObjectPtr pyFaceIdentifier) const
    {
        Expects(_mesh != nullptr);
        auto cppFaceIdOptional = PyFaceIdentifier::pyToCpp(pyFaceIdentifier);

        if (!cppFaceIdOptional) {
            return {};
        }

        return QVector<int>::fromStdVector(_mesh->getNodeIdsForFace(cppFaceIdOptional.get()));
    }

    QVariantList getMeshFaceInfoForFace(PythonQtObjectPtr pyFaceIdentifier) const // TODO: QVector<QVector<int>>
    {
        Expects(_mesh != nullptr);
        auto cppFaceIdOptional = PyFaceIdentifier::pyToCpp(pyFaceIdentifier);

        if (!cppFaceIdOptional) {
            return {};
        }

        auto faceInfos = _mesh->getMeshFaceInfoForFace(cppFaceIdOptional.get());

        auto out = QVariantList{};
        out.reserve(faceInfos.size());
        std::transform(
            faceInfos.begin(), faceInfos.end(), std::back_inserter(out),
            [](const crimson::MeshData::MeshFaceInfo& info) {
                return QVariant::fromValue(QVector<int>{info.elementId, info.globalFaceId, info.nodeIds[0], info.nodeIds[1], info.nodeIds[2]});
            });
        return out;
    }

private:
    const MeshData* _mesh;
};

class   VesselForestDataQtWrapper : public QObject
{
    Q_OBJECT
public:
    using UIDToVesselPathDataMap = std::unordered_map<std::string, const VesselPathAbstractData*>;

    VesselForestDataQtWrapper()
        : _vesselForest(nullptr)
    {
    }
    VesselForestDataQtWrapper(gsl::not_null<const VesselForestData*> data, UIDToVesselPathDataMap uidToVesselPathDataMap)
        : _vesselForest(data)
        , _uidToVesselPathDataMap(std::move(uidToVesselPathDataMap))
    {
    }

    VesselForestDataQtWrapper(const VesselForestDataQtWrapper& other)
        : _vesselForest(other._vesselForest)
        , _uidToVesselPathDataMap(other._uidToVesselPathDataMap)
    {
    }
//    VesselForestDataQtWrapper(const VesselForestDataQtWrapper&) = default;

    const UIDToVesselPathDataMap& getUIDToVesselPathDataMap() const { return _uidToVesselPathDataMap; }

public slots:
    QVariantMap getActiveConnectedComponentsMap() const
    {
        Expects(_vesselForest != nullptr);

        auto result = QVariantMap{};
        for (const auto& uidComponentIndexPair : _vesselForest->computeActiveConnectedComponents()) {
            result[QString::fromStdString(uidComponentIndexPair.first)] = uidComponentIndexPair.second;
        }

        return result;
    }

    QVariantList getClosestPoint(PythonQtObjectPtr pyFaceIdentifier, double x, double y, double z) const
    {
        Expects(_vesselForest != nullptr);

        auto closestPathInfo = _getClosestVesselPath(pyFaceIdentifier, x, y, z);
        if (!closestPathInfo.vesselPath) {
            return {};
        }

        return QVariantList{closestPathInfo.distance, closestPathInfo.t};
    }

    QVector<double> getVesselPathCoordinateFrame(PythonQtObjectPtr pyFaceIdentifier, double x, double y, double z) const
    {
        auto closestPathInfo = _getClosestVesselPath(pyFaceIdentifier, x, y, z);
        if (!closestPathInfo.vesselPath) {
            return{};
        }

        auto pt = closestPathInfo.vesselPath->getPosition(closestPathInfo.t);
        auto tangent = closestPathInfo.vesselPath->getTangentVector(closestPathInfo.t);
        auto normal = closestPathInfo.vesselPath->getNormalVector(closestPathInfo.t);

        return QVector<double>{pt[0], pt[1], pt[2], tangent[0], tangent[1], tangent[2], normal[0], normal[1], normal[2]};
    }

private:
    struct ClosestVesselPathInfo {
        const VesselPathAbstractData* vesselPath = nullptr;
        double distance = 0;
        double t = 0;
    };

    ClosestVesselPathInfo _getClosestVesselPath(PythonQtObjectPtr pyFaceIdentifier, double x, double y, double z) const
    {
        auto cppFaceIdOptional = PyFaceIdentifier::pyToCpp(pyFaceIdentifier);

        if (!cppFaceIdOptional) {
            return {};
        }

        auto&& parentSolidIndices = cppFaceIdOptional->parentSolidIndices;

        auto result = ClosestVesselPathInfo{};

        auto pos = mitk::Point3D{};
        mitk::FillVector3D(pos, x, y, z);
        for (const auto& vesselUID : parentSolidIndices) {
            auto closestPointRequestResult = _uidToVesselPathDataMap.at(vesselUID)->getClosestPoint(pos);

            auto d = closestPointRequestResult.closestPoint.EuclideanDistanceTo(pos);
            if (result.vesselPath == nullptr || d < result.distance) {
                result.vesselPath = _uidToVesselPathDataMap.at(vesselUID);
                result.t = closestPointRequestResult.t;
                result.distance = d;
            }
        }

        return result;
    }

    const VesselForestData* _vesselForest;
    UIDToVesselPathDataMap _uidToVesselPathDataMap;
};

class   PCMRIDataQtWrapper : public QObject
{
	Q_OBJECT
public:
	PCMRIDataQtWrapper()
		: _pcmriData(nullptr)
	{
	}
	PCMRIDataQtWrapper(gsl::not_null<const PCMRIData*> pcmriData)
		: _pcmriData(pcmriData)
	{
	}

	PCMRIDataQtWrapper(const PCMRIDataQtWrapper& other)
		: _pcmriData(other._pcmriData)
	{
	}

public slots:

	QVector<double> getTimepoints() const
	{
		Expects(_pcmriData != nullptr);
		return QVector<double>::fromStdVector(_pcmriData->getTimepoints());
	}

	QVector<double> getFlowWaveform() const
	{
		Expects(_pcmriData != nullptr);
		return QVector<double>::fromStdVector(_pcmriData->getFlowWaveform());
	}

	QVector<double> getSingleMappedPCMRIvector(int pointIndex, int timepointIndex) const
	{
		Expects(_pcmriData != nullptr);
		auto coords = _pcmriData->getSingleMappedPCMRIvector(pointIndex, timepointIndex);
		return QVector<double>{coords[0], coords[1], coords[2]};
	}

private:
	const PCMRIData* _pcmriData;
};
}

Q_DECLARE_METATYPE(crimson::SolidDataQtWrapper)
Q_DECLARE_METATYPE(crimson::MeshDataQtWrapper)
Q_DECLARE_METATYPE(crimson::VesselForestDataQtWrapper)
Q_DECLARE_METATYPE(crimson::PCMRIDataQtWrapper)

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


text:
/******************************************************************************
 *                                                                            *
 * ELECTRONS.C                                                                *
 *                                                                            *
 * ELECTRON THERMODYNAMICS                                                    *
 *                                                                            *
 ******************************************************************************/

#include "decs.h"
#include <gsl/gsl_sf_bessel.h>

/* TODO: Encapsulate electron equation of state in EOS framework.
 *       For now, I'v ebroken encapsulation but made the code
 *	 complain/break if EOS_TYPE is not EOS_TYPE_GAMMA
 *	 when ELECTRONS are active.
 *
 *       The best strategy is probably to make a two-temperature
 *       electron EOS one of the available EOS's.
 *
 *       Until then, I'm worried that we might hit a maintainability
 *       problem becasue if the implementation in ELECTRONS changes
 *       the implementation in EOS_TYPE_GAMMA needs to change too.
 *       ~JMM
 */

#if ELECTRONS
#if EOS == EOS_TYPE_GAMMA
void heat_electrons_zone(int i, int j, int k, double Pi[NVAR], double Ps[NVAR],
    double Pf[NVAR], double Dt);
void fixup_electrons_1zone(double P[NVAR]);

void init_electrons() {
  double uel;

  ZSLOOP(-NG, N1 + NG - 1, -NG, NG + N2 - 1, -NG, NG + N3 - 1) {
    // Set electron internal energy to constant fraction of internal energy
    uel = fel0 * P[i][j][k][UU];

    // Initialize entropies
    P[i][j][k][KTOT] = (gam - 1.) * P[i][j][k][UU] * pow(P[i][j][k][RHO], -gam);
    P[i][j][k][KEL]  = (game - 1.) * uel * pow(P[i][j][k][RHO], -game);
  }

  bound_prim(P);
}

void heat_electrons(
    grid_prim_type Pi, grid_prim_type Ps, grid_prim_type Pf, double Dt) {
  timer_start(TIMER_ELECTRON);
#pragma omp parallel for collapse(3) schedule(dynamic)
  ZLOOP {
    heat_electrons_zone(i, j, k, Pi[i][j][k], Ps[i][j][k], Pf[i][j][k], Dt);
  }
  timer_stop(TIMER_ELECTRON);
}

void heat_electrons_zone(int i, int j, int k, double Pi[NVAR], double Ps[NVAR],
    double Pf[NVAR], double Dt) {
  double ktotharm, ktotadv, fel;

  // Calculated and advected entropy at final time
  ktotharm = (gam - 1.) * Pf[UU] / pow(Pf[RHO], gam);
  ktotadv  = Pf[KTOT];

  // Electron heating fraction
  fel = get_fel(i, j, k, Ps);

  // Update final electron entropy according to Ressler+ 2015 Eqn. 27:
  Pf[KEL] += (game - 1.) / (gam - 1.) * pow(Ps[RHO], gam - game) * fel *
             (ktotharm - ktotadv);

  // Diagnostics
  struct of_geom *geom = &ggeom[i][j][CENT];
  struct of_state q;
  get_state(Ps, geom, &q);
  double uadv = ktotadv / (gam - 1.) * pow(Pf[RHO], gam);
  double Qud  = q.ucon[0] * q.ucov[0] * (Pf[UU] - uadv) *
               pow(Ps[RHO] / Pf[RHO], gam) / Dt;
  Qvisc[i][j][k] = fel * Qud;

  // Reset total entropy
  Pf[KTOT] = ktotharm;
}

double get_fel(int i, int j, int k, double P[NVAR]) {
#if BETA_HEAT == 0
  return fel0;
#endif

  struct of_geom geom;
  double         beta, fel, c1, Trat, Tpr, mbeta, qrat;
  double         pres, bsq;
  double         c2, c3, c22, c32;

  c1 = 0.92;

  Tpr        = (gam - 1.) * P[UU] / P[RHO];
  double uel = 1. / (game - 1.) * P[KEL] * pow(P[RHO], game);
  double Tel = (game - 1.) * uel / P[RHO];
  if (Tel <= 0.)
    Tel = SMALL;
  if (Tpr <= 0.)
    Tpr = SMALL;
  Trat = fabs(Tpr / Tel);

  pres = P[RHO] * Tpr; // Proton pressure

  geom = *get_geometry(i, j, k, CENT);
  bsq  = bsq_calc(P, &geom);

  beta = pres / bsq * 2.;
  if (beta > 1.e20)
    beta = 1.e20;
  mbeta = 2. - 0.2 * log10(Trat);

  if (Trat <= 1.) {
    c2 = 1.6 / Trat;
    c3 = 18. + 5. * log10(Trat);
  } else {
    c2 = 1.2 / Trat;
    c3 = 18.;
  }
  c22 = pow(c2, 2.);
  c32 = pow(c3, 2.);

  qrat = c1 * (c22 + pow(beta, mbeta)) / (c32 + pow(beta, mbeta)) *
         exp(-1. / beta) * pow(MP / ME * Trat, .5);
  fel = 1. / (1. + qrat);

  return fel;
}

// Modified Bessel function of second kind with safe inputs
double safe_Kn(int n, double x) {
  if (x > 100.) {
    return exp(-x) * sqrt(M_PI / (2. * x));
  } else {
    return gsl_sf_bessel_Kn(n, x);
  }
}

#if RADIATION
void coulomb(
    grid_prim_type Pi, grid_prim_type Ps, grid_prim_type Pf, double Dt) {
  timer_start(TIMER_ELECTRON);

#pragma omp parallel for collapse(3)
  ZLOOP {
    double          rho     = Ps[i][j][k][RHO];
    double          thetae  = MP / ME * Ps[i][j][k][KEL] * pow(rho, game - 1.);
    double          ue      = Ps[i][j][k][KEL] * pow(rho, game) / (game - 1.);
    double          up      = Ps[i][j][k][UU] - ue;
    double          n       = rho * Ne_unit;
    double          Ti      = up * U_unit * (gamp - 1.) / (n * KBOL);
    double          thetai  = KBOL * Ti / (MP * CL * CL);
    double          thetam  = 1. / (1. / thetae + 1. / thetai);
    double          logCoul = COULOMB_LOG;
    double          Te      = thetae * ME * CL * CL / KBOL;
    struct of_geom *geom;
    struct of_state q;

    // Sanity checks, although electron fixup routine should catch these
    if (!isnan(Te) && !isnan(Ti) && Te > 0. && Ti > 0.) {
      double Qc, term1, term2;

      // Get Coulomb heating rate.
      // Need to handle cases where Thetai < 1e-2, Thetae < 1e-2, and both
      // Thetae and Thetai < 1e-2 separately due to Bessel functions exploding
      double prefac =
          3. / 2. * ME / MP * n * n * logCoul * CL * KBOL * THOMSON * (Ti - Te);
      double thetaCrit = 1.e-2;
      if (thetae < thetaCrit && thetai < thetaCrit) {
        term1 = sqrt(thetam / (M_PI * thetae * thetai / 2.));
        term2 = sqrt(thetam / (M_PI * thetae * thetai / 2.));
      } else if (thetae < thetaCrit) {
        term1 =
            exp(-1. / thetai) / safe_Kn(2, 1. / thetai) * sqrt(thetam / thetae);
        term2 =
            exp(-1. / thetai) / safe_Kn(2, 1. / thetai) * sqrt(thetam / thetae);
      } else if (thetai < thetaCrit) {
        term1 =
            exp(-1. / thetae) / safe_Kn(2, 1. / thetae) * sqrt(thetam / thetai);
        term2 =
            exp(-1. / thetae) / safe_Kn(2, 1. / thetae) * sqrt(thetam / thetai);
      } else {
        term1 = safe_Kn(1, 1. / thetam) /
                (safe_Kn(2, 1. / thetae) * safe_Kn(2, 1. / thetai));
        term2 = safe_Kn(0, 1. / thetam) /
                (safe_Kn(2, 1. / thetae) * safe_Kn(2, 1. / thetai));
      }
      term1 *= (2. * pow(thetae + thetai, 2) + 1.) / (thetae + thetai);
      term2 *= 2.;
      Qc = prefac * (term1 + term2);

      // Convert to code units
      Qc *= T_unit / U_unit;

      // Update electron internal energy
      geom = &ggeom[i][j][CENT];
      get_state(Ps[i][j][k], geom, &q);

      double ue_f =
          Pf[i][j][k][KEL] * pow(Pf[i][j][k][RHO], game) / (game - 1.);
      ue_f += Qc * Dt / q.ucon[0];

      // Record diagnostic
      Qcoul[i][j][k] = q.ucov[0] * Qc;

      // Update electron entropy
      Pf[i][j][k][KEL] = (game - 1.) * ue_f * pow(Pf[i][j][k][RHO], -game);
    }
  } // ZLOOP

  timer_stop(TIMER_ELECTRON);
}

void apply_rad_force_e(grid_prim_type Prh, grid_prim_type Pr,
    grid_fourvector_type radG, double Dt) {
// Apply only to active zones for this proc -- ghost zone four-force
// depositions already communicated over MPI
#pragma omp parallel for collapse(3) schedule(dynamic)
  ZLOOP {
    struct of_geom *geom = &ggeom[i][j][CENT];
    struct of_state q;
    double          Uel, Urho;
    double          C = 0.;

    // Get fluid state at n + 1/2 where radiation four-force is centered
    get_state(Prh[i][j][k], geom, &q);

    for (int mu = 0; mu < NDIM; mu++) {
      C += -q.ucon[mu] * radG[i][j][k][mu];
    }

    // Get fluid state at n+1 for full update
    get_state(Pr[i][j][k], geom, &q);

    // Remove \sqrt{-g} from radG
    C = C / geom->g;

    Urho = Pr[i][j][k][RHO] * q.ucon[0];
    Uel  = Pr[i][j][k][KEL] * Urho;

    Uel += Dt * (C * (game - 1.) * pow(Prh[i][j][k][RHO], 1. - game));

    // Supercooling diagnostics
    if (Uel < 0.) {
      double          U_1[NVAR], prim_2[NVAR], U_2[NVAR];
      struct of_state q_1, q_2;
      Nsuper[i][j][k]++;

      // (1) Record total energy density after cooling
      get_state(Pr[i][j][k], &ggeom[i][j][CENT], &q_1);
      primtoflux(Pr[i][j][k], &q_1, 0, 0, &ggeom[i][j][CENT], U_1);

      // (2) Calculate total energy density with zero electron energy density
      PLOOP  prim_2[ip] = psupersave[i][j][k][ip];
      double ue         = prim_2[KEL] * pow(prim_2[RHO], game) / (game - 1.);
      prim_2[UU] -= ue;
      get_state(prim_2, &ggeom[i][j][CENT], &q_2);
      primtoflux(prim_2, &q_2, 0, 0, &ggeom[i][j][CENT], U_2);

      // Subtract (2) from (1); integrated over volume, this is fake energy
      Esuper[i][j][k] += fabs((U_1[UU] - U_2[UU]) * dx[1] * dx[2] * dx[3]);
    } // Uel < 0

    Pr[i][j][k][KEL] = Uel / Urho;

    // Reset total entropy
    Pr[i][j][k][KTOT] =
        (gam - 1.) * Pr[i][j][k][UU] * pow(Pr[i][j][k][RHO], -gam);
  } // ZSLOOP
}
#endif // RADIATION

void fixup_electrons(grid_prim_type P) {
  timer_start(TIMER_ELECTRON);
#pragma omp parallel for collapse(3) schedule(dynamic)
  ZLOOP { fixup_electrons_1zone(P[i][j][k]); }
  timer_stop(TIMER_ELECTRON);
}

void fixup_electrons_1zone(double P[NVAR]) {
  double kelmax =
      P[KTOT] * pow(P[RHO], gam - game) / (tptemin + (gam - 1.) / (game - 1.));
  double kelmin =
      P[KTOT] * pow(P[RHO], gam - game) / (tptemax + (gam - 1.) / (game - 1.));

  // Replace NANs with cold electrons
  if (isnan(P[KEL]))
    P[KEL] = kelmin;

  // Enforce maximum Tp/Te
  P[KEL] = MY_MAX(P[KEL], kelmin);

  // Enforce minimum Tp/Te
  P[KEL] = MY_MIN(P[KEL], kelmax);
}
#endif // EOS == EOS_TYPE_GAMMA
#endif // ELECTRONS


text:
#include <math.h>
#include <stdio.h>
//#include "gsl_headers.h"
#include <iostream>
//#include "laplaceEP.h"
//#include "matrix.h"
#include <time.h>

//#include <dai/alldai.h>
//#include <dai/bp.h>
// #include <map>
// #include <boost/shared_ptr.hpp>
// #include <string>
//
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include "gsl/gsl_cdf.h"
#include "gsl/gsl_randist.h"

//Transformations
double fre_1(double x, double func, double mu, double w);
double f_1(double x, double func, double mu, double w);
// double f(double x, double w);

// General functions
double compute_vector_mean(int N, double missing, gsl_vector *v);
double compute_vector_var(int N, double missing, gsl_vector *v);
double compute_vector_max(int N, double missing, gsl_vector *v);
double compute_vector_min(int N, double missing, gsl_vector *v);
int factorial (int N);
int poissrnd(double lambda);
//gsl_matrix *double2gsl(double *Amat, int nRows, int nCols);
void matrix_multiply(gsl_matrix *A,gsl_matrix *B,gsl_matrix *C,double alpha,double beta,CBLAS_TRANSPOSE_t TransA,CBLAS_TRANSPOSE_t TransB);
double *column_to_row_major_order(double *A,int nRows,int nCols);
//double *row_to_column_major_order(double *A,int nRows,int nCols);
double det_get(gsl_matrix *Amat, int Arows, int Acols, int inPlace);
double lndet_get(gsl_matrix *Amat, int Arows, int Acols, int inPlace);
// gsl_matrix *inverse(gsl_matrix *Amat, int Asize);
void inverse(gsl_matrix *Amat, int Asize);
double trace(gsl_matrix *Amat, int Asize);
double logFun(double x);
double expFun(double x);
//Sampling functions
int mnrnd(double *p, int nK);
void mvnrnd(gsl_vector *X, gsl_matrix *Sigma,gsl_vector *Mu, int K, const gsl_rng *seed);
double truncnormrnd(double mu, double sigma, double xlo, double xhi);

max_stars_repo_name: ferjorosa/test-glfm

max_issues_repo_name: ferjorosa/test-glfm

max_forks_repo_name: ferjorosa/test-glfm


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

#ifndef _MCAS_HSTORE_NUPM_H
#define _MCAS_HSTORE_NUPM_H

#include "pool_manager.h"

#include "alloc_key.h" /* AK_FORMAL */
#include "hstore_nupm_types.h"
#include "hstore_open_pool.h"
#include "persister_nupm.h"

#include <common/string_view.h>
#include <gsl/pointers>

#include <cstring> /* strerror */

#include <cinttypes> /* PRIx64 */
#include <cstdlib> /* getenv */

struct dax_manager;

template <typename PersistData, typename Heap>
  struct region;

#pragma GCC diagnostic push
/* Note: making enable_shared_from_this private avoids the non-virtual-dtor error but
 * generates a different error with no error text (G++ 5.4.0)
 */
#pragma GCC diagnostic ignored "-Wnon-virtual-dtor"

/* Region is region<persist_data_t, heap_rc>, Table is hstore::table_t Allocator is table_t::allocator_type, LockType is hstore::locK_type_t */
template <typename Region, typename Table, typename Allocator, typename LockType>
  struct hstore_nupm
    : public pool_manager<::open_pool<non_owner<Region>>>
  {
    using region_type = Region;
  private:
    using table_t = Table;
    using allocator_t = Allocator;
    using lock_type_t = LockType;
    using string_view = common::string_view;
  public:
    using open_pool_handle = ::open_pool<non_owner<region_type>>;
    using base = pool_manager<open_pool_handle>;
  private:
    std::unique_ptr<dax_manager> _dax_manager;
    unsigned _numa_node;

    static unsigned name_to_numa_node(const string_view name);
  public:
    hstore_nupm(unsigned debug_level_, const string_view, const string_view name_, std::unique_ptr<dax_manager> mgr_);

    virtual ~hstore_nupm();

    const std::unique_ptr<dax_manager> & get_dax_manager() const override { return _dax_manager; }
    void pool_create_check(std::size_t) override;

    auto pool_create_1(
      const pool_path &path_
      , std::size_t size_
    ) -> nupm::region_descriptor override;

    auto pool_create_2(
      AK_FORMAL
      const nupm::region_descriptor &rac
      , component::IKVStore::flags_t flags
      , std::size_t expected_obj_count
    ) -> std::unique_ptr<open_pool_handle> override;

    nupm::region_descriptor pool_open_1(
      const pool_path &path_
    ) override;

    auto pool_open_2(
      AK_FORMAL
      const nupm::region_descriptor & v_
      , component::IKVStore::flags_t flags_
    ) -> std::unique_ptr<open_pool_handle> override;

    void pool_close_check(const string_view) override;

    void pool_delete(const pool_path &path_) override;

    /* ERROR: want get_pool_regions(<proper type>, std::vector<::iovec>&) */
    nupm::region_descriptor pool_get_regions(const open_pool_handle &) const override;
  };
#pragma GCC diagnostic pop

#include "hstore_nupm.tcc"

#endif

max_stars_repo_name: fQuinzan/mcas

max_issues_repo_name: fQuinzan/mcas

max_forks_repo_name: fQuinzan/mcas


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code for the initialization of the instrumental noise for LIGO/VIRGO detectors.
 *
 */

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "LLVnoise.h"

/******************************************************************************/
/****** Global variables storing interpolating splines for the noise PSD ******/

gsl_spline* __LLVSimFD_LHONoiseSpline_init = NULL; /* for initialization only */
gsl_spline** const __LLVSimFD_LHONoiseSpline = &__LLVSimFD_LHONoiseSpline_init;
gsl_spline* __LLVSimFD_LLONoiseSpline_init = NULL; /* for initialization only */
gsl_spline** const __LLVSimFD_LLONoiseSpline = &__LLVSimFD_LLONoiseSpline_init;
gsl_spline* __LLVSimFD_VIRGONoiseSpline_init = NULL; /* for initialization only */
gsl_spline** const __LLVSimFD_VIRGONoiseSpline = &__LLVSimFD_VIRGONoiseSpline_init;
gsl_interp_accel* __LLVSimFD_LHONoiseAccel_init = NULL; /* for initialization only */
gsl_interp_accel** const __LLVSimFD_LHONoiseAccel = &__LLVSimFD_LHONoiseAccel_init;
gsl_interp_accel* __LLVSimFD_LLONoiseAccel_init = NULL; /* for initialization only */
gsl_interp_accel** const __LLVSimFD_LLONoiseAccel = &__LLVSimFD_LLONoiseAccel_init;
gsl_interp_accel* __LLVSimFD_VIRGONoiseAccel_init = NULL; /* for initialization only */
gsl_interp_accel** const __LLVSimFD_VIRGONoiseAccel = &__LLVSimFD_VIRGONoiseAccel_init;
double __LLVSimFD_LHONoise_fLow = 0;
double __LLVSimFD_LHONoise_fHigh = 0;
double __LLVSimFD_LLONoise_fLow = 0;
double __LLVSimFD_LLONoise_fHigh = 0;
double __LLVSimFD_VIRGONoise_fLow = 0;
double __LLVSimFD_VIRGONoise_fHigh = 0;
int __LLVSimFD_Noise_setup = FAILURE;

/* The number of points in the noise data files - required as the Read_Vector function needs a gsl vector already initialized to the right length */
//#define noisedata_pts 3000
#define noisedata_pts 16365

/**************************************************************/
/****** Functions loading and evaluating the noise PSD  *******/

/* Function parsing the environment variable $LLV_NOISE_DATA_PATH and trying to run LLVSimFD_Noise_Init in each */
int LLVSimFD_Noise_Init_ParsePath(void)
{
  if (!__LLVSimFD_Noise_setup) return(SUCCESS);

  int ret = FAILURE;
  char *envpath = NULL;
  char path[32768];
  char *brkt, *word;
  envpath = getenv("LLV_NOISE_DATA_PATH");
  if(!envpath) {
    printf("Error: the environment variable LLV_NOISE_DATA_PATH, giving the path to the noise data, seems undefined\n");
    exit(1);
  }
  strncpy(path, envpath, sizeof(path));

  for(word=strtok_r(path,":",&brkt); word; word=strtok_r(NULL,":",&brkt))
  {
    //printf("%s\n", word);
    ret = LLVSimFD_Noise_Init(word);
    if(ret == SUCCESS) break;
  }
  if(ret!=SUCCESS) {
    printf("Error: unable to find LLVSimFD noise data files in $LLV_NOISE_DATA_PATH\n");
    exit(FAILURE);
  }
  __LLVSimFD_Noise_setup = ret;
  return(ret);
}
/* Function loading the noise data from a directory */
int LLVSimFD_Noise_Init(const char dir[]) {
  if(!__LLVSimFD_Noise_setup) {
    printf("Error: LLVSimFD noise was already set up!");
    exit(1);
  }

  /* Loading noise data in gsl_vectors */
  int ret = SUCCESS;
  gsl_matrix* noise_LHO = gsl_matrix_alloc(noisedata_pts, 2);
  gsl_matrix* noise_LLO = gsl_matrix_alloc(noisedata_pts, 2);
  gsl_matrix* noise_VIRGO = gsl_matrix_alloc(noisedata_pts, 2);
  char* file_LIGO = malloc(strlen(dir)+64);
  char* file_VIRGO = malloc(strlen(dir)+64);
  //sprintf(file_LIGO, "%s", "LIGO-P1200087-v18-aLIGO_DESIGN.txt");
  //sprintf(file_VIRGO, "%s", "LIGO-P1200087-v18-AdV_DESIGN.txt");
  sprintf(file_LIGO, "%s", "aLIGO_sensitivity.dat");
  sprintf(file_VIRGO, "%s", "aVirgo_sensitivity.dat");
  ret |= Read_Text_Matrix(dir, file_LIGO, noise_LHO);
  ret |= Read_Text_Matrix(dir, file_LIGO, noise_LLO);
  ret |= Read_Text_Matrix(dir, file_VIRGO, noise_VIRGO);

  if(ret==FAILURE) {
    printf("Error: problem reading LLV noise data.");
    exit(1);
  }

  /* Linear interpolation of the data, after setting the gsl_spline structures */
  else if(ret==SUCCESS) {
    /* Extracting te vectors for the frequencies and data */
    gsl_vector* noise_LHO_freq = gsl_vector_alloc(noisedata_pts);
    gsl_vector* noise_LLO_freq = gsl_vector_alloc(noisedata_pts);
    gsl_vector* noise_VIRGO_freq = gsl_vector_alloc(noisedata_pts);
    gsl_vector* noise_LHO_data = gsl_vector_alloc(noisedata_pts);
    gsl_vector* noise_LLO_data = gsl_vector_alloc(noisedata_pts);
    gsl_vector* noise_VIRGO_data = gsl_vector_alloc(noisedata_pts);
    gsl_matrix_get_col(noise_LHO_freq, noise_LHO, 0);
    gsl_matrix_get_col(noise_LLO_freq, noise_LLO, 0);
    gsl_matrix_get_col(noise_VIRGO_freq, noise_VIRGO, 0);
    gsl_matrix_get_col(noise_LHO_data, noise_LHO, 1);
    gsl_matrix_get_col(noise_LLO_data, noise_LLO, 1);
    gsl_matrix_get_col(noise_VIRGO_data, noise_VIRGO, 1);
    /* Setting the global variables that indicate the range in frequency of these splines */
    __LLVSimFD_LHONoise_fLow = gsl_vector_get(noise_LHO_freq, 0);
    __LLVSimFD_LHONoise_fHigh = gsl_vector_get(noise_LHO_freq, noise_LHO_freq->size - 1);
    __LLVSimFD_LLONoise_fLow = gsl_vector_get(noise_LLO_freq, 0);
    __LLVSimFD_LLONoise_fHigh = gsl_vector_get(noise_LLO_freq, noise_LLO_freq->size - 1);
    __LLVSimFD_VIRGONoise_fLow = gsl_vector_get(noise_VIRGO_freq, 0);
    __LLVSimFD_VIRGONoise_fHigh = gsl_vector_get(noise_VIRGO_freq, noise_VIRGO_freq->size - 1);
    /* Initializing the splines and accelerators */
    *__LLVSimFD_LHONoiseSpline = gsl_spline_alloc(gsl_interp_linear, noisedata_pts);
    *__LLVSimFD_LLONoiseSpline = gsl_spline_alloc(gsl_interp_linear, noisedata_pts);
    *__LLVSimFD_VIRGONoiseSpline = gsl_spline_alloc(gsl_interp_linear, noisedata_pts);
    *__LLVSimFD_LHONoiseAccel = gsl_interp_accel_alloc();
    *__LLVSimFD_LLONoiseAccel = gsl_interp_accel_alloc();
    *__LLVSimFD_VIRGONoiseAccel = gsl_interp_accel_alloc();
    gsl_spline_init(*__LLVSimFD_LHONoiseSpline, gsl_vector_const_ptr(noise_LHO_freq, 0), gsl_vector_const_ptr(noise_LHO_data, 0), noisedata_pts);
    gsl_spline_init(*__LLVSimFD_LLONoiseSpline, gsl_vector_const_ptr(noise_LLO_freq, 0), gsl_vector_const_ptr(noise_LLO_data, 0), noisedata_pts);
    gsl_spline_init(*__LLVSimFD_VIRGONoiseSpline, gsl_vector_const_ptr(noise_VIRGO_freq, 0), gsl_vector_const_ptr(noise_VIRGO_data, 0), noisedata_pts);
    /* Setting the global tag to success and clean up */
    gsl_matrix_free(noise_LHO);
    gsl_matrix_free(noise_LLO);
    gsl_matrix_free(noise_VIRGO);
    gsl_vector_free(noise_LHO_freq);
    gsl_vector_free(noise_LLO_freq);
    gsl_vector_free(noise_VIRGO_freq);
    gsl_vector_free(noise_LHO_data);
    gsl_vector_free(noise_LLO_data);
    gsl_vector_free(noise_VIRGO_data);
    __LLVSimFD_Noise_setup = SUCCESS;
  }

  /* Cleaning and output */
  free(file_LIGO);
  free(file_VIRGO);
  return(ret);
}

/* The noise functions themselves */
double NoiseSnLHO(const double f) {
  if(__LLVSimFD_Noise_setup==FAILURE) {
    printf("Error: noise interpolation has not been set up\n");
    exit(1);
  }
  if ((f < __LLVSimFD_LHONoise_fLow) || (f > __LLVSimFD_LHONoise_fHigh)) {
    return INFINITY;
  }
  else {
    double sqrtSn = gsl_spline_eval(*__LLVSimFD_LHONoiseSpline, f, *__LLVSimFD_LHONoiseAccel);
    return sqrtSn * sqrtSn;
  }
}
double NoiseSnLLO(const double f) {
  if(__LLVSimFD_Noise_setup==FAILURE) {
    printf("Error: noise interpolation has not been set up\n");
    exit(1);
  }
  if ((f < __LLVSimFD_LLONoise_fLow) || (f > __LLVSimFD_LLONoise_fHigh)) {
    return INFINITY;
  }
  else {
    double sqrtSn = gsl_spline_eval(*__LLVSimFD_LLONoiseSpline, f, *__LLVSimFD_LLONoiseAccel);
    return sqrtSn * sqrtSn;
  }
}
double NoiseSnVIRGO(const double f) {
  if(__LLVSimFD_Noise_setup==FAILURE) {
    printf("Error: noise interpolation has not been set up\n");
    exit(1);
  }
  if ((f < __LLVSimFD_VIRGONoise_fLow) || (f > __LLVSimFD_VIRGONoise_fHigh)) {
    return INFINITY;
  }
  else {
    double sqrtSn = gsl_spline_eval(*__LLVSimFD_VIRGONoiseSpline, f, *__LLVSimFD_VIRGONoiseAccel);
    return sqrtSn * sqrtSn;
  }
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#include "sly/defines.h"

#include <gsl/gsl>
#include <memory.h>
#include <optional>
#include <stdio.h>

#ifdef _WIN32
    #include <windows.h>
#endif

#include "sly/debug/logger.h"
#include "sly/enum.h"
#include "sly/errors.h"
#include "sly/macros.h"
#include "sly/retval.h"
#include "sly/types.h"

max_stars_repo_name: Gibbeon/sly

max_issues_repo_name: Gibbeon/sly

max_forks_repo_name: Gibbeon/sly


text:
#include <math.h>
#include <gsl/gsl_math.h>

#include <options/options.h>

#include "../csm/csm_all.h"

int main(int argc, const char * argv[]) {
	sm_set_program_name(argv[0]);

	const char*input_filename;
	const char*output_filename;
	double diff[3] = {0, 0, 0};

	double laser[3] ={0,0,0};

	double omega[2] ={0,0};
	double omega_vel = 0;
	double vel[2] ={0,0};
	double vel_omega = 0;

	struct option* ops = options_allocate(15);
	options_string(ops, "in", &input_filename, "stdin", "input file");
	options_string(ops, "out", &output_filename, "stdout", "output file");
	options_double(ops, "l_x", &(laser[0]), 0.0, "laser x (m)");
	options_double(ops, "l_y", &(laser[1]), 0.0, "laser y (m)");
	options_double(ops, "l_theta", &(laser[2]), 0.0, "laser theta (rad)");

	options_double(ops, "omega0", &(omega[0]), 0.0, "omega (rad)");
	options_double(ops, "omega1", &(omega[1]), 0.0, "omega (linear)");
	options_double(ops, "omega_vel", &(omega_vel), 0.0, "omega x vel");
	options_double(ops, "vel0", &(vel[0]), 0.0, "vel (m)");
	options_double(ops, "vel1", &(vel[1]), 0.0, "vel (linear)");
	options_double(ops, "vel_omega", &(vel_omega), 0.0, "vel x omega");

	if(!options_parse_args(ops, argc, argv)) {
		fprintf(stderr, " Corrects bias in odometry.\n");
		options_print_help(ops, stderr);
		return -1;
	}

	FILE * input_stream = open_file_for_reading(input_filename);
	FILE *output_stream = open_file_for_writing(output_filename);

	if(!input_stream || !output_stream) return -1;

	LDP laser_ref = ld_read_smart(input_stream);
	if(!laser_ref) {
		sm_error("Cannot read first scan.\n");
		return -2;
	}

	copy_d(laser_ref->odometry, 3, laser_ref->estimate);

	double old_odometry[3];
	copy_d(laser_ref->odometry, 3, old_odometry);

	LDP laser_sens;
	ld_write_as_json(laser_ref, output_stream);

	while((laser_sens = ld_read_smart(input_stream))) {
		double guess[3], old_guess[3];
		pose_diff_d(laser_sens->odometry, old_odometry, guess);
		pose_diff_d(laser_sens->odometry, old_odometry, old_guess);
		copy_d(laser_sens->odometry, 3, old_odometry);

		if(fabs(guess[2]) > 1e-7)
			guess[2] = old_guess[2] * omega[1] + omega_vel * guess[0] + omega[0];

		if(fabs(guess[0]) > 1e-7)
			guess[0] = old_guess[0] * vel[1] + vel_omega * guess[2] + vel[0];

		fprintf(stderr, "odo: %f %f %f Corr: %f rad \t%f m  \n", guess[0], guess[1], guess[2], guess[2]-old_guess[2], guess[0]-old_guess[0]);

		oplus_d(laser_ref->odometry, guess, laser_sens->odometry);
		oplus_d(laser_sens->odometry, laser, laser_sens->estimate);

		fprintf(stderr, "ref odo: %s  ref est: %s \n",
			friendly_pose(laser_ref->odometry),
			friendly_pose(laser_ref->estimate));
		fprintf(stderr, "sens odo: %s  sens est: %s \n",
			friendly_pose(laser_sens->odometry),
			friendly_pose(laser_sens->estimate));

			ld_write_as_json(laser_sens, output_stream);

		ld_free(laser_ref);
		laser_ref = laser_sens;
	}

	return 0;
}

max_stars_repo_name: ozaslan/basics

max_issues_repo_name: ozaslan/basics

max_forks_repo_name: ozaslan/basics


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/

#ifndef __DRIVERS_common_utils_h__
#define __DRIVERS_common_utils_h__

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

#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscpc.h>
#include <petscksp.h>
#include <petscsnes.h>

#include <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include "petsc/private/vecimpl.h"
     #include "petsc/private/matimpl.h"
     #include "petsc/private/pcimpl.h"
     #include "petsc/private/kspimpl.h"
     #include "petsc/private/snesimpl.h"
  #else
     #include "petsc-private/vecimpl.h"
     #include "petsc-private/matimpl.h"
     #include "petsc-private/pcimpl.h"
     #include "petsc-private/kspimpl.h"
     #include "petsc-private/snesimpl.h"
  #endif
#else
  #include "private/vecimpl.h"
  #include "private/matimpl.h"
  #include "private/pcimpl.h"
  #include "private/kspimpl.h"
  #include "private/snesimpl.h"
#endif

#include "../../../../StgDomain/Utils/src/PETScCompatibility.h"

#include "pc_GtKG.h"
#include "pc_ScaledGtKG.h"
#include "stokes_block_scaling.h"
#include "stokes_mvblock_scaling.h"
#include "ksp_scale.h"

/* BSSCR_create_execute_script.c */
int BSSCR_create_execute_script( void );

/* create_petsc_header.c */
PetscErrorCode BSSCR_GeneratePetscHeader_for_file( FILE *fd, MPI_Comm comm );
PetscErrorCode BSSCR_GeneratePetscHeader_for_viewer( PetscViewer viewer );

/* MatIdentityOperator.c */
//PetscErrorCode MatCreateIdentityOperator( Mat A, Mat *_I );

/* read_matrices.c */
//PetscErrorCode StokesReadOperators( MPI_Comm comm, Mat *K, Mat *G, Mat *D, Mat *C, Vec *f, Vec *h );

/* stokes_residual.c */
double BSSCR_StokesMomentumResidual( Mat K, Mat G, Vec F, Vec u, Vec p );
double BSSCR_StokesContinuityResidual( Mat G, Mat C, Vec H, Vec u, Vec p );

/* preconditioner.c */
PetscErrorCode BSSCR_StokesReadPCSchurMat( MPI_Comm comm, Mat *S );
PetscErrorCode BSSCR_StokesCreatePCSchur( Mat K, Mat G, PC pc_S );

/* nullspace */
/*
PetscErrorCode BSSCR_MatGtKinvG_ContainsConstantNullSpace(
		Mat K, Mat G, Mat M,Vec t1, Vec ustar,  Vec r, Vec l,
		KSP ksp, PetscTruth *has_cnst_nullsp );
PetscErrorCode BSSCR_VecRemoveConstNullspace(Vec v, Vec nsp_vec);
*/
/* PetscErrorCode BSSCR_MatContainsConstantNullSpace( Mat A, PetscTruth *has_cnst_nullsp ); */
//PetscErrorCode BSSCR_MatContainsConstNullSpace(Mat A, Vec nsp_vec, PetscTruth *has);

/* timed residuals */
//PetscErrorCode KSPLogDestroyMonitor(KSP ksp);
PetscErrorCode BSSCR_KSPLogSetMonitor(KSP ksp,PetscInt len,PetscInt *monitor_index);
PetscErrorCode BSSCR_KSPLogGetTimeHistory(KSP ksp,PetscInt monitor_index,PetscInt *na,PetscLogDouble **log);
PetscErrorCode BSSCR_KSPLogGetResidualHistory(KSP ksp,PetscInt monitor_index,PetscInt *na,PetscReal **rlog);
PetscErrorCode BSSCR_KSPLogGetResidualTimeHistory(KSP ksp,PetscInt monitor_index,PetscInt *na,PetscReal **rlog,PetscLogDouble **tlog);
PetscErrorCode BSSCR_KSPLogSolve(PetscViewer v,PetscInt monitor_index,KSP ksp);
PetscErrorCode BSSCR_BSSCR_KSPLogSolveSummary(PetscViewer v,PetscInt monitor_index,KSP ksp);

/* operator summary */
PetscErrorCode BSSCR_MatInfoLog(PetscViewer v,Mat A,const char name[]);
PetscErrorCode BSSCR_VecInfoLog(PetscViewer v,Vec x,const char name[]);
PetscErrorCode BSSCR_StokesCreateOperatorSummary( Mat K, Mat G, Mat C, Vec f, Vec h, const char filename[] );

/* standard output */
PetscErrorCode BSSCR_solver_output( KSP ksp, PetscInt monitor_index, const char res_file[], const char solver_conf[] );
PetscErrorCode BSSCR_KSPLogConvergenceRate(PetscViewer v,PetscInt monitor_index, KSP ksp);

/* list available operations */
PetscErrorCode BSSCR_MatListOperations( Mat A, PetscViewer v );
PetscErrorCode BSSCR_VecListOperations( Vec x, PetscViewer v );
PetscErrorCode BSSCR_KSPListOperations( KSP ksp, PetscViewer v );
PetscErrorCode BSSCR_PCListOperations( PC pc, PetscViewer v );
PetscErrorCode BSSCR_SNESListOperations( SNES snes, PetscViewer v );

/* Register any petsc objects defined in drivers */
PetscErrorCode BSSCR_PetscExtStokesSolversInitialize( void );
PetscErrorCode BSSCR_PetscExtStokesSolversFinalize( void );

/* stokes block */
/*
PetscErrorCode BSSCR_CreateStokesBlockOperators( MPI_Comm comm,
		Mat K, Mat G, Mat D, Mat C,
		Vec u, Vec p, Vec f, Vec h,
		Mat *A, Vec *x, Vec *b );
*/
/* stokes output */
PetscErrorCode BSSCR_stokes_output( PetscViewer v, Mat stokes_A, Vec stokes_b, Vec stokes_x, KSP ksp, PetscInt monitor_index );

//PetscErrorCode BSSCR_NSPRemoveAll(Vec v, void *_data);
//PetscErrorCode BSSCR_CheckNullspace(KSP ksp_S, Mat S, Vec h_hat, MatStokesBlockScaling BA, Vec * _nsp_vec);
#endif


text:
#ifndef CTETRA_SUBMESH_H
#define CTETRA_SUBMESH_H

#include <math.h>
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "ecache.h"

void submesh_ijk_to_k(int na, int nb, int nc, int i, int j, int k, double k_opt[3]);

int submesh_ijk_index(int na, int nb, int nc, int i, int j, int k);

int **subcells_around_ijk(int na, int nb, int nc, int i, int j, int k, int *subcell_num);

void get_k_orig(double k_opt[3], int G_order[3], int G_neg[3], double k_orig[3]);

void MinMaxVals(EnergyCache *Ecache, double *emin, double *emax);

void OptimizeGs(gsl_matrix *R, int G_order[3], int G_neg[3]);

#endif // CTETRA_SUBMESH_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra