text:
// Copyright © Matt Jones and Contributors. Licensed under the MIT Licence (MIT). See LICENCE.md in the repository root
// for more information.

#ifndef NOVELRT_EXPERIMENTAL_GRAPHICS_H
#define NOVELRT_EXPERIMENTAL_GRAPHICS_H

// Graphics dependencies
#include "../../Graphics/RGBAColour.h"
#include "../../Maths/Maths.h"
#include "../../Utilities/Event.h"
#include "../../Utilities/Misc.h"
#include "../EngineConfig.h"
#include "../Threading/Threading.h"
#include <chrono>
#include <cstdint>
#include <filesystem>
#include <gsl/span>
#include <list>
#include <memory>
#include <mutex>
#include <optional>
#include <string>
#include <typeindex>
#include <utility>
#include <vector>

/**
 * @brief The experimental Graphics plugin API. Comes with built-in support for the ECS.
 */
namespace NovelRT::Experimental::Graphics
{
    enum class ShaderProgramKind : uint32_t;
    class GraphicsDeviceObject;
    enum class GraphicsResourceAccess : uint32_t;
    enum class GraphicsSurfaceKind : uint32_t;
    class IGraphicsSurface;
    class GraphicsAdapter;
    class GraphicsDevice;
    class GraphicsResource;
    class GraphicsBuffer;
    class GraphicsTexture;
    class ShaderProgram;
    class GraphicsPipeline;
    class GraphicsPipelineSignature;
    class GraphicsPipelineInput;
    class GraphicsPipelineResource;
    class GraphicsPipelineInputElement;
    enum class GraphicsPipelineInputElementKind : uint32_t;
    enum class GraphicsPipelineResourceKind : uint32_t;
    enum class ShaderProgramVisibility : uint32_t;
    class GraphicsContext;
    class GraphicsFence;
    class GraphicsPrimitive;
    class GraphicsProvider;
    class GraphicsMemoryAllocator;
    struct GraphicsMemoryAllocatorSettings;
    enum class GraphicsTextureKind : uint32_t;
    class IGraphicsAdapterSelector;
    enum class GraphicsMemoryRegionAllocationFlags : uint32_t;
    class GraphicsMemoryBlockCollection;
    class GraphicsMemoryBlock;
    enum class GraphicsMemoryRegionAllocationFlags : uint32_t;
    class GraphicsMemoryBudget;
    enum class GraphicsBufferKind : uint32_t;
    enum class TexelFormat : uint32_t;
    class GraphicsSurfaceContext;
}

// Graphics types
// clang-format off

#include "ShaderProgramKind.h"
#include "EcsDefaultRenderingComponentTypes.h"
#include "EcsDefaultRenderingSystem.h"
#include "GraphicsAdapter.h"
#include "GraphicsDeviceObject.h"
#include "GraphicsContext.h"
#include "GraphicsFence.h"
#include "GraphicsMemoryAllocatorSettings.h"
#include "GraphicsMemoryRegion.h"
#include "IGraphicsMemoryRegionCollection.h"
#include "GraphicsMemoryRegionAllocationFlags.h"
#include "TexelFormat.h"
#include "GraphicsMemoryAllocator.h"
#include "GraphicsMemoryBlockCollection.h"
#include "GraphicsMemoryBudget.h"
#include "GraphicsMemoryBlock.h"
#include "GraphicsResourceAccess.h"
#include "GraphicsSurfaceKind.h"
#include "GraphicsTextureKind.h"
#include "IGraphicsSurface.h"
#include "GraphicsSurfaceContext.h"
#include "GraphicsDevice.h"
#include "GraphicsResource.h"
#include "GraphicsBufferKind.h"
#include "GraphicsBuffer.h"
#include "GraphicsTexture.h"
#include "IGraphicsAdapterSelector.h"
#include "ShaderProgram.h"
#include "GraphicsPipeline.h"
#include "GraphicsPipelineSignature.h"
#include "GraphicsPrimitive.h"
#include "GraphicsProvider.h"
#include "GraphicsPipelineInput.h"
#include "GraphicsPipelineInputElement.h"
#include "GraphicsPipelineInputElementKind.h"
#include "GraphicsPipelineResource.h"
#include "GraphicsPipelineResourceKind.h"
#include "ShaderProgramVisibility.h"

// clang-format on

#endif // !NOVELRT_EXPERIMENTAL_GRAPHICS_H

max_stars_repo_name: Exadon/NovelRT

max_issues_repo_name: BanalityOfSeeking/NovelRT

max_forks_repo_name: BanalityOfSeeking/NovelRT


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_hpmv (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0627557f, -0.839323f, -0.0877262f, -0.169208f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1118) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1118) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0627557f, -0.839323f, -0.0877262f, -0.169208f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1119) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1119) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0037603f, -0.816761f, -0.0392456f, -0.121154f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1120) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1120) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0037603f, -0.816761f, -0.0392456f, -0.121154f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1121) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1121) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0627557f, -0.839323f, -0.0877262f, -0.169208f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1122) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1122) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0627557f, -0.839323f, -0.0877262f, -0.169208f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1123) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1123) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0037603f, -0.816761f, -0.0392456f, -0.121154f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1124) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1124) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   int N = 2;
   float A[] = { 0.339f, -0.102f, 0.908f, 0.097f, -0.808f, 0.236f };
   float X[] = { 0.993f, -0.502f, -0.653f, 0.796f };
   int incX = -1;
   float Y[] = { -0.35f, 0.339f, -0.269f, -0.122f };
   int incY = -1;
   float y_expected[] = { -0.0037603f, -0.816761f, -0.0392456f, -0.121154f };
   cblas_chpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "chpmv(case 1125) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "chpmv(case 1125) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.745218, -0.60699, -0.37301, -0.983688 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1126) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1126) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.745218, -0.60699, -0.37301, -0.983688 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1127) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1127) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.646956, -0.542012, -0.282168, -0.912668 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1128) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1128) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.646956, -0.542012, -0.282168, -0.912668 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1129) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1129) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.745218, -0.60699, -0.37301, -0.983688 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1130) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1130) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.745218, -0.60699, -0.37301, -0.983688 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1131) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1131) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.646956, -0.542012, -0.282168, -0.912668 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1132) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1132) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   double alpha[2] = {1, 0};
   double beta[2] = {0, 1};
   int N = 2;
   double A[] = { 0.543, -0.737, 0.281, -0.053, -0.098, -0.826 };
   double X[] = { 0.67, -0.857, -0.613, -0.927 };
   int incX = -1;
   double Y[] = { -0.398, -0.934, -0.204, 0.183 };
   int incY = -1;
   double y_expected[] = { 0.646956, -0.542012, -0.282168, -0.912668 };
   cblas_zhpmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zhpmv(case 1133) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zhpmv(case 1133) imag");
     };
   };
  };

}

ext: c


text:
/* WidiPadi is a motif discovery tool for DNA sequences */
/* Copyright (C) 2006 Bioinformatics Centre */
/* Author: Eivind Valen and Ole Winther */

/* 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 <stdlib.h>
#include <stdio.h>
#include <getopt.h>
#include <gsl/gsl_sf_gamma.h>

#include "error.h"
#include "sesa.h"

#define USAGE "Usage: widipadi [OPTIONS] SEQUENCES"

enum bool {false, true};

/* FUNCTION PROTOTYPES */
PSSM hits_to_count_matrix(SESA sesa, struct HitTable *hits, const double threshold, int length, int alphabetSize);
int compare_hit_scores(const void *hit1, const void *hit2);
int parse_options(int argc, char **argv);
void print_help();

/* Acceptance */
int accept_all(SESA sesa, struct HitTable *hits);
int acceptIf_over_300_hits(SESA sesa, struct HitTable *hits);
int acceptIf_required_presence(SESA sesa, struct HitTable *hits);

/* DEFAULTS */
unsigned int length = 8;
unsigned int alphabetSize = 4;
unsigned int mismatches = 1;
int (*acceptanceFunc)(SESA, struct HitTable *) = acceptIf_required_presence;
double required_presence = 0.7;
unsigned int exhaustive = false;

/**
 * Prints a help message to stdout.
 */
void print_help() {
  printf(USAGE);
  printf("\nOptions:\n");
/*   printf("  -a                 Size of alphabet\n"); */
  printf("  -h                 Prints this help and exits\n");
  printf("  -l LENGTH          Length of WM\n");
  printf("  -m MISMATCHES      Number of mismatches\n");
  printf("  -r RATIO           Sets the ratio of sequences where the PSSM must match\n");
  printf("  -x                 Do exhaustive search\n");
  printf("\n");
}

/**
 * Parses command line options
 */
int parse_options(int argc, char **argv) {
  int count = 1;
  int c;

  while ((c = getopt(argc, argv, "hl:m:r:x")) != -1) {
    count++;
    switch (c) {
    case 'h':
      print_help();
      exit(EXIT_SUCCESS);
    case 'l':
      length = (unsigned int) strtol(optarg, (char **) NULL, 0);
      count++;
      break;
    case 'm':
      mismatches = (unsigned int) strtol(optarg, (char **) NULL, 0);
      count++;
      break;
    case 'r':
      acceptanceFunc = acceptIf_required_presence;
      required_presence = strtod(optarg, (char **) NULL);
      count++;
      break;
    case 'x':
      exhaustive = true;
      break;
    default:
      fprintf(stderr, USAGE);
      exit(EXIT_FAILURE);
    }
  }

  return count;
}

/*
 * Converts all hits above scoring threshold to a count matrix
 *
 * Only 0th order
 */
PSSM hits_to_count_matrix(SESA sesa, struct HitTable *hits, const double threshold, int length, int alphabetSize) {
  PSSM pssm = initMatrix(0, length, alphabetSize);
  ESA esa = sesa->esa;
  struct HitEntry *hit;
  int i, j;

  i = pssm->offsets[pssm->length];
  while (i--) pssm->counts[i] = 0;

  for (i = 0; i < hits->nScores; i++) {
    hit = &hits->pScores[i];

    if (hit->score < threshold)
      break;

    for (j = 0; j < length; j++) {
      int ltr = esa->pStr[hit->position + j];
      /*       fprintf(stderr, "%c", ind2chr(ltr)); */
      pssm->counts[ pssm->offsets[j] + ltr ]++;
    }
/*     fprintf(stderr, " (%f)\n", hit->score); */
  }

  return pssm;
}

/**
 * Comparing function for sort. Sorts on scores from low to high.
 */
int compare_hit_scores(const void *hit1, const void *hit2) {
  return ((struct HitEntry *)hit2)->score - ((struct HitEntry *)hit1)->score;
}

/**
 * Accepts all PSSMs (for testing)
 */
