text:
/* linalg/rqr.c
 *
 * Copyright (C) 2019 Patrick Alken, Julien Langou
 *
 * 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 <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/*
 * this module contains routines for the QR factorization of a matrix
 * using the recursive Level 3 BLAS algorithm of Elmroth and Gustavson with
 * additional modifications courtesy of Julien Langou.
 */

static int unpack_Q1(gsl_matrix * Q);
static int unpack_Q2(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q);
static int aux_ULT(const gsl_matrix * L, gsl_matrix * U);
static int aux_mLU(gsl_matrix * A);
static int aux_ApUBT(const gsl_matrix * U, const gsl_matrix * B, gsl_matrix * A);

/*
gsl_linalg_QR_decomp_r()
  QR decomposition using Level 3 BLAS recursive algorithm of:

Elmroth, E. and Gustavson, F.G., 2000. Applying recursion to serial and parallel
  QR factorization leads to better performance. IBM Journal of Research and Development,
  44(4), pp.605-624.

Inputs: A - matrix to be factored, M-by-N with M >= N
        T - N-by-N upper triangular factor of block reflector

Return: success/error

Notes:
1) on output, upper triangle of A contains R; elements below the diagonal
are columns of V. The matrix Q is

Q = I - V T V^T

where T is upper triangular. Note that diag(T) = tau

2) implementation provided by Julien Langou
*/

int
gsl_linalg_QR_decomp_r (gsl_matrix * A, gsl_matrix * T)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M < N)
    {
      GSL_ERROR ("M must be >= N", GSL_EBADLEN);
    }
  else if (T->size1 != T->size2)
    {
      GSL_ERROR ("T matrix must be square", GSL_ENOTSQR);
    }
  else if (T->size1 != N)
    {
      GSL_ERROR ("T matrix does not match dimensions of A", GSL_EBADLEN);
    }
  else if (N == 1)
    {
      /* base case, compute householder transform for single column matrix */

      double * T00 = gsl_matrix_ptr(T, 0, 0);
      gsl_vector_view v = gsl_matrix_column(A, 0);

      *T00 = gsl_linalg_householder_transform(&v.vector);
      return GSL_SUCCESS;
    }
  else
    {
      /*
       * partition matrices:
       *
       *       N1  N2              N1  N2
       * N1 [ A11 A12 ] and  N1 [ T11 T12 ]
       * M2 [ A21 A22 ]      N2 [  0  T22 ]
       */
      int status;
      const size_t N1 = N / 2;
      const size_t N2 = N - N1;
      const size_t M2 = M - N1;

      gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, N1, N1);
      gsl_matrix_view A12 = gsl_matrix_submatrix(A, 0, N1, N1, N2);
      gsl_matrix_view A21 = gsl_matrix_submatrix(A, N1, 0, M2, N1);
      gsl_matrix_view A22 = gsl_matrix_submatrix(A, N1, N1, M2, N2);

      gsl_matrix_view T11 = gsl_matrix_submatrix(T, 0, 0, N1, N1);
      gsl_matrix_view T12 = gsl_matrix_submatrix(T, 0, N1, N1, N2);
      gsl_matrix_view T22 = gsl_matrix_submatrix(T, N1, N1, N2, N2);

      gsl_matrix_view m;

      /*
       * Eq. 2: recursively factor
       *
       * [ A11 ] = Q1 [ R11 ]
       * [ A21 ]      [  0  ]
       *                                          N1
       * Note: Q1 = I - V1 T11 V1^T, where V1 = [ V11 ] N1
       *                                        [ V21 ] M2
       */
      m = gsl_matrix_submatrix(A, 0, 0, M, N1);
      status = gsl_linalg_QR_decomp_r(&m.matrix, &T11.matrix);
      if (status)
        return status;

      /*
       * Eq. 3:
       *
       * [ R12 ] := Q1^T [ A12 ] = [ A12 ] - [ V11 W ]
       * [ A22 ]         [ A22 ]   [ A22 ]   [ V21 W ]
       *
       * where W = T11^T (V11^T A12 + V21^T A22), and using T12 as temporary storage
       */
      gsl_matrix_memcpy(&T12.matrix, &A12.matrix);                                                       /* W := A12 */
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasTrans, CblasUnit, 1.0, &A11.matrix, &T12.matrix);       /* W := V11^T * A12 */
      gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &A21.matrix, &A22.matrix, 1.0, &T12.matrix);         /* W := W + V21^T * A22 */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, 1.0, &T11.matrix, &T12.matrix);    /* W := T11^T * W */
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &A21.matrix, &T12.matrix, 1.0, &A22.matrix);      /* A22 = A22 - V21 * W */
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, 1.0, &A11.matrix, &T12.matrix);     /* tmp = V11 * W */
      gsl_matrix_sub(&A12.matrix, &T12.matrix);                                                          /* R12 := A12 - V11 * W */

      /*
       * Eq. 4: recursively factor
       *
       * A22 = Q2~ R22
       *
       *              N1 M2
       * Note: Q2 = [ I   0  ] N1
       *            [ 0  Q2~ ] M2
       */
      status = gsl_linalg_QR_decomp_r(&A22.matrix, &T22.matrix);
      if (status)
        return status;

      /*
       * Eq. 13: update T12 := -T11 * V1^T * V2 * T22
       *
       * where:
       *
       *        N1                N2
       * V1 = [ V11 ] N1   V2 = [  0  ] N1
       *      [ V21 ] N2        [ V22 ] N2
       *      [ V31 ] M-N       [ V32 ] M-N
       *
       * Note: V1^T V2 = V21^T V22 + V31^T V32
       * Also, V11, V22 are unit lower triangular
       */

      m = gsl_matrix_submatrix(&A21.matrix, 0, 0, N2, N1);                                               /* V21 */
      gsl_matrix_transpose_memcpy(&T12.matrix, &m.matrix);                                               /* T12 := V21^T */

      m = gsl_matrix_submatrix(A, N1, N1, N2, N2);                                                       /* V22 */
      gsl_blas_dtrmm(CblasRight, CblasLower, CblasNoTrans, CblasUnit, 1.0, &m.matrix, &T12.matrix);      /* T12 := V21^T * V22 */

      if (M > N)
        {
          gsl_matrix_view V31 = gsl_matrix_submatrix(A, N, 0, M - N, N1);
          gsl_matrix_view V32 = gsl_matrix_submatrix(A, N, N1, M - N, N2);

          gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &V31.matrix, &V32.matrix, 1.0, &T12.matrix);     /* T12 := T12 + V31^T * V32 */
        }

      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &T11.matrix, &T12.matrix); /* T12 := -T11 * T12 */
      gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, &T22.matrix, &T12.matrix); /* T12 := T12 * T22 */

      return GSL_SUCCESS;
    }
}

/* Solves the square system A x = b for x using the QR factorisation,
 *
 *  R x = Q^T b
 *
 * where Q = I - V T V^T
 */

int
gsl_linalg_QR_solve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x)
{
  const size_t N = QR->size2;

  if (QR->size1 != N)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (T->size1 != QR->size1 || T->size2 != QR->size2)
    {
      GSL_ERROR ("T matrix must be N-by-N", GSL_EBADLEN);
    }
  else if (N != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (N != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute Q^T b = [I - V T^T V^T] b */

      /* x := V^T b */
      gsl_vector_memcpy(x, b);
      gsl_blas_dtrmv(CblasLower, CblasTrans, CblasUnit, QR, x);

      /* x = T^T * x */
      gsl_blas_dtrmv(CblasUpper, CblasTrans, CblasNonUnit, T, x);

      /* x = V * x */
      gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasUnit, QR, x);

      /* x = b - V * x */
      for (i = 0; i < N; ++i)
        {
          double * xi = gsl_vector_ptr(x, i);
          double bi = gsl_vector_get(b, i);
          *xi = bi - (*xi);
        }

      /* Solve R x = Q^T b, storing x in-place */
      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      return GSL_SUCCESS;
    }
}

/* Find the least squares solution to the overdetermined system
 *
 *   A x = b
 *
 * for M >= N using the QR factorization A = Q R.
 *
 * Inputs: QR   - [R; V] matrix, M-by-N
 *         T    - upper triangular block reflector, N-by-N
 *         b    - right hand side, size M
 *         x    - (output) solution, size M
 *                x(1:N) = least squares solution vector
 *                x(N+1:M) = vector whose norm equals ||b - Ax||
 *         work - workspace, size N
 */

