text:
/*
 * The MIT License is a permissive free software license, which permits reuse
 * within both open source and proprietary software. The software is licensed
 * as is, and no warranty is given as to fitness for purpose or absence of
 * infringement of third parties' rights, such as patents. Generally use for
 * research activities is allowed regardless of any third party patents, but
 * commercial use may be subject to a separate license.
 *
 *
 * The MIT License (MIT)
 *
 * Copyright (c) 2015 Tuomo Raitio
 *
 * 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.
 *
 *
 *
 *
 * <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
 *                  GlottHMM Speech Synthesis
 * <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
 *
 * This program reads speech parameters and a glottal pulse/
 * pulse library, and synthesizes speech from them.
 *
 * This program has been written in Aalto University,
 * Department of Signal Processign and Acoustics, Espoo, Finland
 *
 * Main author: Tuomo Raitio
 * Acknowledgements: Antti Suni, Paavo Alku, Martti Vainio
 *
 * File Synthesis.c
 * Version: 1.1
 *
 */

/***********************************************/
/*                 INCLUDE                     */
/***********************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <gsl/gsl_vector.h>		/* GSL, Vector */
#include <gsl/gsl_matrix.h>		/* GSL, Matrix */
#include <gsl/gsl_fft_real.h>		/* GSL, FFT */
#include <libconfig.h>				/* Configuration file */
#include "SynthesisFunctions.h"

/*******************************************************************/
/*                          MAIN                                   */
/*******************************************************************/

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

	/* Start counter */
	double time1 = (double)clock();

	/* Check command line format */
	if(Check_command_line(argc) == EXIT_FAILURE)
		return EXIT_FAILURE;

	/* Read default configuration file and assign parameters */
	struct config_t *conf_def = Read_config(argv[2]);
	if(conf_def == NULL)
		return EXIT_FAILURE;
	PARAM params;
	if(Assign_config_parameters(argv[1],conf_def,&params,DEF_CONF) == EXIT_FAILURE) {
		config_destroy(conf_def);
		free(conf_def);
		return EXIT_FAILURE;
	}

	/* Read user configuration file and assign parameters */
	if(argc == 4) {
		struct config_t *conf_usr = Read_config(argv[3]);
		if(conf_usr == NULL)
			return EXIT_FAILURE;
		if(Assign_config_parameters(argv[1],conf_usr,&params,USR_CONF) == EXIT_FAILURE) {
			config_destroy(conf_usr);
			free(conf_usr);
			return EXIT_FAILURE;
		}
	}

	/* Print synthesis settings */
	Print_synthesis_settings_start(&params);

	/*******************************************************************/
	/*                 READ PULSE LIBRARY / PULSE                      */
	/*******************************************************************/

	/* Allocate space for pulse library and read data */
	gsl_matrix *pulses,*pulses_rs,*plsf,*ptilt,*pharm,*phnr,*pwaveform,*pca_pc,*pca_w_lib;
	gsl_vector *pgain,*pulse_lengths,*ph1h2,*pnaq,*pca_mean,*stoch_env,*stoch_sp;
	if(Read_pulse_library(&params,&pulses,&pulses_rs,&plsf,&ptilt,&pharm,&phnr,&pwaveform,&pca_pc,
			&pca_w_lib,&stoch_env,&stoch_sp,&pgain,&ph1h2,&pnaq,&pca_mean,&pulse_lengths) == EXIT_FAILURE)
		return EXIT_FAILURE;

	/* Read pulse file if pulse library is not used */
	gsl_vector *original_pulse = NULL;
	if(params.use_pulselib == 0) {
		original_pulse = ReadPulseFile(&params);
		if(original_pulse == NULL)
			return EXIT_FAILURE;
	}

	/* Read DNN pulse generation weights and allocate parameters */
	gsl_matrix **DNN_W = (gsl_matrix**)malloc(params.dnn_weight_dims->size/2*sizeof(gsl_matrix*));
	if(Read_DNN_weights(&params,DNN_W) == EXIT_FAILURE)
		return EXIT_FAILURE;
	gsl_vector **input_minmax = (gsl_vector**)malloc(sizeof(gsl_vector*));
	if(Read_input_minmax(&params,input_minmax) == EXIT_FAILURE)
			return EXIT_FAILURE;
	gsl_vector *dnnpulseindices = NULL;
	gsl_vector *dnnpulses = NULL;

	/* Define variables */
	gsl_vector *gain,*fundf,*excitation_voiced,*excitation_unvoiced;
	gsl_vector *gain_new,*pulse_clus_id,*h1h2,*naq,*resynthesis_pulse_index;
	gsl_matrix *LSF,*LSF2,*glflowsp,*glflowsp_new,*hnr,*hnr_new,*harmonics,*waveform,*pulse_clusters,*pca_w,*LSF_interp = NULL;
	int i;

	/****************************************************************************************/
	/*                          START SYNTHESING FILE(S)                                    */
	/****************************************************************************************/

	/* Start loop for synthesis list */
	for(i=0; i<params.synlistlen; i++) {

		/* Initialize params */
		if(Initialize_params(&params,i) == EXIT_FAILURE)
			return EXIT_FAILURE;

		/* Print synthesis settings */
		Print_synthesis_settings_middle(&params);

		/* Compatibility check */
		if(Compatibility_check(&params) == EXIT_FAILURE)
			return EXIT_FAILURE;

		/*****************************************************************************/
		/*                     LOAD AND MODIFY PARAMETERS                            */
		/*****************************************************************************/

		/* Allocate and read synthesis parameters, initialize pulse clustering if used */
		Allocate_params(&excitation_voiced,&excitation_unvoiced,&resynthesis_pulse_index,&gain_new,&glflowsp_new,&hnr_new,&params);
		if(Read_synthesis_parameters(&gain,&fundf,&LSF,&LSF2,&glflowsp,&hnr,&harmonics,&waveform,&h1h2,&naq,&pca_w,&params) == EXIT_FAILURE) continue;
		if(Pulse_clustering(&pulse_clus_id, &pulse_clusters, &params) == EXIT_FAILURE) continue;

		/* Miscellaneous operations */
		Convert_logF0_to_lin(fundf,&params);
		Merge_voiced_unvoiced_spectra(LSF,LSF2,fundf,&params);
		Integrate_LSFs(&LSF,&params);
		if(params.use_tilt == 1)
			Integrate_LSFs(&glflowsp,&params);
		LSF_fix_matrix(LSF);
		LSF_fix_matrix(glflowsp);
		Smooth_matrix(hnr,params.hnr_smooth_len);
		if(params.use_hmm == 0) Smooth_matrix(harmonics,params.harmonics_smooth_len);
		if(params.use_hmm == 0) MA(gain,params.gain_smooth_len);

		/* Modification for noise robust speech */
		Noise_robust_speech2(gain,harmonics,&params);

		/* Formant enhancement for LSFs */
		Postfilter(LSF,&params);

		/* Apply noise reduction */
		Noise_reduction(gain,&params);

		/* Normalize pulse library parameters according to synthesis parameters (normalize mean) */
		Normalize_pulse_library_var(plsf,ptilt,pharm,phnr,pwaveform,pgain,ph1h2,pnaq,LSF,glflowsp,
				harmonics,hnr,waveform,gain,h1h2,naq,fundf,&params);

		/* Adapt synthesis parameters according to pulse library parameters (normalize mean) */
		Adapt_synthesis_parameters_var(plsf,ptilt,pharm,phnr,pwaveform,pgain,ph1h2,pnaq,LSF,glflowsp,
				harmonics,hnr,waveform,gain,h1h2,naq,fundf,pulse_lengths,&params);

		/* Select vocal tract LSFs from pulse library */
		// TODO: Viterbi
		int winsize = 5;
		Select_LSFs_from_pulse_library(LSF,plsf,fundf,winsize,&params);

		/* Smooth and interpolate LSFs to signal length */
		LSF_interp = gsl_matrix_alloc(params.signal_length,params.lpc_order_vt);
		Smooth_interp_lsf(LSF_interp,LSF,params.signal_length,params.use_hmm,params.lsf_smooth_len);

		/* DNN pulse gen */
		if(params.use_dnn_pulsegen == 1) {
			dnnpulseindices = gsl_vector_alloc(params.signal_length);
			dnnpulses = gsl_vector_alloc(2*params.signal_length);
		}

		/*******************************************************************/
		/*                    CREATE EXCITATION                            */
		/*******************************************************************/

		printf("	- Creating excitation...\n");

		/* Create excitation for estimating the HNR of the voiced excitation */
		CreateExcitation(&params,excitation_voiced,excitation_unvoiced,fundf,gain,LSF,glflowsp,hnr,harmonics,waveform,h1h2,naq,
				original_pulse,pulses,pulses_rs,pulse_lengths,pgain,plsf,ptilt,phnr,pharm,pwaveform,ph1h2,pnaq,
				resynthesis_pulse_index,pulse_clus_id,pulse_clusters,gain,glflowsp_new,hnr_new,pca_mean,pca_pc,pca_w,pca_w_lib,
				stoch_env,stoch_sp,DNN_W,input_minmax,dnnpulseindices,dnnpulses);

		/* Compensate HNR etc */
		HNR_compensation(hnr,hnr_new,&params);
		Fill_pulse_indices(resynthesis_pulse_index);
		params.resynth = 1;

		/* Create excitation */
		CreateExcitation(&params,excitation_voiced,excitation_unvoiced,fundf,gain,LSF,glflowsp,hnr,harmonics,waveform,h1h2,naq,
				original_pulse,pulses,pulses_rs,pulse_lengths,pgain,plsf,ptilt,phnr,pharm,pwaveform,ph1h2,pnaq,
				resynthesis_pulse_index,pulse_clus_id,pulse_clusters,gain,glflowsp_new,hnr_new,pca_mean,pca_pc,pca_w,pca_w_lib,
				stoch_env,stoch_sp,DNN_W,input_minmax,dnnpulseindices,dnnpulses);
		Print_elapsed_time(&params);

		/* Compensate for the pre-emphasis during analysis (de-emphasis) */
		if(params.unvoiced_pre_emphasis == 1)
			Integrate(excitation_unvoiced, 0.97); // Slightly less than leaky factor

		/*******************************************************************/
		/*                    SPECTRAL MATCHING                            */
		/*******************************************************************/

		if((params.use_pulselib == 0 && params.use_dnn_pulsegen == 0) || (params.use_pulselib_pca == 1 && params.pca_spectral_matching == 1) ||
			 (params.use_dnn_pulsegen == 1 && params.use_dnn_specmatch == 1)) {
			printf("	- Spectral matching...\n");
			Spectral_match(excitation_voiced, glflowsp, glflowsp_new, &params);
			if(params.use_pulselib_pca == 0 && params.use_dnn_pulsegen == 0 && params.two_pitch_period_diff_pulse == 0)
				LipRadiation(excitation_voiced);
			Print_elapsed_time(&params);
		}

		/*******************************************************************/
		/*                    FILTER EXCITATION                            */
		/*******************************************************************/

		/* Combine excitations and filter */
		printf("	- Filtering excitation...\n");
		gsl_vector_add(excitation_voiced,excitation_unvoiced);
		Filter_excitation(excitation_voiced,LSF_interp,&params);
		Print_elapsed_time(&params);

		/*******************************************************************/
		/*              RESYNTHESIS FOR GAIN NORMALIZATION                 */
		/*******************************************************************/

		/* Evaluate new gain for re-synthesis */
		printf("	- Re-synthesis for gain normalization...\n");
		Evaluate_new_gain(excitation_voiced,gain_new,gain,fundf,&params);

		/* Re-synthesize with new gain */
		CreateExcitation(&params,excitation_voiced,excitation_unvoiced,fundf,gain_new,LSF,glflowsp,hnr,harmonics,waveform,h1h2,naq,
				original_pulse,pulses,pulses_rs,pulse_lengths,pgain,plsf,ptilt,phnr,pharm,pwaveform,ph1h2,pnaq,
				resynthesis_pulse_index,pulse_clus_id,pulse_clusters,gain,glflowsp_new,hnr_new,pca_mean,pca_pc,pca_w,pca_w_lib,
				stoch_env,stoch_sp,DNN_W,input_minmax,dnnpulseindices,dnnpulses);
		if(params.unvoiced_pre_emphasis == 1)
			Integrate(excitation_unvoiced, 0.97); // Slightly less than leaky factor

		/* Spectral matching */
		if((params.use_pulselib == 0 && params.use_dnn_pulsegen == 0) || (params.use_pulselib_pca == 1 && params.pca_spectral_matching == 1) ||
			 (params.use_dnn_pulsegen == 1 && params.use_dnn_specmatch == 1)) {
			Spectral_match(excitation_voiced, glflowsp, glflowsp_new, &params);
			if(params.use_pulselib_pca == 0 && params.use_dnn_pulsegen == 0 && params.two_pitch_period_diff_pulse == 0)
				LipRadiation(excitation_voiced);
		}

		/* Combine excitations and filter */
		gsl_vector_add(excitation_voiced,excitation_unvoiced);

		/* Save excitation to wav file */
		Save_excitation_to_wav(excitation_voiced,&params);

		/* Filter excitation */
		Filter_excitation(excitation_voiced,LSF_interp,&params);

		/*******************************************************************/
		/*              POSTPROCESSING AND SAVE TO FILE                    */
		/*******************************************************************/

		/* Filter out signal below F0 */
		Hp_filt_below_f0(excitation_voiced,fundf,&params);

		/* Compression of speech signal (only in noise robust speech) */
		Compression(excitation_voiced, params.compcoeff);

		/* Scale (if > 1) and save signal to file */
		Scale_signal(excitation_voiced,SCALE_IF_GREATER_THAN_ONE);
		Save_signal_to_file(excitation_voiced,&params,NULL);

		/* Report */
		Print_synthesis_settings_end(&params,time1);
	}

	/*******************************************************************/
	/*                     FREE MEMORY AND FINISH                      */
	/*******************************************************************/

	/* Free memory */
	Free_pulselib_variables(pulses,pulses_rs,pwaveform,pulse_lengths,plsf,ptilt,pharm,phnr,pgain,ph1h2,pnaq,pca_mean,pca_pc,pca_w_lib,stoch_env,stoch_sp,&params);
	Free_variables(original_pulse,excitation_voiced,excitation_unvoiced,fundf,gain,gain_new,LSF,LSF2,LSF_interp,glflowsp,glflowsp_new,
			hnr,hnr_new,harmonics,waveform,h1h2,naq,resynthesis_pulse_index,pulse_clus_id,pulse_clusters,pca_w,
			DNN_W,input_minmax,dnnpulseindices,dnnpulses,&params);

	/* Exit */
  	return EXIT_SUCCESS;
}