int accept_all(SESA sesa, struct HitTable *hits) {
  return true;
}

/**
 * Accepts PSSMs with more than 300 hits (for testing)
 */
int acceptIf_over_300_hits(SESA sesa, struct HitTable *hits) {
  return (hits->nScores > 300 ? true : false);
}

/**
 * Accepts PSSMs with hits present in at least "required_presence" of
 * the sequences
 */
int acceptIf_required_presence(SESA sesa, struct HitTable *hits) {
  int present[sesa->nSeq];
  int i, j;

  i = sesa->nSeq;
  while(i--) present[i] = 0;

  i = hits->nScores;
  while (i--) present[ sesa->seq[ hits->pScores[i].position ] ]++;

  j = 0;
  i = sesa->nSeq;
  while (i--) if (present[i]) j++;

  return ( ((double) j / (double) sesa->nSeq)  > required_presence ? true : false);
}

/*
 * Log likelihood scoring
 */
double log_likelihood(ESA esa, PSSM pssm, int *count_bg_all, double *pcount, double *pcount_bg) {
  int pcount_tot = 0, pcount_bg_tot = 0, count_tot = 0, count_bg_tot = 0;
  int alphlen = pssm->alphabetSize;
  int count_bg[alphlen];
  int ccount, ltr, pos, ord, i;
  double ll = 0.0;

  i = alphlen;
  while (i--) count_bg[i] = count_bg_all[i];

  pos = pssm->length;
  while (pos--) {
    ltr = alphlen;
    while(ltr--) {
      ccount = pssm->counts[pos * alphlen + ltr];
      count_bg[ltr] -= ccount ;
      ll += gsl_sf_lngamma( ccount + pcount[ltr] );
    }
  }

  ltr = alphlen;
  while (ltr--) {
    pcount_tot += pcount[ltr];
    pcount_bg_tot += pcount_bg[ltr];
    count_tot += pssm->counts[ltr];
    count_bg_tot += count_bg[ltr];
    ll += gsl_sf_lngamma( count_bg[ltr] + pcount_bg[ltr] );
  }

  ll -= alphlen * gsl_sf_lngamma( count_tot + pcount_tot ) + gsl_sf_lngamma(count_bg_tot + pcount_bg_tot );

  return ll;
}

/*
 * Log likelihood scoring
 */
double higher_log_likelihood(ESA esa, PSSM pssm, double *pcount, double *pcount_bg) {
  int pos, lcpi, skip, curSuf, i, j, p;
  int numSufs = getSize(esa);
  int curEsaIndex = 0;
  int order = 2;
  int wcount[order+1];

  i = order+1;
  while (i--) wcount[i] = 0;

  /* Iterate all suffixes */
  while(curEsaIndex < numSufs){
    lcpi = getLcp(esa, curEsaIndex);
    curSuf = getSuf(esa, curEsaIndex);
    skip = getSkip(esa, curEsaIndex);

    /* We can skip if the lcp is longer than the order */
    if (lcpi > order) {

      /* FIX: HACK */
      if (skip == 0) {
	curEsaIndex++;
	continue;
      }

      /* Count all the skipped words */
      i = order+1;
      while (i--) wcount[i] += skip - curEsaIndex;

      curEsaIndex = skip;
      continue;

    } else {

      pos = pssm->length;
      while (pos--) { }

      /* Print the hits */
      p = 1;
      for (j = 0; j <= order; j++)
	if (ind2chr(esa->pStr[getSuf(esa, curEsaIndex-1)+j]) == '-')
	  p = 0;

      if (p) {
	for (i = 0; i <= order; i++) fprintf(stderr, "%c", ind2chr(esa->pStr[getSuf(esa, curEsaIndex-1)+i]));
	fprintf(stderr, ": %5i\n", wcount[2]);
      }

      wcount[order] = 1;

      /* For higher order models the first columns have fewer letters */
      i = order;
      while (i--) {
	if (i >= lcpi) {

	  p = 1;
	  for (j = 0; j <= i; j++)
	    if (ind2chr(esa->pStr[getSuf(esa, curEsaIndex-1)+j]) == '-')
		p = 0;

	  if (p) {
	    for (j = 0; j <= i; j++) fprintf(stderr, "%c", ind2chr(esa->pStr[getSuf(esa, curEsaIndex-1)+j]));
	    for (j = i; j < order; j++) fprintf(stderr, " ");
	    fprintf(stderr, ": %5i\n", wcount[i]);
	  }

	  wcount[i] = 1;
	}
      }

      curEsaIndex++;
    }

  }

  exit(0);
}

void print_esa (ESA esa) {
  int numSufs = getSize(esa);
  int i,j;

  for (i = 0; i < numSufs; i++) {
    fprintf(stderr, "%4i: SUF(%3i) LCP(%3i) SKIP(%4i) ", i, getSuf(esa, i), getLcp(esa, i), getSkip(esa, i));

    for (j = 0; j < 8; j++) {
      fprintf(stderr, "%c", ind2chr(esa->pStr[getSuf(esa, i)+j]));
    }
    fprintf(stderr, "\n");
  }
}

int main(int argc, char **argv) {
  PSSMScoreTable results;
  struct HitTable *hits = NULL;
  PSSM pssm, count;
  SESA sesa;
  unsigned int arg_count, j;
  int i;

  arg_count = parse_options(argc, argv);

  printf("Reading file: %s  (%i)\n", argv[arg_count], argc - arg_count);
  if (!argv[arg_count]) {
    fprintf(stderr, "Need an input file\n");
    exit(EXIT_FAILURE);
  }

  sesa = build_SESA(argc - arg_count, &argv[arg_count], "ACGT", "N", 0);

  if (!sesa) {
    fprintf(stderr, "%s\n", getError());
    exit(EXIT_FAILURE);
  }

  printf("Searching %s for motifs of length %i with %i mismatch%s\n", (exhaustive ? "exhaustively" : "exclusively"), length, mismatches, (mismatches > 1 ? "es" : ""));

  if (exhaustive) {
    /* Exhaustive search searches for all words og length n.  */
    results = SESA_exhaustive_search(sesa, length, mismatches, alphabetSize, acceptanceFunc);
  } else {
    /* Exclusive search only searches for PSSMs representing words present in the sequence.  */
    results = SESA_exclusive_search(sesa, length, mismatches, alphabetSize, acceptanceFunc);
  }

  /**
   * Both the above methods return a list of integers that uniquely
   * identifies pssms.  We do not return the hits of the PSSMs even
   * though we have searched with them since this might be very(!)
   * memory demanding, depending on criteria used for acceptance.
   *
   * If we know the criteria is strict we might considering altering
   * this to avoid re-searching and thus speed up the process.
   */

  pssm = initMatrix(0, length, alphabetSize);

  for (i = 0; i < results->nPssm; i++) {
    /* Sets the pssm to the values corresponding to the ID*/
    set_mismatch_scores(pssm, results->pssmNr[i]);
    calcAndSetThresholds(pssm, -((double)mismatches)-0.5);

    /* Search with the PSSM again */
    init_hittable();                 /* Creates a hittable */
    hits = SESA_search(sesa, pssm);  /* Searches again */

    /* Sort the hits */
    qsort((void*)hits->pScores, hits->nScores, sizeof(struct HitEntry), compare_hit_scores);

    printf("\n\n");
    print_pssm(pssm);
    for (j = 0; j <= mismatches; j++) {
      printf("Count matrix for %i mismatches\n", j);
      count = hits_to_count_matrix(sesa, hits, -((double)j)-0.5, pssm->length, pssm->alphabetSize);

      double ost1[] = {1.0, 1.0, 1.0, 1.0};
      double ost2[] = {1.0, 1.0, 1.0, 1.0};

/*       double foo = log_likelihood(count, sesa->bg, ost1, ost2); */

      double foo = higher_log_likelihood(sesa->esa, count, ost1, ost2);

/*       fprintf(stderr, "LL: %g\n", foo); */
/*       print_counts(count);  */
      printf("\n\n");
    }

    release_hits(hits); /* Frees the allocated hittable*/
  }

  exit(EXIT_SUCCESS);
}

max_stars_repo_name: kipkurui/gimmemotifs

max_issues_repo_name: kipkurui/gimmemotifs

max_forks_repo_name: kipkurui/gimmemotifs


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

#ifndef __GSL_VECTOR_DOUBLE_H__
#define __GSL_VECTOR_DOUBLE_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

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

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  double *data;
  gsl_block *block;
  int owner;
}
gsl_vector;

typedef struct
{
  gsl_vector vector;
} _gsl_vector_view;

typedef _gsl_vector_view gsl_vector_view;

typedef struct
{
  gsl_vector vector;
} _gsl_vector_const_view;

typedef const _gsl_vector_const_view gsl_vector_const_view;

/* Allocation */

GSL_FUN gsl_vector *gsl_vector_alloc (const size_t n);
GSL_FUN gsl_vector *gsl_vector_calloc (const size_t n);

GSL_FUN gsl_vector *gsl_vector_alloc_from_block (gsl_block * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector *gsl_vector_alloc_from_vector (gsl_vector * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_free (gsl_vector * v);

/* Views */

GSL_FUN _gsl_vector_view
gsl_vector_view_array (double *v, size_t n);

GSL_FUN _gsl_vector_view
gsl_vector_view_array_with_stride (double *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_const_view
gsl_vector_const_view_array (const double *v, size_t n);

GSL_FUN _gsl_vector_const_view
gsl_vector_const_view_array_with_stride (const double *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_view
gsl_vector_subvector (gsl_vector *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_view
gsl_vector_subvector_with_stride (gsl_vector *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_const_view
gsl_vector_const_subvector (const gsl_vector *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_const_view
gsl_vector_const_subvector_with_stride (const gsl_vector *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_set_zero (gsl_vector * v);
GSL_FUN void gsl_vector_set_all (gsl_vector * v, double x);
GSL_FUN int gsl_vector_set_basis (gsl_vector * v, size_t i);

GSL_FUN int gsl_vector_fread (FILE * stream, gsl_vector * v);
GSL_FUN int gsl_vector_fwrite (FILE * stream, const gsl_vector * v);
GSL_FUN int gsl_vector_fscanf (FILE * stream, gsl_vector * v);
GSL_FUN int gsl_vector_fprintf (FILE * stream, const gsl_vector * v,
                              const char *format);

GSL_FUN int gsl_vector_memcpy (gsl_vector * dest, const gsl_vector * src);

GSL_FUN int gsl_vector_reverse (gsl_vector * v);

GSL_FUN int gsl_vector_swap (gsl_vector * v, gsl_vector * w);
GSL_FUN int gsl_vector_swap_elements (gsl_vector * v, const size_t i, const size_t j);

GSL_FUN double gsl_vector_max (const gsl_vector * v);
GSL_FUN double gsl_vector_min (const gsl_vector * v);
GSL_FUN void gsl_vector_minmax (const gsl_vector * v, double * min_out, double * max_out);

GSL_FUN size_t gsl_vector_max_index (const gsl_vector * v);
GSL_FUN size_t gsl_vector_min_index (const gsl_vector * v);
GSL_FUN void gsl_vector_minmax_index (const gsl_vector * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_add (gsl_vector * a, const gsl_vector * b);
GSL_FUN int gsl_vector_sub (gsl_vector * a, const gsl_vector * b);
GSL_FUN int gsl_vector_mul (gsl_vector * a, const gsl_vector * b);
GSL_FUN int gsl_vector_div (gsl_vector * a, const gsl_vector * b);
GSL_FUN int gsl_vector_scale (gsl_vector * a, const double x);
GSL_FUN int gsl_vector_add_constant (gsl_vector * a, const double x);
GSL_FUN int gsl_vector_axpby (const double alpha, const gsl_vector * x, const double beta, gsl_vector * y);
GSL_FUN double gsl_vector_sum (const gsl_vector * a);

GSL_FUN int gsl_vector_equal (const gsl_vector * u,
                            const gsl_vector * v);

GSL_FUN int gsl_vector_isnull (const gsl_vector * v);
GSL_FUN int gsl_vector_ispos (const gsl_vector * v);
GSL_FUN int gsl_vector_isneg (const gsl_vector * v);
GSL_FUN int gsl_vector_isnonneg (const gsl_vector * v);

GSL_FUN INLINE_DECL double gsl_vector_get (const gsl_vector * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_set (gsl_vector * v, const size_t i, double x);
GSL_FUN INLINE_DECL double * gsl_vector_ptr (gsl_vector * v, const size_t i);
GSL_FUN INLINE_DECL const double * gsl_vector_const_ptr (const gsl_vector * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
double
gsl_vector_get (const gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_set (gsl_vector * v, const size_t i, double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
double *
gsl_vector_ptr (gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (double *) (v->data + i * v->stride);
}

INLINE_FUN
const double *
gsl_vector_const_ptr (const gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const double *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_DOUBLE_H__ */

ext: h

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"

void
cblas_dger (const enum CBLAS_ORDER order, const int M, const int N,
            const double alpha, const double *X, const int incX,
            const double *Y, const int incY, double *A, const int lda)
{
#define BASE double
#include "source_ger.h"
#undef BASE
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
 * Subroutines that deal with the aXe configuration file.
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include "aper_conf.h"
#include "disp_conf.h"
#include "aXe_grism.h"
#include "spc_cfg.h"
#include "disp_conf.h"
#include "aXe_utils.h"

#define MAX(x,y) (((x)>(y))?(x):(y))

/**
 *  Function: get_aperture_descriptor
 *  Read a configuration file and populate an aperture_conf structure
 *  containing a description of the optical set up.
 *  This function populates the geometrical description of all the beams in
 *  all the apertures listed in the configuration file.
 *
 *  Parameters:
 *  @param actinfo  - the header structure with the column names
 *  @param filename - the complete filename to the axe configuration file
 *
 *  Returns:
 *  @return config - the configuration structure created
 */
aperture_conf *
get_aperture_descriptor (char *filename)
{
  char beam[MAXCHAR] = "\0";
  aperture_conf *config;
  int i, ix;
  gsl_vector *v;

  struct CfgStrings AperConfig[] = {
    {"INSTRUMENT", NULL},
    {"CAMERA", NULL},
    {"SCIENCE_EXT", NULL},
    {"ERRORS_EXT", NULL},
    {"DQ_EXT",NULL},
    {"DQMASK",NULL},
    {"EXPTIME",NULL},
    {"GAIN",NULL},
    {"OPTKEY1",NULL},
    {"OPTVAL1",NULL},
    {"FFNAME",NULL},
    {"REFX",NULL},
    {"REFY",NULL},
    /* aXe1.4:*/
    {"DRZRESOLA",NULL},
    {"DRZSCALE",NULL},
    {"DRZLAMB0",NULL},
    {"DRZXINI",NULL},
    {"DRZPFRAC",NULL},
    {"DRZKERNEL",NULL},
    /* aXe1.5: */
    {"PSFCOEFFS",NULL},
    {"PSFRANGE",NULL},
    {"RDNOISE",NULL},
    /* NICMOS HLA */
    {"IPIXFUNCTION",NULL},
    {"POBJSIZE",NULL},
    {"SMFACTOR",NULL},

    {NULL, NULL},
    {NULL, NULL}                /* array terminator. REQUIRED !!! */
  };

  struct CfgStrings BeamConfig[] = {
    {NULL, NULL},
    {NULL, NULL}
  };

  AperConfig[25].name = beam;
  BeamConfig[0].name = beam;

  config = malloc (sizeof (aperture_conf));
  if (config == NULL)
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "Could not allocate memory for aperture configuration");
    }

  CfgRead (filename, AperConfig);

  /* Initialize optional keywords to default values*/
  sprintf(config->instrument,"None");
  sprintf(config->camera,"None");
  sprintf(config->science_ext,"None");
  sprintf(config->errors_ext,"None");
  sprintf(config->dq_ext,"None");
  sprintf(config->exptimekey,"None");
  sprintf(config->gainkey,"None");
  config->dqmask=0;
  sprintf(config->optkey1,"None");
  sprintf(config->optval1,"None");
  sprintf(config->FFname,"None");
  sprintf(config->IPIXfunc,"None");
  sprintf(config->drz_kernel,"square");
  config->refx       = 0;
  config->refy       = 0;
  config->drz_resol  = 0.0;
  config->drz_scale  = 0.0;
  config->drz_lamb0  = 0.0;
  config->drz_xstart = 15.0;
  config->drz_pfrac  = 1.0;
  config->psfcoeffs  = NULL;
  config->psfrange   = NULL;
  config->rdnoise    = 0.0;
  config->pobjsize   = -1.0;
  config->smfactor   = -1.0;

  for (ix = 0; ix < 26; ix++)
    {

      /* Name of the instrument */
      if (!strcmp (AperConfig[ix].name, "INSTRUMENT"))
        {
          if (AperConfig[ix].data == NULL)
            aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                         "INSTRUMENT tag was not read from %s\n",
                         filename);
          sprintf (config->instrument, "%s", AperConfig[ix].data);
        }

      /* Name of the camera */
      if (!strcmp (AperConfig[ix].name, "CAMERA"))
        {
          if (AperConfig[ix].data == NULL)
            aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                         "camera tag was not read from %s\n",
                         filename);
          sprintf (config->camera, "%s", AperConfig[ix].data);
        }

      /* Name of the Science data extension */
      if (!strcmp (AperConfig[ix].name, "SCIENCE_EXT"))
        {
          if (AperConfig[ix].data != NULL)
            sprintf (config->science_ext, "%s", AperConfig[ix].data);
        }

      /* Name of the optional Error FITS extension */
      if (!strcmp (AperConfig[ix].name, "ERRORS_EXT"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->errors_ext, "%s", AperConfig[ix].data);
            }
        }

      /* Name of the optional Data Quality FITS extension */
      if (!strcmp (AperConfig[ix].name, "DQ_EXT"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->dq_ext, "%s", AperConfig[ix].data);
            }
        }

      /* Name of the optional EXPTIME FITS keyword */
      if (!strcmp (AperConfig[ix].name, "EXPTIME"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->exptimekey, "%s", AperConfig[ix].data);
            }
        }

      /* Red in the readout noise */
      if (!strcmp (AperConfig[ix].name, "RDNOISE"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->rdnoise = atof(AperConfig[ix].data);
            }
        }

      /* Name of the optional GAIN FITS keyword */
      if (!strcmp (AperConfig[ix].name, "GAIN"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->gainkey, "%s", AperConfig[ix].data);
            }
        }

      /* Gett the Data Quality Mask */
      if (!strcmp (AperConfig[ix].name, "DQMASK"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->dqmask = atoi(AperConfig[ix].data);
            }
        }

      /* Optional OPTional 1st FITS keyword selection */
      if (!strcmp (AperConfig[ix].name, "OPTKEY1"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->optkey1, "%s", AperConfig[ix].data);
            }
        }

      /* Optional OPTional 1st FITS keyword value */
      if (!strcmp (AperConfig[ix].name, "OPTVAL1"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->optval1, "%s", AperConfig[ix].data);
            }
        }

      /* Optional FFname keyword value */
      if (!strcmp (AperConfig[ix].name, "FFNAME"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->FFname, "%s", AperConfig[ix].data);
            }
        }

      /* Optional intra pixel correction keyword value */
      if (!strcmp (AperConfig[ix].name, "IPIXFUNCTION"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->IPIXfunc, "%s", AperConfig[ix].data);
            }
        }

      /* Optional REFX keyword value */
      if (!strcmp (AperConfig[ix].name, "REFX"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->refx = atoi(AperConfig[ix].data);
            }
        }

      /* Optional REFY keyword value */
      if (!strcmp (AperConfig[ix].name, "REFY"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->refy = atoi(AperConfig[ix].data);
            }
        }
      if (!strcmp (AperConfig[ix].name, "DRZRESOLA"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->drz_resol = atof(AperConfig[ix].data);
            }
        }
      if (!strcmp (AperConfig[ix].name, "DRZSCALE"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->drz_scale = atof(AperConfig[ix].data);
            }
        }
      if (!strcmp (AperConfig[ix].name, "DRZLAMB0"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->drz_lamb0 = atof(AperConfig[ix].data);
            }
        }
      if (!strcmp (AperConfig[ix].name, "DRZXINI"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->drz_xstart = atof(AperConfig[ix].data);
            }
        }
      if (!strcmp (AperConfig[ix].name, "DRZPFRAC"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->drz_pfrac = atof(AperConfig[ix].data);
            }
        }
      if (!strcmp (AperConfig[ix].name, "DRZKERNEL"))
        {
          if (AperConfig[ix].data != NULL)
            {
              sprintf (config->drz_kernel, "%s", AperConfig[ix].data);
            }
        }

      if (!strcmp (AperConfig[ix].name, "PSFCOEFFS"))
        {
          if (AperConfig[ix].data != NULL)
            {
              //              sprintf (config->drz_kernel, "%s", AperConfig[ix].data);
              v = string_to_gsl_array (AperConfig[ix].data);
              config->psfcoeffs = v;
            }
        }

      if (!strcmp (AperConfig[ix].name, "PSFRANGE"))
        {
          if (AperConfig[ix].data != NULL)
            {
              v = string_to_gsl_array (AperConfig[ix].data);
              if (v->size != 2)
                aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                             "configuration file: two items for PSFRANGE, lambda_min and lambda_max must be given, not %i\n", v->size);
              config->psfrange = v;
            }
        }

      // read in the size of point-like objects
      if (!strcmp (AperConfig[ix].name, "POBJSIZE"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->pobjsize = atof(AperConfig[ix].data);
            }
        }

      // read in the adjustment for smoothed flux conversion
      if (!strcmp (AperConfig[ix].name, "SMFACTOR"))
        {
          if (AperConfig[ix].data != NULL)
            {
              config->smfactor = atof(AperConfig[ix].data);
            }
        }

    }

  // release memory
  i=0;
  while(AperConfig[i].name!=NULL)
    free(AperConfig[i++].data);

  /* Looking for up to MAX_BEAMS beams */
  config->nbeams = 0;
  for (i = 0; i < MAX_BEAMS; i++)
    {
      sprintf (beam, "BEAM%c", BEAM (i));
      CfgRead (filename, BeamConfig);
      if (BeamConfig[0].data != NULL)
        {
          v = string_to_gsl_array (BeamConfig[0].data);
          if (v->size != 2)
            {
              aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                           "tag %s does not have the proper format in file %s.\n"
                           "Should be of the form %s [integer] [integer]\n",
                           BeamConfig[0].data, filename,beam);
              //                           beam, filename,beam);
            }
          (config->beam[config->nbeams]).offset.dx0 =
            (int) gsl_vector_get (v, 0);
          (config->beam[config->nbeams]).offset.dx1 =
            (int) gsl_vector_get (v, 1);

          // release the memory
          gsl_vector_free(v);

          (config->beam[config->nbeams]).ID = i;

          // release memory
          ix=0;
          while(BeamConfig[ix].name!=NULL)
            free(BeamConfig[ix++].data);
          BeamConfig[0].data = NULL;

          (config->beam[config->nbeams]).mmag_extract = get_beam_mmag_extract (filename, i);
          (config->beam[config->nbeams]).mmag_mark = get_beam_mmag_mark (filename, i);

          sprintf (beam, "PSF_OFFSET_%c", BEAM (i));
          CfgRead (filename, BeamConfig);
          if (BeamConfig[0].data != NULL)
            {
              (config->beam[config->nbeams]).psf_offset = atof(BeamConfig[0].data);
            }
          else
            {
              (config->beam[config->nbeams]).psf_offset = 0.0;
            }

          // release memory
          ix=0;
          while(BeamConfig[ix].name!=NULL)
            free(BeamConfig[ix++].data);

          BeamConfig[0].data = NULL;

          config->nbeams = config->nbeams + 1;
        }
    }

  if (config->nbeams == 0)
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "No beam was found in %s. At least one beam must be found.\n",
                   filename);
    }
  return config;
}

/**
 * Function: aperture_conf_fprintf
 * Function to output the content of an aperture_conf structure.
 *
 * Parameters:
 * @param file - a pointer to an output stream or file
 * @param conf - a pointer to an existing aperture_conf structure
 */
void
aperture_conf_fprintf (FILE * file, aperture_conf * conf)
{
     int i = 0, id;

     fprintf (file, "Aperture_conf: INSTRUMENT: %s\n", conf->instrument);
     fprintf (file, "Aperture_conf: NBEAMS: %d\n", conf->nbeams);
     for (i = 0; i < conf->nbeams; i++)
       {
         id = conf->beam[i].ID;
         fprintf (file, "Aperture_conf: BEAM%c: %d %d\n", BEAM (id),
                  conf->beam[i].offset.dx0, conf->beam[i].offset.dx1);
       }
     fprintf (file, "FITS SCI extension name: %s, number: %d\n",
              conf->science_ext, conf->science_numext);
     fprintf (file, "FITS ERR extension name: %s, number: %d\n",
              conf->errors_ext, conf->errors_numext);
     fprintf (file, "FITS DQ extension name: %s, number: %d\n",
              conf->dq_ext, conf->dq_numext);
     fprintf (file, "DQ mask: %d\n",conf->dqmask);
     fprintf (file, "FF filename: %s\n",conf->FFname);
     fprintf (file, "Refx: %d\n",conf->refx);
     fprintf (file, "Refy: %d\n",conf->refy);
     fprintf (file, "DRZRESO: %10.3f\n",conf->drz_resol);
     fprintf (file, "DRZSCALE: %10.3f\n",conf->drz_scale);
     fprintf (file, "DRZLAMB0: %10.3f\n",conf->drz_lamb0);
     fprintf (file, "DRZXINI: %10.3f\n",conf->drz_xstart);
     fprintf (file, "DRZPFRAC: %10.3f\n",conf->drz_pfrac);
     fprintf (file, "DRZKERNEL: %s\n",conf->drz_kernel);

}

/**
 * Function: get_extension_numbers
 * This function assigns HDU numbers to an aperture_conf
 * structues after looking for the location of named
 * extension in a given FITS file. Optional values not
 * found are set to -1.
 *
 * Parameters:
 * @param filename - a pointer to a string containing the name
 *                   of an existing FITS file
 * @param conf     - a pointer to a populated aperture_conf structure
 * @param keyword  - a pointer to a string containing the name of
 *                   an optional keyword to read
 * @param keyval   - a pointer to a string containing the value of
 *                   the optional keyword to look for
 */
void
get_extension_numbers(char filename[], aperture_conf * conf,
                      char keyword[],char keyval[])
{
  int extver = -1;

  conf->science_numext = get_hdunum_from_hduname(filename,
                                                 conf->science_ext,keyword
                                                 ,keyval,extver);
  if (conf->science_numext > -1)
    extver = (int)get_float_from_keyword(filename, conf->science_numext, "EXTVER");

  if (strcmp(conf->errors_ext,"None"))
  {
    conf->errors_numext = get_hdunum_from_hduname(filename,
                                                  conf->errors_ext,keyword,
                                                  keyval,extver);
  } else {
    conf->errors_numext = -1;
  }
  if (strcmp(conf->dq_ext,"None")) {
    conf->dq_numext = get_hdunum_from_hduname(filename,conf->dq_ext,
                                              keyword,keyval,extver);
  } else {
    conf->dq_numext = -1;
  }
}

/**
 * Function: free_aperture_conf
 * The function releases all memory allocated
 * in a configuration structure.
 *
 * Parameters:
 * @param conf - the configuration structure
 */
void
free_aperture_conf(aperture_conf * conf)
{
  if (conf->psfcoeffs)
    gsl_vector_free(conf->psfcoeffs);
  if (conf->psfrange)
    gsl_vector_free(conf->psfrange);
  free(conf);
  conf = NULL;
}

/**
 * Function: get_psf_offset
 * The function extracts and returns the psf-offset
 * stored in a given configuration structure and
 * a given beam
 *
 * Parameters:
 * @param conf - the configuration structure
 * @param beam - the beam
 *
 * Returns:
 * @return psf_offset - the psf-offset for the beam
 */
double
get_psf_offset(aperture_conf * conf, const beam actbeam)
{
  // initialize the offset
  double psf_offset=0.0;

  int i=0;

  // go over all beams in the configuration structure
  for (i = 0; i < conf->nbeams; i++)
    {
      // check whether the current beam is the
      // right one
      if (conf->beam[i].ID == actbeam.ID)
        // return the offset
        return conf->beam[i].psf_offset;
    }

  // return the default offset
  return psf_offset;
}

/**
 * Function: get_psf_offset
 * The function extracts and returns the psf-offset
 * stored in a given configuration structure and
 * a given beam
 *
 * Parameters:
 * @param conf - the configuration structure
 * @param beam - the beam
 *
 * Returns:
 * @return psf_offset - the psf-offset for the beam
 */
double
get_max_offset(aperture_conf * conf)
{
  // initialize the offset
  double max_offset=0.0;

  int i=0;

  // go over all beams in the configuration structure
  for (i = 0; i < conf->nbeams; i++)
    {
      // check whether the current beam is the
      // right one
      max_offset = MAX(conf->beam[i].psf_offset, max_offset);
    }

  // return the default offset
  return max_offset;
}

/**
 * Function: extend_config_beams
 * The function extends the size of the beams in a
 * configuration structure.
 *
 * Parameters:
 * @param conf - the configuration structure
 */
void
extend_config_beams(aperture_conf *conf)
{
  int index;

  // go over all beams in the configuration structure
  for (index=0; index < conf->nbeams; index++)
    {
      // lower the left border
      conf->beam[index].offset.dx0 -= BCK_BEAM_EXT;

      // enhance the right border
      conf->beam[index].offset.dx1 += BCK_BEAM_EXT;
    }
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cdf.h>

#include <QVector>

#include <algorithm>
#include <tuple>

class Gaussian2DFit
{
    struct data
    {
        double* x;
        double* y;
        size_t n;

        double* z;
        size_t width;
        size_t height;
    };

    enum
    {
        numOfPara = 7
    };
    /* model function: A *exp(-a*(x-x0)^2 - 2*b*(x-x0)*(y-y0) - c*(y-y0)^2) + d,
    see more detail https://en.wikipedia.org/wiki/Gaussian_function#Two-dimensional_Gaussian_function*/
private:
    gsl_vector* p0; /* initial fitting parameters:amplitude,center(x0,y0),width(a,b,c),offset*/
    gsl_vector* p; /* fitting parameters */
    gsl_vector* f; /* residual yi-y(xi) */
    gsl_matrix* covar; /*covariance matrix, not yet multiplied by sigma, so is not variance-covariance matrix*/

    gsl_multifit_nlinear_fdf         fdf;
    gsl_multifit_nlinear_parameters  fdf_params;
    gsl_multifit_nlinear_workspace*  work;

    data                             fit_data;

    size_t                           max_iter;
    double                           ptol; /* tolerance on fitting parameter p */
    double                           gtol; /* tolerance on gradient */
    double                           ftol;

    int                              info;/*fiting stop reason, for debug*/
    double                           rss; /*mean of sum of residual square, dof corrected*/

    gsl_vector*                      confid95;
    gsl_vector*                      jacEllipseMajorABC;
    gsl_vector*                      jacEllipseMinorABC;
    double                           ellipseMajor95;
    double                           ellipseMinor95;

public:
    Gaussian2DFit() {};
    Gaussian2DFit(size_t n, size_t width, size_t height, double* datax, double* datay, double* dataz,
        double A0, double x00, double y00, double a0, double b0, double c0, double D0,
        size_t max_iter = 200, double ptol = 1.0e-8, double gtol = 1.0e-8, double ftol = 1.0e-8);
    ~Gaussian2DFit();
    void solve_system();
    void set_data(double* datay);
    void set_data(size_t n, size_t width, size_t height, double* datax, double* datay, double* dataz, bool free = true);
    void set_initialP(double A0, double x00, double y00, double a0, double b0, double c0, double D0);

    int getInfo() const { return info; }

    static double gaussian2d(const double A, const double x0, const double y0,
        const double a, const double b, const double c, const double D,
        const double tx, const double ty);
    QVector<double> calcFittedGaussian();

    const QVector<double> fittedPara() const;
    const QVector<double> confidence95Interval() const;
    const std::pair<double, double> MajorMinor95() const;

private:
    static int func_f(const gsl_vector* p, void* datafit, gsl_vector* f);
    static int func_df(const gsl_vector* p, void* datafit, gsl_matrix* J);
    static int func_fvv(const gsl_vector* p, const gsl_vector* v,
        void* datafit, gsl_vector* fvv);

    void inline errorPropaABC(double a, double b, double c, gsl_vector* ellipseJac,
        bool Major = true);/*for calculating jac of the transformed gaussian arg, now is a classical form of ellipse*/
};

max_stars_repo_name: zzpwahaha/VimbaCamJILA


text:
#ifndef OPENDSS_H
#define OPENDSS_H
#include <petsc.h>
#include <complex.h>
extern PetscErrorCode OpenDSSRunCommand(const char*);
extern PetscErrorCode OpenDSSVsourcesSetPU(PetscScalar);
extern PetscErrorCode OpenDSSVsourcesSetAngleDeg(PetscScalar);
extern PetscErrorCode OpenDSSSolutionGetSolve(PetscInt*);
extern PetscErrorCode OpenDSSCircuitGetTotalPower(double complex**);
#endif

max_stars_repo_name: GMLC-TDC/Use-Cases

max_issues_repo_name: GMLC-TDC/Use-Cases

max_forks_repo_name: GMLC-TDC/Use-Cases


text:
#include <petsc.h>

int main(int argc, char **argv) {
  PetscErrorCode ierr;
  PetscMPIInt    rank;
  PetscInt       i;
  PetscReal      localval, globalsum, x;

  ierr = PetscInitialize(&argc,&argv,NULL,
      "Compute exp in parallel with PETSc.\n\n"); if (ierr) return ierr;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);

  PetscOptionsBegin(PETSC_COMM_WORLD, "", "options for expx", "");
  PetscOptionsReal("-x", "input to exp(x) function", "expx.c", x, &x, PETSC_NULL);
  PetscOptionsEnd();

  // compute  exp(x)  where n = (rank of process) + 1
  localval = 1.0;
  for (i = 2; i < rank+1; i++)
      localval /= i;
  localval *= PetscPowReal(x, rank);

  // sum the contributions over all processes
  ierr = MPI_Allreduce(&localval,&globalsum,1,MPIU_REAL,MPIU_SUM,
      PETSC_COMM_WORLD); CHKERRQ(ierr);

  // output estimate of e and report on work from each process
  ierr = PetscPrintf(PETSC_COMM_WORLD,
      "exp(%f) is about %17.15f\n", x, globalsum); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_SELF,
      "rank %d did %d flops\n", rank, (rank > 0) ? rank-1 : 0);
      CHKERRQ(ierr);
  return PetscFinalize();
}


text:
/** @file */

#ifndef __CCL_F1D_H_INCLUDED__
#define __CCL_F1D_H_INCLUDED__

#include <gsl/gsl_spline.h>

CCL_BEGIN_DECLS

typedef enum ccl_f1d_extrap_t {
  ccl_f1d_extrap_0 = 0,  // No extrapolation
  ccl_f1d_extrap_const = 410,  // Constant extrapolation
  ccl_f1d_extrap_linx_liny = 411,  // Linear x, linear y
  ccl_f1d_extrap_linx_logy = 412,  // Linear x, log y
  ccl_f1d_extrap_logx_liny = 413,  // Log x, linear y
  ccl_f1d_extrap_logx_logy = 414,  // Log x, log y
} ccl_f1d_extrap_t;

/*
 * Spline wrapper
 * Used to take care of evaluations outside the supported range
 */
typedef struct {
  gsl_spline *spline;
  double y0,yf; //Constant values to use beyond interpolation limit
  ccl_f1d_extrap_t extrap_lo_type;
  ccl_f1d_extrap_t extrap_hi_type;
  double x_ini, x_end; //Interpolation limits
  double y_ini, y_end;
  double der_lo;
  double der_hi;
} ccl_f1d_t;

ccl_f1d_t *ccl_f1d_t_new(int n,double *x,double *y,double y0,double yf,
			 ccl_f1d_extrap_t extrap_lo_type,
			 ccl_f1d_extrap_t extrap_hi_type, int *status);

double ccl_f1d_t_eval(ccl_f1d_t *spl,double x);

void ccl_f1d_t_free(ccl_f1d_t *spl);

CCL_END_DECLS

#endif

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
#ifndef AFS_H
#define AFS_H
#include <gsl/gsl_matrix.h>
#include <glib.h>

//afsObject represents a two population coalescent
typedef struct afsObject
{
	int npops,nalleles;
	int aCounts[2]; //this keeps track of popn specific allele counts
	gsl_matrix_int **popMats;
	char *matString;

}
afsObject;

typedef struct afsStateSpace
{
	int nstates;
	afsObject **states;
	//below are only to hasten the process of state
	//space enumeration
	//int **nstateArray;
	//afsObject ****tmpStorage;
	GHashTable *stateHash;
}
afsStateSpace;

//for use in programs like cmc_island_mle
struct island_lik_params{
	int snpNumber, maxSampleSize,nstates,n1,n2, hStates,nParams;
	gsl_vector *weights, *sampleSizeVector,*rates;
	//paramVector for HMM == [N2,thetaState1,thetaState2,...,thetaStaten, m1State1,m2State1,m1State2,m2State2,...,m1Staten,m2Staten] <n is number of hmm states>;
	gsl_vector *paramVector,*paramCIUpper,*paramCILower, *mlParams;
	//For straight MLE the paramVector takes the form [N1,N2,m1,m2]
	gsl_matrix *log_afs, *posts, *transMat,*topol,*expAFS, *visitMat;
	gsl_matrix_int *moveType;
	afsStateSpace *stateSpace;
	struct site *data;
	double *cs_work;
	int *dim1, *dim2; //these pointers ease the construction of sparse matrices
	int nzCount;
	double maxLik;
	int optimized;

}island_lik_params;

#define MAXSTATES 100000000
#define MAXSIZE 20

afsObject *afsObjectNew(int n1, int n2);
void afsObjectFree(afsObject *v);
void afsObjectInit(afsObject *v,int n1, int n2);
void afsObjectPrint(afsObject *v);
afsObject *afsObjectNewFrom(afsObject *v);
afsObject *afsObjectDelta(afsObject *a, afsObject *b);
void afsObjectDeltaPre(afsObject *a, afsObject *b,afsObject *tmp);

int afsObjectCount(afsObject *a);

afsStateSpace *afsStateSpaceNew();
void afsStateSpaceFree(afsStateSpace *a);
int afsStateSpaceContainsState(afsStateSpace *S, afsObject *aState);
void afsStateSpaceRemoveAbsorbing(afsStateSpace *S);
int afsStateSpaceContainsStateTmpStorage(afsStateSpace *S, afsObject *aState);
void afsStateSpaceMapPopn(afsStateSpace *S, int *map);
void afsStateSpaceMapAndReducePopn(afsStateSpace *S, int *map, afsStateSpace *reducedSpace, int *reverseMap);

void afsStateSpaceRemovePopulation(afsStateSpace *S, int popn);
afsStateSpace *afsStateSpaceCopy(afsStateSpace *S);

void mcMatsImportFromFile(const char *fileName,int *nnz,double *topol, int *moveA, int *dim1, int *dim2);

int nDescPop1(afsObject *v);
int nDescPop2(afsObject *v);
int compatibleAfsGivenInitial(afsObject *v, int n1, int n2);
int afsObjectsEqual(afsObject *a, afsObject *b);
void acceptOrRejectProposedState(afsStateSpace *S, afsObject *child, int n1, int n2);
int afsObjectQsortCompare(const void *p1, const void *p2);

afsStateSpace *afsStateSpaceImportFromFile(const char *fileName);

void nonZeroEntries(gsl_matrix_int *m, int **coords, int *count);
void entriesGreaterThan(gsl_matrix_int *m, int **coords, int *count, int cut);
void entriesLessThan(gsl_matrix_int *m, int **coords, int *count, int cut);
int countNegativeEntries(gsl_matrix_int *m);

double xchoosey(int n, int k);

int matrixAllGTZ(gsl_matrix_int *m);

int matrixSum(gsl_matrix_int *m);
void gsl_matrix_prettyPrint(gsl_matrix *m);
double matrixSumDouble(gsl_matrix *m);

void preorderEvolveState(afsStateSpace *S, afsObject *aState, int n1, int n2);

void coalMarkovChainTopologyMatrix(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int n1, int n2);
int coalMarkovChainTopologyMatrix_2(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int *dim1, int *dim2);
int coalMarkovChainTopologyMatrix_sparse(afsStateSpace *S,double *topol, int *moveType, int *dim1, int *dim2);

void readTopolFile(const char *fileName, double *topA, int *moveA, int *dim1, int *dim2);

void fillTransitionMatrix(afsStateSpace *S, gsl_matrix *topol, gsl_matrix_int *moveType, gsl_matrix *transMat, gsl_vector *rates,
	double theta1, double theta2, double mig1, double mig2);
struct cs_di_sparse *fillTransitionMatrix_2(afsStateSpace *S, gsl_matrix *topol, gsl_matrix_int *moveType, gsl_vector *rates, int nzCount,
	int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2);

void fillExpectedAFS(afsStateSpace *S, gsl_matrix *transMat,gsl_vector *rates, gsl_matrix *expAFS, int n1, int n2);

void fillExpectedAFS_visits(afsStateSpace *S, gsl_matrix *visitMat,gsl_vector *rates, gsl_matrix *expAFS);
void fillExpectedAFS_visits_2(afsStateSpace *S, double *visitMat,gsl_vector *rates, gsl_matrix *expAFS);

void overlayPopnSplitting(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int n1, int n2);

void maximize_island_lik_const(double *lik, void * p);
void maximize_island_lik_hmm_const(double *lik, void * p);

void maximize_island_lik_CILower(double *lik, void * p);
void maximize_island_lik_CIUpper(double *lik, void * p);
void island_lik_constraint_CI(size_t n, const double *pGuess, void * p, double *result);

void island_lik_constraint(size_t n, const double *pGuess, void * p, double *result);
void island_lik_constraint_hmm(size_t n, const double *pGuess, void * p, double *result);