int
gsl_linalg_QR_lssolve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x, gsl_vector * work)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M < N)
    {
      GSL_ERROR ("QR matrix must have M >= N", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix must be N-by-N", GSL_EBADLEN);
    }
  else if (M != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (M != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (N != work->size)
    {
      GSL_ERROR ("matrix size must match work size", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view R = gsl_matrix_const_submatrix (QR, 0, 0, N, N);
      gsl_vector_view x1 = gsl_vector_subvector(x, 0, N);

      /* compute x = Q^T b */
      gsl_vector_memcpy(x, b);
      gsl_linalg_QR_QTvec_r (QR, T, x, work);

      /* Solve R x = Q^T b */
      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &R.matrix, &x1.vector);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_QR_unpack_r()
  Unpack matrices Q and R

Inputs: QR - packed QR format, M-by-N
        T  - block reflector matrix, N-by-N
        Q  - (output) Q matrix, M-by-M
        R  - (output) R matrix, N-by-N

Return: success/error

Notes:
1) Implementation provided by Julien Langou
*/

int
gsl_linalg_QR_unpack_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q, gsl_matrix * R)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M < N)
    {
      GSL_ERROR ("M must be >= N", GSL_EBADLEN);
    }
  else if (Q->size1 != M || Q->size2 != M)
    {
      GSL_ERROR ("Q matrix must be M-by-M", GSL_EBADLEN);
    }
  else if (R->size1 != N || R->size2 != N)
    {
      GSL_ERROR ("R matrix must be N-by-N", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix must be N-by-N", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view RV = gsl_matrix_const_submatrix(QR, 0, 0, N, N);
      gsl_matrix_view Q1 = gsl_matrix_submatrix(Q, 0, 0, M, N);
      gsl_matrix_view m;

      /*
       * set Q1 = [ T ]
       *          [ V ]
       */
      m = gsl_matrix_submatrix(Q, 0, 0, N, N);
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &m.matrix, T);
      gsl_matrix_tricpy(CblasLower, CblasUnit, &m.matrix, &RV.matrix);

      if (M > N)
        {
          gsl_matrix_const_view tmp = gsl_matrix_const_submatrix(QR, N, 0, M - N, N);
          m = gsl_matrix_submatrix(Q, N, 0, M - N, N);
          gsl_matrix_memcpy(&m.matrix, &tmp.matrix);
        }

      unpack_Q1(&Q1.matrix);

      if (M > N)
        {
          gsl_matrix_view Q2 = gsl_matrix_submatrix(Q, 0, N, M, M - N);
          unpack_Q2(QR, T, &Q2.matrix);
        }

      /* copy R */
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, &RV.matrix);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_QR_QTvec_r()
  Apply M-by-M Q^T to the M-by-1 vector b

Inputs: QR   - [R; V] matrix encoded by gsl_linalg_QR_decomp_r
        T    - block reflector matrix
        b    - M-by-1 vector replaced by Q^T b on output
        work - workspace, length N

Notes:
1) Q^T b = (I - V T^T V^T) b
         = b - V T^T [ V1^T V2^T ] [ b1 ]
                                   [ b2 ]
         = b - V T^T [ V1^T b1 + V2^T b2 ]
         = [ b1 ] - [ V1 w ]
           [ b2 ]   [ V2 w ]

where w = T^T ( V1^T b1 + V2^T b2 )
*/

int
gsl_linalg_QR_QTvec_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_vector * b, gsl_vector * work)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M < N)
    {
      GSL_ERROR ("M must be >= N", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix must be N-by-N", GSL_EBADLEN);
    }
  else if (b->size != M)
    {
      GSL_ERROR ("b vector must have length M", GSL_EBADLEN);
    }
  else if (work->size != N)
    {
      GSL_ERROR ("workspace must be length N", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view V1 = gsl_matrix_const_submatrix(QR, 0, 0, N, N);
      gsl_vector_view b1 = gsl_vector_subvector(b, 0, N);
      gsl_vector_view b2;

      /* work := V1^T b1 */
      gsl_vector_memcpy(work, &b1.vector);
      gsl_blas_dtrmv(CblasLower, CblasTrans, CblasUnit, &V1.matrix, work);

      if (M > N)
        {
          gsl_matrix_const_view V2 = gsl_matrix_const_submatrix(QR, N, 0, M - N, N);

          /* work = work + V2^T b2 */
          b2 = gsl_vector_subvector(b, N, M - N);
          gsl_blas_dgemv(CblasTrans, 1.0, &V2.matrix, &b2.vector, 1.0, work);
        }

      /* work = T^T * work */
      gsl_blas_dtrmv(CblasUpper, CblasTrans, CblasNonUnit, T, work);

      if (M > N)
        {
          /* b2 = b2 - V2 * work */
          gsl_matrix_const_view V2 = gsl_matrix_const_submatrix(QR, N, 0, M - N, N);
          gsl_blas_dgemv(CblasNoTrans, -1.0, &V2.matrix, work, 1.0, &b2.vector);
        }

      /* b1 = b1 - V1 * work */
      gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasUnit, &V1.matrix, work);
      gsl_vector_sub(&b1.vector, work);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_QR_QTmat_r()
  Apply M-by-M Q^T to the M-by-K matrix B

Inputs: QR   - [R; V] matrix encoded by gsl_linalg_QR_decomp_r
        T    - block reflector matrix
        B    - M-by-K matrix replaced by Q^T B on output
        work - N-by-K workspace

Notes:
1) Provided by Julien Langou
*/

int
gsl_linalg_QR_QTmat_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * B, gsl_matrix * work)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;
  const size_t K = B->size2;

  if (M < N)
    {
      GSL_ERROR ("M must be >= N", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix must be N-by-N", GSL_EBADLEN);
    }
  else if (B->size1 != M)
    {
      GSL_ERROR ("B matrix must have M rows", GSL_EBADLEN);
    }
  else if (work->size1 != N || work->size2 != K)
    {
      GSL_ERROR ("workspace must be N-by-K", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view V1 = gsl_matrix_const_submatrix(QR, 0, 0, N, N);
      gsl_matrix_view B1 = gsl_matrix_submatrix(B, 0, 0, N, K);
      gsl_matrix_view B2;

      /* work := V1^T B1 */
      gsl_matrix_memcpy(work, &B1.matrix);
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasTrans, CblasUnit, 1.0, &V1.matrix, work);

      if (M > N)
        {
          gsl_matrix_const_view V2 = gsl_matrix_const_submatrix(QR, N, 0, M - N, N);

          /* work = work + V2^T B2 */
          B2 = gsl_matrix_submatrix(B, N, 0, M - N, K);
          gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &V2.matrix, &B2.matrix, 1.0, work);
        }

      /* work = T^T * work */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, 1.0, T, work);

      if (M > N)
        {
          /* B2 = B2 - V2 * work */
          gsl_matrix_const_view V2 = gsl_matrix_const_submatrix(QR, N, 0, M - N, N);
          gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &V2.matrix, work, 1.0, &B2.matrix);
        }

      /* B1 = B1 - V1 * work */
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, 1.0, &V1.matrix, work);
      gsl_matrix_sub(&B1.matrix, work);

      return GSL_SUCCESS;
    }
}

/*
unpack_Q1()
  Compute Q_1

Inputs: Q  - on input, contains T in upper triangle and V in lower trapezoid
             on output, contains Q_1
             M-by-N

Return: success/error

Notes:
1)                                                         N
Q1 = [ Q1 Q2 ] [ I_n ] = (I - V T V^T) [ I; 0 ] = [ I - V1 T V1^T ] N
               [  0  ]                            [   - V2 T V1^T ] M - N
*/

static int
unpack_Q1(gsl_matrix * Q)
{
  int status;
  const size_t M = Q->size1;
  const size_t N = Q->size2;
  gsl_matrix_view Q1 = gsl_matrix_submatrix(Q, 0, 0, N, N);
  gsl_vector_view diag = gsl_matrix_diagonal(&Q1.matrix);

  /* Q1 := T V1^T */
  status = aux_ULT(&Q1.matrix, &Q1.matrix);
  if (status)
    return status;

  if (M > N)
    {
      /* compute Q2 := - V2 T V1^T */
      gsl_matrix_view V2 = gsl_matrix_submatrix(Q, N, 0, M - N, N);
      gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &Q1.matrix, &V2.matrix);
    }

  /* Q1 := - V1 T V1^T */
  status = aux_mLU(&Q1.matrix);
  if (status)
    return status;

  /* Q1 := I - V1 T V1^T */
  gsl_vector_add_constant(&diag.vector, 1.0);

  return GSL_SUCCESS;
}

/*
unpack_Q2()
  Compute Q_2

Inputs: QR - [R; V] from QR_decomp_r, M-by-N
        T  - upper triangular T factor, N-by-N
        Q  - (output) Q_2 factor, M-by-(M-N)

Return: success/error

Notes:                         N   M-N
1) Since Q = I - V T V^T = M [ Q1  Q2  ], we have

           M-N
Q2 = Q [    0    ] N
       [ I_{M-N} ] M-N

So, Q2 = Q [ 0; I ] = (I - V T V^T) [ 0; I ] = [   - V1 T V2^T ]
                                               [ I - V2 T V2^T ]
*/

static int
unpack_Q2(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M <= N)
    {
      GSL_ERROR ("M must be > N", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix must be N-by-N", GSL_EBADLEN);
    }
  else if (Q->size1 != M || Q->size2 != (M - N))
    {
      GSL_ERROR ("Q matrix must be M-by-(M-N)", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view V1 = gsl_matrix_const_submatrix(QR, 0, 0, N, N);
      gsl_matrix_const_view V2 = gsl_matrix_const_submatrix(QR, N, 0, M - N, N);
      gsl_matrix_view Q1 = gsl_matrix_submatrix(Q, 0, 0, N, M - N);
      gsl_matrix_view Q2 = gsl_matrix_submatrix(Q, N, 0, M - N, M - N);
      gsl_vector_view diag = gsl_matrix_diagonal(&Q2.matrix);

      /* Q1 := V2^T */
      gsl_matrix_transpose_memcpy(&Q1.matrix, &V2.matrix);

      /* Q1 := - T V2^T */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, T, &Q1.matrix);

      /* Q2 := - V2 T V2^T */
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &V2.matrix, &Q1.matrix, 0.0, &Q2.matrix);

      /* Q2 := I - V2 T V2^T */
      gsl_vector_add_constant(&diag.vector, 1.0);

      /* Q1 := - V1 T V2^T */
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, 1.0, &V1.matrix, &Q1.matrix);

      return GSL_SUCCESS;
    }
}

