text:
/* stable/stable_dist.c
 *
 * Main Libstable source file. Definition of the StableDist structures
 * and auxiliary functions to manage alpha-stable distributions.
 *
 * Copyright (C) 2013. Javier Royuela del Val
 *                     Federico Simmross Wattenberg
 *
 * 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; version 3 of the License.
 *
 * 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, see <http://www.gnu.org/licenses/>.
 *
 *
 *  Javier Royuela del Val.
 *  E.T.S.I. Telecomunicación
 *  Universidad de Valladolid
 *  Paseo de Belén 15, 47002 Valladolid, Spain.
 *  jroyval@lpi.tel.uva.es
 */
#include <gsl/gsl_errno.h>

#include "stable.h"

#include <pthread.h>
#include <math.h>
#include <unistd.h>

/*----------------------------------------------------------------------------*/
/*                             Public part                                    */
/*----------------------------------------------------------------------------*/

unsigned int stable_get_THREADS() { return THREADS; }
void stable_set_THREADS(unsigned int value)
{
  if (value == 0) THREADS = sysconf(_SC_NPROCESSORS_ONLN);
  else THREADS = value;
  //printf("\nCPUs = %u\n",THREADS);
}

int stable_get_METHOD1() { return METHOD1; }
void stable_set_METHOD1(int value) { METHOD1 = value; }
int stable_get_METHOD2() { return METHOD2; }
void stable_set_METHOD2(int value) { METHOD2 = value; }
int stable_get_METHOD3() { return METHOD3; }
void stable_set_METHOD3(int value) { METHOD3 = value; }

unsigned int stable_get_IT_MAX() { return IT_MAX; }
void stable_set_IT_MAX(unsigned int value) { IT_MAX = value;}

unsigned int stable_get_INV_MAXITER() { return INV_MAXITER; }
void stable_set_INV_MAXITER(unsigned int value) { INV_MAXITER = value;}

double stable_get_relTOL() { return relTOL; }
void stable_set_relTOL(double value) {
/* Switch between integration methods according to relTOL */
  relTOL = value;
  if (value < 1e-12) {
    stable_set_METHOD1(STABLE_QNG);
    stable_set_METHOD2(STABLE_QAG2);
    stable_set_METHOD3(STABLE_QAG1);
  }
  else if (value < 1e-8) {
    stable_set_METHOD1(STABLE_QNG);
    stable_set_METHOD2(STABLE_QAG2);
    stable_set_METHOD3(STABLE_QAG1);
  }
  else {
    stable_set_METHOD1(STABLE_QNG);
    stable_set_METHOD2(STABLE_QAG1);
    stable_set_METHOD3(STABLE_QAG1);
  }
}

double stable_get_absTOL() { return absTOL; }
void stable_set_absTOL(double value)
{
  absTOL = value;
  //FACTOR = 1e-16;
}

/* Parameter thresholds */

/* When abs(alpha - 1)<ALPHA_TH alpha is set to 1 */
double stable_get_ALPHA_TH() { return ALPHA_TH; }
void stable_set_ALPHA_TH(double value) { ALPHA_TH = value; }

/* When 1-abs(beta)<BETA_TH beta is set to sign(beta)*1.0 */
/* When alpha = 1 and abs(beta)<BETA_TH beta is set to 0.0*/
double stable_get_BETA_TH() { return BETA_TH; }
void stable_set_BETA_TH(double value) { BETA_TH = value; }

/* When abs(x-xxi)<XXI_TH x is set to XXI */
double stable_get_XXI_TH() { return XXI_TH; }
void stable_set_XXI_TH(double value) { XXI_TH = value; }

/* When theta get closer than THETA_TH to integration interval limits theta is set to the limit value */
double stable_get_THETA_TH() { return THETA_TH; }
void stable_set_THETA_TH(double value) { THETA_TH = value; }

/* Debug purposes*/
FILE * stable_get_FINTEG() { return FINTEG; }
FILE * stable_set_FINTEG(char *name)
{
  FINTEG = fopen(name,"wt");
  return FINTEG;
}

/* Log-file configuration */
FILE * stable_get_FLOG() { return FLOG; }
FILE * stable_set_FLOG(char * name)
{
  FLOG = fopen(name,"wt");
  return FLOG;
}

int stable_setparams(StableDist *dist,
                     double alpha, double beta, double sigma, double mu,
                     int parametrization)
{
  int zona;

  if(dist==NULL)
    {
      printf("ERROR");
      return(-1);
    }
  if((zona = stable_checkparams(alpha,beta,sigma,mu,parametrization)) == NOVALID)
    {
    //  printf ("No valid parameters: %lf %lf %lf %lf %d\n",
    //         alpha, beta, sigma, mu, parametrization);
      return zona;
    }

  dist->alpha = alpha;
  dist->beta = beta;
  dist->sigma = sigma;

  switch (zona)
    {
      case STABLE_B1:
        dist->beta = (dist->beta > 0) ? 1.0 : -1.0;
      case STABLE:
        dist->alphainvalpha1 = alpha/(alpha-1.0);
        dist->xi = -beta*tan(0.5*alpha*M_PI);
        dist->theta0 = atan(-dist->xi)/alpha;
        //dist->k1 = pow(1.0+dist->xi*dist->xi,-0.5/(alpha-1.0));
        dist->k1 = -0.5/(alpha-1.0)*log(1.0+dist->xi*dist->xi);
        dist->S = pow(1.0+dist->xi*dist->xi,0.5/alpha);
        dist->Vbeta1 = dist->k1 - dist->alphainvalpha1 * log(dist->alpha)
                                + log(fabs(dist->alpha-1.0));
        dist->stable_pdf_point = &stable_pdf_point_STABLE;
        dist->stable_cdf_point = &stable_cdf_point_STABLE;

        if (alpha < 1.0)
          {
            dist->c1 = 0.5-dist->theta0*M_1_PI;
            dist->c2_part = alpha/((1.0-alpha)*M_PI);
            dist->c3 = M_1_PI;
          }
        else
          {
            dist->c1 = 1.0;
            dist->c2_part = alpha/((alpha-1.0)*M_PI);
            dist->c3 = -M_1_PI;
          }

        if (alpha>1)
         {
          AUX1=log(log(8.5358/(relTOL))/0.9599);/*3.76;*/
          AUX2=log(relTOL);/*-40;*/
         }
        else
         {
          AUX1=log(relTOL);/*-40;*/
          AUX2=log(log(8.5358/(relTOL))/0.9599);/*3.76;*/
         }

        break;

      case ALPHA_1_B1:
        dist->beta = (dist->beta > 0) ? 1.0 : -1.0;
      case ALPHA_1:
        dist->alpha=1;
        dist->c2_part = 0.5/fabs(beta);
        dist->alphainvalpha1 = 0.0;
        dist->xi = 0.0;
        dist->theta0 = M_PI_2;
        dist->k1 = log(2.0*M_1_PI);
        //dist->k1 = 2.0*M_1_PI;
        dist->S = 2.0*M_1_PI;
        dist->c1 = 0.0;
        dist->c3 = M_1_PI;
        dist->Vbeta1=2.0*M_1_PI/M_E;
        dist->stable_pdf_point = &stable_pdf_point_ALPHA_1;
        dist->stable_cdf_point = &stable_cdf_point_ALPHA_1;
        //XXI_TH = 10*EPS;
        if (beta<0)
         {
          AUX1=log(log(8.5358/(relTOL))/0.9599);/*4;*/
          AUX2=log(relTOL);/*-25;*/
         }
        else
         {
          AUX1=log(relTOL);/*-25;*/
          AUX2=log(log(8.5358/(relTOL))/0.9599);/*4;*/
         }
        break;

      case CAUCHY:
        dist->beta=0;
        dist->alpha=1;
        dist->c2_part = 0.0;
        dist->alphainvalpha1 = 0.0;
        dist->xi = 0.0;
        dist->theta0 = M_PI_2;
        dist->k1 = log(2.0*M_1_PI);
        //dist->k1 = 2.0*M_1_PI;
        dist->S = 2.0*M_1_PI;
        dist->c1 = 0.0;
        dist->c3 = M_1_PI;
        dist->Vbeta1=2.0*M_1_PI/M_E;
        dist->stable_pdf_point = &stable_pdf_point_CAUCHY;
        dist->stable_cdf_point = &stable_cdf_point_CAUCHY;
        break;

      case GAUSS:
        dist->alpha=2;
        dist->beta=0.0;
        dist->alphainvalpha1 = 2.0;
        dist->xi = 0.0;
        dist->theta0 = 0.0;
        dist->k1 = log(2.0);
        //dist->k1 = 2.0;
        dist->S = 2.0;
        dist->c1 = 1.0;
        dist->c2_part = 2.0*M_1_PI;
        dist->c3 = -M_1_PI;
        dist->Vbeta1=0.25;
        dist->stable_pdf_point = &stable_pdf_point_GAUSS;
        dist->stable_cdf_point = &stable_cdf_point_GAUSS;
        break;

      case LEVY:
        dist->alpha = 0.5;
        dist->beta = (2.0*(beta>0)-1.0);//  1 ó -1
        dist->alphainvalpha1 = -1.0;
        dist->xi = -dist->beta;
        dist->theta0 = 0.5*M_PI;
        dist->k1 = 0.0;
        dist->S = 1.0;
        dist->c1 = 0.0;
        dist->c2_part = 0.5*M_1_PI;
        dist->c3 = M_1_PI;
        dist->Vbeta1 = dist->k1 - dist->alphainvalpha1 * log(dist->alpha)
                                + log(fabs(dist->alpha-1.0));
        dist->stable_pdf_point = &stable_pdf_point_LEVY;
        dist->stable_cdf_point = &stable_cdf_point_LEVY;
        break;
    }

  if (parametrization == 0)
    {
      dist->mu_0 = mu;
      if (dist->alpha==1)
        dist->mu_1 = mu-dist->beta*2*M_1_PI*dist->sigma*log(dist->sigma);
      else dist->mu_1 = mu+dist->xi*dist->sigma;
    }
  else if (parametrization == 1)
    {
      dist->mu_1 = mu;
      if (dist->alpha==1)
        dist->mu_0 = mu+dist->beta*2*M_1_PI*dist->sigma*log(dist->sigma);
      else dist->mu_0 = mu-dist->xi*dist->sigma;
    }

  dist->theta0_ = dist->theta0;
  dist->beta_ = dist->beta;
  dist->xxipow = 0.0;
  dist->ZONE = zona;

  return zona;
}