/***********/
/*   EOF   */
/***********/

max_stars_repo_name: mjansche/GlottHMM

max_issues_repo_name: mjansche/GlottHMM

max_forks_repo_name: mjansche/GlottHMM


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_errno.h>

#include "avl.c"

int
gsl_spmatrix_memcpy(gsl_spmatrix *dest, const gsl_spmatrix *src)
{
  const size_t M = src->size1;
  const size_t N = src->size2;

  if (M != dest->size1 || N != dest->size2)
    {
      GSL_ERROR("matrix sizes are different", GSL_EBADLEN);
    }
  else if (dest->sptype != src->sptype)
    {
      GSL_ERROR("cannot copy matrices of different storage formats",
                GSL_EINVAL);
    }
  else
    {
      int s = GSL_SUCCESS;
      size_t n;

      if (dest->nzmax < src->nz)
        {
          s = gsl_spmatrix_realloc(src->nz, dest);
          if (s)
            return s;
        }

      /* copy indices and data to dest */
      if (GSL_SPMATRIX_ISTRIPLET(src))
        {
          void *ptr;

          for (n = 0; n < src->nz; ++n)
            {
              dest->i[n] = src->i[n];
              dest->p[n] = src->p[n];
              dest->data[n] = src->data[n];

              /* copy binary tree data */
              ptr = avl_insert(dest->tree_data->tree, &dest->data[n]);
              if (ptr != NULL)
                {
                  GSL_ERROR("detected duplicate entry", GSL_EINVAL);
                }
            }
        }
      else if (GSL_SPMATRIX_ISCCS(src))
        {
          for (n = 0; n < src->nz; ++n)
            {
              dest->i[n] = src->i[n];
              dest->data[n] = src->data[n];
            }

          for (n = 0; n < src->size2 + 1; ++n)
            {
              dest->p[n] = src->p[n];
            }
        }
      else if (GSL_SPMATRIX_ISCRS(src))
        {
          for (n = 0; n < src->nz; ++n)
            {
              dest->i[n] = src->i[n];
              dest->data[n] = src->data[n];
            }

          for (n = 0; n < src->size1 + 1; ++n)
            {
              dest->p[n] = src->p[n];
            }
        }
      else
        {
          GSL_ERROR("invalid matrix type for src", GSL_EINVAL);
        }

      dest->nz = src->nz;

      return s;
    }
} /* gsl_spmatrix_memcpy() */

ext: c

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* stable/stable_fit.c
 *
 * Functions employed by different methods of estimation implemented
 * in Libstable.
 *
 * Copyright (C) 2013. Javier Royuela del Val
 *                     Federico Simmross Wattenberg
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; version 3 of the License.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; If not, see <http://www.gnu.org/licenses/>.
 *
 *
 *  Javier Royuela del Val.
 *  E.T.S.I. Telecomunicación
 *  Universidad de Valladolid
 *  Paseo de Belén 15, 47002 Valladolid, Spain.
 *  jroyval@lpi.tel.uva.es
 */
#include "stable.h"
#include "mcculloch.h"

#include <gsl/gsl_complex.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_fft_real.h>

void stable_fft(double *data, const unsigned int length, double * y)
{
  //int i;

  memcpy ( (void *)y, (const void *) data, length*sizeof(double));

  gsl_fft_real_radix2_transform (y, 1, length);

  return;
}

double stable_loglikelihood(StableDist *dist, double *data, const unsigned int length)
{
  double *pdf=NULL;
  double l=0.0;
  int i;

  pdf=(double*)malloc(sizeof(double)*length);

  stable_pdf(dist,data,length,pdf,NULL);

  for(i=0;i<length;i++)
    {
      if (pdf[i]>0.0) l+=log(pdf[i]);
    }

  free(pdf);
  return l;
}

double stable_loglike_p(stable_like_params *params)
{
  double *pdf;
  double l=0.0;
  int i;

  pdf=(double*)malloc(sizeof(double)*(params->length));

  stable_pdf(params->dist,params->data,params->length,pdf,NULL);

  for(i=0;i<params->length;i++)
    {
      if (pdf[i]>0.0) {l+=log(pdf[i]);}
    }

  free(pdf);

  return l;
}