/* U := U L^T for triangular matrices L and U; L is unit lower triangular */
static int
aux_ULT(const gsl_matrix * L, gsl_matrix * U)
{
  const size_t N = L->size1;

  if (N != L->size2)
    {
      GSL_ERROR ("L matrix must be square", GSL_ENOTSQR);
    }
  else if (U->size1 != N || U->size2 != N)
    {
      GSL_ERROR ("U matrix must be same size as L", GSL_EBADLEN);
    }
  else if (N == 1)
    {
      /* nothing to do */
      return GSL_SUCCESS;
    }
  else
    {
      int status;
      const size_t N1 = N / 2;
      const size_t N2 = N - N1;

      gsl_matrix_const_view L11 = gsl_matrix_const_submatrix(L, 0, 0, N1, N1);
      gsl_matrix_const_view L21 = gsl_matrix_const_submatrix(L, N1, 0, N2, N1);
      gsl_matrix_const_view L22 = gsl_matrix_const_submatrix(L, N1, N1, N2, N2);

      gsl_matrix_view U11 = gsl_matrix_submatrix(U, 0, 0, N1, N1);
      gsl_matrix_view U12 = gsl_matrix_submatrix(U, 0, N1, N1, N2);
      gsl_matrix_view U22 = gsl_matrix_submatrix(U, N1, N1, N2, N2);

      /* U12 = U12 * L22^T */
      gsl_blas_dtrmm(CblasRight, CblasLower, CblasTrans, CblasUnit, 1.0, &L22.matrix, &U12.matrix);

      /* U12 = U12 + U11 * L21^T */
      status = aux_ApUBT(&U11.matrix, &L21.matrix, &U12.matrix);
      if (status)
        return status;

      status = aux_ULT(&L11.matrix, &U11.matrix);
      if (status)
        return status;

      status = aux_ULT(&L22.matrix, &U22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

/* store -L*U in A */
static int
aux_mLU(gsl_matrix * A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (N == 1)
    {
      double *A00 = gsl_matrix_ptr(A, 0, 0);
      *A00 = -(*A00);
      return GSL_SUCCESS;
    }
  else
    {
      int status;
      const size_t N1 = N / 2;
      const size_t N2 = N - N1;

      gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, N1, N1);
      gsl_matrix_view A12 = gsl_matrix_submatrix(A, 0, N1, N1, N2);
      gsl_matrix_view A21 = gsl_matrix_submatrix(A, N1, 0, N2, N1);
      gsl_matrix_view A22 = gsl_matrix_submatrix(A, N1, N1, N2, N2);

      /* A22 = - L22 U22 */
      status = aux_mLU(&A22.matrix);
      if (status)
        return status;

      /* A22 = A22 - L21 U12 */
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &A21.matrix, &A12.matrix, 1.0, &A22.matrix);

      /* A12 - -L11 U12 */
      gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, -1.0, &A11.matrix, &A12.matrix);

      /* A21 = -L21 U11 */
      gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &A11.matrix, &A21.matrix);

      /* A11 = - L11 U11 */
      status = aux_mLU(&A11.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

/* A := A + U B^T where U is upper triangular */
static int
aux_ApUBT(const gsl_matrix * U, const gsl_matrix * B, gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (U->size1 != M || U->size2 != M)
    {
      GSL_ERROR ("U matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (B->size1 != N || B->size2 != M)
    {
      GSL_ERROR ("B matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (M == 1 && N == 1)
    {
      double *aptr = gsl_matrix_ptr(A, 0, 0);
      const double *uptr = gsl_matrix_const_ptr(U, 0, 0);
      const double *bptr = gsl_matrix_const_ptr(B, 0, 0);
      *aptr += (*uptr) * (*bptr);
      return GSL_SUCCESS;
    }
  else if (M == 1)
    {
      double U00 = gsl_matrix_get(U, 0, 0);
      gsl_vector_view v = gsl_matrix_row(A, 0);
      gsl_vector_const_view w = gsl_matrix_const_column(B, 0);
      gsl_blas_daxpy(U00, &w.vector, &v.vector);
      return GSL_SUCCESS;
    }
  else if (N == 1)
    {
      /*
       * partition:
       *
       *         M1  M2
       * B = 1 [ B11 B12 ]
       *
       *       1        M1  M2       1
       * M1 [ A11 ] + [ U11 U12 ] [ B11^T ] M1
       * M2 [ A21 ]   [  0  U22 ] [ B12^T ] M2
       */
      int status;
      const size_t M1 = M / 2;
      const size_t M2 = M - M1;

      gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, M1, 1);
      gsl_matrix_view A21 = gsl_matrix_submatrix(A, M1, 0, M2, 1);
      gsl_vector_view a1 = gsl_matrix_subcolumn(A, 0, 0, M1);

      gsl_matrix_const_view U11 = gsl_matrix_const_submatrix(U, 0, 0, M1, M1);
      gsl_matrix_const_view U12 = gsl_matrix_const_submatrix(U, 0, M1, M1, M2);
      gsl_matrix_const_view U22 = gsl_matrix_const_submatrix(U, M1, M1, M2, M2);

      gsl_matrix_const_view B11 = gsl_matrix_const_submatrix(B, 0, 0, 1, M1);
      gsl_matrix_const_view B12 = gsl_matrix_const_submatrix(B, 0, M1, 1, M2);
      gsl_vector_const_view b2 = gsl_matrix_const_subrow(B, 0, M1, M2);

      /* A(1:M1,1) += U12 * B12^T */
      gsl_blas_dgemv(CblasNoTrans, 1.0, &U12.matrix, &b2.vector, 1.0, &a1.vector);

      /* A11 := A11 + U11 B11^T */
      status = aux_ApUBT(&U11.matrix, &B11.matrix, &A11.matrix);
      if (status)
        return status;

      /* A21 := A21 + U22 B12^T */
      status = aux_ApUBT(&U22.matrix, &B12.matrix, &A21.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
  else
    {
      int status;
      const size_t M1 = M / 2;
      const size_t M2 = M - M1;
      const size_t N1 = N / 2;
      const size_t N2 = N - N1;

      gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, M1, N1);
      gsl_matrix_view A12 = gsl_matrix_submatrix(A, 0, N1, M1, N2);
      gsl_matrix_view A21 = gsl_matrix_submatrix(A, M1, 0, M2, N1);
      gsl_matrix_view A22 = gsl_matrix_submatrix(A, M1, N1, M2, N2);

      gsl_matrix_const_view U11 = gsl_matrix_const_submatrix(U, 0, 0, M1, M1);
      gsl_matrix_const_view U12 = gsl_matrix_const_submatrix(U, 0, M1, M1, M2);
      gsl_matrix_const_view U22 = gsl_matrix_const_submatrix(U, M1, M1, M2, M2);

      gsl_matrix_const_view B11 = gsl_matrix_const_submatrix(B, 0, 0, N1, M1);
      gsl_matrix_const_view B12 = gsl_matrix_const_submatrix(B, 0, M1, N1, M2);
      gsl_matrix_const_view B21 = gsl_matrix_const_submatrix(B, N1, 0, N2, M1);
      gsl_matrix_const_view B22 = gsl_matrix_const_submatrix(B, N1, M1, N2, M2);

      /* A11 := A11 + U11 B11^T */
      status = aux_ApUBT(&U11.matrix, &B11.matrix, &A11.matrix);
      if (status)
        return status;

      /* A11 := A11 + U12 B12^T */
      gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, &U12.matrix, &B12.matrix, 1.0, &A11.matrix);

      /* A12 := A12 + U11 B21^T */
      status = aux_ApUBT(&U11.matrix, &B21.matrix, &A12.matrix);
      if (status)
        return status;

      /* A12 := A12 + U12 B22^T */
      gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, &U12.matrix, &B22.matrix, 1.0, &A12.matrix);

      /* A21 := A21 + U22 B12^T */
      status = aux_ApUBT(&U22.matrix, &B12.matrix, &A21.matrix);
      if (status)
        return status;

      /* A22 := A22 + U22 B22^T */
      status = aux_ApUBT(&U22.matrix, &B22.matrix, &A22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#ifndef HIC_MPLD_POSTDOC_H
#define HIC_MPLD_POSTDOC_H

#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>
#include <map>
#include <math.h>
//#include <mpi.h>
#include <omp.h>
#include <vector>

#include <gsl/gsl_cdf.h>

#include <../IRLS_glm/IRLS.h>

#include <regionDetail2D_pairSite.hpp>
#include <bothEndsMappedFragInfo_withCuttingSite.hpp>

#include <RInside.h>

#define INTTAG 0
#define BOOLTAG 1
#define DIETAG 2

#define NEIGHBDIS 5
#define MORANI 0.001
#define PCUTOFF 0.004
#define PRECISION 100000

//using namespace std;
using std::abs;
using std::basic_string;
using std::cerr;
using std::cin;
using std::cout;
using std::endl;
using std::exception;
using std::flush;
using std::ifstream;
using std::istringstream;
using std::make_pair;
using std::map;
using std::ofstream;
using std::ostringstream;
using std::pair;
using std::sort;
using std::string;
using std::vector;

int postdoc(vector<double>, map<int, int>, const int, const string);
int receiveDomainMap(vector< pair<int, int> >&);
int write_domainFragInteractionMatrix(map< pair< pair< int, pair<int, int> >, pair< int, pair<int, int> > >, vector<fragInfo_bothEndsMapped_withCuttingSite> >&, map<int, int>&, vector< pair<int, int> >&, int, string&);
int getDomainCSinteractionMatrix(map< pair< pair<int, int>, pair<int, int> >, int >&, map< int, vector<int> >&, vector< pair<int, int> >&, int, vector< vector<int> >&, vector< vector< vector<int> > >&, map< int, map<int, int> >&);
inline int convtIndex(int, int, int);
int findIntraDomainInteraction(string, vector< vector<int> >&, vector< vector< vector<int> > >&, map<int, int>&, map< int, map< pair<int, int>, double > >&);
int outputDomainInfo(string, vector< vector<int> >&, vector< vector< vector<int> > >&, map<int, int>&, map< int, map< pair<int, int>, double > >&);
int findInteractingRegion(int, map< int, map< pair<int, int>, double> >&, vector< vector<int> >&, vector<regionDetail_pairSite_2D>&);
int searchNeighbSites(map< pair<int, int>, double >&, map<int, int>&, int, int, int, int, vector< pair< pair<int, int>, double > >&, int&, double&);
int writer_2Dregion(vector<regionDetail_pairSite_2D>&, const string&);
int glm_lasso_fit(vector<double>&, vector<double>&, vector< vector<double> >&);

#endif

ext: h

max_stars_repo_name: Lan-lab/Chrom-Lasso-

max_issues_repo_name: Lan-lab/Chrom-Lasso-

max_forks_repo_name: Lan-lab/Chrom-Lasso-


text:
/* BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 * \file user_blr_model.c
 * \brief allow users to define their own BLR models.
 *
 * User should provide: 1) number of parameters in 1D and 2D;
 *                      2) struct for BLR model;
 *                      3) function for setting parameter range
 *                      4) function for calculating clouds' lag and weight in 1D
 *                      5) function for calculating clouds' lag, velocity, and weight in 2D
 *                      6) function for setting the parameter values used in simulation (if flag_dim<0)
 *
 * The following is an example. The users can make modification with their own BLR models.
 *
 */

#include <gsl/gsl_randist.h>
#include "brains.h"

const int num_params_MyBLRmodel1d = 8;   /*!< number of parameters for 1D model */
const int num_params_MyBLRmodel2d = 17;  /*!< number of parameters for 2D model */

/*!
 *  set parameter range
 */
void set_blr_range_mymodel()
{
  int i;

  i = 0;
  //mu
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 2.0;
  //F
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //inc
  blr_range_model[i][0] = 0.0; //in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in rad
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //gamma
  blr_range_model[i][0] = 1.0;
  blr_range_model[i++][1] = 5.0;
  //xi
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //mbh
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(1.0e3);
  //fellip
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //fflow
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //sigr_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //theta_rot
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 90.0;
  //sig_turb
  blr_range_model[i][0] = log(0.0001);
  blr_range_model[i++][1] = log(0.1);

  return;
}

/*
 * This function caclulate 2d transfer function.
 *
 * The clouds' lag, velocity, and weight are stored in arrays "clouds_tau", "clouds_vel", "clouds_weight",
 * which must be provided.
 */
void gen_cloud_sample_mymodel(const void *pm, int flag_type, int flag_save)
{
  int i, j, nc;
  double r, phi, cos_phi, sin_phi, dis, Lopn_cos;
  double x, y, z, xb, yb, zb, zb0, vx, vy, vz, vxb, vyb, vzb;
  double V, rhoV, theV, Vr, Vph, Vkep, Rs, g;
  double inc, F, beta, mu, k, gam, xi, a, s, sig, rin;
  double mbh, fellip, fflow, sigr_circ, sigthe_circ, sigr_rad, sigthe_rad, theta_rot, sig_turb;
  double Lphi, Lthe, sin_Lphi, cos_Lphi, sin_Lthe, cos_Lthe, sin_inc_cmp, cos_inc_cmp;
  double weight, rnd, rnd_xi;
  MyBLRmodel *model = (MyBLRmodel *)pm;

  Lopn_cos = cos(model->opn*PI/180.0); /* cosine of openning angle */
  inc = acos(model->inc);         /* inclination angle in rad */
  beta = model->beta;
  F = model->F;
  mu = exp(model->mu);                 /* mean radius */
  k = model->k;
  gam = model-> gam;
  xi = model->xi;

  mbh = exp(model->mbh);
  fellip = model->fellip;
  fflow = model->fflow;
  sigr_circ = exp(model->sigr_circ);
  sigthe_circ = exp(model->sigthe_circ);
  sigr_rad = exp(model->sigr_rad);
  sigthe_rad = exp(model->sigthe_rad);
  theta_rot = model->theta_rot*PI/180.0;
  sig_turb = exp(model->sig_turb);

  Rs = 3.0e11*mbh/CM_PER_LD; // Schwarzchild radius in a unit of light-days

  a = 1.0/beta/beta;
  s = mu/a;
  rin=mu*F + Rs;  // include Scharzschild radius
  sig=(1.0-F)*s;

  sin_inc_cmp = cos(inc); //sin(PI/2.0 - inc);
  cos_inc_cmp = sin(inc); //cos(PI/2.0 - inc);

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
// generate a direction of the angular momentum of the orbit
    Lphi = 2.0*PI * gsl_rng_uniform(gsl_r);
    Lthe = acos(Lopn_cos + (1.0-Lopn_cos) * pow(gsl_rng_uniform(gsl_r), gam));
    sin_Lphi = sin(Lphi);
    cos_Lphi = cos(Lphi);
    sin_Lthe = sin(Lthe);
    cos_Lthe = cos(Lthe);

    nc = 0;
    r = rcloud_max_set+1.0;
    while(r>rcloud_max_set || r<rcloud_min_set)
    {
      if(nc > 1000)
      {
        printf("# Error, too many tries in generating ridial location of clouds.\n");
        exit(0);
      }
      rnd = gsl_ran_gamma(gsl_r, a, 1.0);
//    r = mu * F + (1.0-F) * gsl_ran_gamma(gsl_r, 1.0/beta/beta, beta*beta*mu);
      r = rin + sig * rnd;
      nc++;
    }
    phi = 2.0*PI * gsl_rng_uniform(gsl_r);
    cos_phi = cos(phi);
    sin_phi = sin(phi);

    /* Polar coordinates to Cartesian coordinate */
    x = r * cos_phi;
    y = r * sin_phi;
    z = 0.0;

/* right-handed framework
 * first rotate around y axis by an angle of Lthe, then rotate around z axis
 * by an angle of Lphi
 */
  /*xb = cos(Lthe)*cos(Lphi) * x + sin(Lphi) * y - sin(Lthe)*cos(Lphi) * z;
    yb =-cos(Lthe)*sin(Lphi) * x + cos(Lphi) * y + sin(Lthe)*sin(Lphi) * z;
    zb = sin(Lthe) * x + cos(Lthe) * z; */

    xb = cos_Lthe*cos_Lphi * x + sin_Lphi * y;
    yb =-cos_Lthe*sin_Lphi * x + cos_Lphi * y;
    zb = sin_Lthe * x;

    zb0 = zb;
    rnd_xi = gsl_rng_uniform(gsl_r);
    if( (rnd_xi < 1.0 - xi) && zb0 < 0.0)
      zb = -zb;

// counter-rotate around y, LOS is x-axis
    /* x = xb * cos(PI/2.0-inc) + zb * sin(PI/2.0-inc);
       y = yb;
       z =-xb * sin(PI/2.0-inc) + zb * cos(PI/2.0-inc); */

    x = xb * cos_inc_cmp + zb * sin_inc_cmp;
    y = yb;
    z =-xb * sin_inc_cmp + zb * cos_inc_cmp;

    weight = 0.5 + k*(x/r);
    clouds_weight[i] = weight;

#ifndef SA
    dis = r - x;
    clouds_tau[i] = dis;
    if(flag_type == 1)
    {
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
    }
#else
  switch(flag_type)
  {
    case 1:   /* 1D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\n", x, y, z);
      }
      continue;
      break;

    case 2:  /* 2D RM */
      dis = r - x;
      clouds_tau[i] = dis;
      break;

    case 3: /* SA */
      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 4: /* 1D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;

    case 5: /* 2D RM + SA */
      dis = r - x;
      clouds_tau[i] = dis;

      clouds_alpha[i] = y;
      clouds_beta[i] = z;
      break;
  }

#endif

    Vkep = sqrt(mbh/r);

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      rnd = gsl_rng_uniform(gsl_r);

      if(rnd < fellip)
      {
        rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_circ  + 1.0) * Vkep;
        theV =  (gsl_ran_ugaussian(gsl_r) * sigthe_circ + 0.5)*PI;
      }
      else
      {
        if(fflow <= 0.5) /* inflow */
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad   + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad + 1.0) * PI + theta_rot;
        }
        else     /* outflow */
        {
          rhoV = (gsl_ran_ugaussian(gsl_r) * sigr_rad  + 1.0) * Vkep;
          theV = (gsl_ran_ugaussian(gsl_r) * sigthe_rad) * PI + theta_rot;
        }
      }

      Vr = sqrt(2.0) * rhoV * cos(theV);
      Vph = rhoV * sin(theV);

      vx = Vr * cos_phi - Vph * sin_phi;
      vy = Vr * sin_phi + Vph * cos_phi;
      vz = 0.0;

    /*vxb = cos(Lthe)*cos(Lphi) * vx + sin(Lphi) * vy - sin(Lthe)*cos(Lphi) * vz;
      vyb =-cos(Lthe)*sin(Lphi) * vx + cos(Lphi) * vy + sin(Lthe)*sin(Lphi) * vz;
      vzb = sin(Lthe) * vx + cos(Lthe) * vz;*/

      vxb = cos_Lthe*cos_Lphi * vx + sin_Lphi * vy;
      vyb =-cos_Lthe*sin_Lphi * vx + cos_Lphi * vy;
      vzb = sin_Lthe * vx;

      if((rnd_xi < 1.0-xi) && zb0 < 0.0)
        vzb = -vzb;

      /*vx = vxb * cos(PI/2.0-inc) + vzb * sin(PI/2.0-inc);
      vy = vyb;
      vz =-vxb * sin(PI/2.0-inc) + vzb * cos(PI/2.0-inc);*/

      vx = vxb * cos_inc_cmp + vzb * sin_inc_cmp;
      vy = vyb;
      vz =-vxb * sin_inc_cmp + vzb * cos_inc_cmp;

      V = -vx;  //note the definition of the line-of-sight velocity. postive means a receding
                // velocity relative to the observer.

      V += gsl_ran_ugaussian(gsl_r) * sig_turb * Vkep; // add turbulence velocity

      if(fabs(V) >= C_Unit) // make sure that the velocity is smaller than speed of light
        V = 0.9999*C_Unit * (V>0.0?1.0:-1.0);

      g = sqrt( (1.0 + V/C_Unit) / (1.0 - V/C_Unit) ) / sqrt(1.0 - Rs/r); //relativistic effects
      V = (g-1.0)*C_Unit;

      clouds_vel[i*parset.n_vel_per_cloud + j] = V;

      if(flag_save && thistask==roottask)
      {
        if(i%(icr_cloud_save) == 0)
          fprintf(fcloud_out, "%f\t%f\t%f\t%f\t%f\t%f\t%f\n", x, y, z, vx*VelUnit, vy*VelUnit, vz*VelUnit, weight);
      }
    }
  }

  return;
}

/*!
 *  set parameter values used in simulation (if flag_dim < 0).
 */
void set_par_value_mymodel_sim(double *pm)
{
  int i;

	i=0;
  pm[i++] = log(4.0);   // mu
  pm[i++] = 1.0;        // beta
  pm[i++] = 0.25;        // F
  pm[i++] = cos(20.0/180.0*PI); // inc
  pm[i++] = 40.0;       // opn
  pm[i++] = -0.4;        // kappa
  pm[i++] = 5.0;        // gamma
  pm[i++] = 0.5;        // obscuration
  pm[i++] = log(2.0);  //mbh
  pm[i++] = 0.5;       //fellip
  pm[i++] = 0.4;       //fflow
  pm[i++] = log(0.01); //
  pm[i++] = log(0.1);  //
  pm[i++] = log(0.01); //
  pm[i++] = log(0.1);  //
  pm[i++] = 0.0;       // theta_rot
  pm[i++] = log(0.001);  // sig_turb
  pm[i++] = 0.0;       // parameter for spectral broadening
}

#ifdef SA
/*!
 *  set parameter range
 */
void set_sa_blr_range_mymodel()
{
  int i;

  i = 0;
  //mu
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  sa_blr_range_model[i][0] = 0.001;
  sa_blr_range_model[i++][1] = 2.0;
  //F
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //inc
  sa_blr_range_model[i][0] = 0.0; //in cosine
  sa_blr_range_model[i++][1] = 1.0;
  //opn
  sa_blr_range_model[i][0] = 0.0;  // in rad
  sa_blr_range_model[i++][1] = 90.0;
  //k
  sa_blr_range_model[i][0] = -0.5;
  sa_blr_range_model[i++][1] = 0.5;
  //gamma
  sa_blr_range_model[i][0] = 1.0;
  sa_blr_range_model[i++][1] = 5.0;
  //xi
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;

  //mbh
  sa_blr_range_model[i][0] = log(0.1);
  sa_blr_range_model[i++][1] = log(1.0e3);
  //fellip
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //fflow
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 1.0;
  //sigr_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  sa_blr_range_model[i][0] = log(0.001);
  sa_blr_range_model[i++][1] = log(1.0);
  //theta_rot
  sa_blr_range_model[i][0] = 0.0;
  sa_blr_range_model[i++][1] = 90.0;
  //sig_turb
  sa_blr_range_model[i][0] = log(0.0001);
  sa_blr_range_model[i++][1] = log(0.1);

  return;
}

/*!
 *  set parameter values used in simulation (if flag_dim < 0).
 */
void set_sa_par_value_mymodel_sim(double *pm)
{
  int i;

	i=0;
  pm[i++] = log(4.0);   // mu
  pm[i++] = 1.0;        // beta
  pm[i++] = 0.25;        // F
  pm[i++] = cos(20.0/180.0*PI); // inc
  pm[i++] = 40.0;       // opn
  pm[i++] = -0.4;        // kappa
  pm[i++] = 5.0;        // gamma
  pm[i++] = 0.5;        // obscuration
  pm[i++] = log(2.0);  //mbh
  pm[i++] = 0.5;       //fellip
  pm[i++] = 0.4;       //fflow
  pm[i++] = log(0.01); //
  pm[i++] = log(0.1);  //
  pm[i++] = log(0.01); //
  pm[i++] = log(0.1);  //
  pm[i++] = 0.0;       // theta_rot
  pm[i++] = log(0.001);  // sig_turb
  pm[i++] = 0.0;       // parameter for spectral broadening
}
#endif

ext: c

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
#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 <fftw3.h> /* Note: when included AFTER complex.h, fftw_complex type defaults to native double complex */

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

/* Window functions */
double WindowFunction(double x, double xi, double xf, double deltaxi, double deltaxf)
{
  double di =  deltaxi/(20*log(10));
  double df =  deltaxf/(20*log(10));
  if (x <= xi + di) return 0;
  else if (xi +di < x && x < xi + deltaxi - di) {
    return 1./(1 + exp(deltaxi/(x - xi) + deltaxi/(x - (xi + deltaxi))));
  }
  else if (xi + deltaxi - di <= x && x <= xf - deltaxf + df) return 1.;
  else if (xf - deltaxf + df < x && x < xf - df) {
    return 1./(1 + exp(-(deltaxf/(x - (xf - deltaxf))) - deltaxf/(x - xf)));
  }
  else return 0;
}
double WindowFunctionLeft(double x, double xf, double deltaxf)
{
  double df =  deltaxf/(20*log(10));
  if (x <= xf - deltaxf + df) return 1;
  else if (xf - deltaxf + df < x && x < xf - df) {
    return 1./(1 + exp(-deltaxf/(x - (xf - deltaxf)) - deltaxf/(x - xf)));
  }
  else return 0;
}
double WindowFunctionRight(double x, double xi, double deltaxi)
{
  double di =  deltaxi/(20*log(10));
  if (x <= xi + di) return 0;
  else if (xi + di < x && x < xi + deltaxi - di) {
    return 1./(1 + exp(deltaxi/(x - xi) + deltaxi/(x - (xi + deltaxi))));
  }
  else return 1;
}

/* FFT of real time series */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
int FFTRealTimeSeries(ReImFrequencySeries** freqseries, RealTimeSeries* timeseries, double twindowbeg, double twindowend, int nzeropad)
{

  /* deltat of time series */
  /* Warning: assumes linear sampling in time */
  double* times = timeseries->times->data;
  double deltat = times[1] - times[0];
  double tshift = times[0]; /* time shift to be re-applied later */

  /* Initialize vector for windowed, 0-padded FFT input */
  int n = (int) timeseries->times->size;
  int nzeros = (int) pow(2, ((int) ceil(log(n)/log(2))) + nzeropad) - n; /* Here defined with ceil, but with floor in IFFT */
  gsl_vector* hvalues = gsl_vector_alloc(n + nzeros);
  gsl_vector_set_zero(hvalues);

  /* Compute input TD values, with windowing */
  int nbptswindowbeg = (int) ceil(twindowbeg/deltat) + 1;
  int nbptswindowend = (int) ceil(twindowend/deltat) + 1;
  double t1windowbeg = times[0];
  double t2windowbeg = times[nbptswindowbeg-1];
  double t1windowend = times[n-nbptswindowend];
  double t2windowend = times[n-1];
  double deltatwindowbeg = t2windowbeg - t1windowbeg;
  double deltatwindowend = t2windowend - t1windowend;
  double* h = timeseries->h->data;
  double* hval = hvalues->data;

  for (int i=0; i<nbptswindowbeg; i++) {
    hval[i] = WindowFunctionRight(times[i], t1windowbeg, deltatwindowbeg) * h[i];
  }
  for (int i=nbptswindowbeg; i<n-nbptswindowend; i++) {
    hval[i] =  h[i];
  }
  for (int i=n-nbptswindowend; i<n; i++) {
    hval[i] = WindowFunctionLeft(times[i], t2windowend, deltatwindowend) * h[i];
  }

  /* FFT - uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
  int N = (int) hvalues->size;
  fftw_plan p;
  double* in = hvalues->data;
  fftw_complex* out;
  out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2+1)); /* Note: N/2+1 elements */
  p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
  fftw_execute(p);

  /* Initialize output structure */
  ReImFrequencySeries_Init(freqseries, N/2); /* Note: N/2+1 elements as output of fftw, we drop the last one (at Nyquist frequency) */

  /* Extracting and converting data from FFTW output */
  double deltaf = 1./(N*deltat);
  double* freq = (*freqseries)->freq->data;
  double* hreal = (*freqseries)->h_real->data;
  double* himag = (*freqseries)->h_imag->data;
  double f;
  double complex hcomplex;
  double factorshift = 2*PI*tshift; /* Reminder: Flipped sign convention */
  for(int i=0; i<N/2; i++) {
    f = i*deltaf;
    freq[i] = f;
    hcomplex = deltat * conj(out[i]) * cexp(I*factorshift*f); /* Note that we convert here the FT sign convention */
    hreal[i] = creal(hcomplex);
    himag[i] = cimag(hcomplex);
  }

  /* Clean up */
  gsl_vector_free(hvalues);
  fftw_destroy_plan(p);
  fftw_free(out);

  return SUCCESS;
}

/* FFT of Re/Im time series */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
/* NOTE: only keeps positive frequencies in the output */
int FFTTimeSeries(ReImFrequencySeries** freqseries, ReImTimeSeries* timeseries, double twindowbeg, double twindowend, int nzeropad)
{
  /* deltat of time series */
  /* Warning: assumes linear sampling in time */
  double* times = timeseries->times->data;
  double deltat = times[1] - times[0];
  double tshift = times[0]; /* time shift to be re-applied later */

  /* Initialize vector for windowed, 0-padded FFT input */
  int n = (int) timeseries->times->size;
  int nzeros = (int) pow(2, ((int) ceil(log(n)/log(2))) + nzeropad) - n; /* Here defined with ceil, but with floor in IFFT */
  int N = n + nzeros;

  /* Compute input TD values, with windowing */
  int nbptswindowbeg = (int) ceil(twindowbeg/deltat) + 1;
  int nbptswindowend = (int) ceil(twindowend/deltat) + 1;
  double t1windowbeg = times[0];
  double t2windowbeg = times[nbptswindowbeg-1];
  double t1windowend = times[n-nbptswindowend];
  double t2windowend = times[n-1];
  double deltatwindowbeg = t2windowbeg - t1windowbeg;
  double deltatwindowend = t2windowend - t1windowend;
  double* htdreal = timeseries->h_real->data;
  double* htdimag = timeseries->h_imag->data;

  fftw_complex* in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  for (int i=0; i<nbptswindowbeg; i++) {
    in[i] = WindowFunctionRight(times[i], t1windowbeg, deltatwindowbeg) * (htdreal[i] + I*htdimag[i]);
  }
  for (int i=nbptswindowbeg; i<n-nbptswindowend; i++) {
    in[i] =  (htdreal[i] + I*htdimag[i]);
  }
  for (int i=n-nbptswindowend; i<n; i++) {
    in[i] = WindowFunctionLeft(times[i], t2windowend, deltatwindowend) * (htdreal[i] + I*htdimag[i]);
  }

  /* FFT - uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
  /* Represented here by the use of FFTW_BACKWARD (plus sign in the exp) */
  fftw_plan p;
  fftw_complex* out;
  out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); /* Note: N/2+1 elements */
  p = fftw_plan_dft_1d(N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(p);

  /* Initialize output structure */
  ReImFrequencySeries_Init(freqseries, N/2); /* NOTE: N/2 first elements of output of fftw are positive freqs, we eliminate negative frequency (the second half of the series) */

  /* Extracting and converting data from FFTW output */
  double deltaf = 1./(N*deltat);
  double* freq = (*freqseries)->freq->data;
  double* hreal = (*freqseries)->h_real->data;
  double* himag = (*freqseries)->h_imag->data;
  double f;
  double complex hcomplex;
  double factorshift = 2*PI*tshift; /* Reminder: Flipped sign convention */
  for(int i=0; i<N/2; i++) {
    f = i*deltaf;
    freq[i] = f;
    hcomplex = deltat * out[i] * cexp(I*factorshift*f); /* Note the FT sign convention  */
    hreal[i] = creal(hcomplex);
    himag[i] = cimag(hcomplex);
  }

  /* Clean up */
  fftw_destroy_plan(p);
  fftw_free(in);
  fftw_free(out);

  return SUCCESS;
}

/* IFFT of frequency series */
/* Note: assumes frequency series is FT of real data - produces real time series */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
int IFFTFrequencySeriesReal(RealTimeSeries** timeseries, ReImFrequencySeries* freqseries, double f1windowbeg, double f2windowbeg, double f1windowend, double f2windowend, int nzeropad)
{
  /* Checking sanity of windowing frequencies */
  if(!((f1windowbeg<f2windowbeg)&&(f1windowend<f2windowend)&&(f2windowbeg<f1windowend))) {
    printf("Error in IFFTFrequencySeriesReal: inconsistent windowing frequencies.\n");
    printf("(f1windowbeg, f2windowbeg, f1windowend, f2windowend) = (%g, %g, %g, %g)\n", f1windowbeg, f2windowbeg, f1windowend, f2windowend);
  }

  /* deltaf of frequency series */
  /* Warning: assumes linear sampling in frequency */
  double* freq = freqseries->freq->data;
  double deltaf = freq[1] - freq[0];

  /* Initialize vector for windowed, 0-padded FFT input */
  int n = (int) freqseries->freq->size;
  while(freq[n-1] > f2windowend) n--;
  int nzeros = (int) pow(2, ((int) ceil(log(n)/log(2))) + nzeropad) - n; /* Here defined with floor, but with ceil in FFT */
  gsl_vector* hrealvalues = gsl_vector_alloc(n + nzeros);
  gsl_vector* himagvalues = gsl_vector_alloc(n + nzeros);
  gsl_vector_set_zero(hrealvalues);
  gsl_vector_set_zero(himagvalues);

  /* Compute input FD values, with windowing */
  double deltafwindowbeg = f2windowbeg - f1windowbeg;
  double deltafwindowend = f2windowend - f1windowend;
  double* hreal = freqseries->h_real->data;
  double* himag = freqseries->h_imag->data;
  double* hrealval = hrealvalues->data;
  double* himagval = himagvalues->data;
  for (int i=0; i<n; i++) {
    double window = WindowFunction(freq[i], f1windowbeg, f2windowend, deltafwindowbeg, deltafwindowend);
    hrealval[i] = window * hreal[i];
    himagval[i] = window * himag[i];
  }

  /* Input as array of fftw_complex */
  int N = (int) hrealvalues->size;
  fftw_complex* in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  /* NOTE: Restoring the standard sign convention for the FT, used by FFTW - change in convention amounts to f->-f, equivalent to a conjugation only for FFT of a real series  */
  for(int i=0; i<n; i++) {
    in[i] = hrealval[i] - I*himagval[i] ;
  }
  for(int i=n; i<N; i++) {
    in[i] = 0;
  }

  /* FFT - uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
  fftw_plan p;
  double* out = fftw_malloc(sizeof(double) * N);
  p = fftw_plan_dft_c2r_1d(N, in, out, FFTW_ESTIMATE);
  fftw_execute(p);

  /* Initialize output structure */
  RealTimeSeries_Init(timeseries, N);

  /* Extracting and converting data from FFTW output - moving negative times to the left */
  double deltat = 1./(N*deltaf);
  double fac = 1./(N*deltat); /* Additional 1/N to put the FFTW convention in agreement with numpy IFFT */
  double* times = (*timeseries)->times->data;
  double* h = (*timeseries)->h->data;
  for(int i=0; i<N/2; i++) {
    times[i] = (i-N/2)*deltat;
    times[N/2+i] = i*deltat;
    h[i] = fac * out[N/2+i];
    h[N/2+i] = fac * out[i];
  }

  /* Clean up */
  gsl_vector_free(hrealvalues);
  gsl_vector_free(himagvalues);
  fftw_destroy_plan(p);
  fftw_free(in);
  fftw_free(out);

  return SUCCESS;
}

/* IFFT of frequency series */
/* Note: does not assume frequency series is FT of real data - produces complex time series */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
int IFFTFrequencySeries(ReImTimeSeries** timeseries, ReImFrequencySeries* freqseries, double f1windowbeg, double f2windowbeg, double f1windowend, double f2windowend, int nzeropad)
{
  /* Checking sanity of windowing frequencies */
  if(!((f1windowbeg<f2windowbeg)&&(f1windowend<f2windowend)&&(f2windowbeg<f1windowend))) {
    printf("Error in IFFTFrequencySeries: inconsistent windowing frequencies.\n");
    printf("(f1windowbeg, f2windowbeg, f1windowend, f2windowend) = (%g, %g, %g, %g)\n", f1windowbeg, f2windowbeg, f1windowend, f2windowend);
  }

  /* deltaf of frequency series */
  /* Warning: assumes linear sampling in frequency */
  double* freq = freqseries->freq->data;
  double deltaf = freq[1] - freq[0];

  /* Initialize vector for windowed, 0-padded FFT input */
  int n = (int) freqseries->freq->size;
  while(freq[n-1] > f2windowend) n--;
  int nzeros = (int) pow(2, ((int) ceil(log(n)/log(2))) + nzeropad) - n; /* Here defined with floor, but with ceil in FFT */
  gsl_vector* hrealvalues = gsl_vector_alloc(n + nzeros);
  gsl_vector* himagvalues = gsl_vector_alloc(n + nzeros);
  gsl_vector_set_zero(hrealvalues);
  gsl_vector_set_zero(himagvalues);

  /* Compute input FD values, with windowing */
  double deltafwindowbeg = f2windowbeg - f1windowbeg;
  double deltafwindowend = f2windowend - f1windowend;
  double* hreal = freqseries->h_real->data;
  double* himag = freqseries->h_imag->data;
  double* hrealval = hrealvalues->data;
  double* himagval = himagvalues->data;
  for (int i=0; i<n; i++) {
    double window = WindowFunction(freq[i], f1windowbeg, f2windowend, deltafwindowbeg, deltafwindowend);
    hrealval[i] = window * hreal[i];
    himagval[i] = window * himag[i];
  }

  /* Input as array of fftw_complex */
  int N = (int) hrealvalues->size;
  fftw_complex* in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  /* NOTE: the sign convention for the FT used by FFTW is different - change in convention amounts to f->-f, equivalent to a conjugation only for FFT of a real series - for a FFT of a complex series, we do not conjugate and we keep our convention */
  for(int i=0; i<n; i++) {
    in[i] = hrealval[i] + I*himagval[i] ;
  }
  for(int i=n; i<N; i++) {
    in[i] = 0;
  }

  /* FFT - FFTW uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) - our convention is h(f) = int e^(-2ipift)h(t) */
  /* NOTE: due to the difference in convention for the FT, we use the FFTW_FORWARD sign (sign - in the exponential) */
  fftw_plan p;
  double complex* out = fftw_malloc(sizeof(fftw_complex) * N);
  p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_execute(p);

  /* Initialize output structure */
  ReImTimeSeries_Init(timeseries, N);

  /* Extracting and converting data from FFTW output - moving negative times to the left */
  double deltat = 1./(N*deltaf);
  double fac = 1./(N*deltat); /* Additional 1/N to put the FFTW convention in agreement with numpy IFFT */
  double* times = (*timeseries)->times->data;
  double* htdreal = (*timeseries)->h_real->data;
  double* htdimag = (*timeseries)->h_imag->data;
  for(int i=0; i<N/2; i++) {
    times[i] = (i-N/2)*deltat;
    times[N/2+i] = i*deltat;
    htdreal[i] = fac * creal(out[N/2+i]);
    htdimag[i] = fac * cimag(out[N/2+i]);
    htdreal[N/2+i] = fac * creal(out[i]);
    htdimag[N/2+i] = fac * cimag(out[i]);
  }

  /* Clean up */
  gsl_vector_free(hrealvalues);
  gsl_vector_free(himagvalues);
  fftw_destroy_plan(p);
  fftw_free(in);
  fftw_free(out);

  return SUCCESS;
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#ifndef __GSLEXTRA_H__
#define __GSLEXTRA_H__

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>

#endif

void gsl_vector_complex_convert(gsl_vector * source, gsl_vector_complex * target, int length);
void gsl_matrix_complex_convert(gsl_matrix * source, gsl_matrix_complex * target, int rows, int columns);
void gsl_vector_complex_extract(gsl_vector_complex * source, gsl_vector * real, gsl_vector * imag, int length);
void gsl_matrix_complex_extract(gsl_vector_complex * source, gsl_matrix * real,gsl_matrix * imag, int rows, int columns);
void gsl_vector_complex_combine(gsl_vector * real, gsl_vector * imag, gsl_vector_complex * target);
void gsl_matrix_complex_combine(gsl_matrix * real,gsl_matrix * imag, gsl_matrix_complex * target);
void gsl_matrix_diag(gsl_matrix * target, gsl_vector * diag, int length);
void gsl_matrix_complex_diag(gsl_matrix_complex * target, gsl_vector_complex * diag, int length);
void gsl_matrix_mul(gsl_matrix * A, gsl_matrix *B, gsl_matrix * Result,int Acolumn,int Arow,int Bcolumn);
void gsl_matrix_complex_mul(gsl_matrix_complex * A, gsl_matrix_complex *B, gsl_matrix_complex * Result,int Acolumn,int Arow,int Bcolumn);
double gsl_vector_inner_product(gsl_vector * A, gsl_vector * B,int length);
gsl_complex gsl_vector_complex_product(gsl_vector_complex * A, gsl_vector_complex * B, int length);
gsl_complex gsl_vector_complex_inner_product(gsl_vector_complex * A, gsl_vector_complex * B, int length);
void gsl_vector_transform(gsl_vector * vec,gsl_matrix * trf,int length);
void gsl_vector_complex_transform(gsl_vector_complex * vec, gsl_matrix_complex * trf,int length);
void gsl_matrix_unitmatrix(gsl_matrix * m,int length);
void gsl_matrix_complex_unitmatrix(gsl_matrix_complex * m,int length);
void gsl_vector_complex_conjugate(gsl_vector_complex * v, int length);
void gsl_matrix_complex_conjugate(gsl_matrix_complex * m, int rows, int columns);

max_stars_repo_name: Walter-Feng/HomebrewLib

max_issues_repo_name: Walter-Feng/HomebrewLib

max_forks_repo_name: Walter-Feng/HomebrewLib


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "DirectionalLight.h"

namespace Library
{
	class Texture2D;
	class ProxyModel;
}

namespace Rendering
{
	class DisplacementMappingMaterial;

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

		bool RealDisplacementMapEnabled() const;
		void SetRealDisplacementMapEnabled(bool enabled);
		void ToggleRealDisplacementMap();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float DirectionalLightIntensity() const;
		void SetDirectionalLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightDirection() const;
		void RotateDirectionalLight(DirectX::XMFLOAT2 amount);

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

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

	private:
		std::shared_ptr<DisplacementMappingMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		Library::DirectionalLight mDirectionalLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		std::shared_ptr<Library::Texture2D> mRealDisplacementMap;
		std::shared_ptr<Library::Texture2D> mDefaultDisplacementMap;
		bool mUpdateMaterial{ true };
		bool mRealDisplacementMapEnabled{ 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/Exception.h"
#include "ioda/Types/Type_Provider.h"
#include "ioda/defs.h"

namespace ioda {
class Group;
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_
};

/// \brief Data Types can be grouped into a few categories. These are the categories.
/// \note  Not all backends implement all types.
enum class TypeClass {
  Unknown,    ///< Unsupported / unhandled type
  Integer,    ///< All integer types
  Float,      ///< All floating-point types
  String,     ///< All string types (fixed-length, variable, ASCII, UTF-8)
  Bitfield,   ///< All bit fields
  Opaque,     ///< All binary blobs
  Compound,   ///< All compound types (types with member elements)
  Reference,  ///< All object references
  Enum,       ///< All enumerated types
  VlenArray,  ///< All variable-length array types (not strings)
  FixedArray  ///< All fixed-length array types
};

namespace detail {
/// \brief Convenience function to safely copy a string.
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_; }

public:
  virtual ~Type_Base() {}
  std::shared_ptr<Type_Backend> getBackend() const { return backend_; }
  bool isValid() const { return (backend_.use_count() > 0); }

  /// @}
  /// @name General functions
  /// @{

  /// \brief Get the size 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;

  /// \brief Does this type represent a string, an integer, a float, an array, an
  ///   enumeration, a bitset, or any other type?
  virtual TypeClass getClass() const;

  /// \brief Save (commit) the type to a backend
  /// \details From the HDF5 docs:
  ///   Committed datatypes can be used to save space in a file where many
  ///   datasets or attributes use the same datatype or to avoid defining a complex
  ///   compound datatype more than once. Committed datatypes can also be used to ensure
  ///   that multiple instances of the same datatype are truly identical.
  ///
  ///   This is used extensively for enumerated types.
  virtual void commitToBackend(Group& d, const std::string& name) const;

  /// @}
  /// @name Numeric type functions
  /// @{

  /// \brief Is this type signed or unsigned?
  /// \returns true if signed, false if unsigned.
  /// \throws if the type is not a numeric type (i.e. not a simple integer or float).
  virtual bool isTypeSigned() const;

  /// @}
  /// @name String type functions
  /// @{

  /// \brief Is this a variable-length string type?
  /// \returns true if a variable-length string type, false if not.
  ///   False can imply either that the type is a fixed-length string type (if getClass() == TypeClass::String),
  ///   or that the type is not a string type at all.
  virtual bool isVariableLengthStringType() const;

  /// \brief Get the character set of this string type.
  /// \returns Ascii or Unicode.
  /// \throws ioda::Exception on error, or if the type is not a string type.
  /// \note Currently, there is no way to set the character set. Everything is
  ///   assumed to be a UTF-8 string in IODA.
  virtual StringCSet getStringCSet() const;

  /// @}
  /// @name Array type functions
  /// @{

  /// \brief Get the "base" type of an object. For an array, this is the type of
  ///   the array's elements. I.e. an array of int32_t has a base type of int32_t.
  ///   For an enumerated type, this is the type used for the enumeration.
  virtual Type_Implementation getBaseType() const;

  /// \brief Get the dimensions of an array type.
  /// \returns A vector of the dimensions. vector::size() is the rank (dimensionality).
  virtual std::vector<Dimensions_t> getDimensions() 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();
  StringCSet getStringCSet() 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 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
/// \todo Once C++20 support is added, make a distinction between std::string and std::u8string.
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(typeid(DataType), String_Type_Length);
}

// 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: NOAA-EMC/ioda

max_issues_repo_name: NOAA-EMC/ioda

max_forks_repo_name: NOAA-EMC/ioda


text:
/* block/gsl_block_long_double.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_BLOCK_LONG_DOUBLE_H__
#define __GSL_BLOCK_LONG_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_types.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

struct gsl_block_long_double_struct
{
  size_t size;
  long double *data;
};

typedef struct gsl_block_long_double_struct gsl_block_long_double;

GSL_EXPORT gsl_block_long_double *gsl_block_long_double_alloc (const size_t n);
GSL_EXPORT gsl_block_long_double *gsl_block_long_double_calloc (const size_t n);
GSL_EXPORT void gsl_block_long_double_free (gsl_block_long_double * b);

GSL_EXPORT int gsl_block_long_double_fread (FILE * stream, gsl_block_long_double * b);
GSL_EXPORT int gsl_block_long_double_fwrite (FILE * stream, const gsl_block_long_double * b);
GSL_EXPORT int gsl_block_long_double_fscanf (FILE * stream, gsl_block_long_double * b);
GSL_EXPORT int gsl_block_long_double_fprintf (FILE * stream, const gsl_block_long_double * b, const char *format);

GSL_EXPORT int gsl_block_long_double_raw_fread (FILE * stream, long double * b, const size_t n, const size_t stride);
GSL_EXPORT int gsl_block_long_double_raw_fwrite (FILE * stream, const long double * b, const size_t n, const size_t stride);
GSL_EXPORT int gsl_block_long_double_raw_fscanf (FILE * stream, long double * b, const size_t n, const size_t stride);
GSL_EXPORT int gsl_block_long_double_raw_fprintf (FILE * stream, const long double * b, const size_t n, const size_t stride, const char *format);

GSL_EXPORT size_t gsl_block_long_double_size (const gsl_block_long_double * b);
GSL_EXPORT long double * gsl_block_long_double_data (const gsl_block_long_double * b);

__END_DECLS

#endif /* __GSL_BLOCK_LONG_DOUBLE_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/*
    KGSX: Biomolecular Kino-geometric Sampling and Fitting of Experimental Data
    Yao et al, Proteins. 2012 Jan;80(1):25-43
    e-mail: latombe@cs.stanford.edu, vdbedem@slac.stanford.edu, julie.bernauer@inria.fr

        Copyright (C) 2011-2013 Stanford University

        Permission is hereby granted, free of charge, to any person obtaining a copy of
        this software and associated documentation files (the "Software"), to deal in
        the Software without restriction, including without limitation the rights to
        use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
        of the Software, and to permit persons to whom the Software is furnished to do
        so, subject to the following conditions:

        This entire text, including the above copyright notice and this permission notice
        shall be included in all copies or substantial portions of the Software.

        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
        OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
        FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
        IN THE SOFTWARE.
*/

#ifndef KGS_ATOMPAIRDISTANCEDIRECTION_H
#define KGS_ATOMPAIRDISTANCEDIRECTION_H

#include <gsl/gsl_vector.h>

#include <tuple>
#include <vector>

#include "core/graph/KinTree.h"
#include "directions/Direction.h"
#include "Selection.h"

/**
 * A direction where the gradient seeks to drive configurations toward a desired distance between a particular pair of
 * atoms. This class is conceptually similar to `RelativeMSDDirection` and to a lesser degree `LSNrelativeDirection`,
 * but only takes a single pair of atoms and desired distance.
 */
class AtomPairDistanceDirection: public Direction {
 public:
  AtomPairDistanceDirection(const std::tuple<Atom*, Atom*, double> &relativeDistances);

 protected:
  void computeGradient(Configuration* conf, Configuration* target, gsl_vector* ret);

 private:
  const std::tuple<Atom*, Atom*, double> &m_relativeDistance;
};

#endif //KGS_ATOMPAIRDISTANCEDIRECTION_H


text:
// MIT License Copyright (c) 2020 Jarrett Wendt

#pragma once

#define _SILENCE_CLANG_COROUTINE_MESSAGE

// Standard
#include <algorithm>
#include <cassert>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <exception>
#include <filesystem>
#include <fstream>
#include <functional>
#include <future>
#include <initializer_list>
#include <iostream>
#include <new>
#include <random>
#include <sstream>
#include <stdexcept>
#include <string>
#include <thread>
#include <tuple>
#include <utility>
#include <vector>

#ifdef _WIN32
#include <experimental/coroutine>
#else
#include <coroutine>
#endif

// External
#include <gsl/gsl>

max_stars_repo_name: JarrettWendt/FIEAEngine

max_issues_repo_name: JarrettWendt/FIEAEngine

max_forks_repo_name: JarrettWendt/FIEAEngine


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 __BSSCR_MG_h__
#define __BSSCR_MG_h__

#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscpc.h>
#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>
#include <StgFEM/StgFEM.h>
#include <PICellerator/PICellerator.h>
#include <Underworld/Underworld.h>
#include "Solvers/KSPSolvers/KSPSolvers.h"

#include "common-driver-utils.h"
#include "BSSCR.h" /* includes StokesBlockKSPInterface.h */

#define KSPBSSCR         "bsscr"

typedef struct {
  KSP ksp;
  PC pc;

  PetscTruth useAcceleratingSmoothingMG;
  PetscTruth acceleratingSmoothingMGView;

  /* mg_accelerating_smoothing options */

  PetscInt smoothsMax;
  PetscInt smoothsToStartWith;
  PetscInt currentNumberOfSmooths;
  PetscInt smoothingIncrement;
  PetscInt targetCyclesForTenfoldReduction;

  PetscInt smoothingCountThisSolve;
  PetscInt totalSmoothingCount;
  PetscInt totalMgCycleCount;
} MGContext;

PetscErrorCode KSPCycleEffectivenessMonitorAndAdjust(KSP ksp, PetscInt n, PetscReal rnorm, void *_mgctx );
//PetscErrorCode MG_inner_solver_mgContext_initialise(MGContext *mgCtx);
PetscErrorCode MG_inner_solver_pcmg_shutdown( PC pc_MG );
//PetscErrorCode BSSCR_mgPCApply( void *ctx, Vec x, Vec y );
//PetscErrorCode BSSCR_mgPCAccelerating( void *ctx, Vec x, Vec y );
double setupMG( KSP_BSSCR * bsscrp_self, KSP ksp_inner, PC pc_MG, Mat K, MGContext *mgCtx );
#endif


text:
#include <stdlib.h>
#include <assert.h>
#include <cconfigspace.h>
#include <gsl/gsl_rng.h>

static void test_rng_create_with_type() {
	const gsl_rng_type **t1, **t2, *t;
	size_t               type_count = 0;
	int32_t              selected, refcount;
	ccs_rng_t            rng = NULL;
	ccs_result_t         err = CCS_SUCCESS;
	ccs_object_type_t    otype;

	t1 = t2 = gsl_rng_types_setup();
	while (*t1++)
		type_count++;
	selected = rand() % type_count;
	err = ccs_rng_create_with_type(t2[selected], NULL);
	assert( err == -CCS_INVALID_VALUE );
	err = ccs_rng_create_with_type(NULL, &rng);
	assert( err == -CCS_INVALID_VALUE );
	err = ccs_rng_create_with_type(t2[selected], &rng);
	assert( err == CCS_SUCCESS );
	assert( rng );
	err = ccs_rng_get_type(rng, &t);
	assert( err == CCS_SUCCESS );
	assert( t == t2[selected] );
	err = ccs_object_get_type(rng, &otype);
	assert( err == CCS_SUCCESS );
	assert( otype == CCS_RNG );
	err = ccs_object_get_refcount(rng, &refcount);
	assert( err == CCS_SUCCESS );
	assert( refcount == 1 );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_rng_create() {
	ccs_rng_t           rng = NULL;
	ccs_result_t        err = CCS_SUCCESS;
	const gsl_rng_type *t;

	err = ccs_rng_create(NULL);
	assert( err == -CCS_INVALID_VALUE );
	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	assert( rng );
	err = ccs_rng_get_type(rng, &t);
	assert( err == CCS_SUCCESS );
	assert( t == gsl_rng_default );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_rng_min_max() {
	ccs_rng_t         rng = NULL;
	ccs_result_t      err = CCS_SUCCESS;
	unsigned long int imin = 0;
	unsigned long int imax = 0;
	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_rng_min(NULL, &imin);
	assert( err == -CCS_INVALID_OBJECT );
	err = ccs_rng_min(rng, NULL);
	assert( err == -CCS_INVALID_VALUE );
	err = ccs_rng_min(rng, &imin);
	assert( err == CCS_SUCCESS );
	err = ccs_rng_max(NULL, &imax);
	assert( err == -CCS_INVALID_OBJECT );
	err = ccs_rng_max(rng, NULL);
	assert( err == -CCS_INVALID_VALUE );
	err = ccs_rng_max(rng, &imax);
	assert( err == CCS_SUCCESS );
	assert( imin < imax );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_rng_get() {
	ccs_rng_t         rng = NULL;
	ccs_result_t      err = CCS_SUCCESS;
	unsigned long int i = 0;
	unsigned long int imin = 0;
	unsigned long int imax = 0;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_rng_min(rng, &imin);
	assert( err == CCS_SUCCESS );
	err = ccs_rng_max(rng, &imax);
	assert( err == CCS_SUCCESS );
	for (int j = 0; j < 100; j++) {
		err = ccs_rng_get(rng, &i);
		assert( err == CCS_SUCCESS );
		assert( i >= imin );
		assert( i <= imax );
	}
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_rng_uniform() {
	ccs_rng_t    rng = NULL;
	ccs_result_t err = CCS_SUCCESS;
	double       d = -1.0;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	for (int j = 0; j < 100; j++) {
		err = ccs_rng_uniform(rng, &d);
		assert( err == CCS_SUCCESS );
		assert( d >= 0.0 );
		assert( d <  1.0 );
	}
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

int main() {
	ccs_init();
	test_rng_create_with_type();
	test_rng_create();
	test_rng_min_max();
	test_rng_get();
	test_rng_uniform();
	ccs_fini();
	return 0;
}

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS