text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:44:44 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -hc2hc-forward 5 */

/*
 * This function contains 64 FP additions, 40 FP multiplications,
 * (or, 44 additions, 20 multiplications, 20 fused multiply/add),
 * 27 stack variables, and 40 memory accesses
 */
static const fftw_real K250000000 = FFTW_KONST(+0.250000000000000000000000000000000000000000000);
static const fftw_real K559016994 = FFTW_KONST(+0.559016994374947424102293417182819058860154590);
static const fftw_real K587785252 = FFTW_KONST(+0.587785252292473129168705954639072768597652438);
static const fftw_real K951056516 = FFTW_KONST(+0.951056516295153572116439333379382143405698634);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_hc2hc_forward_5(fftw_real *A, const fftw_complex *W, int iostride, int m, int dist)
{
     int i;
     fftw_real *X;
     fftw_real *Y;
     X = A;
     Y = A + (5 * iostride);
     {
	  fftw_real tmp70;
	  fftw_real tmp67;
	  fftw_real tmp68;
	  fftw_real tmp63;
	  fftw_real tmp71;
	  fftw_real tmp66;
	  fftw_real tmp69;
	  fftw_real tmp72;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp70 = X[0];
	  {
	       fftw_real tmp61;
	       fftw_real tmp62;
	       fftw_real tmp64;
	       fftw_real tmp65;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp61 = X[4 * iostride];
	       tmp62 = X[iostride];
	       tmp67 = tmp62 + tmp61;
	       tmp64 = X[2 * iostride];
	       tmp65 = X[3 * iostride];
	       tmp68 = tmp64 + tmp65;
	       tmp63 = tmp61 - tmp62;
	       tmp71 = tmp67 + tmp68;
	       tmp66 = tmp64 - tmp65;
	  }
	  Y[-iostride] = (K951056516 * tmp63) - (K587785252 * tmp66);
	  Y[-2 * iostride] = (K587785252 * tmp63) + (K951056516 * tmp66);
	  X[0] = tmp70 + tmp71;
	  tmp69 = K559016994 * (tmp67 - tmp68);
	  tmp72 = tmp70 - (K250000000 * tmp71);
	  X[iostride] = tmp69 + tmp72;
	  X[2 * iostride] = tmp72 - tmp69;
     }
     X = X + dist;
     Y = Y - dist;
     for (i = 2; i < m; i = i + 2, X = X + dist, Y = Y - dist, W = W + 4) {
	  fftw_real tmp13;
	  fftw_real tmp52;
	  fftw_real tmp42;
	  fftw_real tmp45;
	  fftw_real tmp49;
	  fftw_real tmp50;
	  fftw_real tmp51;
	  fftw_real tmp54;
	  fftw_real tmp53;
	  fftw_real tmp24;
	  fftw_real tmp35;
	  fftw_real tmp36;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp13 = X[0];
	  tmp52 = Y[-4 * iostride];
	  {
	       fftw_real tmp18;
	       fftw_real tmp40;
	       fftw_real tmp34;
	       fftw_real tmp44;
	       fftw_real tmp23;
	       fftw_real tmp41;
	       fftw_real tmp29;
	       fftw_real tmp43;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp15;
		    fftw_real tmp17;
		    fftw_real tmp14;
		    fftw_real tmp16;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp15 = X[iostride];
		    tmp17 = Y[-3 * iostride];
		    tmp14 = c_re(W[0]);
		    tmp16 = c_im(W[0]);
		    tmp18 = (tmp14 * tmp15) - (tmp16 * tmp17);
		    tmp40 = (tmp16 * tmp15) + (tmp14 * tmp17);
	       }
	       {
		    fftw_real tmp31;
		    fftw_real tmp33;
		    fftw_real tmp30;
		    fftw_real tmp32;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp31 = X[3 * iostride];
		    tmp33 = Y[-iostride];
		    tmp30 = c_re(W[2]);
		    tmp32 = c_im(W[2]);
		    tmp34 = (tmp30 * tmp31) - (tmp32 * tmp33);
		    tmp44 = (tmp32 * tmp31) + (tmp30 * tmp33);
	       }
	       {
		    fftw_real tmp20;
		    fftw_real tmp22;
		    fftw_real tmp19;
		    fftw_real tmp21;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp20 = X[4 * iostride];
		    tmp22 = Y[0];
		    tmp19 = c_re(W[3]);
		    tmp21 = c_im(W[3]);
		    tmp23 = (tmp19 * tmp20) - (tmp21 * tmp22);
		    tmp41 = (tmp21 * tmp20) + (tmp19 * tmp22);
	       }
	       {
		    fftw_real tmp26;
		    fftw_real tmp28;
		    fftw_real tmp25;
		    fftw_real tmp27;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp26 = X[2 * iostride];
		    tmp28 = Y[-2 * iostride];
		    tmp25 = c_re(W[1]);
		    tmp27 = c_im(W[1]);
		    tmp29 = (tmp25 * tmp26) - (tmp27 * tmp28);
		    tmp43 = (tmp27 * tmp26) + (tmp25 * tmp28);
	       }
	       tmp42 = tmp40 - tmp41;
	       tmp45 = tmp43 - tmp44;
	       tmp49 = tmp40 + tmp41;
	       tmp50 = tmp43 + tmp44;
	       tmp51 = tmp49 + tmp50;
	       tmp54 = tmp29 - tmp34;
	       tmp53 = tmp18 - tmp23;
	       tmp24 = tmp18 + tmp23;
	       tmp35 = tmp29 + tmp34;
	       tmp36 = tmp24 + tmp35;
	  }
	  X[0] = tmp13 + tmp36;
	  {
	       fftw_real tmp46;
	       fftw_real tmp48;
	       fftw_real tmp39;
	       fftw_real tmp47;
	       fftw_real tmp37;
	       fftw_real tmp38;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp46 = (K951056516 * tmp42) + (K587785252 * tmp45);
	       tmp48 = (K951056516 * tmp45) - (K587785252 * tmp42);
	       tmp37 = K559016994 * (tmp24 - tmp35);
	       tmp38 = tmp13 - (K250000000 * tmp36);
	       tmp39 = tmp37 + tmp38;
	       tmp47 = tmp38 - tmp37;
	       Y[-4 * iostride] = tmp39 - tmp46;
	       X[iostride] = tmp39 + tmp46;
	       X[2 * iostride] = tmp47 - tmp48;
	       Y[-3 * iostride] = tmp47 + tmp48;
	  }
	  Y[0] = tmp51 + tmp52;
	  {
	       fftw_real tmp55;
	       fftw_real tmp60;
	       fftw_real tmp58;
	       fftw_real tmp59;
	       fftw_real tmp56;
	       fftw_real tmp57;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp55 = (K951056516 * tmp53) + (K587785252 * tmp54);
	       tmp60 = (K951056516 * tmp54) - (K587785252 * tmp53);
	       tmp56 = K559016994 * (tmp49 - tmp50);
	       tmp57 = tmp52 - (K250000000 * tmp51);
	       tmp58 = tmp56 + tmp57;
	       tmp59 = tmp57 - tmp56;
	       X[4 * iostride] = -(tmp55 + tmp58);
	       Y[-iostride] = tmp58 - tmp55;
	       X[3 * iostride] = -(tmp59 - tmp60);
	       Y[-2 * iostride] = tmp60 + tmp59;
	  }
     }
     if (i == m) {
	  fftw_real tmp8;
	  fftw_real tmp3;
	  fftw_real tmp6;
	  fftw_real tmp9;
	  fftw_real tmp12;
	  fftw_real tmp11;
	  fftw_real tmp7;
	  fftw_real tmp10;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp8 = X[0];
	  {
	       fftw_real tmp1;
	       fftw_real tmp2;
	       fftw_real tmp4;
	       fftw_real tmp5;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp1 = X[2 * iostride];
	       tmp2 = X[3 * iostride];
	       tmp3 = tmp1 - tmp2;
	       tmp4 = X[4 * iostride];
	       tmp5 = X[iostride];
	       tmp6 = tmp4 - tmp5;
	       tmp9 = tmp3 + tmp6;
	       tmp12 = tmp4 + tmp5;
	       tmp11 = tmp1 + tmp2;
	  }
	  X[2 * iostride] = tmp8 + tmp9;
	  tmp7 = K559016994 * (tmp3 - tmp6);
	  tmp10 = tmp8 - (K250000000 * tmp9);
	  X[0] = tmp7 + tmp10;
	  X[iostride] = tmp10 - tmp7;
	  Y[0] = -((K951056516 * tmp11) + (K587785252 * tmp12));
	  Y[-iostride] = -((K951056516 * tmp12) - (K587785252 * tmp11));
     }
}

static const int twiddle_order[] =
{1, 2, 3, 4};
fftw_codelet_desc fftw_hc2hc_forward_5_desc =
{
     "fftw_hc2hc_forward_5",
     (void (*)()) fftw_hc2hc_forward_5,
     5,
     FFTW_FORWARD,
     FFTW_HC2HC,
     113,
     4,
     twiddle_order,
};


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

/* gcc gsl_mersenne_twister.c -lgsl -lcblas */

int main(int argc, char **argv){
  const gsl_rng_type * T;
  gsl_rng * r;
  int i, n_points;
  double x;

  T = gsl_rng_mt19937;
  r = gsl_rng_alloc (T);

  n_points  = atoi(argv[1]);
  for(i=0;i<n_points;i++){
    x = gsl_rng_uniform(r);
    fprintf(stdout, "%f\n", x);
  }

  return 0;
}

ext: c


text:
#ifndef __SVEC_H__
#define __SVEC_H__

#define HAVE_INLINE
#include <gsl/gsl_vector.h>

typedef gsl_vector_float *svec;

#define svec_alloc gsl_vector_float_alloc
#define svec_calloc gsl_vector_float_calloc
#define svec_free gsl_vector_float_free
#define svec_get gsl_vector_float_get
#define svec_set gsl_vector_float_set
#define svec_scale gsl_vector_float_scale
#define svec_memcpy gsl_vector_float_memcpy
#define svec_set_zero gsl_vector_float_set_zero
#define svec_add gsl_vector_float_add

void svec_print(svec x);
void svec_mul_print(svec x, double k);
void svec_randomize(svec x);
void svec_normalize(svec x);
void svec_zero(svec x);
float svec_pull(svec x, svec y, float d);
float svec_push(svec x, svec y, float d);
float svec_sqdist(svec x, svec y);

#endif

ext: h

max_stars_repo_name: ai-ku/scode

max_issues_repo_name: ai-ku/scode

max_forks_repo_name: ai-ku/scode


text:
#include <assert.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>

#include <asf.h>
#include "asf_tiff.h"

#include <gsl/gsl_math.h>
#include <proj_api.h>

#include "asf_jpeg.h"
#include <png.h>
#include "envi.h"

#include "dateUtil.h"
#include <time.h>
#include "matrix.h"
#include <asf_nan.h>
#include <asf_endian.h>
#include <asf_meta.h>
#include <asf_export.h>
#include <asf_raster.h>
#include <float_image.h>
#include <spheroids.h>
#include <typlim.h>
#include <hdf5.h>

#define RES 16
#define MAX_PTS 256

void h5_att_double(hid_t data, hid_t space, char *name, double value)
{
  hid_t attr = H5Acreate(data, name, H5T_NATIVE_DOUBLE, space, H5P_DEFAULT,
			 H5P_DEFAULT);
  H5Awrite(attr, H5T_NATIVE_DOUBLE, &value);
  H5Aclose(attr);
}

void h5_att_float(hid_t data, hid_t space, char *name, float value)
{
  hid_t attr = H5Acreate(data, name, H5T_NATIVE_FLOAT, space, H5P_DEFAULT,
			 H5P_DEFAULT);
  H5Awrite(attr, H5T_NATIVE_FLOAT, &value);
  H5Aclose(attr);
}

void h5_att_int(hid_t data, hid_t space, char *name, int value)
{
  hid_t attr = H5Acreate(data, name, H5T_NATIVE_INT, space, H5P_DEFAULT,
			 H5P_DEFAULT);
  H5Awrite(attr, H5T_NATIVE_INT, &value);
  H5Aclose(attr);
}

void h5_att_float2(hid_t data, hid_t space, char *name, float *value)
{
  hid_t attr = H5Acreate(data, name, H5T_NATIVE_FLOAT, space, H5P_DEFAULT,
			 H5P_DEFAULT);
  H5Awrite(attr, H5T_NATIVE_FLOAT, value);
  H5Aclose(attr);
}

void h5_att_str(hid_t data, hid_t space, char *name, char *value)
{
  hid_t str = H5Tcopy(H5T_C_S1);
  H5Tset_size(str, strlen(value));
  hid_t attr = H5Acreate(data, name, str, space, H5P_DEFAULT, H5P_DEFAULT);
  H5Awrite(attr, str, value);
  H5Aclose(attr);
}