double stable_minusloglikelihood(const gsl_vector * theta, void * p)
{
/* Cost function to minimize, with the estimation of sigma and mu given by McCulloch at each iteration*/
  double alpha=1, beta=0, sigma=1.0, mu=0.0;
  double minusloglike=0;
  stable_like_params * params = (stable_like_params *) p;

  alpha = gsl_vector_get(theta,0);
  beta = gsl_vector_get(theta,1);

  /* Update sigma and mu with McCulloch. It needs nu_c nu_z*/
  czab(alpha, beta, params->nu_c, params->nu_z, &sigma, &mu);

  /* Check that the parameters are valid */
  if(stable_setparams(params->dist, alpha, beta, sigma, mu, 0) < 0)
    {
      return GSL_NAN;
    }
  else minusloglike = -stable_loglike_p(params);

  if (isinf(minusloglike) || isnan(minusloglike)) minusloglike=GSL_NAN;

  return minusloglike;
}

int compare (const void * a, const void * b)
{
/* qsort compare function */
  return ((*(double *)b < *(double *)a) - (*(double *)a < *(double *)b));
}

inline void get_original(const gsl_vector *s,double *a,double *b,double *c,double *m)
{
  *a = M_2_PI*atan(gsl_vector_get(s,0))+1.0;
  *b = M_2_PI*atan(gsl_vector_get(s,1));
  *c = exp(gsl_vector_get(s,2));
  *m = gsl_vector_get(s,3);
}
inline void set_expanded(gsl_vector *s,const double a,const double b,const double c,const double m)
{
  gsl_vector_set(s,0,tan(M_PI_2*(a-1.0)));
  gsl_vector_set(s,1,tan(M_PI_2*b));
  gsl_vector_set(s,2,log(c));
  gsl_vector_set(s,3,m);
}

double stable_minusloglikelihood_whole(const gsl_vector * theta, void * p)
{
/* Whole cost function to minimize in a 4D parameter space */
  double alpha=1, beta=0, sigma=1.0, mu=0.0;
  double minusloglike=0;
  stable_like_params * params = (stable_like_params *) p;

  get_original(theta,&alpha,&beta,&sigma,&mu);

  /* Check that the parameters are valid */
  if(stable_setparams(params->dist, alpha, beta, sigma, mu, 0) < 0)
    {
      perror("setparams error");
      return GSL_NAN;
    }
  else minusloglike = -stable_loglike_p(params);

  if (isinf(minusloglike) || isnan(minusloglike)) minusloglike=GSL_NAN;

  return minusloglike;
}

void stable_fit_init(StableDist *dist, const double * data, const unsigned int length, double *pnu_c,double *pnu_z)
{
  /* McCulloch estimation */

  double *sorted=NULL;
  double alpha0, beta0, sigma0, mu0;

  /* We need to sort the data to get percentiles */
  sorted = (double*)malloc(length*sizeof(double));
  memcpy ( (void *)sorted, (const void *) data, length*sizeof(double));
  qsort  ( sorted, length, sizeof(double), compare);

  /* Estimate the parameters. */
  stab((const double *) sorted,length,0,&alpha0,&beta0,&sigma0,&mu0);

  /* Set parameters in the distribution */
  if(stable_setparams(dist,alpha0,beta0,sigma0,mu0, 0)<0)
    {
      perror("INITIAL ESTIMATED PARAMETER ARE NOT VALID");
      return;
    }

  /* Get pnu_c and pnu_z needed for mle2d estimation */
  cztab(sorted, length, pnu_c, pnu_z);

  free(sorted);
  return;
}

int stable_fit_iter(StableDist *dist, const double * data, const unsigned int length,const double nu_c,const double nu_z)
{
  const gsl_multimin_fminimizer_type *T;
  gsl_multimin_fminimizer *s;

  gsl_multimin_function likelihood_func;

  gsl_vector *theta, *ss;

  unsigned int iter = 0;
  int status=0;
  double size=0;

  double a=1,b=0.0,c=1,m=0.0;
  stable_like_params par;

  par.dist=dist;
  par.data=(double *)data;
  par.length=length;
  par.nu_c=nu_c;
  par.nu_z=nu_z;

  /* Inicio: Debe haberse inicializado dist con alpha y beta de McCulloch */
  theta=gsl_vector_alloc(2);
  gsl_vector_set (theta, 0, dist->alpha);
  gsl_vector_set (theta, 1, dist->beta);

  #ifdef DEBUG
    printf("%lf, %lf\n",gsl_vector_get (theta, 0),gsl_vector_get (theta, 1));
  #endif

  /* Saltos iniciales */
  ss = gsl_vector_alloc (2);
  gsl_vector_set_all (ss, 0.01);

  /* Funcion a minimizar */
  likelihood_func.n = 2; // Dimension 2 (alpha y beta)
  likelihood_func.f = &stable_minusloglikelihood;
  likelihood_func.params = (void *) (&par);  // Parametros de la funcion

   /* Creacion del minimizer */
  T = gsl_multimin_fminimizer_nmsimplex2rand;

  s = gsl_multimin_fminimizer_alloc (T, 2); /* Dimension 2*/

  /* Poner funcion, estimacion inicial, saltos iniciales */
  gsl_multimin_fminimizer_set (s, &likelihood_func, theta, ss);

  #ifdef DEBUG
    printf("5\n");
  #endif

  /* Iterar */
  do
    {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);
 //     if (status!=GSL_SUCCESS) {
 //       printf("Minimizer warning: %s\n",gsl_strerror(status));
 //       fflush(stdout);
 //      }

      size = gsl_multimin_fminimizer_size (s);
      status = gsl_multimin_test_size (size, 0.02);
/*
      if (status == GSL_SUCCESS)
        {
          printf ("              converged to minimum at\n");
        }

      printf ("%5d %1.5f %1.5f %1.5f %1.5f f() = %1.8e size = %.5f\n",
              (int)iter,
              gsl_vector_get (s->x, 0),
              gsl_vector_get (s->x, 1),
              p->dist->sigma,
              p->dist->mu_1,
              s->fval, size);
        //}
*/
    } while (status == GSL_CONTINUE && iter < 200);

//  if (status!=GSL_SUCCESS)
//    {
//      printf("Minimizer warning: %s\n",gsl_strerror(status));
//      fflush(stdout);
//    }

  /* Se recupera la estimacion alpha y beta */
  gsl_vector_free(theta);
/*
  theta = gsl_multimin_fminimizer_x (s);
  a = gsl_vector_get (theta, 0);
  b = gsl_vector_get (theta, 1);
*/
  a = gsl_vector_get (s->x, 0);
  b = gsl_vector_get (s->x, 1);

  /* Y se estima sigma y mu para esos alpha y beta */
  czab(a, b, nu_c, nu_z, &c, &m);

  //printf("%5d %10.3e %10.3e %10.3e %10.3e\n",(int)iter,a,b,c,m);

  // Se almacena el punto estimado en la distribucion, comprobando que es valido
  if (stable_setparams(dist,a,b,c,m,0)<0)
   {
    perror("FINAL ESTIMATED PARAMETER ARE NOT VALID\n");
   }

  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);

  return status;
}

int stable_fit(StableDist *dist, const double *data, const unsigned int length)
{
  double nu_c=0.0,nu_z=0.0;
  int status = 0;

  stable_fit_init(dist,data,length,&nu_c,&nu_z);
  status=stable_fit_iter(dist,data,length,nu_c,nu_z);

  return status;
}

int stable_fit_iter_whole(StableDist *dist, const double * data, const unsigned int length)
{
  const gsl_multimin_fminimizer_type *T;
  gsl_multimin_fminimizer *s;

  gsl_multimin_function likelihood_func;

  gsl_vector *theta, *ss;

  unsigned int iter = 0;
  int status=0;
  double size=0;

  double a=1,b=0.0,c=1,m=0.0;
  stable_like_params par;

  par.dist=dist;
  par.data=(double *)data;
  par.length=length;
  par.nu_c=0;
  par.nu_z=0;

  /* Inital params (with McCulloch) */
  theta=gsl_vector_alloc(4);
  set_expanded(theta,dist->alpha,dist->beta,dist->sigma,dist->mu_1);

  #ifdef DEBUG
    printf("%lf, %lf, %lf, %lf\n",gsl_vector_get (theta, 0),gsl_vector_get (theta, 1),gsl_vector_get (theta, 2),gsl_vector_get (theta, 3));
  #endif

  /* Initial steps */
  ss = gsl_vector_alloc (4);
  gsl_vector_set_all (ss, 0.01);

  /* Cost function to minimize */
  likelihood_func.n = 4; // 4 Dimensions (alpha, beta, sigma, mu_0)
  likelihood_func.f = &stable_minusloglikelihood_whole;
  likelihood_func.params = (void *) (&par);  // Cost function arguments

  /* Minimizer creation */
  T = gsl_multimin_fminimizer_nmsimplex2rand;

  s = gsl_multimin_fminimizer_alloc (T, 4); /* 4 dimensions */

  /* Set cost function, initial guess and initial steps */
  gsl_multimin_fminimizer_set (s, &likelihood_func, theta, ss);

  #ifdef DEBUG
    printf("5\n");
  #endif

  /* Start iterations */
  do
    {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);
      if (status!=GSL_SUCCESS) {
        perror("Minimizer warning:\n");
      }

      size   = gsl_multimin_fminimizer_size (s);
      status = gsl_multimin_test_size (size, 0.002);
      /*
      printf(" %03d\t size = %f a_ = %f  b_ = %f  c_ = %f  m_ = %f f_ = %f \n",iter,size,gsl_vector_get (s->x, 0),gsl_vector_get (s->x, 1),
                                                           gsl_vector_get (s->x, 2),gsl_vector_get (s->x, 3), gsl_multimin_fminimizer_minimum(s));
      */

    } while (status == GSL_CONTINUE && iter < 200);

  if (status!=GSL_SUCCESS)
    {
      perror("Minimizer warning");
    }

  /* Get last estimation */

  gsl_vector_free(theta);
  theta = gsl_multimin_fminimizer_x (s);
  get_original(theta,&a,&b,&c,&m);

  /* Set estimated parameters to the distribution and check if their have valid values*/
  if (stable_setparams(dist,a,b,c,m,0)<0)
   {
    perror("FINAL ESTIMATED PARAMETER ARE NOT VALID\n");
   }

  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);

  return status;
}