int stable_checkparams(double alpha, double beta, double sigma, double mu,
                       int parametrization)
{
  /*Check parameters*/
  if (0.0 >= alpha || alpha > 2.0)
    {
      //printf("Alpha must lie between 0.0 and 2.0.");
      return NOVALID;
    }
  else if (beta < -1.0 || beta > 1.0)
    {
      //printf("Beta must lie between -1.0 and 1.0.");
      return NOVALID;
    }
  else if (sigma <= 0.0)
    {
      //printf("Sigma must be positive.");
      return NOVALID;
    }
  else if (isnan(mu) || isinf(mu))
    {
      //printf("Mu must be real.");
      return NOVALID;
    }
  else if (parametrization != 0 && parametrization != 1)
    {
      //printf("Only parametrizations 0 and 1 are accepted.");
      return NOVALID;
    }

  /*ZONE determination*/
  if ((2.0 - alpha) <= ALPHA_TH)
    return GAUSS;  //GAUSS
  else if (fabs(alpha-0.5) <= ALPHA_TH && fabs((fabs(beta)-1.0)) <= BETA_TH)
    return LEVY; //LEVY
  else if (fabs(alpha-1.0) <= ALPHA_TH && fabs(beta) <= BETA_TH)
    return CAUCHY;
  else if (fabs(alpha-1.0) <= ALPHA_TH) {
//    if (fabs(fabs(beta)-1) <= BETA_TH) return ALPHA_1_B1;
    return ALPHA_1;
  }
  else {
//    if (fabs(fabs(beta)-1) <= BETA_TH) return STABLE_B1;
    return STABLE;
  }
  /*When alpha=1,beta=1 Landau distribution is obtained,
    but it is calculated as in the general case*/
}

StableDist * stable_create(double alpha, double beta, double sigma, double mu,
                           int parametrization)
{
  /*gsl_error_handler_t * old_handler;
  old_handler = */gsl_set_error_handler (&error_handler);

  StableDist * dist = (StableDist *) malloc(sizeof (StableDist));

  if (dist == NULL)
    {
      perror("Error during distribution creation.");
      return NULL;
    }
  if((stable_setparams(dist,alpha,beta,sigma,mu,parametrization)) == NOVALID)
    {
      perror ("Error during distribution creation.");
      return NULL;
    }
  gsl_rng_env_setup(); //leemos las variables de entorno
  dist->gslworkspace = gsl_integration_workspace_alloc(IT_MAX);
  dist->gslrand = gsl_rng_alloc (gsl_rng_default);

  //Allow the distribution to use THREADS threads.
  stable_set_THREADS(THREADS);

  return dist;
}

StableDist * stable_copy(StableDist *src_dist)
{
  StableDist *dist;

  dist = stable_create(src_dist->alpha, src_dist->beta,
                       src_dist->sigma, src_dist->mu_0, 0);
  return dist;
}

void stable_free(StableDist *dist)
{
  if (dist == NULL)
    return;
  gsl_integration_workspace_free(dist->gslworkspace);
  gsl_rng_free(dist->gslrand);
  free(dist);
}

max_stars_repo_name: TantraLabs/gonum

max_issues_repo_name: TantraLabs/gonum

max_forks_repo_name: TantraLabs/gonum


text:
#pragma once

#include "Rendering/Utility/peBxDF.h"
#include "Type/peColor.h"
#include <gsl.h>

namespace pe {

enum class ToneMapping { Saturate, Linear, Log };

//! \brief Performs tone mapping from the given input image to a displayable
//! 8-bit RGBA output image \param result Result image \param input Input image
//! \param imageWidth Width of the image
//! \param imageHeight Height of the image
//! \param strategy The tone-mapping strategy to use
void ToneMap(gsl::span<RGBA_8Bit> result, const gsl::span<Spectrum_t> input,
             const uint32_t imageWidth, const uint32_t imageHeight,
             const ToneMapping strategy);

} // namespace pe

ext: h

max_stars_repo_name: Mortano/Prismatic

max_issues_repo_name: Mortano/Prismatic

max_forks_repo_name: Mortano/Prismatic


text:
/* integration/qng.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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.
 */

#include <config.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

#include "err.c"
#include "qng.h"

int
gsl_integration_qng (const gsl_function *f,
                     double a, double b,
                     double epsabs, double epsrel,
                     double * result, double * abserr, size_t * neval)
{
  double fv1[5], fv2[5], fv3[5], fv4[5];
  double savfun[21];  /* array of function values which have been computed */
  double res10, res21, res43, res87;    /* 10, 21, 43 and 87 point results */
  double result_kronrod, err ;
  double resabs; /* approximation to the integral of abs(f) */
  double resasc; /* approximation to the integral of abs(f-i/(b-a)) */

  const double half_length =  0.5 * (b - a);
  const double abs_half_length = fabs (half_length);
  const double center = 0.5 * (b + a);
  const double f_center = GSL_FN_EVAL(f, center);

  int k ;

  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
    {
      * result = 0;
      * abserr = 0;
      * neval = 0;
      GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
                 GSL_EBADTOL);
    };

  /* Compute the integral using the 10- and 21-point formula. */

  res10 = 0;
  res21 = w21b[5] * f_center;
  resabs = w21b[5] * fabs (f_center);

  for (k = 0; k < 5; k++)
    {
      const double abscissa = half_length * x1[k];
      const double fval1 = GSL_FN_EVAL(f, center + abscissa);
      const double fval2 = GSL_FN_EVAL(f, center - abscissa);
      const double fval = fval1 + fval2;
      res10 += w10[k] * fval;
      res21 += w21a[k] * fval;
      resabs += w21a[k] * (fabs (fval1) + fabs (fval2));
      savfun[k] = fval;
      fv1[k] = fval1;
      fv2[k] = fval2;
    }

  for (k = 0; k < 5; k++)
    {
      const double abscissa = half_length * x2[k];
      const double fval1 = GSL_FN_EVAL(f, center + abscissa);
      const double fval2 = GSL_FN_EVAL(f, center - abscissa);
      const double fval = fval1 + fval2;
      res21 += w21b[k] * fval;
      resabs += w21b[k] * (fabs (fval1) + fabs (fval2));
      savfun[k + 5] = fval;
      fv3[k] = fval1;
      fv4[k] = fval2;
    }

  resabs *= abs_half_length ;

  {
    const double mean = 0.5 * res21;

    resasc = w21b[5] * fabs (f_center - mean);

    for (k = 0; k < 5; k++)
      {
        resasc +=
          (w21a[k] * (fabs (fv1[k] - mean) + fabs (fv2[k] - mean))
          + w21b[k] * (fabs (fv3[k] - mean) + fabs (fv4[k] - mean)));
      }
    resasc *= abs_half_length ;
  }

  result_kronrod = res21 * half_length;

  err = rescale_error ((res21 - res10) * half_length, resabs, resasc) ;

  /*   test for convergence. */

  if (err < epsabs || err < epsrel * fabs (result_kronrod))
    {
      * result = result_kronrod ;
      * abserr = err ;
      * neval = 21;
      return GSL_SUCCESS;
    }

  /* compute the integral using the 43-point formula. */

  res43 = w43b[11] * f_center;

  for (k = 0; k < 10; k++)
    {
      res43 += savfun[k] * w43a[k];
    }

  for (k = 0; k < 11; k++)
    {
      const double abscissa = half_length * x3[k];
      const double fval = (GSL_FN_EVAL(f, center + abscissa)
                           + GSL_FN_EVAL(f, center - abscissa));
      res43 += fval * w43b[k];
      savfun[k + 10] = fval;
    }

  /*  test for convergence */

  result_kronrod = res43 * half_length;
  err = rescale_error ((res43 - res21) * half_length, resabs, resasc);

  if (err < epsabs || err < epsrel * fabs (result_kronrod))
    {
      * result = result_kronrod ;
      * abserr = err ;
      * neval = 43;
      return GSL_SUCCESS;
    }

  /* compute the integral using the 87-point formula. */

  res87 = w87b[22] * f_center;

  for (k = 0; k < 21; k++)
    {
      res87 += savfun[k] * w87a[k];
    }

  for (k = 0; k < 22; k++)
    {
      const double abscissa = half_length * x4[k];
      res87 += w87b[k] * (GSL_FN_EVAL(f, center + abscissa)
                          + GSL_FN_EVAL(f, center - abscissa));
    }

  /*  test for convergence */

  result_kronrod = res87 * half_length ;

  err = rescale_error ((res87 - res43) * half_length, resabs, resasc);

  if (err < epsabs || err < epsrel * fabs (result_kronrod))
    {
      * result = result_kronrod ;
      * abserr = err ;
      * neval = 87;
      return GSL_SUCCESS;
    }

  /* failed to converge */

  * result = result_kronrod ;
  * abserr = err ;
  * neval = 87;

  GSL_ERROR("failed to reach tolerance with highest-order rule", GSL_ETOL) ;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
// reserved. See files LICENSE and NOTICE for details.
//
// This file is part of CEED, a collection of benchmarks, miniapps, software
// libraries and APIs for efficient high-order finite element and spectral
// element discretizations for exascale applications. For more information and
// source code availability see http://github.com/ceed.
//
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
// a collaborative effort of two U.S. Department of Energy organizations (Office
// of Science and the National Nuclear Security Administration) responsible for
// the planning and preparation of a capable exascale ecosystem, including
// software, applications, hardware, advanced system engineering and early
// testbed platforms, in support of the nation's exascale computing imperative.

//                        libCEED + PETSc Example: CEED BPs
//
// This example demonstrates a simple usage of libCEED with PETSc to solve the
// CEED BP benchmark problems, see http://ceed.exascaleproject.org/bps,
// on a closed surface, such as the one of a discrete sphere.
//
// The code uses higher level communication protocols in DMPlex.
//
// Build with:
//
//     make bpssphere [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>]
//
// Sample runs:
//
//     bpssphere -problem bp1 -degree 3
//     bpssphere -problem bp2 -degree 3
//     bpssphere -problem bp3 -degree 3
//     bpssphere -problem bp4 -degree 3
//     bpssphere -problem bp5 -degree 3 -ceed /cpu/self
//     bpssphere -problem bp6 -degree 3 -ceed /gpu/cuda
//
//TESTARGS -ceed {ceed_resource} -test -problem bp3 -degree 3 -dm_refine 2

/// @file
/// CEED BPs example using PETSc with DMPlex
/// See bps.c for a "raw" implementation using a structured grid.
/// and bpsdmplex.c for an implementation using an unstructured grid.
static const char help[] = "Solve CEED BPs on a sphere using DMPlex in PETSc\n";

#include <stdbool.h>
#include <string.h>
#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscksp.h>

#include "bpssphere.h"
#include "include/sphereproblemdata.h"
#include "include/petscutils.h"
#include "include/petscversion.h"
#include "include/matops.h"
#include "include/libceedsetup.h"

#if PETSC_VERSION_LT(3,12,0)
#ifdef PETSC_HAVE_CUDA
#include <petsccuda.h>
// Note: With PETSc prior to version 3.12.0, providing the source path to
//       include 'cublas_v2.h' will be needed to use 'petsccuda.h'.
#endif
#endif

int main(int argc, char **argv) {
  PetscInt ierr;
  MPI_Comm comm;
  char ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self",
      filename[PETSC_MAX_PATH_LEN];
  double my_rt_start, my_rt, rt_min, rt_max;
  PetscInt degree = 3, q_extra, l_size, g_size, topo_dim = 2, num_comp_x = 3,
           num_comp_u = 1, xl_size;
  PetscScalar *r;
  PetscBool test_mode, benchmark_mode, read_mesh, write_solution, simplex;
  PetscLogStage solve_stage;
  Vec X, X_loc, rhs, rhs_loc;
  Mat mat_O;
  KSP ksp;
  DM  dm;
  UserO user_O;
  Ceed ceed;
  CeedData ceed_data;
  CeedQFunction qf_error;
  CeedOperator op_error;
  CeedVector rhs_ceed, target;
  BPType bp_choice;
  VecType vec_type;
  PetscMemType mem_type;

  ierr = PetscInitialize(&argc, &argv, NULL, help);
  if (ierr) return ierr;
  comm = PETSC_COMM_WORLD;

  // Read command line options
  ierr = PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); CHKERRQ(ierr);
  bp_choice = CEED_BP1;
  ierr = PetscOptionsEnum("-problem",
                          "CEED benchmark problem to solve", NULL,
                          bp_types, (PetscEnum)bp_choice, (PetscEnum *)&bp_choice,
                          NULL); CHKERRQ(ierr);
  num_comp_u = bp_options[bp_choice].num_comp_u;
  test_mode = PETSC_FALSE;
  ierr = PetscOptionsBool("-test",
                          "Testing mode (do not print unless error is large)",
                          NULL, test_mode, &test_mode, NULL); CHKERRQ(ierr);
  benchmark_mode = PETSC_FALSE;
  ierr = PetscOptionsBool("-benchmark",
                          "Benchmarking mode (prints benchmark statistics)",
                          NULL, benchmark_mode, &benchmark_mode, NULL);
  CHKERRQ(ierr);
  write_solution = PETSC_FALSE;
  ierr = PetscOptionsBool("-write_solution",
                          "Write solution for visualization",
                          NULL, write_solution, &write_solution, NULL);
  CHKERRQ(ierr);
  degree = test_mode ? 3 : 2;
  ierr = PetscOptionsInt("-degree", "Polynomial degree of tensor product basis",
                         NULL, degree, &degree, NULL); CHKERRQ(ierr);
  q_extra = bp_options[bp_choice].q_extra;
  ierr = PetscOptionsInt("-q_extra", "Number of extra quadrature points",
                         NULL, q_extra, &q_extra, NULL); CHKERRQ(ierr);
  ierr = PetscOptionsString("-ceed", "CEED resource specifier",
                            NULL, ceed_resource, ceed_resource,
                            sizeof(ceed_resource), NULL); CHKERRQ(ierr);
  read_mesh = PETSC_FALSE;
  ierr = PetscOptionsString("-mesh", "Read mesh from file", NULL,
                            filename, filename, sizeof(filename), &read_mesh);
  CHKERRQ(ierr);
  simplex = PETSC_FALSE;
  ierr = PetscOptionsBool("-simplex", "Use simplices, or tensor product cells",
                          NULL, simplex, &simplex, NULL); CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  // Setup DM
  if (read_mesh) {
    ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, NULL, PETSC_TRUE,
                                &dm);
    CHKERRQ(ierr);
  } else {
    // Create the mesh as a 0-refined sphere. This will create a cubic surface, not a box
    ierr = DMPlexCreateSphereMesh(PETSC_COMM_WORLD, topo_dim, simplex, 1., &dm);
    CHKERRQ(ierr);
    // Set the object name
    ierr = PetscObjectSetName((PetscObject)dm, "Sphere"); CHKERRQ(ierr);
    // Distribute mesh over processes
    {
      DM dm_dist = NULL;
      PetscPartitioner part;

      ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr);
      ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr);
      ierr = DMPlexDistribute(dm, 0, NULL, &dm_dist); CHKERRQ(ierr);
      if (dm_dist) {
        ierr = DMDestroy(&dm); CHKERRQ(ierr);
        dm  = dm_dist;
      }
    }
    // Refine DMPlex with uniform refinement using runtime option -dm_refine
    ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr);
    ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
    ierr = ProjectToUnitSphere(dm); CHKERRQ(ierr);
    // View DMPlex via runtime option
    ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr);
  }

  // Create DM
  ierr = SetupDMByDegree(dm, degree, num_comp_u, topo_dim, false,
                         (BCFunction)NULL);
  CHKERRQ(ierr);

  // Create vectors
  ierr = DMCreateGlobalVector(dm, &X); CHKERRQ(ierr);
  ierr = VecGetLocalSize(X, &l_size); CHKERRQ(ierr);
  ierr = VecGetSize(X, &g_size); CHKERRQ(ierr);
  ierr = DMCreateLocalVector(dm, &X_loc); CHKERRQ(ierr);
  ierr = VecGetSize(X_loc, &xl_size); CHKERRQ(ierr);
  ierr = VecDuplicate(X, &rhs); CHKERRQ(ierr);

  // Operator
  ierr = PetscMalloc1(1, &user_O); CHKERRQ(ierr);
  ierr = MatCreateShell(comm, l_size, l_size, g_size, g_size,
                        user_O, &mat_O); CHKERRQ(ierr);
  ierr = MatShellSetOperation(mat_O, MATOP_MULT,
                              (void(*)(void))MatMult_Ceed); CHKERRQ(ierr);

  // Set up libCEED
  CeedInit(ceed_resource, &ceed);
  CeedMemType mem_type_backend;
  CeedGetPreferredMemType(ceed, &mem_type_backend);

  ierr = DMGetVecType(dm, &vec_type); CHKERRQ(ierr);
  if (!vec_type) { // Not yet set by user -dm_vec_type
    switch (mem_type_backend) {
    case CEED_MEM_HOST: vec_type = VECSTANDARD; break;
    case CEED_MEM_DEVICE: {
      const char *resolved;
      CeedGetResource(ceed, &resolved);
      if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA;
      else if (strstr(resolved, "/gpu/hip/occa"))
        vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678
      else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP;
      else vec_type = VECSTANDARD;
    }
    }
    ierr = DMSetVecType(dm, vec_type); CHKERRQ(ierr);
  }

  // Print summary
  if (!test_mode) {
    PetscInt P = degree + 1, Q = P + q_extra;
    const char *used_resource;
    CeedGetResource(ceed, &used_resource);
    ierr = PetscPrintf(comm,
                       "\n-- CEED Benchmark Problem %d on the Sphere -- libCEED + PETSc --\n"
                       "  libCEED:\n"
                       "    libCEED Backend                    : %s\n"
                       "    libCEED Backend MemType            : %s\n"
                       "  Mesh:\n"
                       "    Number of 1D Basis Nodes (p)       : %d\n"
                       "    Number of 1D Quadrature Points (q) : %d\n"
                       "    Global nodes                       : %D\n",
                       bp_choice+1, ceed_resource, CeedMemTypes[mem_type_backend], P, Q,
                       g_size/num_comp_u); CHKERRQ(ierr);
  }

  // Create RHS vector
  ierr = VecDuplicate(X_loc, &rhs_loc); CHKERRQ(ierr);
  ierr = VecZeroEntries(rhs_loc); CHKERRQ(ierr);
  ierr = VecGetArrayAndMemType(rhs_loc, &r, &mem_type); CHKERRQ(ierr);
  CeedVectorCreate(ceed, xl_size, &rhs_ceed);
  CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r);

  // Setup libCEED's objects
  ierr = PetscMalloc1(1, &ceed_data); CHKERRQ(ierr);
  ierr = SetupLibceedByDegree(dm, ceed, degree, topo_dim, q_extra, num_comp_x,
                              num_comp_u, g_size, xl_size, bp_options[bp_choice],
                              ceed_data, true, rhs_ceed, &target); CHKERRQ(ierr);

  // Gather RHS
  CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL);
  ierr = VecRestoreArrayAndMemType(rhs_loc, &r); CHKERRQ(ierr);
  ierr = VecZeroEntries(rhs); CHKERRQ(ierr);
  ierr = DMLocalToGlobal(dm, rhs_loc, ADD_VALUES, rhs); CHKERRQ(ierr);
  CeedVectorDestroy(&rhs_ceed);

  // Create the error Q-function
  CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].error,
                              bp_options[bp_choice].error_loc, &qf_error);
  CeedQFunctionAddInput(qf_error, "u", num_comp_u, CEED_EVAL_INTERP);
  CeedQFunctionAddInput(qf_error, "true_soln", num_comp_u, CEED_EVAL_NONE);
  CeedQFunctionAddOutput(qf_error, "error", num_comp_u, CEED_EVAL_NONE);

  // Create the error operator
  CeedOperatorCreate(ceed, qf_error, NULL, NULL, &op_error);
  CeedOperatorSetField(op_error, "u", ceed_data->elem_restr_u,
                       ceed_data->basis_u, CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_error, "true_soln", ceed_data->elem_restr_u_i,
                       CEED_BASIS_COLLOCATED, target);
  CeedOperatorSetField(op_error, "error", ceed_data->elem_restr_u_i,
                       CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);

  // Set up Mat
  user_O->comm = comm;
  user_O->dm = dm;
  user_O->X_loc = X_loc;
  ierr = VecDuplicate(X_loc, &user_O->Y_loc); CHKERRQ(ierr);
  user_O->x_ceed = ceed_data->x_ceed;
  user_O->y_ceed = ceed_data->y_ceed;
  user_O->op = ceed_data->op_apply;
  user_O->ceed = ceed;

  // Setup solver
  ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr);
  {
    PC pc;
    ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
    if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) {
      ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr);
      ierr = PCJacobiSetType(pc, PC_JACOBI_ROWSUM); CHKERRQ(ierr);
    } else {
      ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr);
      MatNullSpace nullspace;

      ierr = MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
      CHKERRQ(ierr);
      ierr = MatSetNullSpace(mat_O, nullspace); CHKERRQ(ierr);
      ierr = MatNullSpaceDestroy(&nullspace); CHKERRQ(ierr);
    }
    ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr);
    ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr);
    ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT,
                            PETSC_DEFAULT); CHKERRQ(ierr);
  }
  ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp, mat_O, mat_O); CHKERRQ(ierr);

  // First run, if benchmarking
  if (benchmark_mode) {
    ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1);
    CHKERRQ(ierr);
    my_rt_start = MPI_Wtime();
    ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr);
    my_rt = MPI_Wtime() - my_rt_start;
    ierr = MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, comm);
    CHKERRQ(ierr);
    // Set maxits based on first iteration timing
    if (my_rt > 0.02) {
      ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 5);
      CHKERRQ(ierr);
    } else {
      ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 20);
      CHKERRQ(ierr);
    }
  }

  // Timed solve
  ierr = VecZeroEntries(X); CHKERRQ(ierr);
  ierr = PetscBarrier((PetscObject)ksp); CHKERRQ(ierr);

  // -- Performance logging
  ierr = PetscLogStageRegister("Solve Stage", &solve_stage); CHKERRQ(ierr);
  ierr = PetscLogStagePush(solve_stage); CHKERRQ(ierr);

  // -- Solve
  my_rt_start = MPI_Wtime();
  ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr);
  my_rt = MPI_Wtime() - my_rt_start;

  // -- Performance logging
  ierr = PetscLogStagePop();

  // Output results
  {
    KSPType ksp_type;
    KSPConvergedReason reason;
    PetscReal rnorm;
    PetscInt its;
    ierr = KSPGetType(ksp, &ksp_type); CHKERRQ(ierr);
    ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr);
    ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
    ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr);
    if (!test_mode || reason < 0 || rnorm > 1e-8) {
      ierr = PetscPrintf(comm,
                         "  KSP:\n"
                         "    KSP Type                           : %s\n"
                         "    KSP Convergence                    : %s\n"
                         "    Total KSP Iterations               : %D\n"
                         "    Final rnorm                        : %e\n",
                         ksp_type, KSPConvergedReasons[reason], its,
                         (double)rnorm); CHKERRQ(ierr);
    }
    if (!test_mode) {
      ierr = PetscPrintf(comm,"  Performance:\n"); CHKERRQ(ierr);
    }
    {
      PetscReal max_error;
      ierr = ComputeErrorMax(user_O, op_error, X, target, &max_error);
      CHKERRQ(ierr);
      PetscReal tol = 5e-4;
      if (!test_mode || max_error > tol) {
        ierr = MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, comm);
        CHKERRQ(ierr);
        ierr = MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, comm);
        CHKERRQ(ierr);
        ierr = PetscPrintf(comm,
                           "    Pointwise Error (max)              : %e\n"
                           "    CG Solve Time                      : %g (%g) sec\n",
                           (double)max_error, rt_max, rt_min); CHKERRQ(ierr);
      }
    }
    if (benchmark_mode && (!test_mode)) {
      ierr = PetscPrintf(comm,
                         "    DoFs/Sec in CG                     : %g (%g) million\n",
                         1e-6*g_size*its/rt_max, 1e-6*g_size*its/rt_min); CHKERRQ(ierr);
    }
  }

  // Output solution
  if (write_solution) {
    PetscViewer vtk_viewer_soln;

    ierr = PetscViewerCreate(comm, &vtk_viewer_soln); CHKERRQ(ierr);
    ierr = PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK); CHKERRQ(ierr);
    ierr = PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu"); CHKERRQ(ierr);
    ierr = VecView(X, vtk_viewer_soln); CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&vtk_viewer_soln); CHKERRQ(ierr);
  }

  // Cleanup
  ierr = VecDestroy(&X); CHKERRQ(ierr);
  ierr = VecDestroy(&X_loc); CHKERRQ(ierr);
  ierr = VecDestroy(&user_O->Y_loc); CHKERRQ(ierr);
  ierr = MatDestroy(&mat_O); CHKERRQ(ierr);
  ierr = PetscFree(user_O); CHKERRQ(ierr);
  ierr = CeedDataDestroy(0, ceed_data); CHKERRQ(ierr);
  ierr = DMDestroy(&dm); CHKERRQ(ierr);

  ierr = VecDestroy(&rhs); CHKERRQ(ierr);
  ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
  CeedVectorDestroy(&target);
  CeedQFunctionDestroy(&qf_error);
  CeedOperatorDestroy(&op_error);
  CeedDestroy(&ceed);
  return PetscFinalize();
}

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
/*
* -----------------------------------------------------------------
*  bst_lib.h
*  Binary Search Tree Library
*  Version: 2.0
*  Last Update: Oct 9, 2019
*
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010-2019, Americo Barbosa da Cunha Junior
*  All rights reserved.
* -----------------------------------------------------------------
*  This is the header file for BST_LIB module, a computational
*  library with routines to manipulate a binary search tree.
* -----------------------------------------------------------------
*/

