text:
//STARTWHOLE
static char help[] = "Newton's method for a two-variable system.\n"
    "No analytical Jacobian.  Run with -snes_fd or -snes_mf.\n\n";

#include <petsc.h>

extern PetscErrorCode FormFunction(SNES, Vec, Vec, void*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    SNES  snes;          // nonlinear solver context
    Vec   x, r;          // solution, residual vectors

    PetscInitialize(&argc,&argv,NULL,help);
    ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,2); CHKERRQ(ierr);
    ierr = VecSetFromOptions(x); CHKERRQ(ierr);
    ierr = VecSet(x,1.0); CHKERRQ(ierr);
    ierr = VecDuplicate(x,&r); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,r,FormFunction,NULL); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,x); CHKERRQ(ierr);
    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    VecDestroy(&x);  VecDestroy(&r);  SNESDestroy(&snes);
    return PetscFinalize();
}

PetscErrorCode FormFunction(SNES snes, Vec x, Vec F, void *ctx) {
    PetscErrorCode ierr;
    const double b = 2.0, *ax;
    double       *aF;

    ierr = VecGetArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecGetArray(F,&aF);CHKERRQ(ierr);
    aF[0] = (1.0 / b) * PetscExpReal(b * ax[0]) - ax[1];
    aF[1] = ax[0] * ax[0] + ax[1] * ax[1] - 1.0;
    ierr = VecRestoreArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecRestoreArray(F,&aF);CHKERRQ(ierr);
    return 0;
}
//ENDWHOLE

ext: c


text:
#pragma once

#include "macros.h"

#include <gsl/span>

#include <array>
#include <optional>
#include <string>

namespace linkollector::win {

enum class activity { url, text };

constexpr std::string_view activity_delimiter = "\uedfd";

constexpr std::array<std::byte, activity_delimiter.size()>
    activity_delimiter_bin = []() {
        std::array<std::byte, activity_delimiter.size()> delim = {};
        for (std::size_t i = 0; i < activity_delimiter.size(); ++i) {
            delim.at(i) = static_cast<std::byte>(activity_delimiter[i]);
        }
        return delim;
    }();

[[nodiscard]] std::string activity_to_string(activity activity_) noexcept;

[[nodiscard]] std::wstring activity_to_wstring(activity activity_) noexcept;

[[nodiscard]] std::optional<activity>
activity_from_string(std::string_view activity_) noexcept;

std::optional<std::pair<activity, std::string>>
deserialize(gsl::span<std::byte> msg) noexcept;

} // namespace linkollector::win

ext: h

max_stars_repo_name: Longhanks/linkollector-win

max_issues_repo_name: Longhanks/linkollector-win

max_forks_repo_name: Longhanks/linkollector-win


text:
/*
Performing approximate Bayesian computation sequential Monte Carlo (Toni et al.
2009) for linear regression.

Synthetic data is generated by `ground_truth_and_analysis.ipynb`.

This script writes particles.csv, where each row corresponds to a particle and
each column corresponds to a round of SMC.

Defining #DEBUG_MODE will silence all writing to stdout. One may then add
printf statements in the code, and perhaps write the output to file as:
`./run.sh > output.txt`

Parameter ordering convention:
0 - gradient
1 - intercept
2 - standard deviation

Author: Juvid Aryaman
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_fit.h>

#define N_DATA 30
#define N_PARAMETERS 3

#define N_PARTICLES 2000
#define N_ROUNDS_SMC 25
#define QUANTILE_ACCEPT_DISTANCE 0.8

#define RND gsl_rng_uniform(r)
#define SEED 1
#define DISTANCE_THRESHOLD_INIT_GRADIENT 2
#define DISTANCE_THRESHOLD_INIT_INTERCEPT 50
#define DISTANCE_THRESHOLD_INIT_SIGMA 2

#define X_DATA_FILENAME "x.csv"
#define Y_DATA_FILENAME "y.csv"

//#define DEBUG_MODE

#include "smc.h"
#include "lin_reg.h"

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

/* set up GSL RNG */
gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
/* end of GSL setup */

gsl_rng_set(r, SEED);

/////////////////////////
/*Read data*/
/////////////////////////

FILE *data_pointer_x, *data_pointer_y;

data_pointer_x = fopen(X_DATA_FILENAME, "r");
data_pointer_y = fopen(Y_DATA_FILENAME, "r");

double data_x[N_DATA];
double data_y[N_DATA];
int i, j, read_error_status_x, read_error_status_y;
for (i=0; i < N_DATA; i++){
	read_error_status_x = fscanf(data_pointer_x, "%lf\n", &data_x[i]);
	read_error_status_y = fscanf(data_pointer_y, "%lf\n", &data_y[i]);
}
if (read_error_status_x != 1){printf("Error reading X data\n"); return 0;}
if (read_error_status_y != 1){printf("Error reading Y data\n"); return 0;}

/////////////////////////
/*Initialise variables*/
/////////////////////////

/*Fit a linear model to the data, which will be used as summary statistics of
the data*/
double gradient_fit_data, intercept_fit_data, sigma_fit_data, cov00, cov01, cov11, sumsq;
int gsl_fit_return_value;
gsl_fit_return_value = gsl_fit_linear(data_x, 1, data_y, 1, N_DATA,
																			&intercept_fit_data, &gradient_fit_data,
																			&cov00, &cov01, &cov11, &sumsq);
if (gsl_fit_return_value != 0) {printf("Fit failed.\n"); return -1;}
sigma_fit_data = sqrt(sumsq/(N_DATA-2));

#ifndef DEBUG_MODE
	printf("gradient ML = %.8f\n", gradient_fit_data);
	printf("intercept ML = %.8f\n", intercept_fit_data);
	printf("sigma ML = %.8f\n", sigma_fit_data);
#endif

/*Make a (N_PARAMETERS X N_ROUNDS_SMC X N_PARTICLES) array to store all
particles at all rounds of SMC*/
double*** theta_particle;
theta_particle = (double***) malloc(N_PARAMETERS * sizeof(double**));
for (i = 0; i < N_PARAMETERS; i++){
	theta_particle[i] = (double**) malloc(N_ROUNDS_SMC * sizeof(double*));
}
for (i = 0; i < N_PARAMETERS; i++){
	for (j = 0; j < N_ROUNDS_SMC; j++){
		theta_particle[i][j] = (double*) malloc(N_PARTICLES * sizeof(double));
	}
}

double** distance;
distance = (double**) malloc(N_PARAMETERS * sizeof(double*));
for (i = 0; i < N_PARAMETERS; i++) {
	distance[i] = (double*) malloc(N_PARTICLES * sizeof(double));
}
double **distance_threshold_all = malloc(N_PARAMETERS * sizeof(double*));
for (i = 0; i < N_PARAMETERS; i++) {
	distance_threshold_all[i] = (double*) malloc(N_ROUNDS_SMC * sizeof(double));
}

double distance_threshold[] = {DISTANCE_THRESHOLD_INIT_GRADIENT,
															 DISTANCE_THRESHOLD_INIT_INTERCEPT,
														   DISTANCE_THRESHOLD_INIT_SIGMA};

double *simulated_data = (double*) malloc(N_DATA * sizeof(double));
double *weight = malloc(N_PARTICLES * sizeof(double));

double weight_normalizer = 0.0;

int time_smc=0; // an index of each round of SMC
int param_index_chosen, prior_violated;
int particle_index;
/////////////////////////
/*Perform ABC SMC*/
/////////////////////////

/*For every round of SMC*/
for (time_smc = 0; time_smc < N_ROUNDS_SMC; time_smc++) {
	#ifndef DEBUG_MODE
		printf("Round %d of SMC\n", time_smc);
	#endif
	for (i = 0; i < N_PARAMETERS; i++) {
		distance_threshold_all[i][time_smc] = distance_threshold[i];
	}

	/*Draw or perturb a particle and compute distance*/
	for (particle_index = 0; particle_index < N_PARTICLES; particle_index++) {
		// printf("%d\n", particle_index);
		for (i = 0; i < N_PARAMETERS; i++) {
			// reset distance of particle along each dimension
			distance[i][particle_index] = distance_threshold[i] + 1.0;
		}
		while ((distance[0][particle_index] > distance_threshold[0])||
					 (distance[1][particle_index] > distance_threshold[1])||
				   (distance[2][particle_index] > distance_threshold[2])) {
			if (time_smc == 0) {
				// Sample from the prior
				sample_prior(r, theta_particle, particle_index);
			}
			else{
				/*Sample from the old weights*/
				param_index_chosen = weighted_choice(r, weight);
				if ((param_index_chosen < 0)||(param_index_chosen >= N_PARTICLES)) {
					printf("Error in param_index_chosen\n");
					printf("time_smc = %d\n", time_smc);
					return -1;
				}

				perturb_particle(r, theta_particle, time_smc, param_index_chosen,
												 particle_index);

				// Check if prior support is 0
				prior_violated = check_prior_violated(theta_particle, time_smc,
																							particle_index);
				if(prior_violated == 1) continue;
				}

			simulate_dataset(r, theta_particle, data_x, simulated_data, time_smc,
											 particle_index);

			// Compute distance between data and simulation
			distance_metric_sum_stats(simulated_data,
															 data_x,
															 gradient_fit_data,
															 intercept_fit_data,
															 sigma_fit_data,
														   distance, particle_index);
			//distance[particle_index] = distance_metric_sum_res(simulated_data, data_y);
			// distance[particle_index] = distance_metric_sum_abs_res(simulated_data, data_y);

		}
	}

	#ifndef DEBUG_MODE
		printf("Particles sampled.\n");
	#endif

	/*Compute weights*/
	if (time_smc==0){ for (i = 0; i < N_PARTICLES; i++) weight[i] = 1.0;}
	else{
		weight_normalizer = 0.0;
		for (particle_index = 0; particle_index < N_PARTICLES; particle_index++) {
			weight_normalizer +=  weight[particle_index]*kernel_pdf();
		}
		// print_double_array(weight, N_PARTICLES);
		// printf("\n" );
		// printf("weight_normalizer=%f\n", weight_normalizer);
		// printf("kernel_pdf()=%f\n", kernel_pdf());
		for (particle_index = 0; particle_index < N_PARTICLES; particle_index++) {
			weight[particle_index] = prior_pdf(theta_particle, time_smc,
																				 particle_index)/weight_normalizer;
		}

	}

	/*Normalise weights*/
	weight_normalizer = 0.0;
	for (i = 0; i < N_PARTICLES; i++)	weight_normalizer += weight[i];
	for (i = 0; i < N_PARTICLES; i++)	weight[i] = weight[i]/weight_normalizer;

	/* Resample weights*/
	distance_threshold[0] = update_distance_threshold(distance[0]);
	distance_threshold[1] = update_distance_threshold(distance[1]);
	distance_threshold[2] = update_distance_threshold(distance[2]);

}

#ifndef DEBUG_MODE
	printf("Writing particles to file\n");
#endif
	write_particles_to_csv(theta_particle);

	char *dist_filename = "distances.txt";
	write_2d_double_array_to_csv(distance_threshold_all, N_PARAMETERS, N_ROUNDS_SMC, dist_filename);
#ifndef DEBUG_MODE
	printf("Done!\n");
#endif

return 0; //return from main
} //close main


text:
#ifndef KGS_NULLSPACEQR_H
#define KGS_NULLSPACEQR_H

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

#include "math/QR.h"
#include "math/Nullspace.h"
#include "TransposeQR.h"

/**
 * An implementation of Nullspace backed by a QR decomposition
 */
class NullspaceQR: public Nullspace {
 public:
  /** Will construct a nullspace using a QR (transpose) decomposition */
  NullspaceQR(TransposeQR* qr);

  ~NullspaceQR();

  /** Update the Nullspace (and underlying SVD) to reflect an updated state of the matrix */
  void updateFromMatrix() override;

private:
  TransposeQR* m_qr;                  ///< SVD underlying this nullspace

  friend class Configuration;
};

#endif //KGS_nullptrSPACE_H


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

/* the code in this module performs fixed-point quadrature calculations for
 * integrands and is based on IQPACK */

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

static int fixed_compute(const double a, const double b, const double alpha, const double beta,
                         gsl_integration_fixed_workspace * w);
static int imtqlx ( const int n, double d[], double e[], double z[] );

gsl_integration_fixed_workspace *
gsl_integration_fixed_alloc(const gsl_integration_fixed_type * type, const size_t n,
                            const double a, const double b, const double alpha, const double beta)
{
  int status;
  gsl_integration_fixed_workspace *w;

  /* check inputs */
  if (n < 1)
    {
      GSL_ERROR_VAL ("workspace size n must be at least 1", GSL_EDOM, 0);
    }

  w = calloc(1, sizeof(gsl_integration_fixed_workspace));
  if (w == NULL)
    {
      GSL_ERROR_VAL ("unable to allocate workspace", GSL_ENOMEM, 0);
    }

  w->weights = malloc(n * sizeof(double));
  if (w->weights == NULL)
    {
      gsl_integration_fixed_free(w);
      GSL_ERROR_VAL ("unable to allocate weights", GSL_ENOMEM, 0);
    }

  w->x = malloc(n * sizeof(double));
  if (w->x == NULL)
    {
      gsl_integration_fixed_free(w);
      GSL_ERROR_VAL ("unable to allocate x", GSL_ENOMEM, 0);
    }

  w->diag = malloc(n * sizeof(double));
  if (w->diag == NULL)
    {
      gsl_integration_fixed_free(w);
      GSL_ERROR_VAL ("unable to allocate diag", GSL_ENOMEM, 0);
    }

  w->subdiag = malloc(n * sizeof(double));
  if (w->subdiag == NULL)
    {
      gsl_integration_fixed_free(w);
      GSL_ERROR_VAL ("unable to allocate subdiag", GSL_ENOMEM, 0);
    }

  w->n = n;
  w->type = type;

  /* compute quadrature weights and nodes */
  status = fixed_compute(a, b, alpha, beta, w);
  if (status)
    {
      gsl_integration_fixed_free(w);
      GSL_ERROR_VAL ("error in integration parameters", GSL_EDOM, 0);
    }

  return w;
}

void
gsl_integration_fixed_free(gsl_integration_fixed_workspace * w)
{
  if (w->weights)
    free(w->weights);

  if (w->x)
    free(w->x);

  if (w->diag)
    free(w->diag);

  if (w->subdiag)
    free(w->subdiag);

  free(w);
}

size_t
gsl_integration_fixed_n(const gsl_integration_fixed_workspace * w)
{
  return w->n;
}

double *
gsl_integration_fixed_nodes(const gsl_integration_fixed_workspace * w)
{
  return w->x;
}

double *
gsl_integration_fixed_weights(const gsl_integration_fixed_workspace * w)
{
  return w->weights;
}

int
gsl_integration_fixed(const gsl_function * func, double * result,
                      const gsl_integration_fixed_workspace * w)
{
  const size_t n = w->n;
  size_t i;
  double sum = 0.0;

  for (i = 0; i < n; ++i)
    {
      double fi = GSL_FN_EVAL(func, w->x[i]);
      sum += w->weights[i] * fi;
    }

  *result = sum;

  return GSL_SUCCESS;
}

/*
fixed_compute()
  Compute quadrature weights and nodes
*/

static int
fixed_compute(const double a, const double b, const double alpha, const double beta,
              gsl_integration_fixed_workspace * w)
{
  int s;
  const size_t n = w->n;
  gsl_integration_fixed_params params;
  size_t i;

  params.a = a;
  params.b = b;
  params.alpha = alpha;
  params.beta = beta;

  /* check input parameters */
  s = (w->type->check)(n, &params);
  if (s)
    return s;

  /* initialize Jacobi matrix */
  s = (w->type->init)(n, w->diag, w->subdiag, &params);
  if (s)
    return s;

  if (params.zemu <= 0.0)
    {
      GSL_ERROR("zeroth moment must be positive", GSL_EINVAL);
    }

  for ( i = 0; i < n; i++ )
    {
      w->x[i] = w->diag[i];
    }

  w->weights[0] = sqrt (params.zemu);

  for ( i = 1; i < n; i++ )
    {
      w->weights[i] = 0.0;
    }

  /* diagonalize the Jacobi matrix */
  s = imtqlx (n, w->x, w->subdiag, w->weights);
  if (s)
    return s;

  for (i = 0; i < n; i++)
    {
      w->weights[i] = w->weights[i] * w->weights[i];
    }

  /*
   * The current weights and nodes are valid for a = 0, b = 1.
   * Now scale them for arbitrary a,b
   */
  {
    double p = pow ( params.slp, params.al + params.be + 1.0 );
    size_t k;

    for ( k = 0; k < n; k++ )
      {
        w->x[k] = params.shft + params.slp * w->x[k];
        w->weights[k] = w->weights[k] * p;
      }
  }

  return GSL_SUCCESS;
}

/******************************************************************************/
/*
  Purpose:

    IMTQLX diagonalizes a symmetric tridiagonal matrix.

  Discussion:

    This routine is a slightly modified version of the EISPACK routine to
    perform the implicit QL algorithm on a symmetric tridiagonal matrix.

    The authors thank the authors of EISPACK for permission to use this
    routine.

    It has been modified to produce the product Q' * Z, where Z is an input
    vector and Q is the orthogonal matrix diagonalizing the input matrix.
    The changes consist (essentially) of applying the orthogonal transformations
    directly to Z as they are generated.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 January 2010

  Author:

    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
    C version by John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

    Roger Martin, James Wilkinson,
    The Implicit QL Algorithm,
    Numerische Mathematik,
    Volume 12, Number 5, December 1968, pages 377-383.

  Parameters:

    Input, int N, the order of the matrix.

    Input/output, double D(N), the diagonal entries of the matrix.
    On output, the information in D has been overwritten.

    Input/output, double E(N), the subdiagonal entries of the
    matrix, in entries E(1) through E(N-1).  On output, the information in
    E has been overwritten.

    Input/output, double Z(N).  On input, a vector.  On output,
    the value of Q' * Z, where Q is the matrix that diagonalizes the
    input symmetric tridiagonal matrix.
*/

static int
imtqlx ( const int n, double d[], double e[], double z[] )
{
  double b;
  double c;
  double f;
  double g;
  int i;
  int ii;
  int itn = 30;
  int j;
  int k;
  int l;
  int m;
  int mml;
  double p;
  double r;
  double s;

  if ( n == 1 )
  {
    return GSL_SUCCESS;
  }

  e[n-1] = 0.0;

  for ( l = 1; l <= n; l++ )
  {
    j = 0;
    for ( ; ; )
    {
      for ( m = l; m <= n; m++ )
      {
        if ( m == n )
        {
          break;
        }

        if ( fabs ( e[m-1] ) <= GSL_DBL_EPSILON * ( fabs ( d[m-1] ) + fabs ( d[m] ) ) )
        {
          break;
        }
      }
      p = d[l-1];
      if ( m == l )
      {
        break;
      }
      if ( itn <= j )
      {
        return GSL_EMAXITER;
      }
      j = j + 1;
      g = ( d[l] - p ) / ( 2.0 * e[l-1] );
      r =  sqrt ( g * g + 1.0 );
      g = d[m-1] - p + e[l-1] / ( g + fabs ( r ) * GSL_SIGN ( g ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ( ii = 1; ii <= mml; ii++ )
      {
        i = m - ii;
        f = s * e[i-1];
        b = c * e[i-1];

        if ( fabs ( g ) <= fabs ( f ) )
        {
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e[i] = f * r;
          s = 1.0 / r;
          c = c * s;
        }
        else
        {
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e[i] = g * r;
          c = 1.0 / r;
          s = s * c;
        }
        g = d[i] - p;
        r = ( d[i-1] - g ) * s + 2.0 * c * b;
        p = s * r;
        d[i] = g + p;
        g = c * r - b;
        f = z[i];
        z[i] = s * z[i-1] + c * f;
        z[i-1] = c * z[i-1] - s * f;
      }
      d[l-1] = d[l-1] - p;
      e[l-1] = g;
      e[m-1] = 0.0;
    }
  }
/*
  Sorting.
*/
  for ( ii = 2; ii <= m; ii++ )
  {
    i = ii - 1;
    k = i;
    p = d[i-1];

    for ( j = ii; j <= n; j++ )
    {
      if ( d[j-1] < p )
      {
         k = j;
         p = d[j-1];
      }
    }

    if ( k != i )
    {
      d[k-1] = d[i-1];
      d[i-1] = p;
      p = z[i-1];
      z[i-1] = z[k-1];
      z[k-1] = p;
    }
  }

  return GSL_SUCCESS;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 *
 * @file core_ssbtype1cb.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @date 2011-05-15
 * @generated s Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

#define A(m,n)   (A + LDA * (n) + ((m)-(n)))
#define V(m)     (V + (m))
#define TAU(m)   (TAU + (m))

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_ssbtype1cb is a kernel that will operate on a region (triangle) of data
 *  bounded by st and ed. This kernel eliminate a column by an column-wise
 *  annihiliation, then it apply a left+right update on the hermitian triangle.
 *  Note that the column to be eliminated is located at st-1.
 *
 *  All detail are available on technical report or SC11 paper.
 *  Azzam Haidar, Hatem Ltaief, and Jack Dongarra. 2011.
 *  Parallel reduction to condensed forms for symmetric eigenvalue problems
 *  using aggregated fine-grained and memory-aware kernels. In Proceedings
 *  of 2011 International Conference for High Performance Computing,
 *  Networking, Storage and Analysis (SC '11). ACM, New York, NY, USA, ,
 *  Article 8 , 11 pages.
 *  http://doi.acm.org/10.1145/2063384.2063394
 *
 *******************************************************************************
 *
 * @param[in] N
 *          The order of the matrix A.
 *
 * @param[in] NB
 *          The size of the band.
 *
 * @param[in, out] A
 *          A pointer to the matrix A of size (2*NB+1)-by-N.
 *
 * @param[in] LDA
 *          The leading dimension of the matrix A. LDA >= max(1,2*NB+1)
 *
 * @param[out] V
 *          float array, dimension N if eigenvalue only
 *          requested or (LDV*blkcnt*Vblksiz) if Eigenvectors requested
 *          The Householder reflectors are stored in this array.
 *
 * @param[out] TAU
 *          float array, dimension (N).
 *          The scalar factors of the Householder reflectors are stored
 *          in this array.
 *
 * @param[in] st
 *          A pointer to the start index where this kernel will operate.
 *
 * @param[in] ed
 *          A pointer to the end index where this kernel will operate.
 *
 * @param[in] sweep
 *          The sweep number that is eliminated. it serve to calculate the
 *          pointer to the position where to store the Vs and Ts.
 *
 * @param[in] Vblksiz
 *          constant which correspond to the blocking used when applying the Vs.
 *          it serve to calculate the pointer to the position where to store the
 *          Vs and Ts.
 *
 * @param[in] WANTZ
 *          constant which indicate if Eigenvalue are requested or both
 *          Eigenvalue/Eigenvectors.
 *
 * @param[in] WORK
 *          Workspace of size nb.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
/***************************************************************************
 *          TYPE 1-BAND Lower-columnwise-Householder
 ***************************************************************************/
void
CORE_ssbtype1cb(int N, int NB,
                float *A, int LDA,
                float *V, float *TAU,
                int st, int ed, int sweep, int Vblksiz, int WANTZ,
                float *WORK)
{
    int len, LDX;
    int blkid, vpos, taupos, tpos;

    /* find the pointer to the Vs and Ts as stored by the bulgechasing
     * note that in case no eigenvector required V and T are stored
     * on a vector of size N
     * */
     if( WANTZ == 0 ) {
         vpos   = ((sweep+1)%2)*N + st;
         taupos = ((sweep+1)%2)*N + st;
     } else {
         findVTpos(N, NB, Vblksiz, sweep, st,
                   &vpos, &taupos, &tpos, &blkid);
     }

    LDX = LDA-1;
    len = ed-st+1;
    *V(vpos) = 1.;

    memcpy( V(vpos+1), A(st+1, st-1), (len-1)*sizeof(float) );
    memset( A(st+1, st-1), 0, (len-1)*sizeof(float) );

    /* Eliminate the col  at st-1 */
    LAPACKE_slarfg_work(len, A(st, st-1), V(vpos+1), 1, TAU(taupos) );

    /* Apply left and right on A(st:ed,st:ed) */
    CORE_slarfy(len, A(st,st), LDX, V(vpos), TAU(taupos), WORK);

    return;
}
/***************************************************************************/
#undef A
#undef V
#undef TAU


text:
#ifndef GRAPH_H_
#define GRAPH_H_

#include "common.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/** A graph represented by an adjacency list (made of dynamic arrays). */
typedef struct
{
    int num_v; /** Number of vertices. */

    int *num_e; /** Number of edges per vertex. */

    int *capacity; /** Space allocated to each vertex. */

    int **adj_list; /** Adjacency list. */

    double **w_list; /** Weights. */
}
graph;

/** Initialize a graph with a fixed number of vertices. */
void graph_init(graph *g, int vertices);

/** Number of edges in the entire graph. */
int graph_edges(const graph *g);

/** Number of proper edges in the entire graph. */
int graph_proper_edges(const graph *g);

/** Number of loops in the entire graph. */
int graph_loops(const graph *g);

/** Number of outgoing edges for vertex \f$v\f$. \f$O(1)\f$.*/
int graph_outdegree(const graph *g, int u);

/** Number of ingoing edges for vertex \f$v\f$. \f$O(|V||E|)\f$.*/
int graph_indegree(const graph *g, int u);

/** Return true if the graph is balanced (i.e.: for each vertex \f$v\f$, indegree(v) == outdegree(v)). \f$O(|V||E|)\f$. */
int graph_is_balanced(const graph *g);

/** Add an edge between vertices 'u' and 'v' with weight. O(1), worst-case O(|N|). */
void graph_add_edge(graph *g, int u, int v, double weight);

/** Add an edge between vertices 'u' and 'v' and vice-versa. */
void graph_add_sym_edges(graph *g, int u, int v, double weight);

/** Remove the edge between vertices 'u' and 'v'. */
int graph_rmv_edge(graph *g, int u, int v);

/** Remove the edges between 'u' and 'v' and vice-versa. Return false if only 1 or 0 edges were removed. */
int graph_rmv_sym_edges(graph *g, int u, int v);

/** Return true if the graph has an edge between 'u' and 'v'. O(|E|). */
int graph_has_edge(graph *g, int u, int v);

/** Return true if the graph is strongly connected. */
int graph_strongly_connected(const graph *g);

/** Print as a svg file (assumes all points are in the [0, 1) range).*/
void graph_svg(const graph *g, double *x, double *y, double size, double offset, FILE *out);

/**
 * Print as a svg file with stronger colors for higher abundances. abun must point to an array
 * with one double in the [0, 1) range for each vertex. color = 0 (red), 1 (green), 2 (blue).
 */
void graph_svg_abun(const graph *g, double *x, double *y, double size, double offset, double *abun, int color, FILE *out);

/** Print the graph in GraphML format. */
void graph_graphml(const graph *g, FILE *out, unsigned int id);

/** Print the adjacency list. */
void graph_print(const graph *g, FILE *out);

/** Free the memory of the struct. */
void graph_free(graph *g);

/**
 * Get a random geometric graph in [0,1]^2 with radius 'r'.
 *
 * 'x' and 'y' must be initialized and the correct amount of space must
 * be allocated. The weight of the realized edges is 'r - d', where 'd'
 * is the geometric distance between the two vertices (in short the
 * weight is always between 0 and 1.0).
 */
void graph_get_rgg(graph *g, int vertices, double r, double *x, double *y, gsl_rng *rng);

/**
 * Get a connected random geometric graph in [0,1]^2 with radius 'r'.
 *
 * 'x' and 'y' must be initialized and the correct amount of space must
 * be allocated. The weight of the realized edges is 'r - d', where 'd'
 * is the geometric distance between the two vertices (in short the
 * weight is always between 0 and 1.0).
 */
void graph_get_crgg(graph *g, int vertices, double r, double *x, double *y, gsl_rng *rng);

/**
 * Get a random geometric graph in a rectangle with an area of 1 and
 * radius 'r'.
 *
 * 'x' and 'y' must be initialized and the correct amount of space must
 * be allocated. The weight of the realized edges is 'r - d', where 'd'
 * is the geometric distance between the two vertices (in short the
 * weight is always between 0 and 1.0).
 */
void graph_get_rec_rgg(graph *g, int vertices, double width, double r, double *x, double *y, gsl_rng *rng);

/**
 * Get a connected random geometric graph in a rectangle with an area of
 * 1 and radius 'r'.
 *
 * 'x' and 'y' must be initialized and the correct amount of space must
 * be allocated. The weight of the realized edges is 'r - d', where 'd'
 * is the geometric distance between the two vertices (in short the
 * weight is always between 0 and 1.0).
 */
void graph_get_rec_crgg(graph *g, int vertices, double width, double r, double *x, double *y, gsl_rng *rng);

/**
 * Return a complete graph (all vertices linked to each other). Realized
 * edges have a default weight of 1.0.
 */
void graph_get_complete(graph *g, int vertices);

/**
 * Return a circle. Realized edges have a default weight of 1.0.
 */
void graph_get_circle(graph *g, int vertices);

/**
 * Return a star. Realized edges have a default weight of 1.0.
 */
void graph_get_star(graph *g, int vertices);

///////////////////////////////////////////////////////////////
// 'Private' functions. You shouldn't need those.

/** Recursive function used to test connectivity. */
void graph_test_cc(const graph *g, int *group, int u);

/** Increse the storage of the lists for vertex 'u'. O(|E|). */
void graph_grow_lists(graph *g, int u);

#endif

max_stars_repo_name: PhDP/Origin

max_issues_repo_name: PhDP/Origin

max_forks_repo_name: PhDP/Origin


text:
/*
 *   Copyright (c) 2007, Michael Lehn
 *
 *   All rights reserved.
 *
 *   Redistribution and use in source and binary forms, with or without
 *   modification, are permitted provided that the following conditions
 *   are met:
 *
 *   1) Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *   2) Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions and the following disclaimer in
 *      the documentation and/or other materials provided with the
 *      distribution.
 *   3) Neither the name of the FLENS development group nor the names of
 *      its contributors may be used to endorse or promote products derived
 *      from this software without specific prior written permission.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 *   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 *   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 *   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 *   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 *   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#ifndef FLENS_FLENS_H
#define FLENS_FLENS_H 1

#define ADDRESS(x) reinterpret_cast<const void *>(&x)

#ifndef ASSERT
#define ASSERT(x) assert(x)
#endif //ASSERT

#include <array.h>
#include <bandstorage.h>
#include <blas.h>
#include <blas_flens.h>
#include <cg.h>
#include <complex_helper.h>
#include <crs.h>
#include <densevector.h>
#include <evalclosure.h>
#include <fullstorage.h>
#include <generalmatrix.h>
#include <generic_blas.h>
#include <refcounter.h>
#include <lapack.h>
#include <lapack_flens.h>
#include <listinitializer.h>
#include <matvec.h>
#include <matvecclosures.h>
#include <matvecio.h>
#include <matvecoperations.h>
#include <multigrid.h>
#include <packedstorage.h>
#include <range.h>
#include <snapshot.h>
#include <storage.h>
#include <sparsematrix.h>
#include <sparse_blas.h>
#include <sparse_blas_flens.h>
#include <symmetricmatrix.h>
#include <traits.h>
#include <triangularmatrix.h>
#include <underscore.h>
#include <uplo.h>

#endif // FLENS_FLENS_H

max_stars_repo_name: wmotte/toolkid

max_issues_repo_name: wmotte/toolkid

max_forks_repo_name: wmotte/toolkid


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int
main (void)
{
  int i; double data[2*128];

  for (i = 0; i < 128; i++)
    {
       REAL(data,i) = 0.0; IMAG(data,i) = 0.0;
    }

  REAL(data,0) = 1.0;

  for (i = 1; i <= 10; i++)
    {
       REAL(data,i) = REAL(data,128-i) = 1.0;
    }

  for (i = 0; i < 128; i++)
    {
      printf ("%d %e %e\n", i,
              REAL(data,i), IMAG(data,i));
    }
  printf ("\n");

  gsl_fft_complex_radix2_forward (data, 1, 128);

  for (i = 0; i < 128; i++)
    {
      printf ("%d %e %e\n", i,
              REAL(data,i)/sqrt(128),
              IMAG(data,i)/sqrt(128));
    }

  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>

#ifndef gauge_config
#include "config.h"
#define gauge_config 1

#define sq(x) ((x)*(x))
#define cb(x) ((x)*(x)*(x))
#endif

// SU(2) Matrix
// only the first row carries useful information (Lang p.80)
// note that for the first row must be normalized, i.e. |a|^2+ |b|^2 = 1
// use su2_normalize() for this
typedef struct{
  gsl_complex a;
  gsl_complex b;
} su2_matrix;

// SU(2) matrix methods --------------------------------------

// print matrix
void su2_print(const su2_matrix A){
  printf("%f + i*%f \t", GSL_REAL(A.a), GSL_IMAG(A.a));
  printf("%f + i*%f \n", GSL_REAL(A.b), GSL_IMAG(A.b));
}

// the identity matrix
su2_matrix su2_unit(){
  su2_matrix u;
  u.a=GSL_COMPLEX_ONE;
  u.b=GSL_COMPLEX_ZERO;
  return u;
}

// normalize matrix
su2_matrix su2_norm(const su2_matrix A){
  su2_matrix result;
  double norm = sqrt(gsl_complex_abs2(A.a)+gsl_complex_abs2(A.b));
  result.a = gsl_complex_div(A.a,gsl_complex_rect(norm,0));
  result.b = gsl_complex_div(A.b,gsl_complex_rect(norm,0));
  return result;
}

// multiply matrices
su2_matrix su2_mul(const su2_matrix A, const su2_matrix B){
  su2_matrix result;
  result.a = gsl_complex_sub(gsl_complex_mul(A.a,B.a), gsl_complex_mul(A.b,gsl_complex_conjugate(B.b))); // ac-bd*
  result.b = gsl_complex_add(gsl_complex_mul(A.a,B.b), gsl_complex_mul(A.b,gsl_complex_conjugate(B.a))); // ad+bc*
  return result;
}

// inverse matrix
su2_matrix su2_inv(const su2_matrix A){
  su2_matrix result;
  result.a = gsl_complex_conjugate(A.a);
  result.b = gsl_complex_sub(GSL_COMPLEX_ZERO,A.b);
  return result;
}

// generate random matrix
su2_matrix su2_rand(const gsl_rng* r){
  su2_matrix result;
  result.a = gsl_complex_rect(1,eps*(gsl_rng_uniform(r)*2-1));
  result.b = gsl_complex_rect(0,eps*(gsl_rng_uniform(r)*2-1));
  return su2_norm(result);
}

double su2_trace(const su2_matrix A){
  return 2*GSL_REAL(A.a);
}

// end of SU(2) matrix methods ---------------------------------------------

max_stars_repo_name: uchuutamashi/lqcd

max_issues_repo_name: uchuutamashi/lqcd

max_forks_repo_name: uchuutamashi/lqcd


text:
#if !defined(STRIPWS_H_INCLUDED)
#define STRIPWS_H_INCLUDED

#include "FileEnumerator.h"
#include "Utils.h"

#include <filesystem>
#include <gsl/gsl>
#include <iosfwd>
#include <string>

class StripWS
{
public:
	static int usage(::std::ostream& strm, const ::std::string& progName,
		const char* pMsg);

	StripWS(::gsl::span<const char*const> args);
	int run() const;

	StripWS(const StripWS&) = delete;
	StripWS& operator=(const StripWS&) = delete;
	StripWS(StripWS&&) = delete;
	StripWS& operator=(StripWS&&) = delete;

PRIVATE_EXCEPT_IN_TEST:
	using Path = ::std::filesystem::path;

	void queryFile(const Path& p) const;
	void translateFile(const Path& p) const;

	/// \brief Scans the input file, counts the occurrences of white space at
	/// the end of a line, and optionally strips them into the provided stream.
	///
	/// If pOut is null, then the method simply counts the occurrences of white
	/// space at the end of a line.
	///
	/// If pOut is not null, then the method additionally translates the file
	/// into *pOut.
	static void scanFile(::std::istream& in, /* out */ size_t& numLinesAffected,
		/* out */ size_t& numSpacesStripped, /* out */ size_t& numTabsStripped,
		::std::ostream* pOut = nullptr);

	static void updateCounts(::std::string const& wsRun,
		/* in-out */ size_t& numLinesAffected,
		/* in-out */ size_t& numSpacesStripped,
		/* in-out */ size_t& numTabsStripped);
	static void replaceOriginalFileWithTemp(const Path& originalPath,
		const Path& tempPath);

	bool				m_isInQueryMode;
	FileEnumerator	m_fileEnumerator;
};

#endif // STRIPWS_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


text:
/* specfunc/gsl_sf_coulomb.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#ifndef __GSL_SF_COULOMB_H__
#define __GSL_SF_COULOMB_H__

#include <gsl/gsl_mode.h>
#include <gsl/gsl_sf_result.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

/* Normalized hydrogenic bound states, radial dependence. */

/* R_1 := 2Z sqrt(Z) exp(-Z r)
 */
int gsl_sf_hydrogenicR_1_e(const double Z, const double r, gsl_sf_result * result);
double gsl_sf_hydrogenicR_1(const double Z, const double r);

/* R_n := norm exp(-Z r/n) (2Z/n)^l Laguerre[n-l-1, 2l+1, 2Z/n r]
 *
 * normalization such that psi(n,l,r) = R_n Y_{lm}
 */
int gsl_sf_hydrogenicR_e(const int n, const int l, const double Z, const double r, gsl_sf_result * result);
double gsl_sf_hydrogenicR(const int n, const int l, const double Z, const double r);

/* Coulomb wave functions F_{lam_F}(eta,x), G_{lam_G}(eta,x)
 * and their derivatives; lam_G := lam_F - k_lam_G
 *
 * lam_F, lam_G > -0.5
 * x > 0.0
 *
 * Conventions of Abramowitz+Stegun.
 *
 * Because there can be a large dynamic range of values,
 * overflows are handled gracefully. If an overflow occurs,
 * GSL_EOVRFLW is signalled and exponent(s) are returned
 * through exp_F, exp_G. These are such that
 *
 *   F_L(eta,x)  =  fc[k_L] * exp(exp_F)
 *   G_L(eta,x)  =  gc[k_L] * exp(exp_G)
 *   F_L'(eta,x) = fcp[k_L] * exp(exp_F)
 *   G_L'(eta,x) = gcp[k_L] * exp(exp_G)
 */
int
gsl_sf_coulomb_wave_FG_e(const double eta, const double x,
                            const double lam_F,
			    const int  k_lam_G,
                            gsl_sf_result * F, gsl_sf_result * Fp,
			    gsl_sf_result * G, gsl_sf_result * Gp,
			    double * exp_F, double * exp_G);

/* F_L(eta,x) as array */
int gsl_sf_coulomb_wave_F_array(
  double lam_min, int kmax,
  double eta, double x,
  double * fc_array,
  double * F_exponent
  );

/* F_L(eta,x), G_L(eta,x) as arrays */
int gsl_sf_coulomb_wave_FG_array(double lam_min, int kmax,
                                double eta, double x,
                                double * fc_array, double * gc_array,
				double * F_exponent,
				double * G_exponent
                                );

/* F_L(eta,x), G_L(eta,x), F'_L(eta,x), G'_L(eta,x) as arrays */
int gsl_sf_coulomb_wave_FGp_array(double lam_min, int kmax,
                                double eta, double x,
                                double * fc_array, double * fcp_array,
				double * gc_array, double * gcp_array,
				double * F_exponent,
				double * G_exponent
                                );

/* Coulomb wave function divided by the argument,
 * F(eta, x)/x. This is the function which reduces to
 * spherical Bessel functions in the limit eta->0.
 */
int gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax,
                                        double eta, double x,
	              	                double * fc_array,
			                double * F_exponent
	      	      	                );

/* Coulomb wave function normalization constant.
 * [Abramowitz+Stegun 14.1.8, 14.1.9]
 */
int gsl_sf_coulomb_CL_e(double L, double eta, gsl_sf_result * result);
int gsl_sf_coulomb_CL_array(double Lmin, int kmax, double eta, double * cl);

__END_DECLS

#endif /* __GSL_SF_COULOMB_H__ */


text:
/* multifit_nlinear/scaling.c
 *
 * Copyright (C) 2015, 2016 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.
 */

/*
 * This module handles the updating of the scaling matrix D_k in the
 * trust region subproblem:
 *
 * min m_k (dx), || D_k dx || <= Delta_k
 *
 * where m_k(dx) is a model which approximates the cost function
 * F(x_k + dx) near the current iteration point x_k
 *
 * D_k can be updated according to several different strategies.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>

static int init_diag_levenberg(const gsl_matrix * J, gsl_vector * diag);
static int update_diag_levenberg(const gsl_matrix * J,
                                 gsl_vector * diag);

static int init_diag_marquardt(const gsl_matrix * J, gsl_vector * diag);
static int update_diag_marquardt (const gsl_matrix * J,
                                  gsl_vector * diag);

static int init_diag_more(const gsl_matrix * J, gsl_vector * diag);
static int update_diag_more(const gsl_matrix * J, gsl_vector * diag);

/* Levenberg scaling, D = I */
static int
init_diag_levenberg(const gsl_matrix * J, gsl_vector * diag)
{
  (void)J; /* avoid unused parameter warning */
  gsl_vector_set_all(diag, 1.0);
  return GSL_SUCCESS;
}

static int
update_diag_levenberg(const gsl_matrix * J, gsl_vector * diag)
{
  (void)J;    /* avoid unused parameter warning */
  (void)diag; /* avoid unused parameter warning */

  /* nothing to do */
  return GSL_SUCCESS;
}

/* initialize diagonal scaling matrix D according to Marquardt method */
static int
init_diag_marquardt(const gsl_matrix * J, gsl_vector * diag)
{
  return update_diag_marquardt(J, diag);
}

/* update diagonal scaling matrix D according to Marquardt method */
static int
update_diag_marquardt (const gsl_matrix * J, gsl_vector * diag)
{
  const size_t p = J->size2;
  size_t j;

  for (j = 0; j < p; j++)
    {
      gsl_vector_const_view v = gsl_matrix_const_column(J, j);
      double norm = gsl_blas_dnrm2(&v.vector);

      if (norm == 0.0)
        norm = 1.0;

      gsl_vector_set(diag, j, norm);
    }

  return GSL_SUCCESS;
}

/* initialize diagonal scaling matrix D according to Eq 6.3 of
 * More, 1978 */
static int
init_diag_more(const gsl_matrix * J, gsl_vector * diag)
{
  int status;

  gsl_vector_set_zero(diag);
  status = update_diag_more(J, diag);

  return status;
}

/* update diagonal scaling matrix D according to Eq. 6.3 of
 * More, 1978 */
static int
update_diag_more (const gsl_matrix * J, gsl_vector * diag)
{
  const size_t p = J->size2;
  size_t j;

  for (j = 0; j < p; j++)
    {
      gsl_vector_const_view v = gsl_matrix_const_column(J, j);
      double norm = gsl_blas_dnrm2(&v.vector);
      double *diagj = gsl_vector_ptr(diag, j);

      if (norm == 0.0)
        norm = 1.0;

      *diagj = GSL_MAX(*diagj, norm);
    }

  return GSL_SUCCESS;
}

static const gsl_multifit_nlinear_scale levenberg_type =
{
  "levenberg",
  init_diag_levenberg,
  update_diag_levenberg
};

static const gsl_multifit_nlinear_scale marquardt_type =
{
  "marquardt",
  init_diag_marquardt,
  update_diag_marquardt
};

static const gsl_multifit_nlinear_scale more_type =
{
  "more",
  init_diag_more,
  update_diag_more
};

const gsl_multifit_nlinear_scale *gsl_multifit_nlinear_scale_levenberg = &levenberg_type;
const gsl_multifit_nlinear_scale *gsl_multifit_nlinear_scale_marquardt = &marquardt_type;
const gsl_multifit_nlinear_scale *gsl_multifit_nlinear_scale_more = &more_type;

max_stars_repo_name: ielomariala/Hex-Game


text:
/*******************************************************************************
NAME: extract_rgb_palette

ALGORITHM DESCRIPTION:
  Dumps RGB indexed color palette from graphics file to stdout

ISSUES:
        Only color indexed TIFF files currently supported

*******************************************************************************/
#include "asf.h"
#include "asf_nan.h"
#include <math.h>
#include <ctype.h>
#include "asf_raster.h"
#include "typlim.h"
#include <float_image.h>
#include <uint8_image.h>
#include <geokeys.h>
#include <geo_tiffp.h>
#include <geo_keyp.h>
#include <geotiff.h>
#include <geotiffio.h>
#include <tiff.h>
#include <tiffio.h>
#include <xtiffio.h>
#include <gsl/gsl_math.h>
#include "geotiff_support.h"
#include "extract_rgb_palette_help.h"

/**** TYPES ****/
typedef enum {
  UNKNOWN_GRAPHICS_TYPE=0,
  ASF_IMG,
  JPEG,
  PGM,
  PPM,
  PBM,
  STD_TIFF,
  GEO_TIFF,
  BMP,
  GIF,
  PNG
} graphics_file_t;

#define MISSING_TIFF_DATA -1
typedef struct {
  uint32 width;
  uint32 height;
  short sample_format;
  short bits_per_sample;
  short planar_config;
  data_type_t data_type; // ASF data type
  short num_bands;
  int is_scanline_format;
  int is_palette_color_tiff;
} tiff_data_t;

#define MISSING_GTIF_DATA -1
typedef struct {
  int gtif_data_exists;
  char *GTcitation;
  char *PCScitation;
  int tie_point_elements; // Number of tie point elements (usually a multiple of 6)
  int num_tie_points; // Number of elements divided by 6 since (i,j,k) maps to (x,y,z)
  double *tie_point;  // Usually only 1, but who knows what evil lurks in the hearts of men?
  int pixel_scale_elements; // Usually a multiple of 3
  int num_pixel_scales;
  double *pixel_scale;  // Should always be 3 of these ...for ScaleX, ScaleY, and ScaleZ
  short model_type;
  short raster_type;
  short linear_units;
  double scale_factor;
  datum_type_t datum;
  char hemisphere;
  unsigned long pro_zone; // UTM zone (UTM only)
  short proj_coords_trans;
  short pcs;
  short geodetic_datum;
  short geographic_datum;
  double false_easting;
  double false_northing;
  double natLonOrigin;
  double lonCenter;
  double falseOriginLon;
  double falseOriginLat;
  double natLatOrigin;
  double latCenter;
  double stdParallel1;
  double stdParallel2;
  double lonPole;
} geotiff_data_t;

/**** PROTOTYPES ****/
graphics_file_t getGraphicsFileType (char *file);
void graphicsFileType_toStr (graphics_file_t type, char *type_str);
float get_maxval(data_type_t data_type);
char *data_type2str(data_type_t data_type);
void get_tiff_info_from_file(char *file, tiff_data_t *t);
void get_tiff_info(TIFF *tif, tiff_data_t *t);
void get_geotiff_keys(char *file, geotiff_data_t *g);
void dump_tiff_rgb_palette(char *inFile);

int main(int argc, char **argv)
{
  char *inFile;
  char msg[1024];
  char type_str[255];

  asfSplashScreen(argc, argv);

  if (argc != 2) {
      usage();
      exit(1);
  }
  else {
      check_for_help(argc, argv);
  }

  inFile=argv[1];

  if (!fileExists(inFile)) {
    sprintf(msg, "File not found: %s\n  => Did you forget to use the filename extension?\n", inFile);
    asfPrintError(msg);
  }
  graphics_file_t type;
  type = getGraphicsFileType(inFile);
  if (type != STD_TIFF &&
      type != GEO_TIFF )
  {
    graphicsFileType_toStr(type, type_str);
    sprintf(msg, "Graphics file type %s is not currently supported (Graphics file: %s)\n",
            type_str, inFile);
    asfPrintError(msg);
  }

  switch (type) {
    case STD_TIFF:
    case GEO_TIFF:
      dump_tiff_rgb_palette(inFile);
      break;
    default:
      sprintf(msg, "Unrecognized image file type found.\n");
      asfPrintError(msg);
      break;
  }

  return (0);
}

graphics_file_t getGraphicsFileType (char *file)
{
  FILE *fp = (FILE *)FOPEN(file, "rb");
  uint8 magic[4];
  graphics_file_t file_type = UNKNOWN_GRAPHICS_TYPE;
  TIFF *itif = NULL;
  GTIF *gtif = NULL;

  if (fp != NULL) {
    int i;

    // Read the magic number from the file header
    for (i=0; i<4; i++) {
      magic[i] = (uint8)fgetc(fp);
    }
    if(fp) FCLOSE(fp);

    // Check for a valid magic number combination
    if (magic[0] == 0xFF && magic[1] == 0xD8) {
      file_type = JPEG;
    }
    else if (magic[0] == 'P' && (magic[1] == '4' || magic[1] == '1')) {
      file_type = PBM;
    }
    else if (magic[0] == 'P' && (magic[1] == '5' || magic[1] == '2')) {
      file_type = PGM;
    }
    else if (magic[0] == 'P' && (magic[1] == '6' || magic[1] == '3')) {
      file_type = PPM;
    }
    else if ((magic[0] == 'I' && magic[1] == 'I') || (magic[0] == 'M' && magic[1] == 'M')) {
      file_type = STD_TIFF;

      itif = XTIFFOpen(file, "rb");
      if (itif != NULL) {
        gtif = GTIFNew(itif);
        if (gtif != NULL) {
          double *tie_point = NULL;
          double *pixel_scale = NULL;
          short model_type, raster_type, linear_units;
          int read_count, tie_points, pixel_scales;

          (gtif->gt_methods.get)(gtif->gt_tif, GTIFF_TIEPOINTS, &tie_points, &tie_point);
          (gtif->gt_methods.get)(gtif->gt_tif, GTIFF_PIXELSCALE, &pixel_scales, &pixel_scale);
          read_count = GTIFKeyGet(gtif, GTModelTypeGeoKey, &model_type, 0, 1);
          read_count += GTIFKeyGet(gtif, GTRasterTypeGeoKey, &raster_type, 0, 0);
          read_count += GTIFKeyGet(gtif, ProjLinearUnitsGeoKey, &linear_units, 0, 1);

          if (tie_points == 6 && pixel_scales == 3 && read_count == 3) {
            file_type = GEO_TIFF;
          }

          if (tie_point != NULL) free(tie_point);
          if (pixel_scale != NULL) free(pixel_scale);
          GTIFFree(gtif);
        }
        XTIFFClose(itif);
      }
    }
    else if (magic[0] == 'B' && magic[1] == 'M') {
      file_type = BMP;
    }
    else if (magic[0] == 'G' && magic[1] == 'I' && magic[2] == 'F') {
      file_type = GIF;
    }
    else if (magic[1] == 'P' && magic[2] == 'N' && magic[3] == 'G') {
      file_type = PNG;
    }
    else {
      file_type = UNKNOWN_GRAPHICS_TYPE;
    }
  }

  return file_type;
}

void graphicsFileType_toStr (graphics_file_t type, char *type_str)
{
  switch (type) {
    case ASF_IMG:
      strcpy (type_str, "ASF_IMG");
      break;
    case JPEG:
      strcpy (type_str, "JPEG");
      break;
    case PBM:
      strcpy (type_str, "PBM");
      break;
    case PGM:
      strcpy (type_str, "PGM");
      break;
    case PPM:
      strcpy (type_str, "PPM");
      break;
    case STD_TIFF:
      strcpy (type_str, "TIFF");
      break;
    case GEO_TIFF:
      strcpy (type_str, "GEOTIFF");
      break;
    case BMP:
      strcpy (type_str, "BMP");
      break;
    case GIF:
      strcpy (type_str, "GIF");
      break;
    case PNG:
      strcpy (type_str, "PNG");
      break;
    case UNKNOWN_GRAPHICS_TYPE:
    default:
      strcpy(type_str, "UNRECOGNIZED");
      break;
  }
}

float get_maxval(data_type_t data_type)
{
  float ret;

  // Only non-complex types with 32 bits or less are supported
  switch (data_type) {
    case BYTE:
      ret = powf(2, sizeof(unsigned char)) - 1.0;
      break;
    case INTEGER16:
      ret = powf(2, sizeof(short int)) - 1.0;
      break;
    case INTEGER32:
      ret = powf(2, sizeof(int)) - 1.0;
      break;
    case REAL32:
      ret = MAXREAL;
      break;
    default:
      ret = 0.0;
      break;
  }

  return ret;
}

// User must free the returned string
char *data_type2str(data_type_t data_type)
{
  char *retstr = (char*)CALLOC(64, sizeof(char));

  switch (data_type) {
    case BYTE:
      strcpy(retstr, "BYTE");
      break;
    case INTEGER16:
      strcpy(retstr, "INTEGER16");
      break;
    case INTEGER32:
      strcpy(retstr, "INTEGER32");
      break;
    case REAL32:
      strcpy(retstr, "REAL32");
      break;
    case REAL64:
      strcpy(retstr, "REAL64");
      break;
    case COMPLEX_BYTE:
      strcpy(retstr, "COMPLEX_BYTE");
      break;
    case COMPLEX_INTEGER16:
      strcpy(retstr, "COMPLEX_INTEGER16");
      break;
    case COMPLEX_INTEGER32:
      strcpy(retstr, "COMPLEX_INTEGER32");
      break;
    case COMPLEX_REAL32:
      strcpy(retstr, "COMPLEX_REAL32");
      break;
    case COMPLEX_REAL64:
      strcpy(retstr, "COMPLEX_REAL64");
      break;
    default:
      strcpy(retstr, "UNKNOWN");
      break;
  }

  return retstr;
}

void get_tiff_info_from_file(char *file, tiff_data_t *t)
{
  TIFF *tif;

  t->sample_format = MISSING_TIFF_DATA;
  t->bits_per_sample = MISSING_TIFF_DATA;
  t->planar_config = MISSING_TIFF_DATA;
  t->data_type = 0;
  t->num_bands = 0;
  t->is_scanline_format = 0;
  t->height = 0;
  t->width = 0;

  tif = XTIFFOpen(file, "rb");
  if (tif != NULL) {
    get_tiff_info(tif, t);
  }

  XTIFFClose(tif);
}

void get_tiff_info(TIFF *tif, tiff_data_t *t)
{
  get_tiff_data_config(tif,
                       &t->sample_format,
                       &t->bits_per_sample,
                       &t->planar_config,
                       &t->data_type,
                       &t->num_bands,
                       &t->is_scanline_format,
                       &t->is_palette_color_tiff,
                       REPORT_LEVEL_NONE);
  TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &t->height);
  TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &t->width);
  if (t->planar_config != PLANARCONFIG_CONTIG &&
      t->planar_config != PLANARCONFIG_SEPARATE &&
      t->num_bands == 1)
  {
    t->planar_config = PLANARCONFIG_CONTIG;
  }
}

void get_geotiff_keys(char *file, geotiff_data_t *g)
{
  TIFF *tif = XTIFFOpen(file, "rb");
  GTIF *gtif = NULL;

  // Init values to 'missing'
  g->gtif_data_exists = 0;
  g->GTcitation = NULL;   // Should be unallocated
  g->PCScitation = NULL;  // Should be unallocated

  // Read geotiff info
  if (tif != NULL) {
    gtif = GTIFNew(tif);
    if (gtif != NULL) {
      int count, read_count;
      int citation_length;
      int typeSize;
      tagtype_t citation_type;

      // Get citations
      citation_length = GTIFKeyInfo(gtif, GTCitationGeoKey, &typeSize, &citation_type);
      if (citation_length > 0) {
        g->GTcitation = (char*)MALLOC(citation_length * typeSize);
        GTIFKeyGet(gtif, GTCitationGeoKey, g->GTcitation, 0, citation_length);
      }
      else {
        g->GTcitation = NULL;
      }
      citation_length = GTIFKeyInfo(gtif, PCSCitationGeoKey, &typeSize, &citation_type);
      if (citation_length > 0) {
        g->PCScitation = (char*)MALLOC(citation_length * typeSize);
        GTIFKeyGet(gtif, PCSCitationGeoKey, g->PCScitation, 0, citation_length);
      }
      else {
        g->PCScitation = NULL;
      }
      if ((g->GTcitation != NULL && strlen(g->GTcitation) > 0) ||
           (g->PCScitation != NULL && strlen(g->PCScitation) > 0))
      {
        g->gtif_data_exists = 1;
      }

      // Get tie points and pixel scale
      (gtif->gt_methods.get)(gtif->gt_tif, GTIFF_TIEPOINTS, &count, &g->tie_point);
      if (count >= 6) {
        g->gtif_data_exists = 1;
        g->tie_point_elements = count;
        g->num_tie_points = count / 6;
      }
      (gtif->gt_methods.get)(gtif->gt_tif, GTIFF_PIXELSCALE, &count, &g->pixel_scale);
      if (count >= 3) {
        g->gtif_data_exists = 1;
        g->pixel_scale_elements = count;
        g->num_pixel_scales = count / 3;
      }

      // Get model type, raster type, and linear units
      read_count = GTIFKeyGet (gtif, GTModelTypeGeoKey, &g->model_type, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      read_count = GTIFKeyGet (gtif, GTRasterTypeGeoKey, &g->raster_type, 0, 0);
      if (read_count >= 1) g->gtif_data_exists = 1;
      read_count = GTIFKeyGet (gtif, ProjLinearUnitsGeoKey, &g->linear_units, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;

      // Get UTM related info if it exists
      read_count = GTIFKeyGet(gtif, ProjectedCSTypeGeoKey, &g->pcs, 0, 1);
      if (read_count == 1 && PCS_2_UTM(g->pcs, &g->hemisphere, &g->datum, &g->pro_zone)) {
        g->gtif_data_exists = 1;
      }
      else {
        read_count = GTIFKeyGet(gtif, ProjectionGeoKey, &g->pcs, 0, 1);
        if (read_count == 1 && PCS_2_UTM(g->pcs, &g->hemisphere, &g->datum, &g->pro_zone)) {
          g->gtif_data_exists = 1;
        }
        else {
          g->hemisphere = '\0';
          g->datum = UNKNOWN_DATUM;
          g->pro_zone = MISSING_GTIF_DATA;
        }
      }

      // Get projection type (ProjCoordTransGeoKey) and other projection parameters
      read_count = GTIFKeyGet(gtif, ProjCoordTransGeoKey, &g->proj_coords_trans, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->proj_coords_trans = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet(gtif, GeographicTypeGeoKey, &g->geographic_datum, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->geographic_datum = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet(gtif, GeogGeodeticDatumGeoKey, &g->geodetic_datum, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->geodetic_datum = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet(gtif, ProjScaleAtNatOriginGeoKey, &g->scale_factor, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->scale_factor = MISSING_GTIF_DATA;

      // Get generic projection parameters (Note: projection type is defined by
      // the g->proj_coords_trans value)
      read_count = GTIFKeyGet (gtif, ProjFalseEastingGeoKey, &g->false_easting, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->false_easting = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjFalseNorthingGeoKey, &g->false_northing, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->false_northing = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjNatOriginLongGeoKey, &g->natLonOrigin, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->natLonOrigin = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjNatOriginLatGeoKey, &g->natLatOrigin, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->natLatOrigin = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjStdParallel1GeoKey, &g->stdParallel1, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->stdParallel1 = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjStdParallel2GeoKey, &g->stdParallel2, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->stdParallel2 = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjCenterLongGeoKey, &g->lonCenter, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->lonCenter = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjCenterLatGeoKey, &g->latCenter, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->latCenter = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjFalseOriginLongGeoKey, &g->falseOriginLon, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->falseOriginLon = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjFalseOriginLatGeoKey, &g->falseOriginLat, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->falseOriginLat = MISSING_GTIF_DATA;
      read_count = GTIFKeyGet (gtif, ProjStraightVertPoleLongGeoKey, &g->lonPole, 0, 1);
      if (read_count >= 1) g->gtif_data_exists = 1;
      else g->lonPole = MISSING_GTIF_DATA;
    }
  }
}

void tiff_get_float_line(TIFF *tif, float *buf, int row, int band_no)
{
  tiff_data_t t;

  get_tiff_info(tif, &t);
  // Note: For single-plane (greyscale) images, planar_config may remain unset so
  // we can't use it as a guide for checking TIFF validity...
  if (t.sample_format == MISSING_TIFF_DATA ||
      t.bits_per_sample == MISSING_TIFF_DATA ||
      t.data_type == 0 ||
      t.num_bands == MISSING_TIFF_DATA ||
      t.is_scanline_format == MISSING_TIFF_DATA ||
      t.height == 0 ||
      t.width == 0)
  {
    asfPrintError("Cannot read tif file\n");
  }

  // Read a scanline
  tsize_t scanlineSize = TIFFScanlineSize(tif);
  tdata_t *tif_buf = _TIFFmalloc(scanlineSize);
  if (t.planar_config == PLANARCONFIG_CONTIG || t.num_bands == 1) {
    TIFFReadScanline(tif, tif_buf, row, 0);
  }
  else {
    // Planar configuration is band-sequential
    TIFFReadScanline(tif, tif_buf, row, band_no);
  }

  int col;
  for (col=0; col<t.width; col++) {
    switch(t.bits_per_sample) {
      case 8:
        switch(t.sample_format) {
          case SAMPLEFORMAT_UINT:
            if (t.planar_config == PLANARCONFIG_CONTIG && t.num_bands > 1) {
              buf[col] = (float)(((uint8*)(tif_buf))[(col*t.num_bands)+band_no]);   // Current sample.
            }
            else {
              // Planar configuration is band-sequential or single-banded
              buf[col] = (float)(((uint8*)(tif_buf))[col]);
            }
            break;
          case SAMPLEFORMAT_INT:
            if (t.planar_config == PLANARCONFIG_CONTIG && t.num_bands > 1) {
              buf[col] = (float)(((int8*)(tif_buf))[(col*t.num_bands)+band_no]);   // Current sample.
            }
            else {
              // Planar configuration is band-sequential or single-banded
              buf[col] = (float)(((int8*)(tif_buf))[col]);   // Current sample.
            }
            break;
          default:
            // There is no such thing as an IEEE 8-bit floating point
            if (tif_buf) _TIFFfree(tif_buf);
            if (tif) XTIFFClose(tif);
            asfPrintError("tiff_get_float_line(): Unexpected data type in TIFF file.\n");
            break;
        }
        break;
      case 16:
        switch(t.sample_format) {
          case SAMPLEFORMAT_UINT:
            if (t.planar_config == PLANARCONFIG_CONTIG && t.num_bands > 1) {
              buf[col] = (float)(((uint16*)(tif_buf))[(col*t.num_bands)+band_no]);   // Current sample.
            }
            else {
              // Planar configuration is band-sequential or single-banded
              buf[col] = (float)(((uint16*)(tif_buf))[col]);   // Current sample.
            }
            break;
          case SAMPLEFORMAT_INT:
            if (t.planar_config == PLANARCONFIG_CONTIG && t.num_bands > 1) {
              buf[col] = (float)(((int16*)(tif_buf))[(col*t.num_bands)+band_no]);   // Current sample.
            }
            else {
              // Planar configuration is band-sequential or single-banded
              buf[col] = (float)(((uint16*)(tif_buf))[col]);   // Current sample.
            }
            break;
          default:
            // There is no such thing as an IEEE 16-bit floating point
            if (tif_buf) _TIFFfree(tif_buf);
            if (tif) XTIFFClose(tif);
            asfPrintError("tiff_get_float_line(): Unexpected data type in TIFF file.\n");
            break;
        }
        break;
      case 32:
        switch(t.sample_format) {
          case SAMPLEFORMAT_UINT:
            if (t.planar_config == PLANARCONFIG_CONTIG && t.num_bands > 1) {
              buf[col] = (float)(((uint32*)(tif_buf))[(col*t.num_bands)+band_no]);   // Current sample.
            }
            else {
              // Planar configuration is band-sequential or single-banded
              buf[col] = (float)(((uint32*)(tif_buf))[col]);   // Current sample.
            }
            break;
          case SAMPLEFORMAT_INT:
            if (t.planar_config == PLANARCONFIG_CONTIG && t.num_bands > 1) {
              buf[col] = (float)(((long*)(tif_buf))[(col*t.num_bands)+band_no]);   // Current sample.
            }
            else {
              // Planar configuration is band-sequential or single-banded
              buf[col] = (float)(((long*)(tif_buf))[col]);   // Current sample.
            }
            break;
          case SAMPLEFORMAT_IEEEFP:
            if (t.planar_config == PLANARCONFIG_CONTIG && t.num_bands > 1) {
              buf[col] = (float)(((float*)(tif_buf))[(col*t.num_bands)+band_no]);   // Current sample.
            }
            else {
              // Planar configuration is band-sequential or single-banded
              buf[col] = (float)(((float*)(tif_buf))[col]);   // Current sample.
            }
            break;
          default:
            if (tif_buf) _TIFFfree(tif_buf);
            if (tif) XTIFFClose(tif);
            asfPrintError("tiff_get_float_line(): Unexpected data type in TIFF file.\n");
            break;
        }
        break;
      default:
        if (tif_buf) _TIFFfree(tif_buf);
        if (tif) XTIFFClose(tif);
        asfPrintError("tiff_get_float_line(): Unexpected data type in TIFF file.\n");
        break;
    }
  }
  if (tif_buf) {
    _TIFFfree(tif_buf);
  }
}

float tiff_image_get_float_pixel(TIFF *tif, int row, int col, int band_no)
{
  tiff_data_t t;
  float cs;

  get_tiff_info(tif, &t);
  // Note: For single-plane (greyscale) images, planar_config may remain unset so
  // we can't use it as a guide for checking TIFF validity...
  if (t.sample_format == MISSING_TIFF_DATA ||
      t.bits_per_sample == MISSING_TIFF_DATA ||
      t.data_type == 0 ||
      t.num_bands == MISSING_TIFF_DATA ||
      t.is_scanline_format == MISSING_TIFF_DATA ||
      t.height == 0 ||
      t.width == 0)
  {
    asfPrintError("Cannot read tif file\n");
  }

  // Get a float line from the file
  float *buf = (float*)MALLOC(t.width*sizeof(float));
  tiff_get_float_line(tif, buf, row, band_no);
  cs = buf[col];
  FREE(buf);

  // Return as float
  return cs;
}

void dump_tiff_rgb_palette(char *inFile)
{
    TIFF *tiff = NULL;
    tiff_data_t tiff_data;

    get_tiff_info_from_file(inFile, &tiff_data);

    tiff = XTIFFOpen(inFile, "rb");
    if (!tiff) asfPrintError("Cannot open TIFF file (%s)\n", inFile);

    int color;
    int num_colors=(int)powl(2,tiff_data.bits_per_sample);
    uint16* red=NULL;
    uint16* green=NULL;
    uint16* blue=NULL;
    uint16 rmin, gmin, bmin, min;
    uint16 rmax, gmax, bmax, max;

    TIFFGetField(tiff, TIFFTAG_COLORMAP, &red, &green, &blue);
    if (red == NULL || green == NULL || blue == NULL) asfPrintError("Cannot read palette from TIFF file (%s)\n", inFile);

    rmin=gmin=bmin=min=65535; // Max uint16
    rmax=gmax=bmax=max=0; // Min uint16
    for (color=0; color<num_colors; color++) {
        unsigned short r, g, b;

        r=red[color];
        g=green[color];
        b=blue[color];

        r = (unsigned short)(((float)r/65535)*255 + .5);
        g = (unsigned short)(((float)g/65535)*255 + .5);
        b = (unsigned short)(((float)b/65535)*255 + .5);

        rmin = (r < rmin) ? r : rmin;
        gmin = (g < gmin) ? g : gmin;
        bmin = (b < bmin) ? b : bmin;
        min = (r < min) ? r : min;
        min = (g < min) ? g : min;
        min = (b < min) ? b : min;

        rmax = (r > rmax) ? r : rmax;
        gmax = (g > gmax) ? g : gmax;
        bmax = (b > bmax) ? b : bmax;
        max = (r > max) ? r : max;
        max = (g > max) ? g : max;
        max = (b > max) ? b : max;
    }
    fprintf(stdout, "# red min = %03d, green min = %03d, blue min = %03d, overall min = %03d\n",
            rmin, gmin, bmin, min);
    fprintf(stdout, "# red max = %03d, green max = %03d, blue max = %03d, overall max = %03d\n",
            rmax, gmax, bmax, max);
    fprintf(stdout, "# \n# color\tred\tgreen\tblue\n");
    for (color=0; color<num_colors; color++) {
        unsigned short r, g, b;

        r=red[color];
        g=green[color];
        b=blue[color];

        r = (unsigned short)(((float)r/65535)*255 + .5);
        g = (unsigned short)(((float)g/65535)*255 + .5);
        b = (unsigned short)(((float)b/65535)*255 + .5);
        fprintf(stdout, "%03d\t%03d\t%03d\t%03d\n",
                color, r, g, b);
    }

    XTIFFClose(tiff);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

int main(void)
{
  size_t i, j;

  gsl_matrix *m = gsl_matrix_alloc(10,10);

  for (i=0; i<10; i++)
    for (j=0; j<10; j++)
      gsl_matrix_set(m, i, j, sin(i) + cos(j));

  for (j=0; j<10; j++)
  {
    gsl_vector_view column = gsl_matrix_column(m, j);
    double d;

    d = gsl_blas_dnrm2(&column.vector);

    printf("matrix column %zu, norm = %g\n", j, d);
  }

  gsl_matrix_free(m);

  return 0;
}

max_stars_repo_name: tcrundall/chronostar

max_issues_repo_name: tcrundall/chronostar

max_forks_repo_name: tcrundall/chronostar