int stable_fit_whole(StableDist *dist, const double *data, const unsigned int length)
{
//  double nu_c=0.0,nu_z=0.0;
  int status=0;
//  stable_fit_init(dist,data,length,&nu_c,&nu_z);
//  printf("McCulloch %d sampless: %f %f %f %f\n",length,dist->alpha,dist->beta,dist->sigma,dist->mu_1);

  status = stable_fit_iter_whole(dist,data,length);

  return status;
}

double * load_rand_data(char * filename, int N)
{
  FILE * f_data;
  double * data;
  int i;

  if ((f_data = fopen(filename,"rt")) == NULL)
   {
    perror("Error when opening file with random data");
   }

  data=(double*)malloc(N*sizeof(double));

  for(i=0;i<N;i++)
   {
    if (EOF==fscanf(f_data,"%le\n",data+i))
     {
      perror("Error when reading data");
     }
   }
  return data;
}

int stable_fit_mle(StableDist *dist, const double *data, const unsigned int length) {
  return stable_fit_whole(dist,data,length);
}

int stable_fit_mle2d(StableDist *dist, const double *data, const unsigned int length) {
  return stable_fit(dist,data,length);
}

max_stars_repo_name: TantraLabs/gonum

max_issues_repo_name: TantraLabs/gonum

max_forks_repo_name: TantraLabs/gonum


text:
#ifndef _FUNCTIONS_
#define _FUNCTIONS_

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

#include <cblas.h>
#include <lapacke.h>

#include "structs.c"

#define PATH_LENGTH 32
#define LINE_LENGTH 26000

#define THRESHOLD 0.05

// Declare functions of functions.c file
int is_centroid(int *, int, int);

double rand_gaussian();

int centroids_transposition(conformation **, int *, int *, int, double **);

double **create_distances_array(conformation **, int, int, int);

#endif

ext: h

max_stars_repo_path: Part 3/Proteins/src/functions.h

max_issues_repo_path: Part 3/Proteins/src/functions.h

max_forks_repo_path: Part 3/Proteins/src/functions.h


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

#ifndef __GSL_VECTOR_LONG_DOUBLE_H__
#define __GSL_VECTOR_LONG_DOUBLE_H__

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

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  long double *data;
  gsl_block_long_double *block;
  int owner;
}
gsl_vector_long_double;

typedef struct
{
  gsl_vector_long_double vector;
} _gsl_vector_long_double_view;

typedef _gsl_vector_long_double_view gsl_vector_long_double_view;

typedef struct
{
  gsl_vector_long_double vector;
} _gsl_vector_long_double_const_view;

typedef const _gsl_vector_long_double_const_view gsl_vector_long_double_const_view;

/* Allocation */

GSL_EXPORT gsl_vector_long_double *gsl_vector_long_double_alloc (const size_t n);
GSL_EXPORT gsl_vector_long_double *gsl_vector_long_double_calloc (const size_t n);

GSL_EXPORT gsl_vector_long_double *gsl_vector_long_double_alloc_from_block (gsl_block_long_double * b,
                                                                           const size_t offset,
                                                                           const size_t n,
                                                                           const size_t stride);

GSL_EXPORT gsl_vector_long_double *gsl_vector_long_double_alloc_from_vector (gsl_vector_long_double * v,
                                                                            const size_t offset,
                                                                            const size_t n,
                                                                            const size_t stride);

GSL_EXPORT void gsl_vector_long_double_free (gsl_vector_long_double * v);

/* Views */