#ifndef __BST_LIB_H__
#define __BST_LIB_H__

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

#undef  BST_RIGHT
#define BST_RIGHT 1

#undef  BST_LEFT
#define BST_LEFT 0

/*
*------------------------------------------------------------
*   structure of binary search tree leaf
*
*   last update: Feb 2, 2009
*------------------------------------------------------------
*/

typedef struct leaf
{
    gsl_vector *phi;   /* composition */
    gsl_vector *Rphi;  /* reaction mapping */
    gsl_matrix *A;     /* mapping gradient matrix */
    gsl_matrix *L;     /* ellipsoid Cholesky matrix */
} bst_leaf;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   structure of binary search tree node
*
*   last update: Feb 2, 2009
*------------------------------------------------------------
*/

typedef struct node
{
    gsl_vector *v;          /* cutting plane normal vector */
    double a;               /* scalar a */
    struct node *r_node;    /* right node */
    struct node *l_node;    /* left node */
    struct leaf *r_leaf;    /* right leaf */
    struct leaf *l_leaf;    /* left leaf */
} bst_node;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

bst_leaf *bst_leaf_alloc();

bst_node* bst_node_alloc();

void bst_leaf_free(void **leaf_bl);

void bst_node_free(void **node_bl);

int bst_leaf_set(gsl_vector *phi,
                 gsl_vector *Rphi,
                 gsl_matrix *A,
                 gsl_matrix *L,
                 bst_leaf *leaf);

double bst_cutplane(gsl_vector *phi_l,
                    gsl_vector *phi_r,
                    gsl_vector *v);

int bst_node_set(bst_leaf *l_leaf,
                 bst_leaf *r_leaf,
                 bst_node *node);

int bst_node_add(unsigned int side,
                 bst_node *old_node,
                 bst_node *new_node);

int bst_search(bst_node *root,
               gsl_vector *phi,
               bst_node **end_node,
               bst_leaf **end_leaf);

int bst_height(bst_node *root);

#endif /* __BST_LIB_H__ */

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
// xi2bandpow.c
// 31.08.2018
// converts xi_+/-, xi_g+/x, and xi_gg to bandpowers

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sf_bessel.h>
#include "bjutils.h"

#define xi_pm_balance 0.5 // default weight for xi_+/- in P_E

void xi_to_derived2pt();
double K_bandpow_ee();
double K_bandpow_ne();
double K_bandpow_nn();
double apod_upper();
double apod_lower();

int main(int argc, char *argv[])
{
  int i,j;
  double dr,xip,xim,scale_terms,kernel;
  FILE *dat;
  char name[400],path[200],ident[200],outident[200];

  if (argc!=14) {
    printf("Syntax: %s\n1: <working directory>\n2: <input file identifier>\n3: <output file identifier>\n4: <number of input angular bins>\n5: <min input separation to use in conversion [arcmin] (xi_+ in case of ee)>\n6: <max input separation to use in conversion [arcmin] (xi_+ in case of ee)>\n7: <min input separation to use in conversion [arcmin] (xi_- in case of ee; otherwise unused)>\n8: <max input separation to use in conversion [arcmin] (xi_- in case of ee; otherwise unused)>\n9: <number of output angular bins>\n10: <min output angular frequency>\n11: <max output angular frequency>\n12: <correlation type (1: ee; 2: ne; 3: gg)>\n13: <log width of apodisation window [total width of apodised range is tmax/tmin=exp(width) in arcmin; <0 for no apodisation]>\n\n",argv[0]);
    exit(-1);
  }
  sprintf(path,"%s",argv[1]);
  sprintf(ident,"%s",argv[2]);
  sprintf(outident,"%s",argv[3]);
  const int NR=atoi(argv[4]);
  const double angmin1=atof(argv[5]);
  const double angmax1=atof(argv[6]);
  const double angmin2=atof(argv[7]);
  const double angmax2=atof(argv[8]);
  const int NOUT=atoi(argv[9]);
  const double ellmin=atof(argv[10]);
  const double ellmax=atof(argv[11]);
  const int corrtype=atoi(argv[12]);
  const double logwidth=atof(argv[13]);

  // allocations
  if (corrtype>3) {
    printf("Error: incorrect input for correlation type!\n");
    exit(-1);
  }

  double *tin=bj_alloc(NR);
  double *ell_centre=bj_alloc(NOUT);
  double *ell_bound=bj_alloc(NOUT+1);
  double *weight=bj_alloc(NR);
  double *weight_xim=bj_alloc(NR);
  double *pfrac=bj_alloc(NOUT);

  gsl_vector *xi=gsl_vector_calloc(2*NR);
  gsl_vector *derived2pt=gsl_vector_calloc(2*NOUT);
  gsl_matrix *trans=gsl_matrix_calloc(2*NOUT,2*NR);
  gsl_matrix *xierr=gsl_matrix_calloc(2*NR,2*NR);
  gsl_matrix *derived2pterr=gsl_matrix_calloc(2*NOUT,2*NOUT);

  double logbinwidth=log(ellmax/ellmin)/(1.*NOUT);  // output binning
  for (i=0;i<NOUT;i++) {
    ell_centre[i]=ellmin*exp((i+0.5)*logbinwidth);   // use log centre of bin
    ell_bound[i]=ellmin*exp(i*logbinwidth);
  }
  ell_bound[NOUT]=ellmin*exp(NOUT*logbinwidth);

  // read input data
  sprintf(name,"%s/xi2bandpow_input_%s.dat",path,ident);
  dat=bj_fileopen("r",name);
  for (i=0;i<NR;i++) {
    if (corrtype==3) {
      fscanf(dat,"%lf  %lf",&tin[i],&xip);  //theta (bin centres; in arcmin), w_gg
      xim=0.0;   // use only first half of vector
    }
    else {
      fscanf(dat,"%lf  %lf  %lf",&tin[i],&xip,&xim);  //theta (bin centres), xi_+, xi_- for ee; theta (bin centres), xi_g+, xi_gx for ne
    }

    gsl_vector_set(xi,i,xip);
    gsl_vector_set(xi,NR+i,xim);
    gsl_matrix_set(xierr,i,i,xip*xip);       // treat xi as noise, uncorrelated so diagonal
    gsl_matrix_set(xierr,NR+i,NR+i,xim*xim);

    if (logwidth<0) {   // no apodisation
      if ((tin[i]>angmin1)&&(tin[i]<angmax1)) {
	weight[i]=1.0;
      }
      else {
	weight[i]=0.0;
      }
    }
    else {
      weight[i]=apod_lower(tin[i],angmin1,logwidth)*apod_upper(tin[i],angmax1,logwidth);
    }

    if (corrtype==1) { //ee
      if (logwidth<0) {   // no apodisation
	if ((tin[i]>angmin2)&&(tin[i]<angmax2)) {
	  weight_xim[i]=1.0;
	}
	else {
	  weight_xim[i]=0.0;
	}
      }
      else {
	weight_xim[i]=apod_lower(tin[i],angmin2,logwidth)*apod_upper(tin[i],angmax2,logwidth);
      }
    }

    tin[i]*=PI/10800.;   // convert from arcmin to radian
  }  // end of i loop
  fclose(dat);

  if ((tin[NR-1]<angmax1*PI/10800.*exp(logwidth/2.))||(corrtype==1 && tin[NR-1]<angmax2*PI/10800.*exp(logwidth/2.))) {
    printf("Error: correlation function not available out to max. range requested: %g < %g ; %g\n",tin[NR-1]/PI*10800.,angmax1*exp(logwidth/2.),angmax2*PI/10800.*exp(logwidth/2.));
    exit(-1);
  }
  if ((tin[0]>angmin1*PI/10800.*exp(-logwidth/2.))||(corrtype==1 && tin[0]>angmin2*PI/10800.*exp(-logwidth/2.))) {
    printf("Error: correlation function not available out to min. range requested: %g > %g ; %g\n",tin[0]/PI*10800.,angmin1*exp(-logwidth/2.),angmin2*PI/10800.*exp(-logwidth/2.));
    exit(-1);
  }

  sprintf(name,"%s/xi2bandpow_pmweights_%s.dat",path,ident);
  if ((dat=fopen(name,"r"))==NULL) {
    printf("No file %s - use pmweight of %g instead.\n",name,xi_pm_balance);
    for (i=0;i<NOUT;i++) {
      pfrac[i]=xi_pm_balance;
    }
  }
  else {
    printf("Using pmweights from file %s\n",name);
    for (i=0;i<NOUT;i++) {
      fscanf(dat,"%lf",&pfrac[i]);  // xi_+/- weight in P_E; one value per output ell bin
      printf("K_+ for Bin %i: %g\n",i+1,pfrac[i]);
    }
    fclose(dat);
  }

  // construct transformation matrix
  logbinwidth=log(tin[NR-1]/tin[0])/(1.*NR-1.);
  for (j=0;j<NR;j++) {
    dr=tin[j]*(exp(0.5*logbinwidth)-exp(-0.5*logbinwidth));  // input bin width
    for (i=0;i<NOUT;i++) {
      scale_terms=2.*PI*dr/(tin[j]*log(ell_bound[i+1]/ell_bound[i]));

      kernel=0.0;
      if (corrtype==1) {
	kernel=K_bandpow_ee(tin[j]*ell_bound[i+1],0)-K_bandpow_ee(tin[j]*ell_bound[i],0);
	gsl_matrix_set(trans,i,j,pfrac[i]*weight[j]*kernel*scale_terms);
	gsl_matrix_set(trans,NOUT+i,j,0.5*weight[j]*kernel*scale_terms);

	kernel=K_bandpow_ee(tin[j]*ell_bound[i+1],1)-K_bandpow_ee(tin[j]*ell_bound[i],1);
	gsl_matrix_set(trans,i,NR+j,(1.-pfrac[i])*weight_xim[j]*kernel*scale_terms);
	gsl_matrix_set(trans,NOUT+i,NR+j,-0.5*weight_xim[j]*kernel*scale_terms);
      }
      if (corrtype==2) {
	kernel=K_bandpow_ne(tin[j]*ell_bound[i+1])-K_bandpow_ne(tin[j]*ell_bound[i]);
	gsl_matrix_set(trans,i,j,weight[j]*kernel*scale_terms);
	gsl_matrix_set(trans,NOUT+i,NR+j,weight[j]*kernel*scale_terms);
      }
      if (corrtype==3) {
	kernel=K_bandpow_nn(tin[j]*ell_bound[i+1])-K_bandpow_nn(tin[j]*ell_bound[i]);
	gsl_matrix_set(trans,i,j,weight[j]*kernel*scale_terms);
      }
    }
  }

  // convert signals
  xi_to_derived2pt(xi,xierr,trans,derived2pt,derived2pterr,2*NR,2*NOUT);

  // output band powers
  sprintf(name,"%s/xi2bandpow_output_%s.dat",path,outident);
  dat=bj_fileopen("w",name);
  fprintf(dat,"#   ell         PEE             PEEerr          PBB                PBBerr\n");
  for (i=0;i<NOUT;i++) {
    if (corrtype==3) {
      fprintf(dat,"%15.10g\t%15.10g\t%15.10g\n",ell_centre[i],gsl_vector_get(derived2pt,i),sqrt(gsl_matrix_get(derived2pterr,i,i)));   // ell, bandpow, err
    }
    else {
      fprintf(dat,"%15.10g\t%15.10g\t%15.10g\t%15.10g\t%15.10g\n",ell_centre[i],gsl_vector_get(derived2pt,i),sqrt(gsl_matrix_get(derived2pterr,i,i)),gsl_vector_get(derived2pt,NOUT+i),sqrt(gsl_matrix_get(derived2pterr,NOUT+i,NOUT+i)));   // ell, E/+, err, B/x, err
    }
  }
  fclose(dat);

  sprintf(name,"%s/xi2bandpow_kernels_%s.dat",path,outident);  // write kernels
  dat=bj_fileopen("w",name);
  for (j=0;j<NR;j++) {
    fprintf(dat,"%15.10g\t",tin[j]/PI*10800.);
    for (i=0;i<NOUT;i++) {
      fprintf(dat,"%15.10g\t",gsl_matrix_get(trans,i,j));   // E-mode, xi_+
    }
    for (i=0;i<NOUT;i++) {
      fprintf(dat,"%15.10g\t",gsl_matrix_get(trans,i,NR+j));   // E-mode, xi_-
    }
    for (i=0;i<NOUT;i++) {
      fprintf(dat,"%15.10g\t",gsl_matrix_get(trans,NOUT+i,j));   // B-mode, xi_+
    }
    for (i=0;i<NOUT;i++) {
      fprintf(dat,"%15.10g\t",gsl_matrix_get(trans,NOUT+i,NR+j));   // B-mode, xi_-
    }
    fprintf(dat,"\n");
  }
  fclose(dat);

  // clean up
  gsl_vector_free(xi);
  gsl_vector_free(derived2pt);
  gsl_matrix_free(trans);
  gsl_matrix_free(xierr);
  gsl_matrix_free(derived2pterr);
  free(tin);
  free(ell_centre);
  free(ell_bound);
  free(weight);
  free(weight_xim);
  free(pfrac);
  return 0;
}

void xi_to_derived2pt(gsl_vector *xi,gsl_matrix *xierr,gsl_matrix *trans,gsl_vector *res,gsl_matrix *reserr,int nin,int nout)
{
  gsl_matrix *tmp=gsl_matrix_calloc(nout,nin);
  gsl_matrix *transt=gsl_matrix_calloc(nin,nout);

  // signal conversion
  mattimesvec_gen(trans,xi,res,nout,nin);

  // error conversion
  mattimesmat_gen(trans,xierr,tmp,nout,nin,nin);
  matrix_transpose(trans,transt,nout,nin);
  mattimesmat_gen(tmp,transt,reserr,nout,nin,nout);

  gsl_matrix_free(tmp);
  gsl_matrix_free(transt);
  return;
}

double K_bandpow_ee(double x,int pm)
{
  double res=0.0;

  if (pm==0) {
    res=x*gsl_sf_bessel_Jn(1,x);
  }
  else if (pm==1) {
    res=(x-8./x)*gsl_sf_bessel_Jn(1,x)-8.*gsl_sf_bessel_Jn(2,x);
  }
  else {
    printf("Error in routine K_bandpow_ee: invalid option for 'pm'.\n");
    exit(-1);
  }

  return(res);
}

double K_bandpow_ne(double x)
{
  return((-1.)*x*gsl_sf_bessel_Jn(1,x)-2.*gsl_sf_bessel_Jn(0,x));
}

double K_bandpow_nn(double x)
{
  return(x*gsl_sf_bessel_Jn(1,x));
}

double apod_upper(double x,double scale,double logwidth)
{
  double res=0.0;
  double logx=log(x);
  double logxmin=log(scale)-logwidth/2.;
  double logxmax=log(scale)+logwidth/2.;

  if (logx<=logxmin) {
    res=1.;
  }
  else if (logx>logxmax) {
    res=0.;
  }
  else {
    res=cos(PI/2.*(logx-logxmin)/(logxmax-logxmin));
    res=res*res;  // cos^2
  }
  return(res);
}