void h5_value_double(hid_t file, char *group, char *name,
		     double value, char *long_name, char *units)
{
  char meta[100];
  sprintf(meta, "%s/%s", group, name);
  hid_t h5_space = H5Screate(H5S_SCALAR);
  hid_t h5_data = H5Dcreate(file, meta, H5T_NATIVE_DOUBLE, h5_space,
			    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(h5_data, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
  h5_att_str(h5_data, h5_space, "long_name", long_name);
  if (units && strlen(units) > 0)
    h5_att_str(h5_data, h5_space, "units", units);
  H5Dclose(h5_data);
  H5Sclose(h5_space);
}

void h5_value_float(hid_t file, char *group, char *name,
		    float value, char *long_name, char *units)
{
  char meta[100];
  sprintf(meta, "%s/%s", group, name);
  hid_t h5_space = H5Screate(H5S_SCALAR);
  hid_t h5_data = H5Dcreate(file, meta, H5T_NATIVE_FLOAT, h5_space,
			    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(h5_data, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
  h5_att_str(h5_data, h5_space, "long_name", long_name);
  if (units && strlen(units) > 0)
    h5_att_str(h5_data, h5_space, "units", units);
  H5Dclose(h5_data);
  H5Sclose(h5_space);
}

void h5_value_int(hid_t file, char *group, char *name,
		  int value, char *long_name, char *units)
{
  char meta[100];
  sprintf(meta, "%s/%s", group, name);
  hid_t h5_space = H5Screate(H5S_SCALAR);
  hid_t h5_data = H5Dcreate(file, meta, H5T_NATIVE_INT, h5_space,
			    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(h5_data, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
  h5_att_str(h5_data, h5_space, "long_name", long_name);
  if (units && strlen(units) > 0)
    h5_att_str(h5_data, h5_space, "units", units);
  H5Dclose(h5_data);
  H5Sclose(h5_space);
}

void h5_value_str(hid_t file, char *group, char *name,
		  char *value, char *long_name, char *units)
{
  char meta[100];
  sprintf(meta, "%s/%s", group, name);
  hid_t h5_space = H5Screate(H5S_SCALAR);
  hid_t h5_str = H5Tcopy(H5T_C_S1);
  H5Tset_size(h5_str, strlen(value));
  hid_t h5_data = H5Dcreate(file, meta, h5_str, h5_space,
			 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(h5_data, h5_str, H5S_ALL, H5S_ALL, H5P_DEFAULT, value);
  h5_att_str(h5_data, h5_space, "long_name", long_name);
  if (units && strlen(units) > 0)
    h5_att_str(h5_data, h5_space, "units", units);
  H5Dclose(h5_data);
  H5Sclose(h5_space);
}

h5_t *initialize_h5_file(const char *output_file_name, meta_parameters *md)
{
  hid_t h5_file, h5_datagroup, h5_metagroup, h5_data, h5_proj;
  hid_t h5_array, h5_string, h5_time, h5_lat, h5_lon, h5_xgrid, h5_ygrid;
  int ii, kk, complex=FALSE, projected=FALSE;
  char *spatial_ref=NULL, *datum=NULL, *spheroid=NULL;
  char dataset[50], group[50], band[5], str_attr[50], tmp[50];
  double lfValue;

  // Convenience variables
  meta_general *mg = md->general;
  meta_sar *ms = md->sar;
  meta_state_vectors *mo = md->state_vectors;
  meta_projection *mp = md->projection;

  // Check whether data is map projected
  if (mp && mp->type != SCANSAR_PROJECTION)
    //asfPrintError("Image is map projected. Wrong initialization function!\n");
    projected = TRUE;

  // Initialize the HDF pointer structure
  h5_t *h5 = (h5_t *) MALLOC(sizeof(h5_t));
  int band_count = mg->band_count;
  h5->var_count = band_count;
  h5->var = (hid_t *) MALLOC(sizeof(hid_t)*h5->var_count);

  // Check for complex data
  if (mg->image_data_type == COMPLEX_IMAGE)
    complex = TRUE;

  // Create new HDF5 file
  h5_file =
    H5Fcreate(output_file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
  h5->file = h5_file;

  // Create data space
  int samples = mg->sample_count;
  int lines = mg->line_count;
  hsize_t dims[2] = { lines, samples };
  hsize_t cdims[2] = { 100, samples };
  hsize_t rdims[2] = { 1, 2 };
  h5_array = H5Screate_simple(2, dims, NULL);
  h5->space = h5_array;
  h5_string = H5Screate(H5S_SCALAR);
  hid_t h5_range = H5Screate(H5S_SIMPLE);
  H5Sset_extent_simple(h5_range, 2, rdims, NULL);

  // Create data structure
  char **band_name = extract_band_names(mg->bands, band_count);
  hid_t h5_plist = H5Pcreate(H5P_DATASET_CREATE);
  H5Pset_chunk(h5_plist, 2, cdims);
  H5Pset_deflate(h5_plist, 6);

  // Create a data group
  sprintf(group, "/data");
  h5_datagroup = H5Gcreate(h5_file, group, H5P_DEFAULT, H5P_DEFAULT,
			   H5P_DEFAULT);
  // Projection information
  if (projected) {
    sprintf(group, "/data/projection");
    h5_proj = H5Gcreate(h5_file, group, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

    if (mp->type == UNIVERSAL_TRANSVERSE_MERCATOR) {

      h5_att_str(h5_proj, h5_string, "grid_mapping_name",
		 "transverse_mercator");
      h5_att_double(h5_proj, h5_string, "scale_factor_at_central_meridian",
		    mp->param.utm.scale_factor);
      h5_att_double(h5_proj, h5_string, "longitude_of_central_meridian",
		    mp->param.utm.lon0);
      h5_att_double(h5_proj, h5_string, "latitude_of_projection_origin",
		    mp->param.utm.lat0);
      h5_att_double(h5_proj, h5_string, "false_easting",
		    mp->param.utm.false_easting);
      h5_att_double(h5_proj, h5_string, "false_northing",
		    mp->param.utm.false_northing);
      h5_att_str(h5_proj, h5_string, "projection_x_coordinate", "xgrid");
      h5_att_str(h5_proj, h5_string, "projection_y_coordinate", "ygrid");
      h5_att_str(h5_proj, h5_string, "units", "meters");
      h5_att_double(h5_proj, h5_string, "grid_boundary_top_projected_y",
		    mp->startY);
      lfValue = mp->startY + mg->line_count * mp->perY;
      h5_att_double(h5_proj, h5_string, "grid_boundary_bottom_projected_y",
		    lfValue);
      lfValue = mp->startX + mg->sample_count * mp->perX;
      h5_att_double(h5_proj, h5_string, "grid_boundary_right_projected_x",
		    lfValue);
      h5_att_double(h5_proj, h5_string, "grid_boundary_left_projected_x",
		    mp->startX);
      spatial_ref = (char *) MALLOC(sizeof(char)*1024);
      datum = (char *) datum_toString(mp->datum);
      spheroid = (char *) spheroid_toString(mp->spheroid);
      double flat = mp->re_major/(mp->re_major - mp->re_minor);
      sprintf(spatial_ref, "PROJCS[\"%s_UTM_Zone_%d%c\",GEOGCS[\"GCS_%s\",DATUM[\"D_%s\",SPHEROID[\"%s\",%.1lf,%-16.11g]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",%.1lf],PARAMETER[\"False_Northing\",%.1lf],PARAMETER[\"Central_Meridian\",%.1lf],PARAMETER[\"Scale_Factor\",%.4lf],PARAMETER[\"Latitude_Of_Origin\",%.1lf],UNIT[\"Meter\",1]]",
	      spheroid, mp->param.utm.zone, mp->hem, spheroid, datum,
	      spheroid, mp->re_major, flat, mp->param.utm.false_easting,
	      mp->param.utm.false_northing, mp->param.utm.lon0,
	      mp->param.utm.scale_factor, mp->param.utm.lat0);
      h5_att_str(h5_proj, h5_string, "spatial_ref", spatial_ref);
      sprintf(str_attr, "+proj=utm +zone=%d", mp->param.utm.zone);
      if (mg->center_latitude < 0)
	strcat(str_attr, " +south");
      h5_att_str(h5_proj, h5_string, "proj4text", str_attr);
      h5_att_int(h5_proj, h5_string, "zone", mp->param.utm.zone);
      h5_att_double(h5_proj, h5_string, "semimajor_radius", mp->re_major);
      h5_att_double(h5_proj, h5_string, "semiminor_radius", mp->re_minor);
      sprintf(str_attr, "%.6lf %.6lf 0 %.6lf 0 %.6lf", mp->startX, mp->perX,
	      mp->startY, mp->perY);
      h5_att_str(h5_proj, h5_string, "GeoTransform", str_attr);
    }
    else if (mp->type == POLAR_STEREOGRAPHIC) {

      h5_att_str(h5_proj, h5_string, "grid_mapping_name",
		 "polar_stereographic");
      h5_att_double(h5_proj, h5_string, "straight_vertical_longitude_from_pole",
		    mp->param.ps.slon);
      h5_att_double(h5_proj, h5_string, "latitude_of_projection_origin",
		    mp->param.ps.slat);
      h5_att_double(h5_proj, h5_string, "scale_factor_at_projection_origin",
		    1.0);
      h5_att_double(h5_proj, h5_string, "false_easting",
		    mp->param.ps.false_easting);
      h5_att_double(h5_proj, h5_string, "false_northing",
		    mp->param.ps.false_northing);
      h5_att_str(h5_proj, h5_string, "projection_x_coordinate", "xgrid");
      h5_att_str(h5_proj, h5_string, "projection_y_coordinate", "ygrid");
      h5_att_str(h5_proj, h5_string, "units", "meters");
      h5_att_double(h5_proj, h5_string, "grid_boundary_top_projected_y",
		    mp->startY);
      lfValue = mp->startY + mg->line_count * mp->perY;
      h5_att_double(h5_proj, h5_string, "grid_boundary_bottom_projected_y",
		    lfValue);
      lfValue = mp->startX + mg->sample_count * mp->perX;
      h5_att_double(h5_proj, h5_string, "grid_boundary_right_projected_x",
		    lfValue);
      h5_att_double(h5_proj, h5_string, "grid_boundary_left_projected_x",
		    mp->startX);
      spatial_ref = (char *) MALLOC(sizeof(char)*1024);
      double flat = mp->re_major/(mp->re_major - mp->re_minor);
      sprintf(spatial_ref, "PROJCS[\"Stereographic_North_Pole\",GEOGCS[\"unnamed ellipse\",DATUM[\"D_unknown\",SPHEROID[\"Unknown\",%.3lf,%-16.11g]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0002247191011236]],PROJECTION[\"Stereographic_North_Pole\"],PARAMETER[\"standard_parallel_1\",%.4lf],PARAMETER[\"central_meridian\",%.4lf],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",%.1lf],PARAMETER[\"false_northing\",%.1lf],UNIT[\"Meter\",1, AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"3411\"]]",
	      mp->re_major, flat, mp->param.ps.slat, mp->param.ps.slon,
	      mp->param.ps.false_easting, mp->param.ps.false_northing);
      h5_att_str(h5_proj, h5_string, "spatial_ref", spatial_ref);
      if (mp->param.ps.is_north_pole)
	sprintf(str_attr, "+proj=stere +lat_0=90.0000 +lat_ts=%.4lf "
		"+lon_0=%.4lf +k=1 +x_0=%.3lf +y_0=%.3lf +a=%.3lf +b=%.3lf "
		"+units=m +no_defs", mp->param.ps.slat, mp->param.ps.slon,
		mp->param.ps.false_easting, mp->param.ps.false_northing,
		mp->re_major, mp->re_minor);
      else
	sprintf(str_attr, "+proj=stere +lat_0=-90.0000 +lat_ts=%.4lf "
		"+lon_0=%.4lf +k=1 +x_0=%.3lf +y_0=%.3lf +a=%.3lf +b=%.3lf "
		"+units=m +no_defs", mp->param.ps.slat, mp->param.ps.slon,
		mp->param.ps.false_easting, mp->param.ps.false_northing,
		mp->re_major, mp->re_minor);
      h5_att_str(h5_proj, h5_string, "proj4text", str_attr);
      h5_att_double(h5_proj, h5_string, "semimajor_radius", mp->re_major);
      h5_att_double(h5_proj, h5_string, "semiminor_radius", mp->re_minor);
      sprintf(str_attr, "%.6lf %.6lf 0 %.6lf 0 %.6lf", mp->startX, mp->perX,
	      mp->startY, mp->perY);
      h5_att_str(h5_proj, h5_string, "GeoTransform", str_attr);
    }
    else if (mp->type == ALBERS_EQUAL_AREA) {

      h5_att_str(h5_proj, h5_string, "grid_mapping_name",
		 "albers_conical_equal_area");
      h5_att_double(h5_proj, h5_string, "standard_parallel_1",
		    mp->param.albers.std_parallel1);
      h5_att_double(h5_proj, h5_string, "standard_parallel_2",
		    mp->param.albers.std_parallel2);
      h5_att_double(h5_proj, h5_string, "longitude_of_central_meridian",
		    mp->param.albers.center_meridian);
      h5_att_double(h5_proj, h5_string, "latitude_of_projection_origin",
		    mp->param.albers.orig_latitude);
      h5_att_double(h5_proj, h5_string, "false_easting",
		    mp->param.albers.false_easting);
      h5_att_double(h5_proj, h5_string, "false_northing",
		    mp->param.albers.false_northing);
      h5_att_str(h5_proj, h5_string, "projection_x_coordinate", "xgrid");
      h5_att_str(h5_proj, h5_string, "projection_y_coordinate", "ygrid");
      h5_att_str(h5_proj, h5_string, "units", "meters");
      h5_att_double(h5_proj, h5_string, "grid_boundary_top_projected_y",
		    mp->startY);
      lfValue = mp->startY + mg->line_count * mp->perY;
      h5_att_double(h5_proj, h5_string, "grid_boundary_bottom_projected_y",
		    lfValue);
      lfValue = mp->startX + mg->sample_count * mp->perX;
      h5_att_double(h5_proj, h5_string, "grid_boundary_right_projected_x",
		    lfValue);
      h5_att_double(h5_proj, h5_string, "grid_boundary_left_projected_x",
		    mp->startX);
      spatial_ref = (char *) MALLOC(sizeof(char)*1024);
      datum = (char *) datum_toString(mp->datum);
      spheroid = (char *) spheroid_toString(mp->spheroid);
      double flat = mp->re_major/(mp->re_major - mp->re_minor);
      sprintf(spatial_ref, "PROJCS[\"Albers_Equal_Area_Conic\",GEOGCS[\"GCS_%s\",DATUM[\"D_%s\",SPHEROID[\"%s\",%.3lf,%-16.11g]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199432955]],PROJECTION[\"Albers\"],PARAMETER[\"False_Easting\",%.3lf],PARAMETER[\"False_Northing\",%.3lf],PARAMETER[\"Central_Meridian\",%.4lf],PARAMETER[\"Standard_Parallel_1\",%.4lf],PARAMETER[\"Standard_Parallel_2\",%.4lf],PARAMETER[\"Latitude_Of_Origin\",%.4lf],UNIT[\"Meter\",1]]",
	      datum, datum, spheroid, mp->re_major, flat,
	      mp->param.albers.false_easting, mp->param.albers.false_northing,
	      mp->param.albers.center_meridian, mp->param.albers.std_parallel1,
	      mp->param.albers.std_parallel2, mp->param.albers.orig_latitude);
      h5_att_str(h5_proj, h5_string, "spatial_ref", spatial_ref);
      sprintf(str_attr, "+proj=aea +lat_1=%.4lf +lat_2=%.4lf +lat_0=%.4lf "
	      "+lon_0=%.4lf +x_0=%.3lf +y_0=%.3lf",
	      mp->param.albers.std_parallel1, mp->param.albers.std_parallel2,
	      mp->param.albers.orig_latitude, mp->param.albers.center_meridian,
	      mp->param.albers.false_easting, mp->param.albers.false_northing);
      h5_att_str(h5_proj, h5_string, "proj4text", str_attr);
      h5_att_double(h5_proj, h5_string, "semimajor_radius", mp->re_major);
      h5_att_double(h5_proj, h5_string, "semiminor_radius", mp->re_minor);
      sprintf(str_attr, "%.6lf %.6lf 0 %.6lf 0 %.6lf", mp->startX, mp->perX,
	      mp->startY, mp->perY);
      h5_att_str(h5_proj, h5_string, "GeoTransform", str_attr);
    }
    else if (mp->type == LAMBERT_CONFORMAL_CONIC) {

      h5_att_str(h5_proj, h5_string, "grid_mapping_name",
		 "lambert_conformal_conic");
      h5_att_double(h5_proj, h5_string, "standard_parallel_1",
		    mp->param.lamcc.plat1);
      h5_att_double(h5_proj, h5_string, "standard_parallel_2",
		    mp->param.lamcc.plat2);
      h5_att_double(h5_proj, h5_string, "longitude_of_central_meridian",
		    mp->param.lamcc.lon0);
      h5_att_double(h5_proj, h5_string, "latitude_of_projection_origin",
		    mp->param.lamcc.lat0);
      h5_att_double(h5_proj, h5_string, "false_easting",
		    mp->param.lamcc.false_easting);
      h5_att_double(h5_proj, h5_string, "false_northing",
		    mp->param.lamcc.false_northing);
      h5_att_str(h5_proj, h5_string, "projection_x_coordinate", "xgrid");
      h5_att_str(h5_proj, h5_string, "projection_y_coordinate", "ygrid");
      h5_att_str(h5_proj, h5_string, "units", "meters");
      h5_att_double(h5_proj, h5_string, "grid_boundary_top_projected_y",
		    mp->startY);
      lfValue = mp->startY + mg->line_count * mp->perY;
      h5_att_double(h5_proj, h5_string, "grid_boundary_bottom_projected_y",
		    lfValue);
      lfValue = mp->startX + mg->sample_count * mp->perX;
      h5_att_double(h5_proj, h5_string, "grid_boundary_right_projected_x",
		    lfValue);
      h5_att_double(h5_proj, h5_string, "grid_boundary_left_projected_x",
		    mp->startX);
      spatial_ref = (char *) MALLOC(sizeof(char)*1024);
      datum = (char *) datum_toString(mp->datum);
      spheroid = (char *) spheroid_toString(mp->spheroid);
      double flat = mp->re_major/(mp->re_major - mp->re_minor);
      sprintf(spatial_ref, "PROJCS[\"Lambert_Conformal_Conic\",GEOGCS[\"GCS_%s\",DATUM[\"D_%s\",SPHEROID[\"%s\",%.3lf,%-16.11g]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199432955]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",%.3lf],PARAMETER[\"False_Northing\",%.3lf],PARAMETER[\"Central_Meridian\",%.4lf],PARAMETER[\"Standard_Parallel_1\",%.4lf],PARAMETER[\"Standard_Parallel_2\",%.4lf],PARAMETER[\"Latitude_Of_Origin\",%.4lf],UNIT[\"Meter\",1]]",
	      datum, datum, spheroid, mp->re_major, flat,
	      mp->param.lamcc.false_easting, mp->param.lamcc.false_northing,
	      mp->param.lamcc.lon0, mp->param.lamcc.plat1,
	      mp->param.lamcc.plat2, mp->param.lamcc.lat0);
      h5_att_str(h5_proj, h5_string, "spatial_ref", spatial_ref);
      sprintf(str_attr, "+proj=lcc +lat_1=%.4lf +lat_2=%.4lf +lat_0=%.4lf "
	      "+lon_0=%.4lf +x_0=%.3lf +y_0=%.3lf",
	      mp->param.lamcc.plat1, mp->param.lamcc.plat2,
	      mp->param.lamcc.lat0, mp->param.lamcc.lon0,
	      mp->param.lamcc.false_easting, mp->param.lamcc.false_northing);
      h5_att_str(h5_proj, h5_string, "proj4text", str_attr);
      h5_att_double(h5_proj, h5_string, "semimajor_radius", mp->re_major);
      h5_att_double(h5_proj, h5_string, "semiminor_radius", mp->re_minor);
      sprintf(str_attr, "%.6lf %.6lf 0 %.6lf 0 %.6lf", mp->startX, mp->perX,
	      mp->startY, mp->perY);
      h5_att_str(h5_proj, h5_string, "GeoTransform", str_attr);
    }
    else if (mp->type == LAMBERT_AZIMUTHAL_EQUAL_AREA) {

      h5_att_str(h5_proj, h5_string, "grid_mapping_name",
		 "lambert_azimuthal_equal_area");
      h5_att_double(h5_proj, h5_string, "longitude_of_projection_origin",
		    mp->param.lamaz.center_lon);
      h5_att_double(h5_proj, h5_string, "latitude_of_projection_origin",
		    mp->param.lamaz.center_lat);
      h5_att_double(h5_proj, h5_string, "false_easting",
		    mp->param.lamaz.false_easting);
      h5_att_double(h5_proj, h5_string, "false_northing",
		    mp->param.lamaz.false_northing);
      h5_att_str(h5_proj, h5_string, "projection_x_coordinate", "xgrid");
      h5_att_str(h5_proj, h5_string, "projection_y_coordinate", "ygrid");
      h5_att_str(h5_proj, h5_string, "units", "meters");
      h5_att_double(h5_proj, h5_string, "grid_boundary_top_projected_y",
		    mp->startY);
      lfValue = mp->startY + mg->line_count * mp->perY;
      h5_att_double(h5_proj, h5_string, "grid_boundary_bottom_projected_y",
		    lfValue);
      lfValue = mp->startX + mg->sample_count * mp->perX;
      h5_att_double(h5_proj, h5_string, "grid_boundary_right_projected_x",
		    lfValue);
      h5_att_double(h5_proj, h5_string, "grid_boundary_left_projected_x",
		    mp->startX);
      spatial_ref = (char *) MALLOC(sizeof(char)*1024);
      datum = (char *) datum_toString(mp->datum);
      spheroid = (char *) spheroid_toString(mp->spheroid);
      double flat = mp->re_major/(mp->re_major - mp->re_minor);
      sprintf(spatial_ref, "PROJCS[\"Lambert_Azimuthal_Equal_Area\",GEOGCS[\"GCS_%s\",DATUM[\"D_%s\",SPHEROID[\"%s\",%.3lf,%-16.11g]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199432955]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",%.3lf],PARAMETER[\"False_Northing\",%.3lf],PARAMETER[\"Central_Meridian\",%.4lf],PARAMETER[\"Latitude_Of_Origin\",%.4lf],UNIT[\"Meter\",1]]",
	      datum, datum, spheroid, mp->re_major, flat,
	      mp->param.lamaz.false_easting, mp->param.lamaz.false_northing,
	      mp->param.lamaz.center_lon, mp->param.lamaz.center_lat);
      h5_att_str(h5_proj, h5_string, "spatial_ref", spatial_ref);
      sprintf(str_attr, "+proj=laea +lat_0=%.4lf +lon_0=%.4lf +x_0=%.3lf "
	      "+y_0=%.3lf",
	      mp->param.lamaz.center_lat, mp->param.lamaz.center_lon,
	      mp->param.lamaz.false_easting, mp->param.lamaz.false_northing);
      h5_att_str(h5_proj, h5_string, "proj4text", str_attr);
      h5_att_double(h5_proj, h5_string, "semimajor_radius", mp->re_major);
      h5_att_double(h5_proj, h5_string, "semiminor_radius", mp->re_minor);
      sprintf(str_attr, "%.6lf %.6lf 0 %.6lf 0 %.6lf", mp->startX, mp->perX,
	      mp->startY, mp->perY);
      h5_att_str(h5_proj, h5_string, "GeoTransform", str_attr);
    }
  }

  for (ii=0; ii<band_count; ii++) {

    // Create data set
    strncpy(band, band_name[ii], 2);
    band[2] = '\0';
    sprintf(dataset, "/data/%s_AMPLITUDE_IMAGE", band);
    h5_data = H5Dcreate(h5_file, dataset, H5T_NATIVE_FLOAT, h5_array,
			H5P_DEFAULT, h5_plist, H5P_DEFAULT);
    h5->var[ii] = h5_data;

    // Add attributes (from CF convention)
    sprintf(str_attr, "%s", mg->sensor);
    if (mg->image_data_type < 9)
      strcat(str_attr, " radar backscatter");
    if (mg->radiometry >= r_SIGMA_DB && mg->radiometry <= r_GAMMA_DB)
      strcat(str_attr, " in dB");
    h5_att_str(h5_data, h5_string, "long_name", str_attr);
    h5_att_str(h5_data, h5_string, "cell_methods", "area: backscatter value");
    h5_att_str(h5_data, h5_string, "units", "1");
    h5_att_str(h5_data, h5_string, "units_description",
	       "unitless normalized radar cross-section");
    if (mg->radiometry >= r_SIGMA && mg->radiometry <= r_GAMMA)
      strcat(str_attr, " stored as powerscale");
    else if (mg->radiometry >= r_SIGMA_DB && mg->radiometry <= r_GAMMA_DB)
      strcat(str_attr, " stored as dB=10*log10(*)");
    h5_att_float(h5_data, h5_string, "_FillValue", -999);
    h5_att_str(h5_data, h5_string, "coordinates", "longitude latitude");
    if (projected)
      h5_att_str(h5_data, h5_string, "grid_mapping", "projection");

    // Close up
    H5Dclose(h5_data);
  }

  // Extra bands - Time
  asfPrintStatus("Storing band 'time' ...\n");
  float serial_date = seconds_from_str(mg->acquisition_date);
  sprintf(dataset, "/data/time");
  h5_time = H5Dcreate(h5_file, dataset, H5T_NATIVE_FLOAT, h5_string,
		      H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  H5Dwrite(h5_time, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
	   &serial_date);
  h5_att_str(h5_time, h5_string, "units", "seconds since 1900-01-01T00:00:00Z");
  h5_att_str(h5_time, h5_string, "references", "scene center time");
  h5_att_str(h5_time, h5_string, "standard_name", "time");
  h5_att_str(h5_time, h5_string, "axis", "T");
  h5_att_str(h5_time, h5_string, "long_name", "serial date");
  H5Dclose(h5_time);

  // Extra bands - Longitude
  int nl = mg->line_count;
  int ns = mg->sample_count;
  long pixel_count = mg->line_count * mg->sample_count;
  double *value = (double *) MALLOC(sizeof(double)*MAX_PTS);
  double *l = (double *) MALLOC(sizeof(double)*MAX_PTS);
  double *s = (double *) MALLOC(sizeof(double)*MAX_PTS);
  double line, sample, lat, lon, first_value;
  float *lons = (float *) MALLOC(sizeof(float)*pixel_count);
  asfPrintStatus("Generating band 'longitude' ...\n");
  meta_get_latLon(md, 0, 0, 0.0, &lat, &lon);
  if (lon < 0.0)
    first_value = lon + 360.0;
  else
    first_value = lon;
  asfPrintStatus("Calculating grid for quadratic fit ...\n");
  for (ii=0; ii<RES; ii++) {
    for (kk=0; kk<RES; kk++) {
      line = ii * nl / RES;
      sample = kk * ns / RES;
      meta_get_latLon(md, line, sample, 0.0, &lat, &lon);
      l[ii*RES+kk] = line;
      s[ii*RES+kk] = sample;
      if (lon < 0.0)
	value[ii*RES+kk] = lon + 360.0;
      else
	value[ii*RES+kk] = lon;
    }
    asfLineMeter(ii, nl);
  }
  quadratic_2d q = find_quadratic(value, l, s, MAX_PTS);
  q.A = first_value;
  for (ii=0; ii<nl; ii++) {
    for (kk=0; kk<ns; kk++) {
      lons[ii*ns+kk] = (float)
	(q.A + q.B*ii + q.C*kk + q.D*ii*ii + q.E*ii*kk + q.F*kk*kk +
	 q.G*ii*ii*kk + q.H*ii*kk*kk + q.I*ii*ii*kk*kk + q.J*ii*ii*ii +
	 q.K*kk*kk*kk) - 360.0;
      if (lons[ii*ns+kk] < -180.0)
	lons[ii*ns+kk] += 360.0;
    }
    asfLineMeter(ii, nl);
  }
  asfPrintStatus("Storing band 'longitude' ...\n");
  sprintf(dataset, "/data/longitude");
  h5_lon = H5Dcreate(h5_file, dataset, H5T_NATIVE_FLOAT, h5_array,
		     H5P_DEFAULT, h5_plist, H5P_DEFAULT);
  H5Dwrite(h5_lon, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, lons);
  h5_att_str(h5_lon, h5_string, "units", "degrees_east");
  h5_att_str(h5_lon, h5_string, "long_name", "longitude");
  h5_att_str(h5_lon, h5_string, "standard_name", "longitude");
  float valid_range[2] = { -180.0, 180.0 };
  h5_att_float2(h5_lon, h5_range, "valid_range", valid_range);
  h5_att_float(h5_lon, h5_string, "_FillValue", -999);
  H5Dclose(h5_lon);
  FREE(lons);

  // Extra bands - Latitude
  float *lats = (float *) MALLOC(sizeof(float)*pixel_count);
  asfPrintStatus("Generating band 'latitude' ...\n");
  meta_get_latLon(md, 0, 0, 0.0, &lat, &lon);
  first_value = lat + 180.0;
  asfPrintStatus("Calculating grid for quadratic fit ...\n");
  for (ii=0; ii<RES; ii++) {
    for (kk=0; kk<RES; kk++) {
      line = ii * nl / RES;
      sample = kk * ns / RES;
      meta_get_latLon(md, line, sample, 0.0, &lat, &lon);
      l[ii*RES+kk] = line;
      s[ii*RES+kk] = sample;
      value[ii*RES+kk] = lat + 180.0;
    }
    asfLineMeter(ii, nl);
  }
  q = find_quadratic(value, l, s, MAX_PTS);
  q.A = first_value;
  for (ii=0; ii<nl; ii++) {
    for (kk=0; kk<ns; kk++) {
      if (mg->orbit_direction == 'A')
	lats[(nl-ii-1)*ns+kk] = (float)
	  (q.A + q.B*ii + q.C*kk + q.D*ii*ii + q.E*ii*kk + q.F*kk*kk +
	   q.G*ii*ii*kk + q.H*ii*kk*kk + q.I*ii*ii*kk*kk + q.J*ii*ii*ii +
	   q.K*kk*kk*kk) - 180.0;
      else
	lats[ii*ns+kk] = (float)
	  (q.A + q.B*ii + q.C*kk + q.D*ii*ii + q.E*ii*kk + q.F*kk*kk +
	   q.G*ii*ii*kk + q.H*ii*kk*kk + q.I*ii*ii*kk*kk + q.J*ii*ii*ii +
	   q.K*kk*kk*kk) - 180.0;
    }
    asfLineMeter(ii, nl);
  }
  asfPrintStatus("Storing band 'latitude' ...\n");
  sprintf(dataset, "/data/latitude");
  h5_lat = H5Dcreate(h5_file, dataset, H5T_NATIVE_FLOAT, h5_array,
		     H5P_DEFAULT, h5_plist, H5P_DEFAULT);
  H5Dwrite(h5_lat, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, lats);
  h5_att_str(h5_lat, h5_string, "units", "degrees_north");
  h5_att_str(h5_lat, h5_string, "long_name", "latitude");
  h5_att_str(h5_lat, h5_string, "standard_name", "latitude");
  valid_range[0] = -90.0;
  valid_range[1] = 90.0;
  h5_att_float2(h5_lat, h5_range, "valid_range", valid_range);
  h5_att_float(h5_lat, h5_string, "_FillValue", -999);
  H5Dclose(h5_lat);
  FREE(lats);

  if (projected) {
    // Extra bands - ygrid
    float *ygrids = (float *) MALLOC(sizeof(float)*pixel_count);
    for (ii=0; ii<nl; ii++) {
      for (kk=0; kk<ns; kk++)
	ygrids[ii*ns+kk] = mp->startY + kk*mp->perY;
      asfLineMeter(ii, nl);
    }
    asfPrintStatus("Storing band 'ygrid' ...\n");
    sprintf(dataset, "/data/ygrid");
    h5_ygrid = H5Dcreate(h5_file, dataset, H5T_NATIVE_FLOAT, h5_array,
			 H5P_DEFAULT, h5_plist, H5P_DEFAULT);
    H5Dwrite(h5_ygrid, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, ygrids);
    h5_att_str(h5_ygrid, h5_string, "units", "meters");
    h5_att_str(h5_ygrid, h5_string, "long_name",
	       "projection_grid_y_coordinates");
    h5_att_str(h5_ygrid, h5_string, "standard_name",
	       "projection_y_coordinates");
    h5_att_str(h5_ygrid, h5_string, "axis", "Y");
    H5Dclose(h5_ygrid);
    FREE(ygrids);

    // Extra bands - xgrid
    float *xgrids = (float *) MALLOC(sizeof(float)*pixel_count);
    for (ii=0; ii<nl; ii++) {
      for (kk=0; kk<ns; kk++)
	xgrids[ii*ns+kk] = mp->startX + kk*mp->perX;
      asfLineMeter(ii, nl);
    }
    asfPrintStatus("Storing band 'xgrid' ...\n");
    sprintf(dataset, "/data/xgrid");
    h5_xgrid = H5Dcreate(h5_file, dataset, H5T_NATIVE_FLOAT, h5_array,
			 H5P_DEFAULT, h5_plist, H5P_DEFAULT);
    H5Dwrite(h5_xgrid, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, xgrids);
    h5_att_str(h5_xgrid, h5_string, "units", "meters");
    h5_att_str(h5_xgrid, h5_string, "long_name",
	       "projection_grid_x_coordinates");
    h5_att_str(h5_xgrid, h5_string, "standard_name",
	       "projection_x_coordinates");
    h5_att_str(h5_xgrid, h5_string, "axis", "X");
    H5Dclose(h5_xgrid);
    FREE(xgrids);
  }

  H5Gclose(h5_datagroup);

  // Adding global attributes
  hid_t h5_global = H5Gopen(h5_file, "/", H5P_DEFAULT);
  h5_att_str(h5_global, h5_string, "institution", "Alaska Satellite Facility");
  sprintf(str_attr, "%s %s %s image", mg->sensor, mg->sensor_name, mg->mode);
  h5_att_str(h5_global, h5_string, "title", str_attr);
  if (mg->image_data_type == AMPLITUDE_IMAGE)
    strcpy(str_attr, "SAR backcatter image");
  h5_att_str(h5_global, h5_string, "source", str_attr);
  h5_att_str(h5_global, h5_string, "original_file", mg->basename);
  ymd_date ymd;
  hms_time hms;
  parse_date(mg->acquisition_date, &ymd, &hms);
  if (strcmp_case(mg->sensor, "RSAT-1") == 0)
    sprintf(str_attr, "Copyright Canadian Space Agency, %d", ymd.year);
  else if (strncmp_case(mg->sensor, "ERS", 3) == 0)
    sprintf(str_attr, "Copyright European Space Agency, %d", ymd.year);
  else if (strcmp_case(mg->sensor, "JERS-1") == 0 ||
	   strcmp_case(mg->sensor, "ALOS") == 0)
    sprintf(str_attr, "Copyright Japan Aerospace Exploration Agency , %d",
	    ymd.year);
  h5_att_str(h5_global, h5_string, "comment", str_attr);
  h5_att_str(h5_global, h5_string, "reference",
	     "Documentation available at: www.asf.alaska.edu");
  time_t t;
  struct tm *timeinfo;
  time(&t);
  timeinfo = gmtime(&t);
  sprintf(str_attr, "%s", asctime(timeinfo));
  chomp(str_attr);
  strcat(str_attr, ", UTC: H5 File created.");
  h5_att_str(h5_global, h5_string, "history", str_attr);
  H5Gclose(h5_global);

  // Metadata
  sprintf(group, "/metadata");
  h5_metagroup = H5Gcreate(h5_file, group, H5P_DEFAULT, H5P_DEFAULT,
			   H5P_DEFAULT);

  // Metadata - General block
  h5_value_str(h5_file, group, "general_name", mg->basename, "file_name", NULL);
  h5_value_str(h5_file, group, "general_sensor", mg->sensor,
	       "imaging satellite", NULL);
  h5_value_str(h5_file, group, "general_sensor_name", mg->sensor_name,
	       "imaging sensor", NULL);
  h5_value_str(h5_file, group, "general_mode", mg->mode, "imaging mode", NULL);
  h5_value_str(h5_file, group, "general_processor", mg->processor,
	       "name and version of processor", NULL);
  h5_value_str(h5_file, group, "general_data_type",
	       data_type2str(mg->data_type), "type of samples (e.g. REAL64)",
	       NULL);
  h5_value_str(h5_file, group, "general_image_data_type",
	       image_data_type2str(mg->image_data_type),
	       "image data type (e.g. AMPLITUDE_IMAGE)", NULL);
  h5_value_str(h5_file, group, "general_radiometry",
	       radiometry2str(mg->radiometry), "radiometry (e.g. SIGMA)", NULL);
  // FIXME: UDUNITS seconds since ...
  h5_value_str(h5_file, group, "general_acquisition_date", mg->acquisition_date,
	       "acquisition date of image", NULL);
  h5_value_int(h5_file, group, "general_orbit", mg->orbit,
	       "orbit number of image", NULL);
  if (mg->orbit_direction == 'A')
    strcpy(str_attr, "Ascending");
  else
    strcpy(str_attr, "Descending");
  h5_value_str(h5_file, group, "general_orbit_direction", str_attr,
	       "orbit direction", NULL);
  h5_value_int(h5_file, group, "general_frame", mg->frame,
	       "frame number of image", NULL);
  h5_value_int(h5_file, group, "general_band_count", mg->band_count,
	       "number of bands in image", NULL);
  h5_value_str(h5_file, group, "general_bands", mg->bands,
	       "bands of the sensor", NULL);
  h5_value_int(h5_file, group, "general_line_count", mg->line_count,
	       "number of lines in image", NULL);
  h5_value_int(h5_file, group, "general_sample_count", mg->sample_count,
	       "number of samples in image", NULL);
  h5_value_int(h5_file, group, "general_start_line", mg->start_line,
	       "first line relative to original image", NULL);
  h5_value_int(h5_file, group, "general_start_sample", mg->start_sample,
	       "first sample relative to original image", NULL);
  h5_value_double(h5_file, group, "general_x_pixel_size", mg->x_pixel_size,
		  "range pixel size", "m");
  h5_value_double(h5_file, group, "general_y_pixel_size", mg->y_pixel_size,
		  "azimuth pixel size", "m");
  h5_value_double(h5_file, group, "general_center_latitude",
		  mg->center_latitude, "approximate image center latitude",
		  "degrees_north");
  h5_value_double(h5_file, group, "general_center_longitude",
		  mg->center_longitude, "approximate image center longitude",
		  "degrees_east");
  h5_value_double(h5_file, group, "general_re_major", mg->re_major,
		  "major (equator) axis of earth", "m");
  h5_value_double(h5_file, group, "general_re_minor", mg->re_minor,
		  "minor (polar) axis of earth", "m");
  h5_value_double(h5_file, group, "general_bit_error_rate", mg->bit_error_rate,
		  "fraction of bits which are in error", NULL);
  h5_value_int(h5_file, group, "general_missing_lines", mg->missing_lines,
	       "number of missing lines in image", NULL);
  h5_value_float(h5_file, group, "general_no_data", mg->no_data,
		 "value indicating no data for a pixel", NULL);

  if (ms) {
    if (ms->image_type == 'S')
      strcpy(str_attr, "slant range");
    else if (ms->image_type == 'G')
      strcpy(str_attr, "ground range");
    else if (ms->image_type == 'P')
      strcpy(str_attr, "projected");
    else if (ms->image_type == 'R')
      strcpy(str_attr, "georeferenced");
    h5_value_str(h5_file, group, "sar_image_type", str_attr, "image type",
		 NULL);
    if (ms->look_direction == 'R')
      strcpy(str_attr, "right");
    else if (ms->look_direction == 'L')
      strcpy(str_attr, "left");
    h5_value_str(h5_file, group, "sar_look_direction", str_attr,
		 "SAR satellite look direction", NULL);
    h5_value_int(h5_file, group, "sar_look_count", ms->look_count,
		 "number of looks to take from SLC", NULL);
    h5_value_int(h5_file, group, "sar_multilook", ms->multilook,
		 "multilooking flag", NULL);
    h5_value_int(h5_file, group, "sar_deskewed", ms->deskewed,
		 "zero doppler deskew flag", NULL);
    h5_value_int(h5_file, group, "sar_original_line_count",
		 ms->original_line_count, "number of lines in original image",
		 NULL);
    h5_value_int(h5_file, group, "sar_original_sample_count",
		 ms->original_sample_count,
		 "number of samples in original image", NULL);
    h5_value_double(h5_file, group, "sar_line_increment",
		    ms->line_increment, "line increment for sampling", NULL);
    h5_value_double(h5_file, group, "sar_sample_increment",
		    ms->sample_increment, "sample increment for sampling",
		    NULL);
    h5_value_double(h5_file, group, "sar_range_time_per_pixel",
		    ms->range_time_per_pixel, "time per pixel in range", "s");
    h5_value_double(h5_file, group, "sar_azimuth_time_per_pixel",
		    ms->azimuth_time_per_pixel, "time per pixel in azimuth",
		    "s");
    h5_value_double(h5_file, group, "sar_slant_range_first_pixel",
		    ms->slant_range_first_pixel, "slant range to first pixel",
		    "m");
    h5_value_double(h5_file, group, "sar_slant_shift", ms->slant_shift,
		    "error correction factor in slant range", "m");
    h5_value_double(h5_file, group, "sar_time_shift", ms->time_shift,
		    "error correction factor in time", "s");
    h5_value_double(h5_file, group, "sar_wavelength", ms->wavelength,
		    "SAR carrier wavelength", "m");
    h5_value_double(h5_file, group, "sar_pulse_repetition_frequency",
		    ms->prf, "pulse repetition frequency", "Hz");
    h5_value_double(h5_file, group, "sar_earth_radius", ms->earth_radius,
		    "earth radius at image center", "m");
    h5_value_double(h5_file, group, "sar_satellite_height",
		    ms->satellite_height,
		    "satellite height from earth's center", "m");
    h5_value_double(h5_file, group, "sar_range_doppler_centroid",
		    ms->range_doppler_coefficients[0], "range doppler centroid",
		    "Hz");
    // FIXME: UDUNITS unit_description
    h5_value_double(h5_file, group, "sar_range_doppler_linear",
		    ms->range_doppler_coefficients[1],
		    "range doppler per range pixel", "Hz/pixel");
    // FIXME: UDUNITS unit description
    h5_value_double(h5_file, group, "sar_range_doppler_quadratic",
		    ms->range_doppler_coefficients[2],
		    "range doppler per range pixel square", "Hz/pixel^2");
    h5_value_double(h5_file, group, "sar_azimuth_doppler_centroid",
		    ms->azimuth_doppler_coefficients[0],
		    "azimuth doppler centroid", "Hz");
    h5_value_double(h5_file, group, "sar_azimuth_doppler_linear",
		    ms->azimuth_doppler_coefficients[1],
		    "azimuth doppler per azimuth pixel", "Hz/pixel");
    h5_value_double(h5_file, group, "sar_azimuth_doppler_quadratic",
		    ms->azimuth_doppler_coefficients[2],
		    "azimuth doppler per azimuth pixel square", "Hz/pixel^2");
  }

  if (mo) {
    int vector_count = mo->vector_count;
    h5_value_int(h5_file, group, "orbit_year", mo->year,
		 "year of image start", NULL);
    h5_value_int(h5_file, group, "orbit_day_of_year", mo->julDay,
		 "day of year at image start", NULL);
    h5_value_double(h5_file, group, "orbit_second_of_day", mo->second,
		    "second of day at image start", "s");
    h5_value_int(h5_file, group, "orbit_vector_count", vector_count,
		 "number of state vectors", NULL);
    for (ii=0; ii<vector_count; ii++) {
      sprintf(tmp, "orbit_vector[%d]_time", ii+1);
      h5_value_double(h5_file, group, tmp, mo->vecs[ii].time,
		      "time relative to image start", "s");
      sprintf(tmp, "orbit_vector[%d]_position_x", ii+1);
      h5_value_double(h5_file, group, tmp, mo->vecs[ii].vec.pos.x,
		      "x coordinate, earth-fixed", "m");
      sprintf(tmp, "orbit_vector[%d]_position_y", ii+1);
      h5_value_double(h5_file, group, tmp, mo->vecs[ii].vec.pos.y,
		      "y coordinate, earth-fixed", "m");
      sprintf(tmp, "orbit_vector[%d]_position_z", ii+1);
      h5_value_double(h5_file, group, tmp, mo->vecs[ii].vec.pos.z,
		      "z coordinate, earth-fixed", "m");
      sprintf(tmp, "orbit_vector[%d]_velocity_x", ii+1);
      h5_value_double(h5_file, group, tmp, mo->vecs[ii].vec.vel.x,
		      "x velocity, earth-fixed", "m/s");
      sprintf(tmp, "orbit_vector[%d]_velocity_y", ii+1);
      h5_value_double(h5_file, group, tmp, mo->vecs[ii].vec.vel.y,
		      "y velocity, earth-fixed", "m/s");
      sprintf(tmp, "orbit_vector[%d]_velocity_z", ii+1);
      h5_value_double(h5_file, group, tmp, mo->vecs[ii].vec.vel.z,
		      "z velocity, earth-fixed", "m/s");
    }
  }
  H5Gclose(h5_metagroup);
  H5Sclose(h5_string);
  H5Sclose(h5_array);

  // Write ASF metadata to XML file
  char *output_file =
    (char *) MALLOC(sizeof(char)*(strlen(output_file_name)+5));
  sprintf(output_file, "%s.xml", output_file_name);
  meta_write_xml(md, output_file);
  FREE(output_file);

  return h5;
}

void finalize_h5_file(h5_t *hdf)
{
  H5Fclose(hdf->file);

  // Clean up
  FREE(hdf->var);
  FREE(hdf);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


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

/*!
 *  \file transfun.c
 *  \brief calculate transfer functions and emission lines
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <mpi.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>

#include "brains.h"

/*
 * calculate (Fcon + ftrend)^Ag
 *
 */
void calculate_con_rm(const void *pm)
{
  int i, m;
  double fcon, A, Ag, ftrend, a0=0.0, tmp;
  double *pmodel = (double *)pm;

  A=exp(pmodel[idx_resp]);
  Ag=pmodel[idx_resp + 1];

  /* add different trend in continuum and emission */
  if(parset.flag_trend_diff > 0)
  {
    tmp = 0.0;
    for(m=1; m<num_params_difftrend+1; m++)
    {
      tmp += pmodel[idx_difftrend + m-1] * pow_Tcon_data[m-1];
    }
    a0 = -tmp;

    for(i=0; i<parset.n_con_recon; i++)
    {
      ftrend = a0;
      tmp = 1.0;
      for(m=1; m<num_params_difftrend+1; m++)
      {
        tmp *= (Tcon[i] - Tmed_data);
        ftrend += pmodel[idx_difftrend + m-1] * tmp;
      }

      fcon = Fcon[i] + ftrend;
      if(fcon > 0.0)
      {
        Fcon_rm[i] = A * pow(fcon, 1.0 + Ag);
      }
      else
      {
        Fcon_rm[i] = 0.0;
      }
    }
  }
  else
  {
    for(i=0; i<parset.n_con_recon; i++)
    {
      Fcon_rm[i] = A * pow(Fcon[i], 1.0+Ag);
    }
  }

  return;
}

/*!
 * This function calculate 1d line from a given BLR model.
 *
 * Note that continuum light curve has been obtained in advance
 *
 * Input:
 *   pm: model parameters
 *   Tl: epochs of line
 *   nl: number of points
 * Output:
 *   Fl: line fluxes
 *
 */
void calculate_line_from_blrmodel(const void *pm, double *Tl, double *Fl, int nl)
{
  int i, j;
  double fline, fcon_rm, tl, tc, tau, dTransTau;

  dTransTau = TransTau[1] - TransTau[0];

  for(i=0;i<nl;i++)
  {
    tl = Tl[i];
    fline = 0.0;
    for(j=0; j<parset.n_tau; j++)
    {
      tau = TransTau[j];
      tc = tl - tau;
      fcon_rm = gsl_interp_eval(gsl_linear, Tcon, Fcon_rm, tc, gsl_acc); /* interpolation */

      fline += Trans1D[j] * fcon_rm;     /*  line response */
    }
    Fl[i] = fline * dTransTau;
  }

  return;
}

/*!
 * This function caclulate 2d line from obtained transfer function.
 */
void calculate_line2d_from_blrmodel(const void *pm, const double *Tl, const double *transv, const double *trans2d,
                                              double *fl2d, int nl, int nv)
{
  int i, j, k;
  double tau, tl, tc, fcon_rm, fnarrow, dTransTau;
  double *pmodel = (double *)pm;

  dTransTau = TransTau[1] - TransTau[0];

  for(j=0;j<nl; j++)
  {
    for(i=0; i<nv; i++)
      fl2d[j*nv + i] = 0.0;

    tl = Tl[j];
    for(k=0; k<parset.n_tau; k++)
    {
      tau = TransTau[k];
      tc = tl - tau;
      fcon_rm = gsl_interp_eval(gsl_linear, Tcon, Fcon_rm, tc, gsl_acc);

      for(i=0; i<nv; i++)
      {
        fl2d[j*nv + i] += trans2d[k*nv+i] * fcon_rm;
      }
    }
    for(i=0; i<nv; i++)
    {
      fl2d[j*nv + i] *= dTransTau;
    }
  }

  /* add intrinsic narrow line */
  if(parset.flag_narrowline != 0)
  {
    double flux, width, shift;
    if(parset.flag_narrowline == 1)  /* fixed narrow line */
    {
      flux = parset.flux_narrowline;
      width = parset.width_narrowline;
      shift = parset.shift_narrowline;
    }
    else if(parset.flag_narrowline == 2) /* narrow line with Gaussian priors */
    {
      flux =  parset.flux_narrowline  + pmodel[num_params_blr_model] * parset.flux_narrowline_err;
      width = parset.width_narrowline + pmodel[num_params_blr_model+1] * parset.width_narrowline_err;
      shift = parset.shift_narrowline + pmodel[num_params_blr_model+2] * parset.shift_narrowline_err;
    }
    else  /* narrow line with logrithmic prior of flux */
    {
      flux =  exp(pmodel[num_params_blr_model]);
      width = parset.width_narrowline + pmodel[num_params_blr_model+1] * parset.width_narrowline_err;
      shift = parset.shift_narrowline + pmodel[num_params_blr_model+2] * parset.shift_narrowline_err;
    }

    width = fmax(1.0e-10, width); /* make sure thant width is not zero */

    for(i=0; i<nv; i++)
    {
      fnarrow = flux * exp( -0.5 * pow( (transv[i] - shift)/(width), 2.0) );
      for(j = 0; j<nl; j++)
      {
        fl2d[j*nv + i] += fnarrow;
      }
    }
  }

  /* smooth the line profile */
  line_gaussian_smooth_2D_FFT(transv, fl2d, nl, nv, pm);
}

/*!
 * This function calculates 1d transfer function.
 */
void transfun_1d_cal_cloud(const void *pm, int flag_save)
{
  /* generate cloud sample and calculate the corresponding time lags and weights */
  gen_cloud_sample(pm, 1, flag_save);
  transfun_1d_cal_with_sample();
  return;
}

/*!
 * This function calculates 1d transfer function.
 */
void transfun_1d_cal_with_sample()
{
  int i, idt;
  double tau_min, tau_max, dTransTau;
  double Anorm, dis;

  tau_min = clouds_tau[0];
  tau_max = clouds_tau[0];
  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    if(tau_min > clouds_tau[i])
      tau_min = clouds_tau[i];
    if(tau_max < clouds_tau[i])
      tau_max = clouds_tau[i];
  }

  dTransTau = (tau_max - tau_min)/(parset.n_tau - 1);
  for(i=0; i<parset.n_tau; i++)
  {
    TransTau[i] = tau_min + dTransTau * i;
    Trans1D[i] = 0.0;
  }

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    dis = clouds_tau[i];
    idt = (dis - tau_min)/dTransTau;
    //Trans1D[idt] += pow(1.0/r, 2.0*(1 + gam)) * weight;
    Trans1D[idt] += clouds_weight[i];
  }

  /* normalize transfer function */
  Anorm = 0.0;
  for(i=0;i<parset.n_tau;i++)
  {
    Anorm += Trans1D[i];
  }
  Anorm *= dTransTau;

  Anorm += EPS;

  for(i=0; i<parset.n_tau; i++)
  {
    Trans1D[i] /= Anorm;
  }
  return;
}

/*!
 * This function calculates 2d transfer function.
 */
void transfun_2d_cal_cloud(const void *pm, double *transv, double *trans2d, int n_vel, int flag_save)
{
  /* generate cloud sample and calculate the corresponding time lags, LOS velocity, and weights */
  gen_cloud_sample(pm, 2, flag_save);
  transfun_2d_cal_with_sample(transv, trans2d, n_vel);
  return;
}

/*!
 * This function calculates 2d transfer function.
 */
void transfun_2d_cal_with_sample(double *transv, double *trans2d, int n_vel)
{
  int i, j, idt, idV;
  double tau_min, tau_max, dTransTau;
  double Anorm, dis, V, dV;

  tau_min = clouds_tau[0];
  tau_max = clouds_tau[0];
  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    if(tau_min > clouds_tau[i])
      tau_min = clouds_tau[i];
    if(tau_max < clouds_tau[i])
      tau_max = clouds_tau[i];
  }

  dTransTau = (tau_max - tau_min)/(parset.n_tau - 1);
  for(i=0; i<parset.n_tau; i++)
  {
    TransTau[i] = tau_min + dTransTau * i;
  }

  dV =(transv[1] - transv[0]); /* velocity grid width */
  for(i=0; i<parset.n_tau; i++)
    for(j=0;j<n_vel;j++)
      trans2d[i*n_vel+j]=0.0;   /* cleanup of transfer function */

  for(i=0; i<parset.n_cloud_per_task; i++)
  {
    dis = clouds_tau[i];
    idt = (dis - tau_min)/dTransTau;

    for(j=0; j<parset.n_vel_per_cloud; j++)
    {
      V = clouds_vel[i*parset.n_vel_per_cloud + j];
      if(V<transv[0] || V >= transv[n_vel-1]+dV)
        continue;
      idV = (V - transv[0])/dV;
      //trans2d[idt*n_vel + idV] += pow(1.0/r, 2.0*(1 + gam)) * weight;
      trans2d[idt*n_vel + idV] += clouds_weight[i];
    }
  }

  /* normalize transfer function */
  Anorm = 0.0;
  for(i=0; i<parset.n_tau; i++)
    for(j=0; j<n_vel; j++)
    {
      Anorm += trans2d[i*n_vel+j];
    }
  Anorm *= (dV * dTransTau);

  Anorm += EPS;
  for(i=0; i<parset.n_tau; i++)
  {
    for(j=0; j<n_vel; j++)
    {
      trans2d[i*n_vel+j] /= Anorm;
    }
  }
  return;
}

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
/*
** compute correlation matrix
**
** G.Lohmann, Feb 2011
*/

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_histogram.h>

#include "viaio/Vlib.h"
#include "viaio/VImage.h"
#include "viaio/mu.h"

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

#define SQR(x) ((x) * (x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

extern void VNormalize(float *data,int nt,VBoolean stddev);
extern long GetAddr(VImage mapimage,int bi,int ri,int ci,int m,int k,int l);

int VNumNeigbours(size_t id,VImage map,VImage mapimage,int adjdef)
{
  int nslices = VImageNBands(mapimage);
  int nrows = VImageNRows(mapimage);
  int ncols = VImageNColumns(mapimage);

  int bi = VPixel(map,0,0,id,VShort);
  int ri = VPixel(map,0,1,id,VShort);
  int ci = VPixel(map,0,2,id,VShort);
  if (bi < 1 || bi >= nslices-2) return 0;
  if (ri < 1 || ri >= nrows-2) return 0;
  if (ci < 1 || ci >= ncols-2) return 0;

  int wn=1,rad2=0;
  if (adjdef == 3) {
    wn = 2;
    rad2 = 2*2;
  }

  int n=0;
  int k,l,m;
  for (m=-wn; m<=wn; m++) {
    for (k=-wn; k<=wn; k++) {
      for (l=-wn; l<=wn; l++) {
      	long jj = GetAddr(mapimage,bi,ri,ci,m,k,l);
      	if (jj < 0) continue;  /* outside of brain mask */

      	if (adjdef == 0) {     /* 6 adjacency */
      	  if (ABS(m)+ABS(k)+ABS(l) > 1) continue;
      	}
      	if (adjdef == 1) {     /* 18-adjacency */
      	  if (ABS(m) > 0 && ABS(k) > 0 && ABS(l) > 0) continue;
      	}
	if (adjdef == 3) {     /* sphere */
	  if (m*m + k*k + l*l > rad2) continue;
	}
      	n++;
      }
    }
  }
  return n;
}

/* convert to ranks */
void GetRank(float *data,gsl_vector *v,gsl_permutation *perm,gsl_permutation *rank)
{
  size_t i;
  size_t n = v->size;
  for (i=0; i<n; i++) gsl_vector_set(v,i,(double)data[i]);
  gsl_sort_vector_index (perm, v);
  gsl_permutation_inverse (rank, perm);
  for (i=0; i<n; i++) data[i] = (float)rank->data[i];
}

double Spearman(const float *data1,const float *data2,int n)
{
  int i;
  double nx = (double)n;
  double kx = nx*(nx*nx-1.0);
  double sxy=0.0;
  for (i=0; i<n; i++) {
    const double u = (double)data1[i];
    const double v = (double)data2[i];
    const double d = (u-v);
    sxy += d*d;
  }
  double rho = 1.0 - 6.0*sxy/kx;
  return rho;
}

double Correlation(const float *data1,const float *data2,int n)
{
  int i;
  double corr=0;
  for (i=0; i<n; i++) {
    const double u = (double)data1[i];
    const double v = (double)data2[i];
    corr += u*v;
  }
  corr /= (double)n;
  return corr;
}

float EdgeCorr(const float *data1,const float *data2,int n,int metric)
{
  float corr=0.0;
  if (metric == 0) corr = Correlation(data1,data2,n);
  if (metric == 1) corr = Spearman(data1,data2,n);
  if (corr > 0) return corr;
  else return 0.0;
}

void GetSNR(gsl_matrix_float **X1,gsl_matrix_float **X2,int *table,int n,gsl_matrix_float *SNR,int metric)
{
  long j,nvox=X1[0]->size1;
  long k,s,nt=X1[0]->size2;

  if (n < 3) VError(" n: %d",n);
  double nx = (double)n;
  double *sum1 = (double *)VCalloc(nt,sizeof(double));
  double *sum2 = (double *)VCalloc(nt,sizeof(double));

  gsl_vector *vec = NULL;
  gsl_permutation *perm = NULL;
  gsl_permutation *rank = NULL;
  if (metric == 1) {  /* only needed for spearman correlation */
    vec = gsl_vector_calloc(nt);
    perm = gsl_permutation_alloc(nt);
    rank = gsl_permutation_alloc(nt);
  }

  double ave=0,var=0,snr=0;

  for (j=0; j<nvox; j++) {
    memset(sum1,0,nt*sizeof(double));
    memset(sum2,0,nt*sizeof(double));
    const float *pp=NULL;
    for (s=0; s<n; s++) {
      if (table[s] == 0)
	pp = gsl_matrix_float_const_ptr(X1[s],j,0);
      else
	pp = gsl_matrix_float_const_ptr(X2[s],j,0);
      for (k=0; k<nt; k++) {
	const double x = (double)(*pp++);
	sum1[k] += x;
	sum2[k] += x*x;
      }
    }
    for (k=0; k<nt; k++) {
      ave = (sum1[k]/nx);
      var = (sum2[k] - nx*ave*ave) / (nx - 1.0);
      snr = 0;
      if (var > 0) snr = ave / sqrt(var);
      gsl_matrix_float_set(SNR,j,k,snr);
    }

    if (metric == 0) {  /* pearson correlation */
      float *qq = gsl_matrix_float_ptr(SNR,j,0);
      VNormalize(qq,nt,TRUE);
    }
    if (metric == 1) {  /* spearman ranks      */
      float *qq = gsl_matrix_float_ptr(SNR,j,0);
      GetRank(qq,vec,perm,rank);
    }
  }
  VFree(sum1);
  VFree(sum2);
}

float ZMatrix(gsl_histogram *histogram,gsl_matrix_float *SNR1,gsl_matrix_float *SNR2,int n1,int n2,
	      VImage roi,VImage map,VImage mapimage,int adjdef,float elength,float quantile,int step,int metric)
{
  size_t i;
  size_t nvox=SNR1->size1;
  int nt=SNR1->size2;
  size_t progress=0;
  int rad2 = (int)(elength*elength);
  double tiny=1.0e-8;
  gsl_set_error_handler_off ();

  fprintf(stderr," Computing matrix...\n");

  gsl_histogram_reset(histogram);
  size_t nbins = gsl_histogram_bins (histogram);
  double hmax = gsl_histogram_max (histogram);
  double hmin = gsl_histogram_min (histogram);
  double zmin = 99999.0;
  double zmax = -99999.0;
  int minadj = 1;

#pragma omp parallel for shared(progress,histogram) schedule(dynamic) firstprivate(SNR1,SNR2)
  for (i=0; i<nvox; i+=step) {
    if (i%1000 == 0) fprintf(stderr," %ld000\r",(long)progress++);
    int bi = (int)VPixel(map,0,0,i,VShort);
    int ri = (int)VPixel(map,0,1,i,VShort);
    int ci = (int)VPixel(map,0,2,i,VShort);
    int nadjx = VNumNeigbours(i,map,mapimage,adjdef);
    if (nadjx < minadj) continue;

    int roiflagi = 0;
    if (roi != NULL) {
      if (VGetPixel(roi,bi,ri,ci) > 0.5) roiflagi = 1;
    }

    gsl_histogram *tmphist = gsl_histogram_alloc (nbins);
    gsl_histogram_set_ranges_uniform (tmphist,hmin,hmax);
    gsl_histogram_reset(tmphist);

    const float *datax1 = gsl_matrix_float_const_ptr(SNR1,i,0);
    const float *datax2 = gsl_matrix_float_const_ptr(SNR2,i,0);

    size_t j=0;
    for (j=0; j<i; j+=step) {
      int bj = (int)VPixel(map,0,0,j,VShort);
      int rj = (int)VPixel(map,0,1,j,VShort);
      int cj = (int)VPixel(map,0,2,j,VShort);
      int d = SQR(bi-bj) + SQR(ri-rj) + SQR(ci-cj);
      if (d < rad2) continue;
      int nadjy = VNumNeigbours(j,map,mapimage,adjdef);
      if (nadjy < minadj) continue;

      int roiflagj = 0;
      if (roi != NULL) {
	if (VGetPixel(roi,bj,rj,cj) > 0.5) roiflagj = 1;
	if (roiflagi + roiflagj != 1) continue;
      }
      const float *datay1 = gsl_matrix_float_const_ptr(SNR1,j,0);
      const float *datay2 = gsl_matrix_float_const_ptr(SNR2,j,0);

      /* edge z-value */
      double z1 = EdgeCorr(datax1,datay1,nt,metric);
      double z2 = EdgeCorr(datax2,datay2,nt,metric);
      double z = (z1-z2);
      if (z < zmin) zmin = z;
      if (z > zmax) zmax = z;

      if (z < hmin) z = hmin;
      if (z > hmax-tiny) z = hmax-tiny;
      gsl_histogram_increment (tmphist,z);
    }

#pragma omp critical
    {
      gsl_histogram_add (histogram,tmphist);
    }
    gsl_histogram_free (tmphist);
  }

  /* get quantile cutoff */
  gsl_histogram_pdf *pdf = gsl_histogram_pdf_alloc(nbins);
  gsl_histogram_pdf_init(pdf,histogram);
  double lower=0,upper=0;

  size_t i0=0;
  for (i=nbins-1; i>=0; i--) {
    gsl_histogram_get_range (histogram,i,&lower,&upper);
    if (pdf->sum[i] < quantile) {
      if (gsl_histogram_get_range (histogram,i,&lower,&upper) == GSL_EDOM) VError(" err hist");
      i0 = i;
      break;
    }
  }
  double sum=0;
  for (i=nbins-1; i>i0; i--) {
    sum += gsl_histogram_get(histogram,i);
  }
  return (float)upper;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


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

#pragma once

#ifndef loopblinn_875a5488_30b0_4449_b97a_e55c8d8a0db1_h
#define loopblinn_875a5488_30b0_4449_b97a_e55c8d8a0db1_h

#include <gslib/std.h>
#include <gslib/rtree.h>
#include <ariel/painterpath.h>
#include <ariel/delaunay.h>

__ariel_begin__

class lb_joint;
class lb_end_joint;
class lb_control_joint;
class lb_line;
class lb_polygon;
class lb_shrink_line;
class lb_shrink;

typedef vector<lb_joint*> lb_joint_list;
typedef vector<lb_line*> lb_line_list;
typedef vector<lb_polygon*> lb_polygon_list;
typedef vector<lb_shrink_line*> lb_shrink_lines;
typedef stack<lb_polygon*> lb_polygon_stack;

enum lb_joint_type
{
    lbt_end_joint,
    lbt_control_joint,
};

/*
 * The reason why the original point and the NDC(normalized device coordinate) point should
 * be separated was that we need the original point for the convenience of the error handling
 * works, and we need the NDC point to do the loopblinn calculations so that the float
 * point won't overflow.
 */
class __gs_novtable lb_joint abstract
{
protected:
    lb_line*            _prev;
    lb_line*            _next;
    vec2                _point;
    vec2                _ndcpoint;
    void*               _binding;

public:
    lb_joint()
    {
        _prev = _next = nullptr;
        _binding = nullptr;
    }
    virtual ~lb_joint() {}
    virtual lb_joint_type get_type() const = 0;
    virtual const vec2& get_point() const { return _point; }
    virtual const vec2& get_ndc_point() const { return _ndcpoint; }

public:
    void set_point(const vec2& p) { _point = p; }
    void set_ndc_point(const vec2& p) { _ndcpoint = p; }
    void set_prev_line(lb_line* p) { _prev = p; }
    void set_next_line(lb_line* p) { _next = p; }
    lb_line* get_prev_line() const { return _prev; }
    lb_line* get_next_line() const { return _next; }
    lb_joint* get_prev_joint() const;
    lb_joint* get_next_joint() const;
    bool is_adjacent_joint(const lb_joint* p) const { return p == get_prev_joint() || p == get_next_joint(); }
    void set_binding(void* p) { _binding = p; }
    void* get_binding() const { return _binding; }
};

class lb_end_joint:
    public lb_joint
{
protected:
    vec3                _klm[2];

public:
    lb_end_joint() {}
    lb_joint_type get_type() const override { return lbt_end_joint; }
    bool prev_is_curve() const;
    bool next_is_curve() const;
    void set_klm(int i, const vec3& p) { _klm[i] = p; }
    const vec3& get_klm(int i) const { return _klm[i]; }
};

class lb_control_joint:
    public lb_joint
{
protected:
    vec3                _klm;

public:
    lb_control_joint() {}
    lb_joint_type get_type() const override { return lbt_control_joint; }
    void set_klm(const vec3& p) { _klm = p; }
    const vec3& get_klm() const { return _klm; }
};

class lb_line
{
protected:
    lb_joint*           _joint[2];
    bool                _opened;

public:
    lb_line()
    {
        _joint[0] = _joint[1] = nullptr;
        _opened = false;
    }
    void set_opened(bool b) { _opened = b; }
    bool is_opened() const { return _opened; }
    void set_prev_joint(lb_joint* p) { _joint[0] = p; }
    void set_next_joint(lb_joint* p) { _joint[1] = p; }
    lb_joint* get_prev_joint() const { return _joint[0]; }
    lb_joint* get_next_joint() const { return _joint[1]; }
    lb_line* get_prev_line() const;
    lb_line* get_next_line() const;
    const vec2& get_prev_point() const { return _joint[0]->get_point(); }
    const vec2& get_next_point() const { return _joint[1]->get_point(); }
    static lb_line* get_line_between(lb_joint* j1, lb_joint* j2);
};

enum lb_span_type
{
    lst_linear,
    lst_quad,
    lst_cubic,
};

class __gs_novtable lb_span abstract
{
protected:
    rectf               _rc;

public:
    virtual ~lb_span() {}
    virtual lb_span_type get_type() const = 0;
    virtual bool can_split() const = 0;
    virtual bool is_overlapped(const lb_span* span) const = 0;

public:
    const rectf& get_rect() const { return _rc; }
    float get_area() const { return _rc.width() * _rc.height(); }
};

typedef vector<lb_span*> lb_span_list;
typedef rtree_entity<lb_span*> lb_rtree_entity;
typedef rtree_node<lb_rtree_entity> lb_rtree_node;
typedef _tree_allocator<lb_rtree_node> lb_rtree_alloc;
typedef tree<lb_rtree_entity, lb_rtree_node, lb_rtree_alloc> lb_tree;
typedef rtree<lb_rtree_entity, quadratic_split_alg<8, 3, lb_tree>, lb_rtree_node, lb_rtree_alloc> lb_rtree;
typedef delaunay_triangulation lb_triangulator;

/*
 * The polygon should be decomposed to the form like boundary - holes,
 * this procedure could also be called flattening.
 * After we retrieve the boundary we need to create a shrink for the boundary
 * for including test, if necessary, so that we could decide if a sub path
 * (case : cw - ccw - {cw}?) was a new path that ends the previous boundary.
 */
class lb_polygon
{
protected:
    lb_line*            _boundary;
    lb_line_list        _holes;
    lb_rtree            _mytree;
    lb_triangulator     _cdt;
    dt_input_joints     _dtjoints;

public:
    lb_polygon() { _boundary = nullptr; }
    void set_boundary(lb_line* p) { _boundary = p; }
    void add_hole(lb_line* p) { _holes.push_back(p); }
    lb_line* get_boundary() const { return _boundary; }
    const lb_line_list& get_holes() const { return _holes; }
    bool is_inside(const vec2& p) const;
    lb_rtree& get_rtree() { return _mytree; }
    void convert_to_ndc(const mat3& m);
    void create_dt_joints();
    void pack_constraints();
    void build_cdt();
    lb_triangulator& get_cdt_result() { return _cdt; }
    void tracing() const;
    void trace_boundary() const;
    void trace_last_hole() const;
    void trace_holes() const;
    void trace_rtree() const { _mytree.tracing(); }
};

/*
 * Notice that the path put into the loopblinn categorizer MUST be a simple polygon,
 * which means you should run a xor clip process before if you can't tell whether a
 * path was simple or complex.
 * Another point was that the path MUST be Winding rule.
 * You can also convert a path of OddEven rule to Winding by clipping.
 */
class loop_blinn_processor
{
public:
    loop_blinn_processor(float w, float h) { _width = w, _height = h; }
    ~loop_blinn_processor();
    void proceed(const painter_path& path);
    lb_polygon_list& get_polygons() { return _polygons; }
    lb_joint_list& get_joints() { return _joint_holdings; }
    lb_line_list& get_lines() { return _line_holdings; }
    void trace_polygons() const;
    void trace_rtree() const;

protected:
    float               _width;
    float               _height;
    lb_line_list        _line_holdings;
    lb_joint_list       _joint_holdings;
    lb_polygon_list     _polygons;
    lb_span_list        _span_holdings;

protected:
    template<class _joint>
    lb_joint* create_joint(const vec2& p);
    lb_line* create_line();
    lb_polygon* create_polygon();
    void flattening(const painter_path& path);
    int flattening(const painter_path& path, int start, lb_polygon* parent, lb_polygon_stack& st);
    int create_patch(lb_line*& line, const painter_path& path, int start);
    lb_joint* create_segment(lb_joint* prev, const painter_path& path, int i);
    void check_boundary(lb_polygon* poly);
    void check_holes(lb_polygon* poly);
    void check_span(lb_polygon* poly, lb_control_joint* joint);
    void check_rtree(lb_polygon* poly);
    void check_rtree_span(lb_polygon* poly, lb_span* span);
    void split_quadratic(lb_line* line1, lb_line* line2, lb_joint* sp[5]);
    void split_cubic(lb_line* line1, lb_line* line2, lb_line* line3, lb_joint* sp[7], float t);
    int try_split_cubic(lb_line* line1, lb_line* line2, lb_line* line3, lb_joint* sp[7], float t);
    lb_span* split_rtree_span(lb_span* span);
    void split_span_recursively(lb_polygon* poly, lb_span* span);
    void calc_klm_coords();
    void calc_klm_coords(lb_polygon* poly);
    void calc_klm_coords(lb_polygon* poly, lb_line* start);
    lb_line* calc_klm_span(lb_polygon* poly, lb_line* line);
};

extern lb_line* lb_get_span(lb_line_list& span, lb_line* start);
extern void lb_get_current_span(lb_line_list& span, lb_control_joint* joint);

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/*
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "aXe_errors.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_spline.h>
#include "fringe_conf.h"
#include "fringe_model.h"

#define AXE_IMAGE_PATH  "AXE_IMAGE_PATH"
#define AXE_OUTPUT_PATH "AXE_OUTPUT_PATH"
#define AXE_CONFIG_PATH "AXE_CONFIG_PATH"

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

  char        fconf_file[MAXCHAR];
  char        fconf_file_path[MAXCHAR];

  char        fimage_file[MAXCHAR];
  char        fimage_file_path[MAXCHAR];

  fringe_conf  *fconf;
  //interpolator *filter_through;
  gsl_matrix   *fringe_image;

  if ((argc < 3) || (opt = get_online_option("help", argc, argv))) {
    fprintf(stdout,
	    "aXe_FRIGEN Version %s:\n"
	    "\n", RELEASE);
    exit(1);
  }
  fprintf(stdout, "aXe_FRIGEN: Starting...\n");

  /* Get the name of the fringe configuration file*/
  strcpy(fconf_file, argv[1]);
  build_path(AXE_CONFIG_PATH, fconf_file, fconf_file_path);

  /* Get the configuration file name */
  strcpy(fimage_file, argv[2]);
  build_path(AXE_OUTPUT_PATH, fimage_file, fimage_file_path);

  /* report on the input and output that will be used */
  fprintf(stdout,
	  "aXe_FRIGEN: Input fringe configuration file:   %s\n",
	  fconf_file_path);
  fprintf(stdout,
	  "aXe_FRIGEN: Output fringe image name:          %s\n",
	  fimage_file_path);

  // load the fringe configuration file
  fconf = load_fringe_conf(fconf_file_path);

  // check whether all necessary information
  // is in place, provide defaults
  check_fringe_conf(fconf);

  // check whether a filter throughput file is set
  if (!fconf->filter_through)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "aXe_FRIGEN: No filter throughput table\n"
		 "is set in the fringe configuration file!\n");

  // comute the fringe amplitude
  fringe_image = compute_fringe_amplitude(fconf);

  // save the matrix to a fits-image
  gsl_to_FITSimage (fringe_image, fimage_file_path, 1, NULL);

  // release the space in the fringe image
  gsl_matrix_free(fringe_image);

  // release the allocated memory
  free_fringe_conf(fconf);

  // exit the program
  fprintf(stdout, "aXe_FRIGEN: Done...\n");
  exit(0);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/* popGenTools.h
/
/
/ Andrew Kern
/
/
*/

#define MAXSNPS 1000000
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

struct snp{ /*SNP struct definition- meant to hold position, derived freq, sample size */
  int pos, i, n, ascSize;
};

struct my_F_params{
  int i, n, ascSize;
};

struct my_snpProb_params{
  int i, n, ascSize;
  double denom;
};

struct my_f_params{
  int i,n, ascSize;
  double beta;
};

struct my_lik_params{
  int snpNumber, maxSampleSize, ascSize;
  gsl_vector *weights, *sampleSizeVector, *betas;
  struct snp *data;
  gsl_matrix *sfsBools;

};
struct my_likCI_params{
  int snpNumber, maxSampleSize;
  gsl_vector *weights, *sampleSizeVector;
  struct snp *data;
  gsl_matrix *sfsBools;
  double lMax, logUnits, beta_hat;
};

struct my_jointLik_params{
  int snpNumber, maxSampleSize, ascSize, nstates;
  gsl_vector *weights, *sampleSizeVector, *betas, *startingPoint;
  struct snp *data;
  gsl_matrix *sfsBools, *posts;

};

double my_f(double q, void * params);
double my_F(double beta, void * params);
double my_F2(double beta, void * p);
double my_lik(double beta, void * params);
double sfsLikBetaVector(gsl_vector *beta, void * p);
double weightedLik(double beta, void * params);
double weightedLikLook(double beta, void * params);
double weightedLikLookCI(double beta, void * params);
double likWrap(double beta);
double snpProb(double beta, void * params);
double snpProbDenomLookup(double beta, void * params);
double snpProbDenom(double beta, int sampleSize);
gsl_vector *makeSnpProbDenomVector(double beta, int maxSampleSize, gsl_vector *sampleSizeVec);
double ml_est(double * lik_beta_hat);
double weighted_ml_est(double * lik_beta_hat, void * p);
double weighted_ml_est_lookup(double * lik_beta_hat, void * p);
double weighted_ml_est_lookup_errReport(double * lik_beta_hat, void * p);
double weighted_ml_CILower_lookup(double * lik_beta_hat, void * p);
double weighted_ml_CIUpper_lookup(double * lik_beta_hat, void * p);
double wlikWrap(double beta, void * p);
double  wlikWrapLook(double beta, void * p);
double  wlikCIWrapLook(double beta, void * p);
gsl_vector *sampleSizeVector(struct snp data[], int snpNumber, int maxSampleSize);
gsl_matrix *summarizeSFS(int maxSampleSize, struct snp data[], int snpNumber);
gsl_matrix *summarizeSFSBool(int maxSampleSize, struct snp data[], int snpNumber);
gsl_matrix *snpProbMatrix(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools);
gsl_matrix *snpProbMatrixNotLog(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools);
gsl_matrix *snpProbMatrixNotLogFull(double beta,  int maxSampleSize,gsl_vector *sampleSizeVector);
gsl_vector *snpProbVectorNotLog(double beta,  int sampleSize);
int simulateSiteFreq(int sampleSize, double alpha, void *r);
double simpleAscertain(double sampleFreq, void *p);
double probAscertainmentGivenModel(double beta, void *p);
gsl_matrix *snpAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools, int ascSize);
gsl_matrix *estimateAscSFS(int maxSampleSize, gsl_vector *sampleSizeVector,  gsl_matrix *sfsBools, int ascSize, gsl_matrix *sfsSummary);
double weightedLikLookAsc(double beta, void * p);
double  wlikWrapLookAsc(double beta, void * p);
double weighted_ml_est_lookup_asc(double * lik_beta_hat, void * p);
double probAscertainmentGivenModelLookup(double beta, int sampSize, gsl_matrix *snpProbs, gsl_matrix *ascProbs);
double probAscertainmentGivenModelHemiLookup(double beta, int sampSize, gsl_vector *snpProbs, gsl_matrix *ascProbs);
double outgroupAscertain(double sampleFreq, void * p);
gsl_matrix *snpOutgroupAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, int ascSize);
double weightedLikLookOutgroupAsc(double beta, void * p);
double  wlikWrapLookOutgroupAsc(double beta, void * p);
double weighted_ml_est_lookup_outgroup_asc(double * lik_beta_hat, void * p);
double probOutgroupAscertainmentGivenModel(double beta, void *p);
double sfsLikBetaVectorOutgroupAsc(gsl_vector *betas, void * p);

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


text:
#include "stdlib.h"
#include "stdio.h"
#include "/home/lillian/work/install_fpdebug/valgrind-3.7.0/fpdebug/fpdebug.h"
#include <gsl/gsl_sf.h>
int main(int argc, const char * argv[]) {
unsigned long int hexdouble;
int a;
a = atoi(argv[1]);
int b;
b = atoi(argv[2]);
int c;
c = atoi(argv[3]);
int d;
d = atoi(argv[4]);
int e;
e = atoi(argv[5]);
int f;
f = atoi(argv[6]);
double result = gsl_sf_coupling_6j(a, b, c, d, e, f);
//printf("%.15f\n", result);
VALGRIND_PRINT_ERROR("result", &result);
return 0;
}

max_stars_repo_name: floatfeather/FpGenetic

max_issues_repo_name: floatfeather/FpGenetic

max_forks_repo_name: floatfeather/FpGenetic


text:
#include <string.h>
#include <sys/types.h>
#include <stdint.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_wavelet.h>
#include <fftw3.h>
#include "common.h"
#include "utils.h"
#include "filters.h"

filters_t *filters_init(ANALYSIS_WAVEFORM_BASE_TYPE *inWav, size_t n)
{
    filters_t *fHdl;
    fHdl = (filters_t*)malloc(sizeof(filters_t));

    fHdl->wavLen = n;
    fHdl->respLen = 0;
    fHdl->malloced = 0;
    fHdl->fftUsed = 0;
    fHdl->fftwNThreads = FFTW_NTHREADS_DEFAULT;
    fHdl->fftwFlags = FFTW_FLAGS_DEFAULT;

    if(inWav == NULL) {
        fHdl->inWav = (ANALYSIS_WAVEFORM_BASE_TYPE*)
            calloc(fHdl->wavLen, sizeof(ANALYSIS_WAVEFORM_BASE_TYPE));
        fHdl->malloced = 1;
    } else {
        fHdl->inWav = inWav;
    }

    fHdl->outWav = (ANALYSIS_WAVEFORM_BASE_TYPE*)
        calloc(fHdl->wavLen, sizeof(ANALYSIS_WAVEFORM_BASE_TYPE));
    fHdl->waveletWav = (WAVELET_BASE_TYPE*)
        calloc(fHdl->wavLen, sizeof(WAVELET_BASE_TYPE));

    fHdl->gslDWT = gsl_wavelet_alloc(gsl_wavelet_daubechies_centered, 10);
    fHdl->gslDWTWork = gsl_wavelet_workspace_alloc(fHdl->wavLen);

    return fHdl;
}

/* for convolution */
filters_t *filters_init_for_convolution(ANALYSIS_WAVEFORM_BASE_TYPE *inWav, size_t n, size_t np)
{
    filters_t *fHdl;

    if((np > 0) && (np % 2 == 0)) {
        error_printf("%s(): np = %zd is not odd!\n", __FUNCTION__, np);
        return NULL;
    }

    fHdl = filters_init(inWav, n);
    fHdl->respLen = np;
    fHdl->fftUsed = 1;
    fHdl->respWav = (ANALYSIS_WAVEFORM_BASE_TYPE*)
        calloc(fHdl->respLen, sizeof(ANALYSIS_WAVEFORM_BASE_TYPE));
    if(np > 0) {
        fHdl->fftLen = (fHdl->wavLen + fHdl->respLen+1); /* zero padding */
        if(fHdl->fftLen % 2) fHdl->fftLen++; /* ensure fHdl->fftLen is even */
    } else {
        fHdl->fftLen = fHdl->wavLen; /* for spectrum calculation */
    }

    if(fHdl->fftwNThreads > 0) {
        if(FFTW(init_threads)() == 0) {
            error_printf("fftw_init_threads error!\n");
        }
        FFTW(plan_with_nthreads)(fHdl->fftwNThreads);
    }

    fHdl->fftwWork = (FFT_BASE_TYPE*) FFTW(malloc)(sizeof(FFT_BASE_TYPE) * fHdl->fftLen);
    fHdl->fftwWork1 = (FFT_BASE_TYPE*) FFTW(malloc)(sizeof(FFT_BASE_TYPE) * fHdl->fftLen);
    fHdl->fftwWin = (FFT_BASE_TYPE*) FFTW(malloc)(sizeof(FFT_BASE_TYPE) * fHdl->fftLen);
    filters_hanning_window(fHdl); /* Hanning window as the default */
    fHdl->dt = 1.0;

    fHdl->fftwPlan = FFTW(plan_r2r_1d)(fHdl->fftLen, fHdl->fftwWork, fHdl->fftwWork,
                                       FFTW_R2HC, fHdl->fftwFlags);
    fHdl->fftwPlan1 = FFTW(plan_r2r_1d)(fHdl->fftLen, fHdl->fftwWork1, fHdl->fftwWork1,
                                        FFTW_R2HC, fHdl->fftwFlags);
    fHdl->fftwPlan2 = FFTW(plan_r2r_1d)(fHdl->fftLen, fHdl->fftwWork, fHdl->fftwWork,
                                        FFTW_HC2R, fHdl->fftwFlags);

    return fHdl;
}

int filters_close(filters_t *fHdl)
{
    if(fHdl->malloced) {
        if(fHdl->inWav)
            free(fHdl->inWav);
    }
    if(fHdl->outWav)
        free(fHdl->outWav);
    if(fHdl->fftUsed) {
        if(fHdl->respWav)
            free(fHdl->respWav);
        FFTW(destroy_plan)(fHdl->fftwPlan);
        FFTW(destroy_plan)(fHdl->fftwPlan1);
        FFTW(destroy_plan)(fHdl->fftwPlan2);
        FFTW(free)(fHdl->fftwWork);
        FFTW(free)(fHdl->fftwWork1);
        FFTW(free)(fHdl->fftwWin);
        if(fHdl->fftwNThreads > 0) {
            FFTW(cleanup_threads)();
            FFTW(cleanup)();
        }
    }
    gsl_wavelet_free(fHdl->gslDWT);
    gsl_wavelet_workspace_free(fHdl->gslDWTWork);
    if(fHdl->waveletWav)
        free(fHdl->waveletWav);
    return 0;
}

/* Filters should directly write into fHdl->respWav the real space
 * response waveform in wrapped around order */

int filters_SavitzkyGolay(filters_t *fHdl, int m, int ld)
/* m: order of polynomial, np: number of points, ld: degree of derivative*/
{
    int np;
    ANALYSIS_WAVEFORM_BASE_TYPE *c;
    int ipj, imj, mm, j, k, nl, nr;
    double fac, sum;
    gsl_permutation * p;
    gsl_vector *b;
    gsl_matrix *a;

    np = fHdl->respLen;
    if(np<1 || np<m-1 || np%2==0 || ld>m || np!=fHdl->respLen) {
        error_printf("%s(): improper arguments, returning...\n", __FUNCTION__);
        return 1;
    }

    c = calloc(np, sizeof(ANALYSIS_WAVEFORM_BASE_TYPE));

    p = gsl_permutation_alloc (m+1);
    b = gsl_vector_alloc(m+1);
    a = gsl_matrix_alloc(m+1, m+1);

    nl = np/2;
    nr = nl;

    for(ipj=0;ipj<=(m << 1);ipj++) {
        sum=(ipj ? 0.0 : 1.0);
        for(k=1;k<=nr;k++) sum += pow((double)(k),(double)(ipj));
        for(k=1;k<=nl;k++) sum += pow((double)(-k),(double)(ipj));
        mm=MIN(ipj,2*m-ipj);
        for(imj=-mm;imj<=mm;imj+=2) gsl_matrix_set(a,(ipj+imj)/2,(ipj-imj)/2,sum);
    }

    gsl_linalg_LU_decomp(a, p, &k);
    for (j=0;j<m+1;j++) gsl_vector_set(b,j,0.0);
    gsl_vector_set(b,ld,1.0);

    gsl_linalg_LU_solve (a, p, b, b);

    for(k = -nl;k<=nr;k++) {
        sum = gsl_vector_get(b,0);
        fac = 1.0;
        for (mm=1;mm<=m;mm++) sum += gsl_vector_get(b,mm)*(fac *= k);

        j=(np-k) % np;
        c[j]=sum; // c is in wraparound order, convenient for fft convolute

        // c[nl + k] = sum;
    }
    memcpy(fHdl->respWav, c, np * sizeof(ANALYSIS_WAVEFORM_BASE_TYPE));
/*
    for(j=0; j<np; j++) {
        fprintf(stderr, "%g\n", c[j]);
    }
*/
    gsl_vector_free(b);
    gsl_matrix_free(a);
    gsl_permutation_free(p);

/*
    for(k=nl; k<wavlen - nr; k++) {
        sum = 0.0;
        for(j=0; j<np; j++) {
            sum += c[j] * inwav[k+j-nl];
            outwav[k] = sum;
        }
    }
*/
    free(c);
    return 0;
}

int filters_raisedCosine(filters_t *fHdl)
{
    ssize_t i;
    ANALYSIS_WAVEFORM_BASE_TYPE x;

    for(i=0; i<(fHdl->respLen+1)/2; i++) { /* positive side */
        x = 2.0*M_PI/(ANALYSIS_WAVEFORM_BASE_TYPE)(fHdl->respLen-1)*(ANALYSIS_WAVEFORM_BASE_TYPE)i;
        fHdl->respWav[i] = (1.0 + cos(x))/(ANALYSIS_WAVEFORM_BASE_TYPE)fHdl->respLen;
    }
    for(i=-((fHdl->respLen-1)/2); i<0; i++) { /* negative side */
        /* i=-(fHdl->respLen-1)/2 cast to wrong value */
        x = 2.0*M_PI/(ANALYSIS_WAVEFORM_BASE_TYPE)(fHdl->respLen-1)*(ANALYSIS_WAVEFORM_BASE_TYPE)i;
        fHdl->respWav[fHdl->respLen+i] = (1.0 + cos(x))/(ANALYSIS_WAVEFORM_BASE_TYPE)fHdl->respLen;
    }
    return 0;
}

int filters_convolute(filters_t *fHdl)
{
    size_t i;
    ANALYSIS_WAVEFORM_BASE_TYPE re, im;

    for(i=0; i<fHdl->wavLen; i++) {
        fHdl->fftwWork[i] = fHdl->inWav[i];
    }
    for(i=fHdl->wavLen; i<fHdl->fftLen; i++) {
        fHdl->fftwWork[i] = 0.0;
    }

    // fill in with respwav in wrap-around order
    fHdl->fftwWork1[0] = fHdl->respWav[0];
    for(i=1; i<(fHdl->respLen+1)/2; i++) {
        fHdl->fftwWork1[i] = fHdl->respWav[i];
        fHdl->fftwWork1[fHdl->fftLen-i] = fHdl->respWav[fHdl->respLen-i];
    }
    for(i=(fHdl->respLen+1)/2; i<(fHdl->fftLen-(fHdl->respLen+1)/2); i++) {
        fHdl->fftwWork1[i] = 0.0;
    }

    // do fft
    FFTW(execute)(fHdl->fftwPlan);
    FFTW(execute)(fHdl->fftwPlan1);

    // multiply in complex fourier space, half-complex format
    fHdl->fftwWork[0] = fHdl->fftwWork[0] * fHdl->fftwWork1[0]
        / (ANALYSIS_WAVEFORM_BASE_TYPE)fHdl->fftLen;
    for(i=1; i<fHdl->fftLen/2; i++) {
        re = fHdl->fftwWork[i] * fHdl->fftwWork1[i]
            - fHdl->fftwWork[fHdl->fftLen-i] * fHdl->fftwWork1[fHdl->fftLen-i];
        im = fHdl->fftwWork[i] * fHdl->fftwWork1[fHdl->fftLen-i]
            + fHdl->fftwWork[fHdl->fftLen-i] * fHdl->fftwWork1[i];
        fHdl->fftwWork[i] = re / (ANALYSIS_WAVEFORM_BASE_TYPE)fHdl->fftLen;
        fHdl->fftwWork[fHdl->fftLen-i] = im / (ANALYSIS_WAVEFORM_BASE_TYPE)fHdl->fftLen;
    }
    fHdl->fftwWork[fHdl->fftLen/2] = fHdl->fftwWork[fHdl->fftLen/2] * fHdl->fftwWork1[fHdl->fftLen/2]
        / (ANALYSIS_WAVEFORM_BASE_TYPE)fHdl->fftLen;

    // ifft
    FFTW(execute)(fHdl->fftwPlan2);

    // copy the output to outwav
    memcpy(fHdl->outWav, fHdl->fftwWork, fHdl->wavLen * sizeof(ANALYSIS_WAVEFORM_BASE_TYPE));
    return 0;
}

int filters_hanning_window(filters_t *fHdl)
{
    size_t i;

    fHdl->fftwS1 = 0.0; fHdl->fftwS2 = 0.0;
    for(i=0; i<fHdl->fftLen; i++) {
        fHdl->fftwWin[i] = 0.5 * (1.0 - cos(2*M_PI*i/(double)fHdl->fftLen));
        fHdl->fftwS1 += fHdl->fftwWin[i];
        fHdl->fftwS2 += fHdl->fftwWin[i] * fHdl->fftwWin[i];
    }

    return 0;
}

int filters_fft_spectrum(filters_t *fHdl)
{
    size_t i;

    for(i=0; i<fHdl->fftLen; i++) {
        fHdl->fftwWork[i] = fHdl->inWav[i] * fHdl->fftwWin[i];
    }
    FFTW(execute)(fHdl->fftwPlan);

    /* Compute linearized power spectrum into fftwWork1, in [V] for example, normalized.
     * Total length should be (int)(n/2)+1 */
    fHdl->fftwWork1[0] = fHdl->fftwWork[0] / fHdl->fftwS1;
    for(i=1; i<(fHdl->fftLen+1)/2; i++) {
        fHdl->fftwWork1[i] = hypot(fHdl->fftwWork[i],
                                   fHdl->fftwWork[fHdl->fftLen - i]) * sqrt(2.0) / fHdl->fftwS1;
    }
    if(fHdl->fftLen % 2 == 0) { /* even number */
        fHdl->fftwWork1[fHdl->fftLen/2] = fHdl->fftwWork[fHdl->fftLen/2] / fHdl->fftwS1;
    }
    /* For spectra density, normalization should be * sqrt(2.0 * dt / fftwS2) */
    for(i=0; i<=fHdl->fftLen/2; i++) {
        fHdl->fftwWork[i] = fHdl->fftwWork1[i] * fHdl->fftwS1 * sqrt(fHdl->dt / fHdl->fftwS2);
    }
    return 0;
}

int filters_DWT(filters_t *fHdl) /* discrete wavelet transform */
{
    gsl_wavelet_transform_forward(fHdl->gslDWT, fHdl->waveletWav, 1, fHdl->wavLen,
                                  fHdl->gslDWTWork);
    return 0;
}

int filters_median(filters_t *fHdl, size_t n) /* median filter with moving window size n */
{
    size_t i, mid=n/2;

    for(i=0; i<mid;i++) /* overhang at the beginning */
        fHdl->outWav[i] = quickselect(fHdl->inWav, mid+i+1, (mid+i+1)/2);
    for(i=0; i<fHdl->wavLen-1 - mid; i++)
        fHdl->outWav[mid+i] = quickselect(fHdl->inWav + i, n, n/2);
    for(i=0; i<mid; i++) /* overhang at the end */
        fHdl->outWav[fHdl->wavLen-1 - i] = quickselect(fHdl->inWav + fHdl->wavLen - (mid-i+1),
                                                       mid-i+1, (mid-i+1)/2);
    return 0;
}

int filters_trapezoidal(filters_t *fHdl, size_t k, size_t l, double M)
{
    double s, pp;
    ssize_t i, j, jk, jl, jkl;
    double vj, vjk, vjl, vjkl, dkl;

    s = 0.0; pp = 0.0;

    for(i=0; i<fHdl->wavLen; i++) {
        j=i; jk = j-k; jl = j-l; jkl = j-k-l;
        vj   = j>=0   ? fHdl->inWav[j]   : fHdl->inWav[0];
        vjk  = jk>=0  ? fHdl->inWav[jk]  : fHdl->inWav[0];
        vjl  = jl>=0  ? fHdl->inWav[jl]  : fHdl->inWav[0];
        vjkl = jkl>=0 ? fHdl->inWav[jkl] : fHdl->inWav[0];

        dkl = vj - vjk - vjl + vjkl;
        pp = pp + dkl;
        if(M>=0.0) {
            s = s + pp + dkl * M;
        } else { /* infinit decay time, so the input is a step function */
            s = s + dkl;
        }
        fHdl->outWav[i] = s / (fabs(M) * (double)k);
    }
    return 0;
}

#ifdef FILTERS_DEBUG_ENABLEMAIN
int main(int argc, char **argv)
{
    #define PLEN 1503
    ANALYSIS_WAVEFORM_BASE_TYPE pulse[PLEN] = {0.0};
    size_t i;
    filters_t *fHdl;

    i = 2;
    pulse[i++]=0.0; pulse[i++]=1.0; pulse[i++]=10.0; pulse[i++]=8.0; pulse[i++]=6.0;
    pulse[i++]=4.0; pulse[i++]=2.0; pulse[i++]=1.0; pulse[i++]=0.5; pulse[i++]=0.2;
    i = 100;
    pulse[i--]=0.0; pulse[i--]=1.0; pulse[i--]=10.0; pulse[i--]=8.0; pulse[i--]=6.0;
    pulse[i--]=4.0; pulse[i--]=2.0; pulse[i--]=1.0; pulse[i--]=0.5; pulse[i--]=0.2;

    fHdl = filters_init(pulse, PLEN);
    filters_median(fHdl, 11);
    for(i=0; i<fHdl->wavLen; i++) {
        printf("%g %g\n", fHdl->inWav[i], fHdl->outWav[i]);
    }
    printf("\n\n");

#if 0
    for(i=0; i<PLEN; i++) {
        pulse[i] += 10.0 * cos(2.0 * M_PI/3.0 * i);
    }

    fHdl = filters_init_for_convolution(pulse, PLEN, 0);
    filters_fft_spectrum(fHdl);
    for(i=0; i<fHdl->fftLen; i++) {
        printf("%g %g\n", fHdl->fftwWork[i], fHdl->fftwWork1[i]);
    }
    printf("\n\n");
#endif
#if 0
    fHdl = filters_init_for_convolution(NULL, PLEN, 31);
    memcpy(fHdl->inWav, pulse, PLEN * sizeof(ANALYSIS_WAVEFORM_BASE_TYPE));

    filters_raisedCosine(fHdl);
    filters_convolute(fHdl);

    for(i=0; i<fHdl->fftLen/*PLEN*/; i++) {
        // printf("%g %g\n", pulse[i], fHdl->outWav[i]);
        printf("%g %g\n", fHdl->fftwWork[i], fHdl->fftwWork1[i]);
    }
    printf("\n\n");

    filters_SavitzkyGolay(fHdl, 5, 0);
    filters_convolute(fHdl);

    for(i=0; i<fHdl->fftLen/*PLEN*/; i++) {
        // printf("%g %g\n", pulse[i], fHdl->outWav[i]);
        printf("%g %g\n", fHdl->fftwWork[i], fHdl->fftwWork1[i]);
    }
#endif

    filters_close(fHdl);
    return EXIT_SUCCESS;
}
#endif