GSL_EXPORT
_gsl_vector_long_double_view
gsl_vector_long_double_view_array (long double *v, size_t n);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_vector_long_double_view_array_with_stride (long double *base,
                                         size_t stride,
                                         size_t n);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_vector_long_double_const_view_array (const long double *v, size_t n);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_vector_long_double_const_view_array_with_stride (const long double *base,
                                               size_t stride,
                                               size_t n);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_vector_long_double_subvector (gsl_vector_long_double *v,
                            size_t i,
                            size_t n);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_vector_long_double_subvector_with_stride (gsl_vector_long_double *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_vector_long_double_const_subvector (const gsl_vector_long_double *v,
                                  size_t i,
                                  size_t n);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_vector_long_double_const_subvector_with_stride (const gsl_vector_long_double *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_EXPORT long double gsl_vector_long_double_get (const gsl_vector_long_double * v, const size_t i);
GSL_EXPORT void gsl_vector_long_double_set (gsl_vector_long_double * v, const size_t i, long double x);

GSL_EXPORT long double *gsl_vector_long_double_ptr (gsl_vector_long_double * v, const size_t i);
GSL_EXPORT const long double *gsl_vector_long_double_const_ptr (const gsl_vector_long_double * v, const size_t i);

GSL_EXPORT void gsl_vector_long_double_set_zero (gsl_vector_long_double * v);
GSL_EXPORT void gsl_vector_long_double_set_all (gsl_vector_long_double * v, long double x);
GSL_EXPORT int gsl_vector_long_double_set_basis (gsl_vector_long_double * v, size_t i);

GSL_EXPORT int gsl_vector_long_double_fread (FILE * stream, gsl_vector_long_double * v);
GSL_EXPORT int gsl_vector_long_double_fwrite (FILE * stream, const gsl_vector_long_double * v);
GSL_EXPORT int gsl_vector_long_double_fscanf (FILE * stream, gsl_vector_long_double * v);
GSL_EXPORT int gsl_vector_long_double_fprintf (FILE * stream, const gsl_vector_long_double * v,
                                               const char *format);

GSL_EXPORT int gsl_vector_long_double_memcpy (gsl_vector_long_double * dest, const gsl_vector_long_double * src);

GSL_EXPORT int gsl_vector_long_double_reverse (gsl_vector_long_double * v);

GSL_EXPORT int gsl_vector_long_double_swap (gsl_vector_long_double * v, gsl_vector_long_double * w);
GSL_EXPORT int gsl_vector_long_double_swap_elements (gsl_vector_long_double * v, const size_t i, const size_t j);

GSL_EXPORT long double gsl_vector_long_double_max (const gsl_vector_long_double * v);
GSL_EXPORT long double gsl_vector_long_double_min (const gsl_vector_long_double * v);
GSL_EXPORT void gsl_vector_long_double_minmax (const gsl_vector_long_double * v, long double * min_out, long double * max_out);

GSL_EXPORT size_t gsl_vector_long_double_max_index (const gsl_vector_long_double * v);
GSL_EXPORT size_t gsl_vector_long_double_min_index (const gsl_vector_long_double * v);
GSL_EXPORT void gsl_vector_long_double_minmax_index (const gsl_vector_long_double * v, size_t * imin, size_t * imax);

GSL_EXPORT int gsl_vector_long_double_add (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_EXPORT int gsl_vector_long_double_sub (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_EXPORT int gsl_vector_long_double_mul (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_EXPORT int gsl_vector_long_double_div (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_EXPORT int gsl_vector_long_double_scale (gsl_vector_long_double * a, const double x);
GSL_EXPORT int gsl_vector_long_double_add_constant (gsl_vector_long_double * a, const double x);

GSL_EXPORT int gsl_vector_long_double_isnull (const gsl_vector_long_double * v);

#ifdef HAVE_INLINE

extern inline
long double
gsl_vector_long_double_get (const gsl_vector_long_double * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_long_double_set (gsl_vector_long_double * v, const size_t i, long double x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
long double *
gsl_vector_long_double_ptr (gsl_vector_long_double * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (long double *) (v->data + i * v->stride);
}

extern inline
const long double *
gsl_vector_long_double_const_ptr (const gsl_vector_long_double * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const long double *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_LONG_DOUBLE_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#ifndef _BLAS_UTIL_H_
#define _BLAS_UTIL_H_

#ifdef __BLAS_LEGACY__
#include <math.h>
#include "cblas.h"
//#include <lapacke.h> //! hasn't used
#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))

#elif defined __APPLE__
#include <Accelerate/Accelerate.h>

#elif defined __USE_MKL__
#include <mkl.h>

#endif

#endif /* _BLAS_UTIL_H_ */


text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_BLAS3_MATRIX_PRODUCT_H
#define SCILIB_LINALG_BLAS3_MATRIX_PRODUCT_H

#ifdef USE_MKL
#include <mkl.h>
#else
#include <cblas.h>
#endif

#include <scilib/mdarray.h>
#include <scilib/linalg_impl/lapack_types.h>
#include <experimental/mdspan>
#include <complex>
#include <type_traits>

namespace Sci {
namespace Linalg {

namespace stdex = std::experimental;

template <class T_a,
          stdex::extents<>::size_type nrows_a,
          stdex::extents<>::size_type ncols_a,
          class Layout_a,
          class Accessor_a,
          class T_b,
          stdex::extents<>::size_type nrows_b,
          stdex::extents<>::size_type ncols_b,
          class Layout_b,
          class Accessor_b,
          class T_c,
          stdex::extents<>::size_type nrows_c,
          stdex::extents<>::size_type ncols_c,
          class Layout_c,
          class Accessor_c>
    requires(!std::is_const_v<T_c>)
inline void matrix_product(
    stdex::mdspan<T_a, stdex::extents<nrows_a, ncols_a>, Layout_a, Accessor_a>
        a,
    stdex::mdspan<T_b, stdex::extents<nrows_b, ncols_b>, Layout_b, Accessor_b>
        b,
    stdex::mdspan<T_c, stdex::extents<nrows_c, ncols_c>, Layout_c, Accessor_c>
        c)
{
    static_assert(a.static_extent(1) == b.static_extent(0));

    using size_type = stdex::extents<>::size_type;

    const size_type n = a.extent(0);
    const size_type m = a.extent(1);
    const size_type p = b.extent(1);

    for (size_type i = 0; i < n; ++i) {
        for (size_type j = 0; j < p; ++j) {
            c(i, j) = T_c{0};
            for (size_type k = 0; k < m; ++k) {
                c(i, j) += a(i, k) * b(k, j);
            }
        }
    }
}

template <class Layout>
inline void matrix_product(Sci::Matrix_view<double, Layout> a,
                           Sci::Matrix_view<double, Layout> b,
                           Sci::Matrix_view<double, Layout> c)
{
    constexpr double alpha = 1.0;
    constexpr double beta = 0.0;

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(b.extent(1));
    const BLAS_INT k = static_cast<BLAS_INT>(a.extent(1));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = k;
    BLAS_INT ldb = n;
    BLAS_INT ldc = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
        ldb = k;
        ldc = m;
    }
    cblas_dgemm(matrix_layout, CblasNoTrans, CblasNoTrans, m, n, k, alpha,
                a.data(), lda, b.data(), ldb, beta, c.data(), ldc);
}

template <class Layout>
inline void matrix_product(Sci::Matrix_view<const double, Layout> a,
                           Sci::Matrix_view<const double, Layout> b,
                           Sci::Matrix_view<double, Layout> c)
{
    constexpr double alpha = 1.0;
    constexpr double beta = 0.0;

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(b.extent(1));
    const BLAS_INT k = static_cast<BLAS_INT>(a.extent(1));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = k;
    BLAS_INT ldb = n;
    BLAS_INT ldc = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
        ldb = k;
        ldc = m;
    }
    cblas_dgemm(matrix_layout, CblasNoTrans, CblasNoTrans, m, n, k, alpha,
                a.data(), lda, b.data(), ldb, beta, c.data(), ldc);
}

#ifdef USE_MKL
template <class Layout>
inline void matrix_product(Sci::Matrix_view<std::complex<double>, Layout> a,
                           Sci::Matrix_view<std::complex<double>, Layout> b,
                           Sci::Matrix_view<std::complex<double>, Layout> c)
{
    constexpr std::complex<double> alpha = {1.0, 0.0};
    constexpr std::complex<double> beta = {0.0, 0.0};

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(b.extent(1));
    const BLAS_INT k = static_cast<BLAS_INT>(a.extent(1));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = k;
    BLAS_INT ldb = n;
    BLAS_INT ldc = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
        ldb = k;
        ldc = m;
    }
    cblas_zgemm(matrix_layout, CblasNoTrans, CblasNoTrans, m, n, k, &alpha,
                a.data(), lda, b.data(), ldb, &beta, c.data(), ldc);
}

template <class Layout>
inline void
matrix_product(Sci::Matrix_view<const std::complex<double>, Layout> a,
               Sci::Matrix_view<const std::complex<double>, Layout> b,
               Sci::Matrix_view<std::complex<double>, Layout> c)
{
    constexpr std::complex<double> alpha = {1.0, 0.0};
    constexpr std::complex<double> beta = {0.0, 0.0};

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(b.extent(1));
    const BLAS_INT k = static_cast<BLAS_INT>(a.extent(1));

    auto matrix_layout = CblasRowMajor;
    BLAS_INT lda = k;
    BLAS_INT ldb = n;
    BLAS_INT ldc = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = CblasColMajor;
        lda = m;
        ldb = k;
        ldc = m;
    }
    cblas_zgemm(matrix_layout, CblasNoTrans, CblasNoTrans, m, n, k, &alpha,
                a.data(), lda, b.data(), ldb, &beta, c.data(), ldc);
}
#endif

template <class T, class Layout, class Allocator>
inline Sci::Matrix<T, Layout, Allocator>
matrix_product(const Sci::Matrix<T, Layout, Allocator>& a,
               const Sci::Matrix<T, Layout, Allocator>& b)
{
    using size_type = stdex::extents<>::size_type;

    const size_type n = a.extent(0);
    const size_type p = b.extent(1);

    Sci::Matrix<T, Layout, Allocator> res(n, p);
    matrix_product(a.view(), b.view(), res.view());
    return res;
}

template <class T, class Layout, class Allocator>
inline void matrix_product(const Sci::Matrix<T, Layout, Allocator>& a,
                           const Sci::Matrix<T, Layout, Allocator>& b,
                           Sci::Matrix<T, Layout, Allocator>& c)
{
    matrix_product(a.view(), b.view(), c.view());
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_BLAS3_MATRIX_PRODUCT_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


text:
/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>

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

#include "os-features.h"
#include "sip-utils.h"
#include "gslutils.h"
#include "starutil.h"
#include "mathutil.h"
#include "errors.h"
#include "log.h"

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

void tan_rotate(const tan_t* tanin, tan_t* tanout, double angle) {
    double s,c;
    double newcd[4];
    memmove(tanout, tanin, sizeof(tan_t));
    s = sin(deg2rad(angle));
    c = cos(deg2rad(angle));
    newcd[0] = c*tanin->cd[0][0] + s*tanin->cd[1][0];
    newcd[1] = c*tanin->cd[0][1] + s*tanin->cd[1][1];
    newcd[2] = -s*tanin->cd[0][0] + c*tanin->cd[1][0];
    newcd[3] = -s*tanin->cd[0][1] + c*tanin->cd[1][1];
    tanout->cd[0][0] = newcd[0];
    tanout->cd[0][1] = newcd[1];
    tanout->cd[1][0] = newcd[2];
    tanout->cd[1][1] = newcd[3];
}

int sip_ensure_inverse_polynomials(sip_t* sip) {
    if ((sip->a_order == 0 && sip->b_order == 0) ||
        (sip->ap_order > 0  && sip->bp_order > 0)) {
        return 0;
    }
    sip->ap_order = sip->bp_order = MAX(sip->a_order, sip->b_order) + 1;
    return sip_compute_inverse_polynomials(sip, 0, 0, 0, 0, 0, 0);
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    return 0;
}

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

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

anbool sip_is_inside_image(const sip_t* wcs, double ra, double dec) {
    double x,y;
    if (!sip_radec2pixelxy(wcs, ra, dec, &x, &y))
        return FALSE;
    return sip_pixel_is_inside_image(wcs, x, y);
}

anbool tan_is_inside_image(const tan_t* wcs, double ra, double dec) {
    double x,y;
    if (!tan_radec2pixelxy(wcs, ra, dec, &x, &y))
        return FALSE;
    return tan_pixel_is_inside_image(wcs, x, y);
}

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

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

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

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

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

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

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

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

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

    *p_Ngood = Ngood;

    return inds;
}

void sip_get_radec_center(const sip_t* wcs,
                          double* p_ra, double* p_dec) {
    double px = wcs_pixel_center_for_size(wcs->wcstan.imagew);
    double py = wcs_pixel_center_for_size(wcs->wcstan.imageh);
    sip_pixelxy2radec(wcs, px, py, p_ra, p_dec);
}

void tan_get_radec_center(const tan_t* wcs,
                          double* p_ra, double* p_dec) {
    double px = wcs_pixel_center_for_size(wcs->imagew);
    double py = wcs_pixel_center_for_size(wcs->imageh);
    tan_pixelxy2radec(wcs, px, py, p_ra, p_dec);
}

double sip_get_radius_deg(const sip_t* wcs) {
    return arcsec2deg(sip_pixel_scale(wcs) * hypot(wcs->wcstan.imagew, wcs->wcstan.imageh)/2.0);
}

double tan_get_radius_deg(const tan_t* wcs) {
    return arcsec2deg(tan_pixel_scale(wcs) * hypot(wcs->imagew, wcs->imageh)/2.0);
}

void sip_get_radec_center_hms(const sip_t* wcs,
                              int* rah, int* ram, double* ras,
                              int* decsign, int* decd, int* decm, double* decs) {
    double ra, dec;
    sip_get_radec_center(wcs, &ra, &dec);
    ra2hms(ra, rah, ram, ras);
    dec2dms(dec, decsign, decd, decm, decs);
}

void sip_get_radec_center_hms_string(const sip_t* wcs,
                                     char* rastr, char* decstr) {
    double ra, dec;
    sip_get_radec_center(wcs, &ra, &dec);
    ra2hmsstring(ra, rastr);
    dec2dmsstring(dec, decstr);
}

void sip_get_field_size(const sip_t* wcs,
                        double* pw, double* ph,
                        char** units) {
    double minx = 0.5;
    double maxx = (wcs->wcstan.imagew + 0.5);
    double midx = (minx + maxx) / 2.0;
    double miny = 0.5;
    double maxy = (wcs->wcstan.imageh + 0.5);
    double midy = (miny + maxy) / 2.0;
    double ra1, dec1, ra2, dec2, ra3, dec3;
    double w, h;

    // measure width through the middle
    sip_pixelxy2radec(wcs, minx, midy, &ra1, &dec1);
    sip_pixelxy2radec(wcs, midx, midy, &ra2, &dec2);
    sip_pixelxy2radec(wcs, maxx, midy, &ra3, &dec3);
    w = arcsec_between_radecdeg(ra1, dec1, ra2, dec2) +
        arcsec_between_radecdeg(ra2, dec2, ra3, dec3);
    // measure height through the middle
    sip_pixelxy2radec(wcs, midx, miny, &ra1, &dec1);
    sip_pixelxy2radec(wcs, midx, midy, &ra2, &dec2);
    sip_pixelxy2radec(wcs, midx, maxy, &ra3, &dec3);
    h = arcsec_between_radecdeg(ra1, dec1, ra2, dec2) +
        arcsec_between_radecdeg(ra2, dec2, ra3, dec3);

    if (MIN(w, h) < 60.0) {
        *units = "arcseconds";
        *pw = w;
        *ph = h;
    } else if (MIN(w, h) < 3600.0) {
        *units = "arcminutes";
        *pw = w / 60.0;
        *ph = h / 60.0;
    } else {
        *units = "degrees";
        *pw = w / 3600.0;
        *ph = h / 3600.0;
    }
}

void sip_walk_image_boundary(const sip_t* wcs, double stepsize,
                             void (*callback)(const sip_t* wcs, double x, double y, double ra, double dec, void* token),
                             void* token) {
    int i, side;
    // Walk the perimeter of the image in steps of stepsize pixels
    double W = wcs->wcstan.imagew;
    double H = wcs->wcstan.imageh;
    {
        double Xmin = 0.5;
        double Xmax = W + 0.5;
        double Ymin = 0.5;
        double Ymax = H + 0.5;
        double offsetx[] = { Xmin, Xmax, Xmax, Xmin };
        double offsety[] = { Ymin, Ymin, Ymax, Ymax };
        double stepx[] = { +stepsize, 0, -stepsize, 0 };
        double stepy[] = { 0, +stepsize, 0, -stepsize };
        int Nsteps[] = { ceil(W/stepsize), ceil(H/stepsize), ceil(W/stepsize), ceil(H/stepsize) };

        for (side=0; side<4; side++) {
            for (i=0; i<Nsteps[side]; i++) {
                double ra, dec;
                double x, y;
                x = MIN(Xmax, MAX(Xmin, offsetx[side] + i * stepx[side]));
                y = MIN(Ymax, MAX(Ymin, offsety[side] + i * stepy[side]));
                sip_pixelxy2radec(wcs, x, y, &ra, &dec);
                callback(wcs, x, y, ra, dec, token);
            }
        }
    }
}

struct radecbounds {
    double rac, decc;
    double ramin, ramax, decmin, decmax;
};

static void radec_bounds_callback(const sip_t* wcs, double x, double y, double ra, double dec, void* token) {
    struct radecbounds* b = token;
    b->decmin = MIN(b->decmin, dec);
    b->decmax = MAX(b->decmax, dec);
    if (ra - b->rac > 180)
        // wrap-around: racenter < 180, ra has gone < 0 but been wrapped around to > 180
        ra -= 360;
    if (b->rac - ra > 180)
        // wrap-around: racenter > 180, ra has gone > 360 but wrapped around to > 0.
        ra += 360;

    b->ramin = MIN(b->ramin, ra);
    b->ramax = MAX(b->ramax, ra);
}

void sip_get_radec_bounds(const sip_t* wcs, int stepsize,
                          double* pramin, double* pramax,
                          double* pdecmin, double* pdecmax) {
    struct radecbounds b;

    sip_get_radec_center(wcs, &(b.rac), &(b.decc));
    b.ramin  = b.ramax = b.rac;
    b.decmin = b.decmax = b.decc;
    sip_walk_image_boundary(wcs, stepsize, radec_bounds_callback, &b);

    // Check for poles...
    // north pole
    if (sip_is_inside_image(wcs, 0, 90)) {
        b.ramin = 0;
        b.ramax = 360;
        b.decmax = 90;
    }
    if (sip_is_inside_image(wcs, 0, -90)) {
        b.ramin = 0;
        b.ramax = 360;
        b.decmin = -90;
    }

    if (pramin) *pramin = b.ramin;
    if (pramax) *pramax = b.ramax;
    if (pdecmin) *pdecmin = b.decmin;
    if (pdecmax) *pdecmax = b.decmax;
}

void sip_shift(const sip_t* sipin, sip_t* sipout,
               double xlo, double xhi, double ylo, double yhi) {
    memmove(sipout, sipin, sizeof(sip_t));
    tan_transform(&(sipin->wcstan), &(sipout->wcstan),
                  xlo, xhi, ylo, yhi, 1.0);
}

void tan_transform(const tan_t* tanin, tan_t* tanout,
                   double xlo, double xhi, double ylo, double yhi,
                   double scale) {
    memmove(tanout, tanin, sizeof(tan_t));
    tanout->imagew = (xhi - xlo + 1) * scale;
    tanout->imageh = (yhi - ylo + 1) * scale;
    tanout->crpix[0] = (tanout->crpix[0] - (xlo - 1)) * scale;
    tanout->crpix[1] = (tanout->crpix[1] - (ylo - 1)) * scale;
    tanout->cd[0][0] /= scale;
    tanout->cd[0][1] /= scale;
    tanout->cd[1][0] /= scale;
    tanout->cd[1][1] /= scale;
}

void tan_scale(const tan_t* tanin, tan_t* tanout,
               double scale) {
    memmove(tanout, tanin, sizeof(tan_t));
    tanout->imagew *= scale;
    tanout->imageh *= scale;

    tanout->crpix[0] = 0.5 + scale * (tanin->crpix[0] - 0.5);
    tanout->crpix[1] = 0.5 + scale * (tanin->crpix[1] - 0.5);
    tanout->cd[0][0] /= scale;
    tanout->cd[0][1] /= scale;
    tanout->cd[1][0] /= scale;
    tanout->cd[1][1] /= scale;
}

void sip_scale(const sip_t* wcsin, sip_t* wcsout,
               double scale) {
    int i, j;
    memmove(wcsout, wcsin, sizeof(sip_t));
    tan_scale(&(wcsin->wcstan), &(wcsout->wcstan), scale);
    for (i=0; i<=wcsin->a_order; i++) {
        for (j=0; j<=wcsin->a_order; j++) {
            if (i + j > wcsin->a_order)
                continue;
            wcsout->a[i][j] *= pow(scale, 1 - (i+j));
        }
    }
    for (i=0; i<=wcsin->b_order; i++) {
        for (j=0; j<=wcsin->b_order; j++) {
            if (i + j > wcsin->b_order)
                continue;
            wcsout->b[i][j] *= pow(scale, 1 - (i+j));
        }
    }
    for (i=0; i<=wcsin->ap_order; i++) {
        for (j=0; j<=wcsin->ap_order; j++) {
            if (i + j > wcsin->ap_order)
                continue;
            wcsout->ap[i][j] *= pow(scale, 1 - (i+j));
        }
    }
    for (i=0; i<=wcsin->bp_order; i++) {
        for (j=0; j<=wcsin->bp_order; j++) {
            if (i + j > wcsin->bp_order)
                continue;
            wcsout->bp[i][j] *= pow(scale, 1 - (i+j));
        }
    }
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
#include "tools.h"

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

gsl_rng **_random = NULL;

double tool_get_time()
{
    struct timeval t;
    gettimeofday(&t, NULL);
    return (double)t.tv_sec + (double)t.tv_usec * 1e-6;
}

double tool_cut(double x)
{
    if (x <= 0.0)
        return 0.0;
    if (x >= 1.0)
        return 1.0;
    return x;
}

void print_array(double c[], int l)
{
    for (int i = 0; i < l - 1; i++)
    {
        printf("%.3f ", c[i]);
    }
    printf("%.3f\n", c[l - 1]);
}

void print_array_b(bool c[], int l)
{
    for (int i = 0; i < l - 1; i++)
    {
        printf("%d,", c[i] ? 1 : 0);
    }
    printf("%d\n", c[l - 1] ? 1 : 0);
}

void print_array_i(int c[], int l)
{
    for (int i = 0; i < l - 1; i++)
    {
        printf("%d,", c[i]);
    }
    printf("%d\n", c[l - 1]);
}

void print_array_ui(unsigned int c[], int l)
{
    for (int i = 0; i < l - 1; i++)
    {
        printf("%u,", c[i]);
    }
    printf("%u\n", c[l - 1]);
}

void print_array_2d(bool **c, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
    {
        for (int y = 0; y < ylen; y++)
        {
            printf("%d,", c[x][y]);
        }
        puts("");
    }
}

void print_array_2d_d(double **c, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
    {
        for (int y = 0; y < ylen; y++)
        {
            printf("%.2f,", c[x][y]);
        }
        puts("");
    }
}

void tool_free_array_d(double **a, int xlen)
{
    for (int x = 0; x < xlen; x++)
    {
        free(a[x]);
    }
    free(a);
}

void tool_free_array_b(bool **a, int xlen)
{
    for (int x = 0; x < xlen; x++)
    {
        free(a[x]);
    }
    free(a);
}

void tool_free_array_ui(unsigned int **a, int xlen)
{
    for (int x = 0; x < xlen; x++)
    {
        free(a[x]);
    }
    free(a);
}

bool **tool_alloc_array_b(int xlen, int ylen)
{

    bool **a = (bool **)malloc(xlen * sizeof(bool *));
    for (int x = 0; x < xlen; x++)
    {
        a[x] = (bool *)calloc((size_t)ylen, sizeof(bool));
    }
    return a;
}

unsigned int **tool_alloc_array_ui(int xlen, int ylen)
{

    unsigned int **a = (unsigned int **)malloc(xlen * sizeof(unsigned int *));
    for (int x = 0; x < xlen; x++)
    {
        a[x] = (unsigned int *)calloc((size_t)ylen, sizeof(unsigned int));
    }
    return a;
}

double **tool_alloc_array_d(int xlen, int ylen)
{
    double **a = (double **)malloc(xlen * sizeof(double *));
    for (int x = 0; x < xlen; x++)
    {
        a[x] = (double *)calloc((size_t)ylen, sizeof(double));
    }
    return a;
}

void tool_array_bool_to_double(double *res, bool *src, int len)
{
    for (int i = 0; i < len; i++)
    {
        res[i] = src[i] ? 1.0 : 0.0;
    }
}

void tool_round(double **a, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
        for (int y = 0; y < ylen; y++)
            a[x][y] = round(a[x][y]);
}

static int tool_current_thread()
{
#ifdef _USE_OPENMP
    return omp_get_thread_num();
#else
    return 0;
#endif
}

double tool_rand()
{
    return gsl_rng_uniform(_random[tool_current_thread()]);
}

/*
 * Returns random int from [min, max-1]
 */
int tool_rand_i(int min, int max)
{
    return min + gsl_rng_uniform_int(_random[tool_current_thread()], max - min);
}

void tool_fill_random1(double *a, int xlen, bool deterministic)
{
    for (int x = 0; x < xlen; x++)
    {
        a[x] = tool_rand();
        if (deterministic)
        {
            a[x] = a[x] > 0.5 ? 1 : 0;
        }
    }
}

void tool_fill_random1ui(unsigned int *a, int xlen, int min, int max)
{
    for (int x = 0; x < xlen; x++)
    {
        a[x] = (unsigned int)tool_rand_i(min, max);
    }
}

void tool_fill_random1i(int *a, int xlen, int min, int max)
{
    for (int x = 0; x < xlen; x++)
    {
        a[x] = tool_rand_i(min, max);
    }
}

void tool_fill_random1b(bool *a, int xlen)
{
    for (int x = 0; x < xlen; x++)
    {
        a[x] = tool_rand() < 0.5;
    }
}

void tool_fill_random(double **a, int xlen, int ylen, bool deterministic)
{
    for (int x = 0; x < xlen; x++)
    {
        for (int y = 0; y < ylen; y++)
        {
            a[x][y] = tool_rand();
            if (deterministic)
            {
                a[x][y] = a[x][y] > 0.5 ? 1 : 0;
            }
        }
    }
}

void tool_fill_random_d(double **a, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
    {
        for (int y = 0; y < ylen; y++)
        {
            a[x][y] = tool_rand(); // < 0.5 ? 1 : 0;
        }
    }
}

void tool_fill_random_b(bool **a, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
    {
        double p = tool_rand();
        for (int y = 0; y < ylen; y++)
        {
            a[x][y] = tool_rand() <= p;
        }
    }
}

void tool_fill_random_ui(unsigned int **a, int xlen, int ylen, int min, int max)
{
    for (int x = 0; x < xlen; x++)
    {
        for (int y = 0; y < ylen; y++)
        {
            a[x][y] = (unsigned int)tool_rand_i(min, max);
        }
    }
}

double tool_sum(double **a, int xlen, int ylen)
{
    double res = 0;
    for (int x = 0; x < xlen; x++)
        for (int y = 0; y < ylen; y++)
            res += a[x][y];
    return res;
}

int tool_foldedSize(int s, int r)
{
    return s + 2 * r;
}

void tool_fold(bool *f_res, bool *u_src, int u_len, int r)
{
    memcpy(f_res + r, u_src, (size_t)u_len * sizeof(bool));
    memcpy(f_res, u_src + u_len - r, (size_t)r * sizeof(bool));
    memcpy(f_res + r + u_len, u_src, (size_t)r * sizeof(bool));
}

void tool_fold_d(double *f_res, double *u_src, int u_len, int r)
{
    memcpy(f_res + r, u_src, (size_t)u_len * sizeof(double));
    memcpy(f_res, u_src + u_len - r, (size_t)r * sizeof(double));
    memcpy(f_res + r + u_len, u_src, (size_t)r * sizeof(double));
}

void tool_unfold(bool u_res[], bool f_src[], int f_len, int r)
{
    memcpy(u_res, f_src + r, (size_t)(f_len - 2 * r) * sizeof(bool));
}

void tool_unfold_d(double u_res[], double f_src[], int f_len, int r)
{
    memcpy(u_res, f_src + r, (size_t)(f_len - 2 * r) * sizeof(double));
}

void tool_add(double **res, bool **src, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
        for (int y = 0; y < ylen; y++)
            res[x][y] += src[x][y] ? 1.0 : 0.0;
}

void tool_mult(double **a, double c, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
        for (int y = 0; y < ylen; y++)
            a[x][y] *= c;
}

void tool_abs_subst(double **res, bool **a, double **b, int xlen, int ylen)
{
    for (int x = 0; x < xlen; x++)
        for (int y = 0; y < ylen; y++)
        {
            res[x][y] = fabs((a[x][y] ? 1.0 : 0.0) - b[x][y]);
        }
}

double tool_1d_max(double *a, int len)
{
    double max = 0.0;
    for (int i = 0; i < len; i++)
    {
        if (a[i] > max)
            max = a[i];
    }
    return max;
}

int tool_differences(bool **a, bool **b, int xlen, int ylen)
{
    int result = 0;
    for (int x = 0; x < xlen; x++)
    {
        for (int y = 0; y < ylen; y++)
        {
            if (a[x][y] && !b[x][y])
                result++;
            if (!a[x][y] && b[x][y])
                result++;
        }
    }
    return result;
}

static int tool_num_threads()
{
#ifdef _USE_OPENMP
    return omp_get_num_threads();
#else
    return 1;
#endif
}

void tool_fake_random()
{
    int threads = 0;
#ifdef _USE_OPENMP
#pragma omp parallel
#pragma omp master
#endif
    {
        threads = tool_num_threads();
        const gsl_rng_type *T;
        gsl_rng_env_setup();
        T = gsl_rng_default;
        _random = (gsl_rng **)malloc(threads * sizeof(gsl_rng *));
        for (int i = 0; i < threads; i++)
        {
            _random[i] = gsl_rng_alloc(T);
            gsl_rng_set(_random[i], 1234567890);
        }
    }
}

void tool_init_random()
{
    int threads = 0;
#ifdef _USE_OPENMP
#pragma omp parallel
#pragma omp master
#endif
    {
        if (_random == NULL)
        {
            threads = tool_num_threads();

#ifndef _USE_RAND
#ifdef __APPLE__
            srandomdev();
#else
            srandom((unsigned int)tool_get_time());
#endif
#else
            srand((unsigned int)tool_get_time());
#endif

            const gsl_rng_type *T;

            gsl_rng_env_setup();

            T = gsl_rng_default;

            _random = (gsl_rng **)malloc(threads * sizeof(gsl_rng *));
            for (int i = 0; i < threads; i++)
            {
                _random[i] = gsl_rng_alloc(T);
                unsigned long int seed;
#ifdef _USE_RAND
                seed = (unsigned long int)rand();
#else
                seed = (unsigned long int)random();
#endif
                gsl_rng_set(_random[i], seed);
            }
        }
    }
}

double tool_rand_gauss(double sigma)
{
    return gsl_ran_gaussian_ziggurat(_random[tool_current_thread()], sigma);
}

int tool_cap(int a, int min, int max)
{
    if (a > max)
        return max;
    if (a < min)
        return min;
    return a;
}

void tool_range(int *result, int start, int len)
{
    for (int i = start; i < start + len; i++)
    {
        result[i - start] = i;
    }
}

int *tool_alloc_range(int start, int len)
{
    int *a = (int *)malloc(sizeof(int) * len);
    tool_range(a, start, len);
    return a;
}

#ifdef __APPLE__
int _compr(void *f, const void *a, const void *b)
{
#else
int _compr(const void *a, const void *b, void *f)
{
#endif
    double *fitness = (double *)f;
    int i1 = (int)(*(int *)a);
    int i2 = (int)(*(int *)b);
    double diff = fitness[i1] - fitness[i2];
    if (diff == 0.0)
        return 0;
    if (diff > 0.0)
        return 1;
    return -1;
}

int *tool_sort_fitness_index(int *index, double *fitness, int count)
{
#ifdef __APPLE__
    qsort_r(index, count, sizeof(int), fitness, _compr);
#else
    qsort_r(index, count, sizeof(int), _compr, fitness);
#endif
    return index;
}

#ifdef __APPLE__
int _compr_int(void *f, const void *a, const void *b)
{
#else
int _compr_int(const void *a, const void *b, void *f)
{
#endif
    int *fitness = (int *)f;
    int i1 = (int)(*(int *)a);
    int i2 = (int)(*(int *)b);
    int diff = fitness[i1] - fitness[i2];
    if (diff == 0)
        return 0;
    if (diff > 0)
        return 1;
    return -1;
}

int *tool_sort_index_int(int *index, int *array, int count)
{
#ifdef __APPLE__
    qsort_r(index, count, sizeof(int), array, _compr_int);
#else
    qsort_r(index, count, sizeof(int), _compr_int, array);
#endif
    return index;
}

bool tool_row_cmp(bool *r1, bool *r2, int len)
{
    for (int i = 0; i < len; i++)
    {
        if (r1[i] != r2[i])
            return false;
    }
    return true;
}

double tool_avg(double *a, int len)
{
    double result = 0;
    for (int i = 0; i < len; i++)
        result += a[i];
    return result / len;
}

double tool_max(double *a, int len)
{
    double result = a[0];
    for (int i = 1; i < len; i++)
    {
        if (result < a[i])
            result = a[i];
    }
    return result;
}

int tool_max_i(int *a, int len)
{
    int result = a[0];
    for (int i = 1; i < len; i++)
    {
        if (result < a[i])
            result = a[i];
    }
    return result;
}

double tool_min(double *a, int len)
{
    double result = a[0];
    for (int i = 1; i < len; i++)
    {
        if (result > a[i])
            result = a[i];
    }
    return result;
}

int tool_min_i(int *a, int len)
{
    int result = a[0];
    for (int i = 1; i < len; i++)
    {
        if (result > a[i])
            result = a[i];
    }
    return result;
}

double tool_max_ix(double *a, int len, int *id)
{
    double result = a[0];
    (*id) = 0;
    for (int i = 1; i < len; i++)
    {
        if (result < a[i])
        {
            result = a[i];
            (*id) = i;
        }
    }
    return result;
}

double tool_min_ix(double *a, int len, int *id)
{
    double result = a[0];
    (*id) = 0;
    for (int i = 1; i < len; i++)
    {
        if (result > a[i])
        {
            result = a[i];
            (*id) = i;
        }
    }
    return result;
}

double tool_1d_sum(double *a, int len)
{
    double result = 0;
    for (int i = 0; i < len; i++)
    {
        result += a[i];
    }
    return result;
}

bool tool_compar_b(bool *a, bool *b, int len)
{
    for (int i = 0; i < len; i++)
    {
        if (a[i] != b[i])
            return false;
    }
    return true;
}

bool tool_compar_d(double *a, double *b, int len)
{
    for (int i = 0; i < len; i++)
    {
        if (a[i] != b[i])
            return false;
    }
    return true;
}

max_stars_repo_name: houp/identify

max_issues_repo_name: houp/identify

max_forks_repo_name: houp/identify


text:
// xll_roots.h - GSL 1-dim root finding
// http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html#One-dimensional-Root_002dFinding
#pragma once
#include <functional>
#include <memory>
#include <stdexcept>
#include <tuple>
#include "gsl/gsl_errno.h"
#include "gsl/gsl_roots.h"
#include "xll_math.h"

namespace gsl {

namespace root {

	// root bracketing solvers
	class fsolver {
		gsl_root_fsolver* s;
		gsl::function f;
	public:
		explicit fsolver(const gsl_root_fsolver_type * type)
			: s{gsl_root_fsolver_alloc(type)}
		{ }
		fsolver(const fsolver&) = delete;
		fsolver& operator=(const fsolver&) = delete;
		~fsolver()
		{
			gsl_root_fsolver_free(s);
		}

		// needed for gsl_root_fsolver_* routines
		gsl_root_fsolver* get() const
		{
			return s;
		}
		// syntactic sugar
		operator gsl_root_fsolver*() const
		{
			return get();
		}

		int set(const std::function<double(double)>& f_, double lo, double hi)
		{
			f = f_;

			return gsl_root_fsolver_set(s, &f, lo, hi);
		}

		// forward to gsl_root_fsolver_* functions
		int iterate()
		{
			return gsl_root_fsolver_iterate(s);
		}
		double x_lower() const
		{
			return gsl_root_fsolver_x_lower(s);
		}
		double x_upper() const
		{
			return gsl_root_fsolver_x_upper(s);
		}
		double root() const
		{
			return gsl_root_fsolver_root(s);
		}

		// specify convergence condition
		double solve(const std::function<bool(const fsolver&)>& done)
		{
			while (GSL_SUCCESS == iterate()) {
				if (done(*this))
					break;
			}

			return root();
		}
	};

	// convergence helper functions
	inline auto test_interval(double epsabs, double epsrel)
	{
		return [epsabs,epsrel](const fsolver& s) {
			return GSL_SUCCESS == gsl_root_test_interval(s.x_lower(), s.x_upper(), epsabs, epsrel);
		};
	}

	// root finding using derivatives
	class fdfsolver {
		gsl_root_fdfsolver* s;

		// function and its derivative
		gsl::function_fdf FdF;
	public:
		fdfsolver(const gsl_root_fdfsolver_type * type)
			: s{gsl_root_fdfsolver_alloc(type)}
		{ }
		fdfsolver(const fdfsolver&) = delete;
		fdfsolver& operator=(const fdfsolver&) = delete;
		~fdfsolver()
		{
			gsl_root_fdfsolver_free(s);
		}
		// needed for gsl_root_fdfsolver_* routines
		gsl_root_fdfsolver* get() const
		{
			return s;
		}
		// syntactic sugar
		operator gsl_root_fdfsolver*() const
		{
			return s;
		}

		int set(const std::function<double(double)>& f, const std::function<double(double)>& df, double x0)
		{
			FdF = gsl::function_fdf(f, df);

			return gsl_root_fdfsolver_set(s, &FdF, x0);
		}

		// forward to gsl_root_fdfsolver_* functions
		int iterate()
		{
			return gsl_root_fdfsolver_iterate(s);
		}
		double root() const
		{
			return gsl_root_fdfsolver_root(s);
		}

		// specify convergence condition given previous root
		double solve(const std::function<bool(const fdfsolver&,double)>& done)
		{
			double x = root();
			while (GSL_SUCCESS == iterate()) {
				if (done(*this, x))
					break;
				x = root();
			}

			return root();
		}
	};

	// convergence helper functions
	inline auto test_delta(double epsabs, double epsrel)
	{
		return [epsabs,epsrel](const fdfsolver& s, double x) -> bool {
			return GSL_SUCCESS == gsl_root_test_delta(x, s.root(), epsabs, epsrel);
		};
	}

	} // root
} // gsl

#ifdef _DEBUG
#include <cassert>
#include <vector>
//#include <gsl/gsl_poly.h>

// http://www.gnu.org/software/gsl/manual/html_node/Root-Finding-Examples.html#Root-Finding-Examples

// a x^2 + b x + c
struct quadratic_params {
	double a, b, c;
};

// This is bad. Why???
inline double quadratic(double x, void *params)
{
	struct quadratic_params *p
		= (struct quadratic_params *) params;

	double a = p->a;
	double b = p->b;
	double c = p->c;

	return (a * x + b) * x + c;
}

inline double quadratic_deriv(double x, void *params)
{
	struct quadratic_params *p
		= (struct quadratic_params *) params;

	double a = p->a;
	double b = p->b;
//	double c = p->c;

	return 2.0 * a * x + b;
}

inline void quadratic_fdf(double x, void *params, double *y, double *dy)
{
	struct quadratic_params *p
		= (struct quadratic_params *) params;

	double a = p->a;
	double b = p->b;
	double c = p->c;

	*y = (a * x + b) * x + c;
	*dy = 2.0 * a * x + b;
}

inline void test_gsl_root_fsolver()
{
	// GSL example
	{
		gsl_root_fsolver* s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);

		quadratic_params params{1.,0.,-5.}; // x^2 - 5
		gsl_function F;
		F.function = quadratic;
		F.params = &params;
		double x_lo = 0.0, x_hi = 5.0;
		double epsabs = 0, epsrel = 1e-5;

		gsl_root_fsolver_set(s, &F, x_lo, x_hi);
		while (GSL_SUCCESS == gsl_root_fsolver_iterate(s)) {
			x_lo = gsl_root_fsolver_x_lower(s);
			x_hi = gsl_root_fsolver_x_upper(s);

			// |x_lo - x_hi| < epsabs + epsrel * min(|x_lo|, |x_hi|)
			if (GSL_SUCCESS == gsl_root_test_interval(x_lo, x_hi, epsabs, epsrel))
				break;
		}
		double root = gsl_root_fsolver_root(s);
		double sqrt5 = sqrt(5.);
		assert (fabs(root - sqrt5) < sqrt5*epsrel);

		gsl_root_fsolver_free(s);
	}
	// root_fsolver class with function<double(double)>
	{
		gsl::root::fsolver s(gsl_root_fsolver_brent);

		std::vector<double> params{-5,0,1}; // -5 + x^2
		auto F = [params](double x) { return params[0] + x*(params[1] + x*params[2]); };
		gsl::function F_(F);
		assert (F_(0) == -5);
		assert (F_.call(0) == -5);
		double x_lo = 0.0, x_hi = 5.0, epsrel = 1e-6;
		s.set(F, x_lo, x_hi);

		double root = s.solve(gsl::root::test_interval(0, epsrel));
		double sqrt5 = sqrt(5.);
		assert (fabs(root - sqrt5) < sqrt5*epsrel);
	}
}
inline void test_gsl_root_fdfsolver()
{
	// root_fdfsolver class with function<tuple<double,double>(double)>
	{
		gsl::root::fdfsolver s(gsl_root_fdfsolver_newton);

		// F(x) = x^2 - 5
		auto F = [](double x) {
			return x*x - 5;
		};
		auto dF = [](double x) {
			return 2*x;
		};
		double x = 5.0, epsrel = 1e-6;
		s.set(F, dF, x);

		double root = s.solve(gsl::root::test_delta(0, epsrel));
		double sqrt5 = sqrt(5.);
		assert (fabs(root - sqrt5) < sqrt5*epsrel);
	}
}

#endif // _DEBUG