double apod_lower(double x,double scale,double logwidth)
{
  double res=0.0;
  double logx=log(x);
  double logxmin=log(scale)-logwidth/2.;
  double logxmax=log(scale)+logwidth/2.;

  if (logx<=logxmin) {
    res=0.;
  }
  else if (logx>logxmax) {
    res=1.;
  }
  else {
    res=cos(PI/2.*(logx-logxmax)/(logxmax-logxmin));
    res=res*res;  // cos^2
  }
  return(res);
}


text:
/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */

#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>

#include "compat.h"
#include "sip-utils.h"
#include "gslutils.h"
#include "mathutil.h"
#include "sip.h"

double wcs_pixel_center_for_size(double size) {
    return 0.5 + 0.5 * size;
}

int sip_compute_inverse_polynomials(sip_t* sip, int NX, int NY,
                                    double xlo, double xhi,
                                    double ylo, double yhi) {
    int inv_sip_order;
    int M, N;
    int i, j, p, q, gu, gv;
    double maxu, maxv, minu, minv;
    double u, v, U, V;
    gsl_matrix *mA;
    gsl_vector *b1, *b2, *x1, *x2;
    tan_t* tan;

    assert(sip->a_order == sip->b_order);
    assert(sip->ap_order == sip->bp_order);
    tan = &(sip->wcstan);

    logverb("sip_compute-inverse_polynomials: A %i, AP %i\n",
            sip->a_order, sip->ap_order);

    /*
     basic idea: lay down a grid in image, for each gridpoint, push
     through the polynomial to get yourself into warped image
     coordinate (but not yet lifted onto the sky).  Then, using the
     set of warped gridpoints as inputs, fit back to their original
     grid locations as targets.
     */
    inv_sip_order = sip->ap_order;

    // Number of grid points to use:
    if (NX == 0)
        NX = 10 * (inv_sip_order + 1);
    if (NY == 0)
        NY = 10 * (inv_sip_order + 1);
    if (xhi == 0)
        xhi = tan->imagew;
    if (yhi == 0)
        yhi = tan->imageh;

    logverb("NX,NY %i,%i, x range [%f, %f], y range [%f, %f]\n",
            NX,NY, xlo, xhi, ylo, yhi);

    // Number of coefficients to solve for:
    // We only compute the upper triangle polynomial terms
    N = (inv_sip_order + 1) * (inv_sip_order + 2) / 2;

    // Number of samples to fit.
    M = NX * NY;

    mA = gsl_matrix_alloc(M, N);
    b1 = gsl_vector_alloc(M);
    b2 = gsl_vector_alloc(M);
    assert(mA);
    assert(b1);
    assert(b2);

    /*
     *  Rearranging formula (4), (5), and (6) from the SIP paper gives the
     *  following equations:
     *
     *    +----------------------- Linear pixel coordinates in PIXELS
     *    |                        before SIP correction
     *    |                   +--- Intermediate world coordinates in DEGREES
     *    |                   |
     *    v                   v
     *                   -1
     *    U = [CD11 CD12]   * x
     *    V   [CD21 CD22]     y
     *
     *    +---------------- PIXEL distortion delta from telescope to
     *    |                 linear coordinates
     *    |    +----------- Linear PIXEL coordinates before SIP correction
     *    |    |       +--- Polynomial U,V terms in powers of PIXELS
     *    v    v       v
     *
     *    -f(u1,v1) =  p11 p12 p13 p14 p15 ... * ap1
     *    -f(u2,v2) =  p21 p22 p23 p24 p25 ...   ap2
     *    ...
     *
     *    -g(u1,v1) =  p11 p12 p13 p14 p15 ... * bp1
     *    -g(u2,v2) =  p21 p22 p23 p24 p25 ...   bp2
     *    ...
     *
     *  which recovers the A and B's.
     */

    minu = xlo - tan->crpix[0];
    maxu = xhi - tan->crpix[0];
    minv = ylo - tan->crpix[1];
    maxv = yhi - tan->crpix[1];

    // Sample grid locations.
    i = 0;
    for (gu=0; gu<NX; gu++) {
        for (gv=0; gv<NY; gv++) {
            double fuv, guv;
            // Calculate grid position in original image pixels
            u = (gu * (maxu - minu) / (NX-1)) + minu;
            v = (gv * (maxv - minv) / (NY-1)) + minv;
            // compute U=u+f(u,v) and V=v+g(u,v)
            sip_calc_distortion(sip, u, v, &U, &V);
            fuv = U - u;
            guv = V - v;
            // Polynomial terms...
            j = 0;
            for (p = 0; p <= inv_sip_order; p++)
                for (q = 0; q <= inv_sip_order; q++) {
                    if (p + q > inv_sip_order)
                        continue;
                    assert(j < N);
                    gsl_matrix_set(mA, i, j,
                                   pow(U, (double)p) * pow(V, (double)q));
                    j++;
                }
            assert(j == N);
            gsl_vector_set(b1, i, -fuv);
            gsl_vector_set(b2, i, -guv);
            i++;
        }
    }
    assert(i == M);

    // Solve the linear equation.
    if (gslutils_solve_leastsquares_v(mA, 2, b1, &x1, NULL, b2, &x2, NULL)) {
        ANERROR("Failed to solve SIP inverse matrix equation!");
        return -1;
    }

    // Extract the coefficients
    j = 0;
    for (p = 0; p <= inv_sip_order; p++)
        for (q = 0; q <= inv_sip_order; q++) {
            if ((p + q > inv_sip_order))
                continue;
            assert(j < N);
            sip->ap[p][q] = gsl_vector_get(x1, j);
            sip->bp[p][q] = gsl_vector_get(x2, j);
            j++;
        }
    assert(j == N);

    // Check that we found values that actually invert the polynomial.
    // The error should be particularly small at the grid points.
    if (true) {
        // rms error accumulators:
        double sumdu = 0;
        double sumdv = 0;
        int Z;
        for (gu = 0; gu < NX; gu++) {
            for (gv = 0; gv < NY; gv++) {
                double newu, newv;
                // Calculate grid position in original image pixels
                u = (gu * (maxu - minu) / (NX-1)) + minu;
                v = (gv * (maxv - minv) / (NY-1)) + minv;
                sip_calc_distortion(sip, u, v, &U, &V);
                sip_calc_inv_distortion(sip, U, V, &newu, &newv);
                sumdu += square(u - newu);
                sumdv += square(v - newv);
            }
        }
        sumdu /= (NX*NY);
        sumdv /= (NX*NY);
        debug("RMS error of inverting a distortion (at the grid points, in pixels):\n");
        debug("  du: %g\n", sqrt(sumdu));
        debug("  dv: %g\n", sqrt(sumdu));
        debug("  dist: %g\n", sqrt(sumdu + sumdv));

        sumdu = 0;
        sumdv = 0;
        Z = 1000;
        for (i=0; i<Z; i++) {
            double newu, newv;
            u = uniform_sample(minu, maxu);
            v = uniform_sample(minv, maxv);
            sip_calc_distortion(sip, u, v, &U, &V);
            sip_calc_inv_distortion(sip, U, V, &newu, &newv);
            sumdu += square(u - newu);
            sumdv += square(v - newv);
        }
        sumdu /= Z;
        sumdv /= Z;
        debug("RMS error of inverting a distortion (at random points, in pixels):\n");
        debug("  du: %g\n", sqrt(sumdu));
        debug("  dv: %g\n", sqrt(sumdu));
        debug("  dist: %g\n", sqrt(sumdu + sumdv));
    }

    gsl_matrix_free(mA);
    gsl_vector_free(b1);
    gsl_vector_free(b2);
    gsl_vector_free(x1);
    gsl_vector_free(x2);

    return 0;
}

bool tan_pixel_is_inside_image(const tan_t* wcs, double x, double y) {
    return (x >= 1 && x <= wcs->imagew && y >= 1 && y <= wcs->imageh);
}

bool sip_pixel_is_inside_image(const sip_t* wcs, double x, double y) {
    return tan_pixel_is_inside_image(&(wcs->wcstan), x, y);
}