void fillLogAFS_ptr_sp(void * p);
void fillLogAFS_ptr_sp2(void * p);
void fillLogAFS_ptr_sp_hmm(int state,void * p);

void itoa(int n, char s[]);
void reverse(char s[]);

void NumericalGradient_island_lik(size_t n, const double *v, void *params, double *df);
void island_lik_fdf(const size_t n, const double *x,void *fparams,double *fval,double *grad);

#endif

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
/* gsl_multilarge.h
 *
 * Copyright (C) 2015 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MULTILARGE_H__
#define __GSL_MULTILARGE_H__

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.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

/* iteration solver type */
typedef struct
{
  const char *name;
  void * (*alloc) (const size_t p);
  int (*reset) (void *);
  int (*accumulate) (gsl_matrix * X, gsl_vector * y,
                     void *);
  int (*solve) (const double lambda, gsl_vector * c,
                double * rnorm, double * snorm, void *);
  int (*rcond) (double * rcond, void *);
  int (*lcurve) (gsl_vector * reg_param, gsl_vector * rho,
                 gsl_vector * eta, void *);
  void (*free) (void *);
} gsl_multilarge_linear_type;

typedef struct
{
  const gsl_multilarge_linear_type * type;
  void * state;
  size_t p;
} gsl_multilarge_linear_workspace;

/* available types */
GSL_VAR const gsl_multilarge_linear_type * gsl_multilarge_linear_normal;
GSL_VAR const gsl_multilarge_linear_type * gsl_multilarge_linear_tsqr;

/*
 * Prototypes
 */
gsl_multilarge_linear_workspace *
gsl_multilarge_linear_alloc(const gsl_multilarge_linear_type * T,
                            const size_t p);

void gsl_multilarge_linear_free(gsl_multilarge_linear_workspace * w);

const char *gsl_multilarge_linear_name(const gsl_multilarge_linear_workspace * w);

int gsl_multilarge_linear_reset(gsl_multilarge_linear_workspace * w);

int gsl_multilarge_linear_accumulate(gsl_matrix * X,
                                     gsl_vector * y,
                                     gsl_multilarge_linear_workspace * w);

int gsl_multilarge_linear_solve(const double lambda, gsl_vector * c,
                                double * rnorm, double * snorm,
                                gsl_multilarge_linear_workspace * w);

int gsl_multilarge_linear_rcond(double *rcond, gsl_multilarge_linear_workspace * w);

int gsl_multilarge_linear_lcurve(gsl_vector * reg_param, gsl_vector * rho,
                                 gsl_vector * eta,
                                 gsl_multilarge_linear_workspace * w);

int gsl_multilarge_linear_wstdform1 (const gsl_vector * L,
                                     const gsl_matrix * X,
                                     const gsl_vector * w,
                                     const gsl_vector * y,
                                     gsl_matrix * Xs,
                                     gsl_vector * ys,
                                     gsl_multilarge_linear_workspace * work);

int gsl_multilarge_linear_stdform1 (const gsl_vector * L,
                                    const gsl_matrix * X,
                                    const gsl_vector * y,
                                    gsl_matrix * Xs,
                                    gsl_vector * ys,
                                    gsl_multilarge_linear_workspace * work);

int gsl_multilarge_linear_L_decomp (gsl_matrix * L, gsl_vector * tau);

int gsl_multilarge_linear_wstdform2 (const gsl_matrix * LQR,
                                     const gsl_vector * Ltau,
                                     const gsl_matrix * X,
                                     const gsl_vector * w,
                                     const gsl_vector * y,
                                     gsl_matrix * Xs,
                                     gsl_vector * ys,
                                     gsl_multilarge_linear_workspace * work);

int gsl_multilarge_linear_stdform2 (const gsl_matrix * LQR,
                                    const gsl_vector * Ltau,
                                    const gsl_matrix * X,
                                    const gsl_vector * y,
                                    gsl_matrix * Xs,
                                    gsl_vector * ys,
                                    gsl_multilarge_linear_workspace * work);

int gsl_multilarge_linear_genform1 (const gsl_vector * L,
                                    const gsl_vector * cs,
                                    gsl_vector * c,
                                    gsl_multilarge_linear_workspace * work);

int gsl_multilarge_linear_genform2 (const gsl_matrix * LQR,
                                    const gsl_vector * Ltau,
                                    const gsl_vector * cs,
                                    gsl_vector * c,
                                    gsl_multilarge_linear_workspace * work);

__END_DECLS

#endif /* __GSL_MULTILARGE_H__ */

max_stars_repo_name: shi-bash-cmd/qtTest

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <float.h>
#include <inttypes.h>
#include <limits.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <valgrind/callgrind.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

#include "betree.h"
#include "debug.h"
#include "hashmap.h"
#include "helper.h"
#include "printer.h"
#include "tree.h"
#include "utils.h"

#define MAX_EXPRS 10000
#define MAX_EVENTS 5000
#define DEFAULT_SEARCH_COUNT 10

// wc -L filename
#define MAX_EVENT_CHARACTERS 20000
#define MAX_EXPR_CHARACTERS 17000
#define MAX_CONSTANT_CHARACTERS 20

const char* EXPRS_FILE = "data/betree_exprs";
const char* CONSTANTS_FILE = "data/betree_constants";
const char* EVENTS_FILE = "data/betree_events";

struct betree_events {
    size_t count;
    char** events;
};

void add_event(char* event, struct betree_events* events)
{
    if(events->count == 0) {
        events->events = calloc(1, sizeof(*events->events));
        if(events == NULL) {
            fprintf(stderr, "%s calloc failed", __func__);
            abort();
        }
    }
    else {
        char** new_events
            = realloc(events->events, sizeof(*new_events) * ((events->count) + 1));
        if(new_events == NULL) {
            fprintf(stderr, "%s realloc failed", __func__);
            abort();
        }
        events->events = new_events;
    }
    events->events[events->count] = strdup(event);
    events->count++;
}

char* strip_chars(const char* string, const char* chars)
{
    char* new_string = malloc(strlen(string) + 1);
    int counter = 0;

    for(; *string; string++) {
        if(!strchr(chars, *string)) {
            new_string[counter] = *string;
            ++counter;
        }
    }
    new_string[counter] = 0;
    return new_string;
}

int event_parse(const char* text, struct betree_event** event);

size_t read_betree_events(struct betree_events* events)
{
    FILE* f = fopen(EVENTS_FILE, "r");
    size_t count = 0;

    char line[MAX_EVENT_CHARACTERS]; // Arbitrary from what I've seen
    while(fgets(line, sizeof(line), f)) {
        if(MAX_EVENTS != 0 && events->count == MAX_EVENTS) {
            break;
        }

        add_event(line, events);
        count++;
    }
    fclose(f);
    return count;
}

size_t read_betree_exprs(struct betree* tree)
{

    FILE* f = fopen(EXPRS_FILE, "r");
    FILE* constants_f = fopen(CONSTANTS_FILE, "r");

    //char* lines[MAX_EXPRS];
    char line[MAX_EXPR_CHARACTERS]; // Arbitrary from what I've seen
    char constants_line[MAX_CONSTANT_CHARACTERS];
    size_t count = 0;
    const struct betree_sub* subs[MAX_EXPRS];

    enum e { constant_count = 3 };
    while(fgets(line, sizeof(line), f)) {
        char* ignore = fgets(constants_line, sizeof(constants_line), constants_f);
        (void)ignore;
        char* copy = strdup(constants_line);
        char* rest = copy;
        int64_t campaign_id = strtoll(strtok_r(rest, ",", &rest), NULL, 10);
        int64_t advertiser_id = strtoll(strtok_r(rest, ",", &rest), NULL, 10);
        int64_t flight_id = strtoll(strtok_r(rest, "\n", &rest), NULL, 10);
        const struct betree_constant* constants[constant_count] = {
            betree_make_integer_constant("campaign_id", campaign_id),
            betree_make_integer_constant("advertiser_id", advertiser_id),
            betree_make_integer_constant("flight_id", flight_id),
        };

        const struct betree_sub* sub = betree_make_sub(tree, count, constant_count, constants, line);
        subs[count] = sub;
        count++;
        betree_free_constants(constant_count, (struct betree_constant**) constants);
        free(copy);
        if(MAX_EXPRS != 0 && count == MAX_EXPRS) {
            break;
        }
    }
    /*for(size_t i = 0; i < tree->config->attr_domain_count; i++) {*/
        /*const struct attr_domain* attr_domain = tree->config->attr_domains[i];*/
        /*print_attr_domain(attr_domain);*/
    /*}*/
    fclose(f);
    fclose(constants_f);

    for(size_t i = 0; i < count; i++) {
        const struct betree_sub* sub = subs[i];
        if(!betree_insert_sub(tree, sub)) {
            printf("Can't insert expr %zu\n", i);
            abort();
        }
    }

    return count;
}

void read_betree_defs(struct betree* tree)
{
    FILE* f = fopen("data/betree_defs", "r");

    char line[LINE_MAX];
    while(fgets(line, sizeof(line), f)) {
        add_variable_from_string(tree, line);
    }

    fclose(f);
}

int compare_int( const void* a, const void* b )
{
    if( *(int*)a == *(int*)b ) return 0;
    return *(int*)a < *(int*)b ? -1 : 1;
}