int* sip_filter_stars_in_field(const sip_t* sip, const tan_t* tan,
                               const double* xyz, const double* radec,
                               int N, double** p_xy, int* inds, int* p_Ngood) {
    int i, Ngood;
    int W, H;
    double* xy = NULL;
    bool allocd = false;

    assert(sip || tan);
    assert(xyz || radec);
    assert(p_Ngood);

    Ngood = 0;
    if (!inds) {
        inds = malloc(N * sizeof(int));
        allocd = true;
    }

    if (p_xy)
        xy = malloc(N * 2 * sizeof(double));

    if (sip) {
        W = (int)sip->wcstan.imagew;
        H = (int)sip->wcstan.imageh;
    } else {
        W = (int)tan->imagew;
        H = (int)tan->imageh;
    }

    for (i=0; i<N; i++) {
        double x, y;
        if (xyz) {
            if (sip) {
                if (!sip_xyzarr2pixelxy(sip, xyz + i*3, &x, &y))
                    continue;
            } else {
                if (!tan_xyzarr2pixelxy(tan, xyz + i*3, &x, &y))
                    continue;
            }
        } else {
            if (sip) {
                if (!sip_radec2pixelxy(sip, radec[i*2], radec[i*2+1], &x, &y))
                    continue;
            } else {
                if (!tan_radec2pixelxy(tan, radec[i*2], radec[i*2+1], &x, &y))
                    continue;
            }
        }
        // FIXME -- check half- and one-pixel FITS issues.
        if ((x < 0) || (y < 0) || (x >= W) || (y >= H))
            continue;

        inds[Ngood] = i;
        if (xy) {
            xy[Ngood * 2 + 0] = x;
            xy[Ngood * 2 + 1] = y;
        }
        Ngood++;
    }

    if (allocd)
        inds = realloc(inds, Ngood * sizeof(int));

    if (xy)
        xy = realloc(xy, Ngood * 2 * sizeof(double));
    if (p_xy)
        *p_xy = xy;

    *p_Ngood = Ngood;

    return inds;
}

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: mgreter/astrometrylib

max_forks_repo_name: mgreter/astrometrylib


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

#pragma once

// IWYU pragma: begin_exports
#include <units/bits/external/fixed_string.h>
#include <compare>
#include <cstdint>
#include <cstddef>
// IWYU pragma: end_exports

#include <gsl/gsl-lite.hpp>

namespace units {

namespace detail {

constexpr void validate_ascii_char([[maybe_unused]] char c) noexcept { gsl_Expects((c & 0x80) == 0); }

template<std::size_t N>
constexpr void validate_ascii_string([[maybe_unused]] const char (&s)[N + 1]) noexcept
{
#ifndef NDEBUG
  if constexpr (N != 0)
    for (size_t i = 0; i < N; ++i)
      validate_ascii_char(s[i]);
#endif
}

}  // namespace detail

/**
 * @brief A symbol text representation
 *
 * This class template is responsible for definition and handling of a symbol text
 * representation. In the libary it is used to define symbols of units and prefixes.
 * Each symbol can have two versions: Unicode and ASCI-only.
 *
 * @tparam StandardCharT Character type to be used for a Unicode representation
 * @tparam N The size of a Unicode symbol
 * @tparam M The size of the ASCII-only symbol
 */
template<typename StandardCharT, std::size_t N, std::size_t M>
struct basic_symbol_text {
  basic_fixed_string<StandardCharT, N> standard_;
  basic_fixed_string<char, M> ascii_;

  constexpr basic_symbol_text(char std) noexcept: standard_(std), ascii_(std) { detail::validate_ascii_char(std); }
  constexpr basic_symbol_text(StandardCharT std, char a) noexcept: standard_(std), ascii_(a) { detail::validate_ascii_char(a); }
  constexpr basic_symbol_text(const char (&std)[N + 1]) noexcept: standard_(std), ascii_(std) { detail::validate_ascii_string<N>(std); }
  constexpr basic_symbol_text(const basic_fixed_string<char, N>& std) noexcept: standard_(std), ascii_(std) { detail::validate_ascii_string<N>(std.data_); }
  constexpr basic_symbol_text(const StandardCharT (&std)[N + 1], const char (&a)[M + 1]) noexcept: standard_(std), ascii_(a) { detail::validate_ascii_string<M>(a); }
  constexpr basic_symbol_text(const basic_fixed_string<StandardCharT, N>& std, const basic_fixed_string<char, M>& a) noexcept: standard_(std), ascii_(a) { detail::validate_ascii_string<M>(a.data_); }

  [[nodiscard]] constexpr auto& standard() { return standard_; }
  [[nodiscard]] constexpr const auto& standard() const { return standard_; }
  [[nodiscard]] constexpr auto& ascii() { return ascii_; }
  [[nodiscard]] constexpr const auto& ascii() const { return ascii_; }

  template<std::size_t N2, std::size_t M2>
  [[nodiscard]] constexpr friend basic_symbol_text<StandardCharT, N + N2, M + M2> operator+(
      const basic_symbol_text& lhs, const basic_symbol_text<StandardCharT, N2, M2>& rhs) noexcept
  {
    return basic_symbol_text<StandardCharT, N + N2, M + M2>(
        lhs.standard() + rhs.standard(), lhs.ascii() + rhs.ascii());
  }

  template<std::size_t N2>
  [[nodiscard]] constexpr friend basic_symbol_text<StandardCharT, N + N2, M + N2> operator+(
      const basic_symbol_text& lhs, const basic_fixed_string<StandardCharT, N2>& rhs) noexcept
  {
    return lhs + basic_symbol_text<StandardCharT, N2, N2>(rhs);
  }

  template<std::size_t N2>
  [[nodiscard]] constexpr friend basic_symbol_text<StandardCharT, N + N2, M + N2> operator+(
      const basic_fixed_string<StandardCharT, N2>& lhs, const basic_symbol_text& rhs) noexcept
  {
    return basic_symbol_text<StandardCharT, N2, N2>(lhs) + rhs;
  }

  template<std::size_t N2>
  [[nodiscard]] constexpr friend basic_symbol_text<StandardCharT, N + N2 - 1, M + N2 - 1> operator+(
      const basic_symbol_text& lhs, const StandardCharT (&rhs)[N2]) noexcept
  {
    return lhs + basic_symbol_text<StandardCharT, N2 - 1, N2 - 1>(rhs);
  }

  template<std::size_t N2>
  [[nodiscard]] constexpr friend basic_symbol_text<StandardCharT, N + N2 - 1, M + N2 - 1> operator+(
      const StandardCharT (&lhs)[N2], const basic_symbol_text& rhs) noexcept
  {
    return basic_symbol_text<StandardCharT, N2 - 1, N2 - 1>(lhs) + rhs;
  }

  [[nodiscard]] constexpr friend basic_symbol_text<StandardCharT, N + 1, M + 1> operator+(
      const basic_symbol_text& lhs, StandardCharT rhs) noexcept
  {
    return lhs + basic_symbol_text<StandardCharT, 1, 1>(rhs);
  }

  [[nodiscard]] constexpr friend basic_symbol_text<StandardCharT, N + 1, M + 1> operator+(
      StandardCharT lhs, const basic_symbol_text& rhs) noexcept
  {
    return basic_symbol_text<StandardCharT, 1, 1>(lhs) + rhs;
  }

  template<typename StandardCharT2, std::size_t N2, std::size_t M2>
  [[nodiscard]] friend constexpr auto operator<=>(const basic_symbol_text& lhs,
                                                  const basic_symbol_text<StandardCharT2, N2, M2>& rhs) noexcept
  {
    if (const auto cmp = lhs.standard() <=> rhs.standard(); cmp != 0) return cmp;
    return lhs.ascii() <=> rhs.ascii();
  }

  template<typename StandardCharT2, std::size_t N2, std::size_t M2>
  [[nodiscard]] friend constexpr bool operator==(const basic_symbol_text& lhs,
                                                 const basic_symbol_text<StandardCharT2, N2, M2>& rhs) noexcept
  {
    return lhs.standard() == rhs.standard() && lhs.ascii() == rhs.ascii();
  }
};

basic_symbol_text(char) -> basic_symbol_text<char, 1, 1>;

template<typename StandardCharT>
basic_symbol_text(StandardCharT, char) -> basic_symbol_text<StandardCharT, 1, 1>;

template<std::size_t N>
basic_symbol_text(const char (&)[N]) -> basic_symbol_text<char, N - 1, N - 1>;

template<std::size_t N>
basic_symbol_text(const basic_fixed_string<char, N>&) -> basic_symbol_text<char, N, N>;

template<typename StandardCharT, std::size_t N, std::size_t M>
basic_symbol_text(const StandardCharT (&)[N], const char (&)[M]) -> basic_symbol_text<StandardCharT, N - 1, M - 1>;

template<typename StandardCharT, std::size_t N, std::size_t M>
basic_symbol_text(const basic_fixed_string<StandardCharT, N>&,
                  const basic_fixed_string<char, M>&) -> basic_symbol_text<StandardCharT, N, M>;

}  // namespace units

max_stars_repo_name: bebuch/units

max_issues_repo_name: bebuch/units

max_forks_repo_name: bebuch/units