int main(int argc, char** argv)
{
    size_t search_count = DEFAULT_SEARCH_COUNT;
    if(argc > 1) {
        search_count = atoi(argv[1]);
    }
    if(access("data/betree_defs", F_OK) == -1 || access("data/betree_events", F_OK) == -1
        || access("data/betree_exprs", F_OK) == -1 || access("data/betree_constants", F_OK) == -1) {
        fprintf(stderr, "Missing files, skipping the tests");
        return 0;
    }
    struct timespec start, insert_done, gen_event_done, search_done;

    // Init
    struct betree* tree = betree_make();
    read_betree_defs(tree);

    clock_gettime(CLOCK_MONOTONIC_RAW, &start);

    // Insert
    size_t expr_count = read_betree_exprs(tree);

    clock_gettime(CLOCK_MONOTONIC_RAW, &insert_done);
    uint64_t insert_us = (insert_done.tv_sec - start.tv_sec) * 1000000
        + (insert_done.tv_nsec - start.tv_nsec) / 1000;
    printf("    Insert took %" PRIu64 "\n", insert_us);

    struct betree_events events = { .count = 0, .events = NULL };
    size_t event_count = read_betree_events(&events);

    uint64_t evaluated_sum = 0;
    uint64_t matched_sum = 0;
    uint64_t memoized_sum = 0;
    uint64_t shorted_sum = 0;

    const size_t search_us_count = search_count * events.count;
    double search_us_data[search_us_count];

    CALLGRIND_START_INSTRUMENTATION;

    size_t search_us_i = 0;

    for(size_t j = 0; j < search_count; j++) {
        for(size_t i = 0; i < events.count; i++) {
            clock_gettime(CLOCK_MONOTONIC_RAW, &gen_event_done);

            char* event = events.events[i];
            struct report* report = make_report();
            if(betree_search(tree, event, report) == false) {
                fprintf(stderr, "Failed to search with event\n");
                abort();
            }

            clock_gettime(CLOCK_MONOTONIC_RAW, &search_done);

            uint64_t search_us = (search_done.tv_sec - gen_event_done.tv_sec) * 1000000
                + (search_done.tv_nsec - gen_event_done.tv_nsec) / 1000;

            search_us_data[search_us_i] = (double)search_us;

            evaluated_sum += report->evaluated;
            matched_sum += report->matched;
            memoized_sum += report->memoized;
            shorted_sum += report->shorted;
            free_report(report);
            search_us_i++;
        }
        printf("Finished run %zu/%zu\n", j, search_count);
    }

    CALLGRIND_STOP_INSTRUMENTATION;
    CALLGRIND_DUMP_STATS;

    for(size_t i = 0; i < events.count; i++) {
        free(events.events[i]);
    }

    double evaluated_average = (double)evaluated_sum / (double)MAX_EVENTS;
    double matched_average = (double)matched_sum / (double)MAX_EVENTS;
    double memoized_average = (double)memoized_sum / (double)MAX_EVENTS;
    double shorted_average = (double)shorted_sum / (double)MAX_EVENTS;
    printf("%zu searches, %zu expressions, %zu events, %zu preds, %zu memoize preds. Evaluated %.2f, matched %.2f, memoized %.2f, shorted %.2f\n",
        search_us_count,
        expr_count,
        event_count,
        tree->config->pred_map->pred_count,
        tree->config->pred_map->memoize_count,
        evaluated_average / search_count,
        matched_average / search_count,
        memoized_average / search_count,
        shorted_average / search_count);

    double search_us_min = gsl_stats_min(search_us_data, 1, search_us_count);
    double search_us_max = gsl_stats_max(search_us_data, 1, search_us_count);
    double search_us_mean = gsl_stats_mean(search_us_data, 1, search_us_count);
    gsl_sort(search_us_data, 1, search_us_count);
    double search_us_90 = gsl_stats_quantile_from_sorted_data(search_us_data, 1, search_us_count, 0.90);
    double search_us_95 = gsl_stats_quantile_from_sorted_data(search_us_data, 1, search_us_count, 0.95);
    double search_us_99 = gsl_stats_quantile_from_sorted_data(search_us_data, 1, search_us_count, 0.99);

    printf("Min: %.1f, Mean: %.1f, Max: %.1f, 90: %.1f, 95: %.1f, 99: %.1f\n", search_us_min, search_us_mean, search_us_max, search_us_90, search_us_95, search_us_99);

    printf("| %lu | %.1f | %.1f | %.1f | %.1f | %.1f | %.1f | |\n", insert_us, search_us_min, search_us_mean, search_us_max, search_us_90, search_us_95, search_us_99);

    // DEBUG
    write_dot_file(tree);
    // DEBUG

    free(events.events);
    betree_free(tree);
    return 0;
}

max_stars_repo_name: jonahharris/be-tree

max_issues_repo_name: jonahharris/be-tree

max_forks_repo_name: jonahharris/be-tree


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>

#include <gsl/gsl_errno.h>

#include "ccl.h"

// Error handling policy: whether to exit on error (C default) or continue
// (Python or other binding default)
static CCLErrorPolicy _ccl_error_policy = CCL_ERROR_POLICY_EXIT;

// Debug mode policy: whether to print error messages as they are raised. This
// is useful for the Python wrapper, which normally allows errors to be
// overwritten by the C code until control returns to Python. If debug mode is
// switched on, the errors are always printed by the C code when they occur.
// Setting the debug mode to warning allows ccl_raise_warning to print warnings
// but keeps the behavior of ccl_raise_exception as if debug mode is set to off.
static CCLDebugModePolicy _ccl_debug_mode_policy = CCL_DEBUG_MODE_WARNING;

// Set error policy
void ccl_set_error_policy(CCLErrorPolicy error_policy)
{
  _ccl_error_policy = error_policy;
}

// Set debug mode policy
void ccl_set_debug_policy(CCLDebugModePolicy debug_policy)
{
    _ccl_debug_mode_policy = debug_policy;
}

// Convenience function to raise exceptions in an appropriate way
void ccl_raise_exception(int err, const char* msg, ...)
{
  char message[256];

  va_list va;
  va_start(va, msg);
  vsnprintf(message, 250, msg, va);
  va_end(va);

  // Print error message and exit if fatal errors are enabled
  if ((_ccl_error_policy == CCL_ERROR_POLICY_EXIT) && (err)) {
    fprintf(stderr, "ERROR %d: %s\n", err, message);
    exit(1);
  }
  // Print error message and exit if debug output is enabled
  else if ((_ccl_debug_mode_policy == CCL_DEBUG_MODE_ON) && (err)){
    fprintf(stderr, "ERROR %d: %s\n", err, message);
  }
}

// Convenience function to handle warnings
void ccl_raise_warning(int err, const char* msg, ...)
{
  char message[256];

  va_list va;
  va_start(va, msg);
  vsnprintf(message, 250, msg, va);
  va_end(va);

  // For now just print warning to stderr if debug is enabled.
  // TODO: Implement some kind of error stack that can be passed on to, e.g.,
  // the python binding.
  if( (_ccl_debug_mode_policy == CCL_DEBUG_MODE_ON)
      || (_ccl_debug_mode_policy == CCL_DEBUG_MODE_WARNING) ) {
    fprintf(stderr, "WARNING: %s\n", message);
  }
}

// Convenience function to handle warnings
void ccl_raise_gsl_warning(int gslstatus, const char* msg, ...)
{
  char message[256];

  va_list va;
  va_start(va, msg);
  vsnprintf(message, 250, msg, va);
  va_end(va);

  ccl_raise_warning(gslstatus, "%s GSL error: %s", message, gsl_strerror(gslstatus));
  return;
}

void ccl_check_status(ccl_cosmology *cosmo, int * status)
{

  switch (*status) {
  case 0: // all good, nothing to do
    return;
  case CCL_ERROR_LINSPACE:	// spacing allocation error, always terminate
    ccl_raise_exception(*status, cosmo->status_message);
  case CCL_ERROR_SPLINE:	// spline allocation error, always terminate
    ccl_raise_exception(*status, cosmo->status_message);
  case CCL_ERROR_COMPUTECHI:	// compute_chi error //RH
    ccl_raise_exception(*status, cosmo->status_message);
  case CCL_ERROR_HMF_INTERP: // terminate if hmf definition not supported
    ccl_raise_exception(*status, cosmo->status_message);
  case CCL_ERROR_NU_INT: // error in getting the neutrino integral spline: exit. No status_message in cosmo because can't pass cosmology to the function.
    ccl_raise_exception(*status, "Error, in ccl_neutrinos.c. ccl_calculate_nu_phasespace_spline(): Error in setting neutrino phasespace spline.");
  case CCL_ERROR_NU_SOLVE: // error in converting Omeganuh2-> Mnu: exit. No status_message in cosmo because can't pass cosmology to the function.
    ccl_raise_exception(*status, "Error, in ccl_neutrinos.c. Omeganuh2_to_Mnu(): Root finding did not converge.");
    // TODO: Implement softer error handling, e.g. for integral convergence here
  default:
    ccl_raise_exception(*status, cosmo->status_message);
  }
}

/* ------- ROUTINE: ccl_check_status_nocosmo ------
   INPUTS: pointer to a status integer
   TASK: Perform a check on status for the case where it is not possible to have a cosmology object where the status check is required.
*/
void ccl_check_status_nocosmo(int * status)
{
  switch (*status) {
  case 0: // Nothing to do
    return;
  case CCL_ERROR_LINSPACE:
    // Spacing allocation error, always terminate
    ccl_raise_exception(*status, "CCL_ERROR_LINSPACE: Spacing allocation error.");
  case CCL_ERROR_SPLINE:
    // Spline allocation error, always terminate
    ccl_raise_exception(*status, "CCL_ERROR_SPLINE: Spline allocation error.");
  case CCL_ERROR_COMPUTECHI:
    // Compute_chi error
    ccl_raise_exception(*status,
             "CCL_ERROR_COMPUTECHI: Comoving distance chi computation failed.");
  case CCL_ERROR_HMF_INTERP:
    // Terminate if hmf definition not supported
    ccl_raise_exception(*status,
          "CCL_ERROR_HMF_INTERP: Halo mass function definition not supported.");
  case CCL_ERROR_NU_INT:
    // Error in getting the neutrino integral spline: exit. No status_message
    // in cosmo because can't pass cosmology to the function.
    ccl_raise_exception(*status,
      "CCL_ERROR_NU_INT: Error getting the neutrino phase-space integral spline.");
  case CCL_ERROR_NU_SOLVE:
    // Error in converting Omeganuh2-> Mnu: exit. No status_message in cosmo
    // because can't pass cosmology to the function.
    ccl_raise_exception(*status,
                      "CCL_ERROR_NU_SOLVE: Error converting Omeganuh2 -> Mnu.");
  case CCL_ERROR_MNU_UNPHYSICAL:
    // Error in the sum of mnu or Omeganu passed for the hierarchy requested.
	  ccl_raise_exception(*status,
      "CCL_ERROR_MNU_UNPHYSICAL: Sum of neutrinos masses for this Omeganu value is incompatible with the requested mass hierarchy.");
  case CCL_ERROR_NOT_IMPLEMENTED:
    ccl_raise_exception(*status,
      "CCL_ERROR_NOT_IMPLEMENTED: the type of m_nu specified is not supported.");
  default:
    ccl_raise_exception(*status,
             "Unrecognized error code (see gsl_errno.h for error codes 1-32).");
  }
}

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
#ifndef __GSL_STATISTICS_H__
#define __GSL_STATISTICS_H__

#include <gsl/gsl_statistics_long_double.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_statistics_float.h>

#include <gsl/gsl_statistics_ulong.h>
#include <gsl/gsl_statistics_long.h>

#include <gsl/gsl_statistics_uint.h>
#include <gsl/gsl_statistics_int.h>

#include <gsl/gsl_statistics_ushort.h>
#include <gsl/gsl_statistics_short.h>

#include <gsl/gsl_statistics_uchar.h>
#include <gsl/gsl_statistics_char.h>

#endif /* __GSL_STATISTICS_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
#pragma once

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

#define SIM_SUCCESS 0
#define SIM_FAILURE -1

#define SIM_INFO(...) do {fprintf(stdout,##__VA_ARGS__);} while(0)
#define SIM_ERROR(...) do {fprintf(stdout,##__VA_ARGS__);} while(0)

#define SHUTDOWN_MESSAGE "shutdown\n"

typedef struct {
  gsl_vector *yd;
  gsl_vector *y;
} VehicleState;

typedef struct {
  gsl_vector *(*feedback) (gsl_vector *yd, gsl_vector *y);
} Controller;

max_stars_repo_name: umd-agrc/SimpleControlSim

max_issues_repo_name: umd-agrc/SimpleControlSim

max_forks_repo_name: umd-agrc/SimpleControlSim


text:
/* randist/cauchy.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, 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 <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The Cauchy probability distribution is

   p(x) dx = (1/(pi a)) (1 + (x/a)^2)^(-1) dx

   It is also known as the Lorentzian probability distribution */

double
gsl_ran_cauchy (const gsl_rng * r, const double a)
{
  double u;
  do
    {
      u = gsl_rng_uniform (r);
    }
  while (u == 0.5);

  return a * tan (M_PI * u);
}

double
gsl_ran_cauchy_pdf (const double x, const double a)
{
  double u = x / a;
  double p = (1 / (M_PI * a)) / (1 + u * u);
  return p;
}


text:
#ifndef BRATU3D_H
#define BRATU3D_H

#include <petsc.h>

typedef struct Params {
  double lambda_;
} Params;

PetscErrorCode FormInitGuess(DM da, Vec x, Params *p);
PetscErrorCode FormFunction(DM da, Vec x, Vec F, Params *p);
PetscErrorCode FormJacobian(DM da, Vec x, Mat J, Params *p);

#endif/*BRATU3D_H*/


text:
/**
 * Subroutines used in the "aXe_AF2PET.c" for the task aXe_AF2PET
 */

#include <stdarg.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_nan.h>
#include <gsl/gsl_sys.h>

#include "spc_wl_calib.h"
#include "aXe_grism.h"
#include "aXe_errors.h"
#include "disp_conf.h"

/**
 * The function performs a wavelength calibration on an allready
 * existing ap_pixel table. It fills in the lambda and dlambda fields.
 *
 * @param ap_p  - the ap_pixel table to work on, with the pathlength field
 *                filled out
 * @param calib - the calibration function
 */
void
wl_calib (ap_pixel * ap_p, const calib_function * const calib)
{
  double l1,l2;

  while (ap_p->p_x != -1)
    {
      ap_p->lambda = calib->func (ap_p->xi, calib->order, calib->coeffs);

      l1 = calib->func (ap_p->xi - .5, calib->order, calib->coeffs);
      l2 = calib->func (ap_p->xi + .5, calib->order, calib->coeffs);
      ap_p->dlambda = fabs(l2-l1);
      ap_p++;
    }
}

/*
 * Function: pwise_wl_calib
 * Make for every PET pixel an individual wavelength
 * calibration. The function was developed for FORS2 MXU spectra
 * and makes especially sense for long spectra and/or fast
 * variations of the 2D calibration solution.
 */
void
pwise_wl_calib(const global_disp *gdisp, const d_point pixel,
	       const beam actbeam, const int for_grism, ap_pixel *ap_p, const calib_function *old_calib)
{

  calib_function *calib;
  d_point actpix;

  double l1,l2;
  double tmp;

  int i=0;

  while (ap_p[i].p_x != -1)
    {
      // get the current xy-position of the pixel
      actpix.x = pixel.x + ap_p[i].dist*cos(actbeam.orient);
      actpix.y = pixel.y + ap_p[i].dist*sin(actbeam.orient);

      // create a calibration for that position
      calib = create_calib_from_gsl_vector(for_grism, get_calvector_from_gdisp(gdisp, actpix));

      // determine the wavelength value
      ap_p[i].lambda = calib->func (ap_p[i].xi, calib->order, calib->coeffs);

      // determine the dlambda value
      l1 = calib->func (ap_p[i].xi - .5, calib->order, calib->coeffs);
      l2 = calib->func (ap_p[i].xi + .5, calib->order, calib->coeffs);
      ap_p[i].dlambda = fabs(l2-l1);

      // determine the trace length
      tmp =  old_calib->func (ap_p[i].xi, old_calib->order, old_calib->coeffs);
      ap_p[i].xi = ap_p[i].xi + (ap_p[i].lambda - tmp) / ap_p[i].dlambda;

      // free the calibration
      // structure
      free_calib(calib);

      // enhance the counter
      i++;
    }
}

/**
 *  The function trims down an ap_pixel table to only allowed values
 *  for the trace path length or trace distance. The input ap_pixel
 *  is only changed. Valid entries are moved to the front side,
 *  and a new endmark is set after all valid entries
 *
 *  @param  ap_pixel       the input ap_pixel table
 *  @param  prange         the valid trace range
 *  @param  wl_calibration the calibration structure, needed to get the
 *                         first dispersion coefficient
 *
 *  @return ret            ap_pixel table output
 */
ap_pixel *
prange_cut(ap_pixel * in_p, const gsl_vector * prange,
           const calib_function * wl_calibration)
{

  //long nelem;
  ap_pixel * ret, *cur_ap, *ap_p;
  double prel;
  double trmin=1.0E+06, trmax=-1.0E6;

  // set all PET list pointers
  // to the beginning of the array
  ret    = in_p;
  cur_ap = in_p;
  ap_p   = in_p;

  while (ap_p->p_x != -1)
    {
      // search for new minimum
      // and amximum in tracelength
      if (ap_p->xi > trmax)
        trmax = ap_p->xi;
      if (ap_p->xi < trmin)
        trmin = ap_p->xi;

      // compute the tracelength with
      // respect to the zeropoint
      prel = ap_p->xi - wl_calibration->coeffs[0];

      // check whether the realitive zeropoint is within the allowed range
      if (prel > gsl_vector_get(prange,0) && prel < gsl_vector_get(prange,1))
        {
          // transfer the properties
          // of the good pixel to the next
          // PET element in the list of good PET's
          cur_ap->p_x    = ap_p->p_x;
          cur_ap->p_y    = ap_p->p_y;
          cur_ap->x      = ap_p->x;
          cur_ap->y      = ap_p->y;
          cur_ap->dist   = ap_p->dist;
          cur_ap->xs     = ap_p->xs;
          cur_ap->ys     = ap_p->ys;
          cur_ap->dxs    = ap_p->dxs;
          cur_ap->xi     = ap_p->xi;
          cur_ap->count  = ap_p->count;
          cur_ap->error  = ap_p->error;
          cur_ap->weight = ap_p->weight;
          cur_ap->contam = ap_p->contam;
          cur_ap->dq     = ap_p->dq;

          // forward the pointer
          // of good pixels
          cur_ap++;
        }

      // always forward the
      // general pointer
      ap_p++;
    }

  // make a new endmark
  cur_ap->p_x = -1;
  cur_ap->p_y = -1;
  cur_ap->count = -1;

  // return the pointer
  return ret;
}

/**
 *  Calibration function for the grism: simple polynomial
 *
 * @param xi     the path length of the sample point
 * @param order  the order of the polynomial, i.e. the number of coefficients
 *               minus one
 * @param coeffs a pointer to order+1 doubles, starting with the zeroth
 *               coefficent
 *
 * @return res   wavelength at xi
 */
static double
grism_calib_func (const double xi, const int order,
		  const double *const coeffs)
{
  int i;
  double res = 0;

  for (i = 0; i <= order; i++)
    {
      res += coeffs[i] * pow(xi,i*1.0);
    }

  return res;
}

/**
 * Calibration function for the prism: reciprocal polynomial in the
 * form: lambda = a1 + a2/(xi-a0) + a3/(xi-a0)**2 + a4/(xi-a0)**3
 *
 * @param  xi     the path length of the sample point
 * @param  order  the order of the polynomial, i.e. the number of coefficients
 *                minus one
 * @param  coeffs a pointer to order+1 doubles, starting with the zeroth
 *                coefficent
 * @return res    wavelength at xi
 */
static double
prism_calib_func (const double xi, const int order,
		  const double *const coeffs)
{
  int i;
  double res = 0;

  if (xi==0.) return GSL_NAN;
  for (i = 0; i <= order-1; i++)
    {
      //	    res += coeffs[i] / pow(xi,i*1.0);
      res += coeffs[i+1] / pow((xi-coeffs[0]),i*1.0);
    }

  return res;
}

/**
 * Create a prism or grism calibration function, i.e. a function of one of
 * the forms "sum(a_i 1/x^i, i=0..order)" or "sum(a_i x^i, i=0..order)"
 * for transforming pathlength along the trace to lambda.
 *
 *@param for_grism if true, a the second form of the calibration function
 *                 is used (normal polynom), otherwise the first one.
 *@param order     order of the polynom
 *@param a         a GSL vector of length n containing the n coefficients
 *                 describing the wavelength calibration polynomial.
 *@return calib    an allocated structure containing an allocated array
 *                 for the coefficients.  Free using calib_free.
 */
calib_function *
create_calib_from_gsl_vector (const int for_grism, const gsl_vector * a)
{
  calib_function *calib;
  int i, order;

  // set the polynomial order
  order = a->size - 1;

  // allocate memory
  calib = malloc(sizeof (calib_function));

  // complain if memory
  // allocation failed
  if (!calib)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__, "Out of memory");

  // store the order of the
  // calibration
  calib->order = order;

  // transfer the calibration functions
  if (for_grism)
    calib->func = &grism_calib_func;
  else
    calib->func = &prism_calib_func;

  // allocate memory; complain if problems occur
  calib->coeffs = malloc ((order + 1) * sizeof (double));
  if (!calib->coeffs)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__, "Out of memory");

  // transfer the coefficients
  calib->coeffs[0] = gsl_vector_get (a, 0);
  for (i = 1; i <= order; i++)
      calib->coeffs[i] = gsl_vector_get (a, i);

  // set the allowed range for
  // the trace lengths to NULL
  calib->pr_range = NULL;

  // return the calibration
  // structure
  return calib;
}

/**
 * Function: free_calib
 * Free a calibration function and its assoicated data.
 *
 * @param  calib the calibration structure to be freed
 *
 */
void
free_calib(calib_function * calib)
{
  free (calib->coeffs);
  calib->coeffs = NULL;
  if (calib->pr_range != NULL)
    {
      gsl_vector_free(calib->pr_range);
      calib->pr_range = NULL;
    }
  free (calib);
  calib = NULL;
}