text:
/**
 * \file IIRFilter.h
 */

#ifndef ATK_EQ_IIRFILTER_H
#define ATK_EQ_IIRFILTER_H

#include <algorithm>
#include <cassert>
#include <vector>

#include <gsl/gsl>

#include <ATK/config.h>
#include <ATK/Core/TypeTraits.h>
#include <ATK/EQ/config.h>

namespace ATK
{
  /// IIR filter template class (Direct Form I)
  template<class Coefficients >
  class IIRFilter final : public Coefficients
  {
  protected:
    /// Simplify parent calls
    using Parent = Coefficients;
    using typename Parent::DataType;
    using typename Parent::AlignedScalarVector;
    using Parent::converted_inputs;
    using Parent::outputs;
    using Parent::coefficients_in;
    using Parent::coefficients_out;
    using Parent::input_sampling_rate;
    using Parent::output_sampling_rate;
    using Parent::nb_input_ports;
    using Parent::nb_output_ports;

    using Parent::in_order;
    using Parent::out_order;
    using Parent::input_delay;
    using Parent::output_delay;
    using Parent::setup;

  public:
    /*!
     * @brief Constructor
     * @param nb_channels is the number of input and output channels
     */
    explicit IIRFilter(gsl::index nb_channels = 1)
      :Parent(nb_channels)
    {
    }

    /// Move constructor
    IIRFilter(IIRFilter&& other)
    :Parent(std::move(other))
    {
    }

    void setup() final
    {
      Parent::setup();
      input_delay = in_order;
      output_delay = out_order;

      if (out_order > 0)
      {
        coefficients_out_2.resize(out_order, 0);
        for (unsigned int i = 1; i < out_order; ++i)
        {
          coefficients_out_2[i] = coefficients_out[out_order - 1] * coefficients_out[i] + coefficients_out[i - 1];
        }
        coefficients_out_2[0] = coefficients_out[out_order - 1] * coefficients_out[0];
      }
      if (out_order > 1)
      {
        coefficients_out_3.resize(out_order, 0);
        for (unsigned int i = 0; i < 2; ++i)
        {
          coefficients_out_3[i] = coefficients_out[out_order - 2]  * coefficients_out[i] + coefficients_out[out_order - 1] * coefficients_out_2[i];
        }
        for (unsigned int i = 2; i < out_order; ++i)
        {
          coefficients_out_3[i] = coefficients_out[out_order - 2]  * coefficients_out[i] + coefficients_out[out_order - 1] * coefficients_out_2[i] + coefficients_out[i - 2];
        }
        if (out_order > 2)
        {
          coefficients_out_4.resize(out_order, 0);
          for (unsigned int i = 0; i < 3; ++i)
          {
            coefficients_out_4[i] = coefficients_out[out_order - 3]  * coefficients_out[i] + coefficients_out[out_order - 2] * coefficients_out_2[i] + coefficients_out[out_order - 1] * coefficients_out_3[i];
          }
          for (unsigned int i = 3; i < out_order; ++i)
          {
            coefficients_out_4[i] = coefficients_out[out_order - 3]  * coefficients_out[i] + coefficients_out[out_order - 2] * coefficients_out_2[i] + coefficients_out[out_order - 1] * coefficients_out_3[i] + coefficients_out[i - 3];
          }
        }
        else // out_order = 2
        {
          coefficients_out_4.resize(out_order, 0);
          for (unsigned int i = 0; i < 2; ++i)
          {
            coefficients_out_4[i] = coefficients_out[out_order - 2] * coefficients_out_2[i] + coefficients_out[out_order - 1] * coefficients_out_3[i];
          }
        }
      }
    }

    template<typename T>
    void handle_recursive_iir(const T* ATK_RESTRICT coefficients_out_ptr, const T* ATK_RESTRICT coefficients_out_2_ptr, const T* ATK_RESTRICT coefficients_out_3_ptr, const T* ATK_RESTRICT coefficients_out_4_ptr, DataType* ATK_RESTRICT output, gsl::index size) const
    {
      gsl::index i = 0;
      if (out_order > 2)
      {
        for (i = 0; i < std::min(size - 3, size); i += 4)
        {
          DataType tempout = output[i];
          DataType tempout2 = output[i] * coefficients_out_ptr[out_order - 1] + output[i + 1];
          DataType tempout3 = output[i] * coefficients_out_ptr[out_order - 2] + tempout2 * coefficients_out_ptr[out_order - 1] + output[i + 2];
          DataType tempout4 = output[i] * coefficients_out_ptr[out_order - 3] + tempout2 * coefficients_out_ptr[out_order - 2] + tempout3 * coefficients_out_ptr[out_order - 1] + output[i + 3];

          ATK_VECTORIZE_REMAINDER for (unsigned int j = 0; j < out_order; ++j)
          {
            tempout += coefficients_out_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
            tempout2 += coefficients_out_2_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
            tempout3 += coefficients_out_3_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
            tempout4 += coefficients_out_4_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
          }
          output[i] = tempout;
          output[i + 1] = tempout2;
          output[i + 2] = tempout3;
          output[i + 3] = tempout4;
        }
      }
      else if(out_order == 2)
      {
        for (i = 0; i < std::min(size - 3, size); i += 4)
        {
          DataType tempout = output[i];
          DataType tempout2 = output[i] * coefficients_out_ptr[out_order - 1] + output[i + 1];
          DataType tempout3 = output[i] * coefficients_out_ptr[out_order - 2] + tempout2 * coefficients_out_ptr[out_order - 1] + output[i + 2];
          DataType tempout4 = tempout2 * coefficients_out_ptr[out_order - 2] + tempout3 * coefficients_out_ptr[out_order - 1] + output[i + 3];

          ATK_VECTORIZE_REMAINDER for (unsigned int j = 0; j < out_order; ++j)
          {
            tempout += coefficients_out_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
            tempout2 += coefficients_out_2_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
            tempout3 += coefficients_out_3_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
            tempout4 += coefficients_out_4_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
          }
          output[i] = tempout;
          output[i + 1] = tempout2;
          output[i + 2] = tempout3;
          output[i + 3] = tempout4;
        }
      }
      for (; i < size; ++i)
      {
        DataType tempout = output[i];
        for (unsigned int j = 0; j < out_order; ++j)
        {
          tempout += coefficients_out_ptr[j] * output[static_cast<int64_t>(i) - out_order + j];
        }
        output[i] = tempout;
      }
    }

    void process_impl(gsl::index size) const final
    {
      assert(input_sampling_rate == output_sampling_rate);
      assert(nb_input_ports == nb_output_ports);
      assert(coefficients_in.data());
      assert(out_order == 0 || coefficients_out.data() != nullptr);

      const auto* ATK_RESTRICT coefficients_in_ptr = coefficients_in.data();
      const auto* ATK_RESTRICT coefficients_out_ptr = coefficients_out.data();
      const auto* ATK_RESTRICT coefficients_out_2_ptr = coefficients_out_2.data();
      const auto* ATK_RESTRICT coefficients_out_3_ptr = coefficients_out_3.data();
      const auto* ATK_RESTRICT coefficients_out_4_ptr = coefficients_out_4.data();

      for(gsl::index channel = 0; channel < nb_input_ports; ++channel)
      {
        const DataType* ATK_RESTRICT input = converted_inputs[channel] - static_cast<int64_t>(in_order);
        DataType* ATK_RESTRICT output = outputs[channel];

        for(gsl::index i = 0; i < size; ++i)
        {
          output[i] = 0;
        }

        for (gsl::index j = 0; j < in_order + 1; ++j)
        {
          for (gsl::index i = 0; i < size; ++i)
          {
            output[i] += coefficients_in_ptr[j] * input[i + j];
          }
        }

        handle_recursive_iir(coefficients_out_ptr, coefficients_out_2_ptr, coefficients_out_3_ptr, coefficients_out_4_ptr, output, size);
      }
    }

    /// Returns the vector of internal coefficients for the MA section
    const AlignedScalarVector& get_coefficients_in() const
    {
      return coefficients_in;
    }

    /// Returns the vector of internal coefficients for the AR section, without degree 0 implicitely set to -1
    const AlignedScalarVector& get_coefficients_out() const
    {
      return coefficients_out;
    }

  protected:
    AlignedScalarVector coefficients_out_2;
    AlignedScalarVector coefficients_out_3;
    AlignedScalarVector coefficients_out_4;
  };

  /// IIR filter template class. Transposed Direct Form II implementation
  template<class Coefficients >
  class IIRTDF2Filter final : public Coefficients
  {
  public:
    /// Simplify parent calls
    using Parent = Coefficients;
    using typename Parent::DataType;
    using typename Parent::AlignedScalarVector;
    using Parent::converted_inputs;
    using Parent::outputs;
    using Parent::coefficients_in;
    using Parent::coefficients_out;
    using Parent::input_sampling_rate;
    using Parent::output_sampling_rate;
    using Parent::nb_input_ports;
    using Parent::nb_output_ports;

    using Parent::in_order;
    using Parent::out_order;
    using Parent::input_delay;
    using Parent::output_delay;
    using Parent::setup;
  protected:
    mutable typename Parent::AlignedVector state;
  public:
    explicit IIRTDF2Filter(gsl::index nb_channels = 1)
      :Parent(nb_channels)
    {
    }

    /// Move constructor
    IIRTDF2Filter(IIRTDF2Filter&& other)
      :Parent(std::move(other))
    {
    }

    void setup() final
    {
      Parent::setup();
      input_delay = in_order;
      output_delay = out_order;
      state.assign(nb_input_ports * (std::max(input_delay, output_delay) + 1), TypeTraits<DataType>::Zero());
    }

    void process_impl(gsl::index size) const final
    {
      assert(input_sampling_rate == output_sampling_rate);

      for(gsl::index channel = 0; channel < nb_input_ports; ++channel)
      {
        const DataType* ATK_RESTRICT input = converted_inputs[channel];
        DataType* ATK_RESTRICT output = outputs[channel];
        DataType* ATK_RESTRICT current_state = &state[channel * (std::max(input_delay, output_delay) + 1)];

        for(gsl::index i = 0; i < size; ++i)
        {
          output[i] = coefficients_in[in_order] * input[i] + current_state[0];
          auto min_order = std::min(input_delay, output_delay);

          for(gsl::index j = 0; j < min_order; ++j)
          {
            current_state[j] = current_state[j + 1] + input[i] * coefficients_in[in_order - static_cast<int64_t>(j) - 1] + output[i] * coefficients_out[out_order - static_cast<int64_t>(j) - 1];
          }
          for(gsl::index j = min_order; j < input_delay; ++j)
          {
            current_state[j] = current_state[j + 1] + input[i] * coefficients_in[in_order - static_cast<int64_t>(j) - 1];
          }
          for(gsl::index j = min_order; j < output_delay; ++j)
          {
            current_state[j] = current_state[j + 1] + output[i] * coefficients_out[out_order - static_cast<int64_t>(j) - 1];
          }
        }
      }
    }

    /// Returns the vector of internal coefficients for the MA section
    const AlignedScalarVector& get_coefficients_in() const
    {
      return coefficients_in;
    }

    /// Returns the vector of internal coefficients for the AR section, without degree 0 implicitely set to -1
    const AlignedScalarVector& get_coefficients_out() const
    {
      return coefficients_out;
    }
  };

}

#endif

max_stars_repo_name: D-J-Roberts/AudioTK

max_issues_repo_name: D-J-Roberts/AudioTK

max_forks_repo_name: D-J-Roberts/AudioTK


text:
#ifndef TOY_VN_GENERATOR_H
#define TOY_VN_GENERATOR_H

#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <complex>
#include <vector>
#include <iomanip>
#include <stdio.h>
#include <algorithm>
#include <chrono>
#include <random>
#include <omp.h>
#include "sampled_distribution.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

using namespace std;

namespace flow_generator
{
	// ------------------------------------
	// Initialize default parameters here.
	// ------------------------------------

	int n_pT                               = 6;
	double max_pT                          = 3000000.0;
	double min_pT				           = 0.0;
	int order                              = 2;
	//double eta_low                         = 2.0;
	//double eta_high                        = 10.0;

	int N_events_to_generate               = 1000;
	int N_particles_per_event              = 100;
	long N_total_events                    = 0;
	long N_running_total_events            = 0;

	bool use_seed                          = true;
	int seed_index                         = -1;

	bool print_randomly_generated_data     = false;

	bool include_psi2_PTslope_fluctuations = false;
	double v2_magnitude                    = 0.1;

	string resultsDirectory                = "./results/";
	string output_name                     = resultsDirectory + "/data.dat";

	// ------------------------------------
	// Define structs and vectors.
	// ------------------------------------

	// struct for event information
	typedef struct
	{
		int eventID;
		int Nparticles;

		vector<vector<double> > eta_vec;
		vector<vector<double> > pphi_vec;
		vector<vector<double> > pT_vec;
	} event_data;

	// declare needed vectors
	vector<event_data>                         all_events;

	vector<vector<double> >                    pT_vec;
	vector<vector<double> >                    eta_vec;
	vector<vector<double> >                    phi_vec;
	vector<vector<int> >                       event_vec;

	//---------------------------------------------------------
	// Generate data randomly to mimic (non-factorizable) flow.
	void generate_random_data(int dataset=-1)
	{
		double delta_pT = (max_pT-min_pT)/(n_pT);

		pT_vec.clear();
		eta_vec.clear();
		phi_vec.clear();
		event_vec.clear();
		pT_vec = vector<vector<double> >(n_pT, vector<double> (0));
		eta_vec = vector<vector<double> >(n_pT, vector<double> (0));
		phi_vec = vector<vector<double> >(n_pT, vector<double> (0));
		event_vec = vector<vector<int> >(n_pT, vector<int> (0));

		// Set up distributions
		unsigned seed = chrono::system_clock::now().time_since_epoch().count();
		default_random_engine generator;
		if ( use_seed )
			generator = default_random_engine(seed);
		else
			generator = default_random_engine(seed_index);

		uniform_real_distribution<double> mean_pT_distribution(0.90, 1.10);
		//const double mean_pT = mean_pT_distribution(generator);
		const double mean_pT = 1.0;
		exponential_distribution<double> pT_distribution( mean_pT );
		normal_distribution<double> eta_distribution(0.0, 2.0);
		uniform_real_distribution<double> v2_distribution(0.05, 0.45);
		uniform_real_distribution<double> psi2_distribution(0.0, 2.0*M_PI);
		uniform_real_distribution<double> pphi_distribution(0.0, 2.0*M_PI);
		normal_distribution<double> psi2_pTslope_distribution(0.0,0.1);

		auto phiFunc = [](double x, double v2) { return x + v2*sin(2.0*x); };
	    std::mt19937 gen;
		if ( use_seed )
			gen = std::mt19937(seed);
		else
			gen = std::mt19937(seed_index);

		double fluctuation_switch_factor
				= double(include_psi2_PTslope_fluctuations);

		for (int iEvent = 0; iEvent < N_events_to_generate; iEvent++)
		{

			//double v2 = v2_distribution(generator);
			double v2 = v2_magnitude;
			double psi2 = 0.0*psi2_distribution(generator);
			double psi2_pTslope = fluctuation_switch_factor
									* psi2_pTslope_distribution(generator);

			for (int iParticle = 0; iParticle < N_particles_per_event; iParticle++)
			{
				// ---------------------------------------------
				// Try to mimic different types of flow effects
				// ---------------------------------------------
				const double pT_value = pT_distribution(generator);
			    //Sampled_distribution<> dist(phiFunc, 0.0, 2.0*M_PI, v2*pT_value / (mean_pT+pT_value));
			    Sampled_distribution<> dist(phiFunc, 0.0, 2.0*M_PI, v2);
				const double pphi_value = dist(gen) + psi2 + psi2_pTslope*pT_value;
				//const double pphi_value = dist(gen) + psi2 + 0.5*M_PI*pT_value / (mean_pT+pT_value);

				// Sort into bins.
				double pT = pT_value;
				double eta = 0.0*eta_distribution(generator);
				double phi = pphi_value;
				int event = iEvent;

				if(pT > max_pT)
				{
					//cout << "Warning: pT = " << pT << " larger than max_pT = " << max_pT << " in event = " << event << endl;
					continue;
				}
				if(pT < min_pT)
				{
					//cout << "Warning: pT = " << pT << " smaller than min_pT = " << min_pT << " in event = " << event << endl;
					continue;
				}
				else
				{
					int bin = int((pT-min_pT)/delta_pT);
					eta_vec[bin].push_back(eta);
					phi_vec[bin].push_back(phi);
					pT_vec[bin].push_back(pT);
					event_vec[bin].push_back(event);
					//cout << "Check: " << pT << "   " << eta << "   " << phi << "   " << event << " " << bin << endl;
				}
			}
		}

		//----------------------------------------------
		all_events.clear();
		all_events = vector<event_data>(N_events_to_generate);
		for ( auto & this_event : all_events )
		{
			this_event.eta_vec = vector<vector<double> >( n_pT, vector<double>(0) );
			this_event.pphi_vec = vector<vector<double> >( n_pT, vector<double>(0) );
			this_event.pT_vec = vector<vector<double> >( n_pT, vector<double>(0) );
		}
		for (int bin = 0; bin < n_pT; bin++)
		{
			const int n_particles_in_this_bin = event_vec[bin].size();
			for (int iParticle = 0; iParticle < n_particles_in_this_bin; iParticle++)
			{
				int event_ID_of_this_particle = event_vec[bin][iParticle];

				all_events[event_ID_of_this_particle].eta_vec[bin].push_back( eta_vec[bin][iParticle] );
				all_events[event_ID_of_this_particle].pphi_vec[bin].push_back( phi_vec[bin][iParticle] );
				all_events[event_ID_of_this_particle].pT_vec[bin].push_back( pT_vec[bin][iParticle] );
				all_events[event_ID_of_this_particle].eventID = event_ID_of_this_particle;
			}
		}
		//----------------------------------------------

		//----------------------------------------------
		// Tally number of completed events
		N_running_total_events += N_events_to_generate;
		cout << "Generated " << N_particles_per_event * N_events_to_generate
				<< " particles belonging to "
				<< N_events_to_generate << " events."
				<< endl;
		cout << "Currently analyzed " << N_running_total_events
				<< " / " << N_total_events << " events." << endl;

		//----------------------------------------------
		// Dump simulated data to file if desired
		if ( print_randomly_generated_data )
		{
			int n_zero = static_cast<int>( log10( N_particles_per_event*N_total_events ) )-1;
			string dataset_stem = ( dataset >= 0 ) ?
									std::string(n_zero - std::to_string(dataset).length(), '0')
									+ std::to_string( dataset ) : "";
			string RNG_filename = resultsDirectory + "/event_" + dataset_stem + ".dat";
			ofstream RNG_output( RNG_filename.c_str() );
			for ( auto & this_event : all_events )
			for (int bin = 0; bin < n_pT; bin++)
			{
				const int n_particles_in_this_bin = this_event.eta_vec[bin].size();
				const auto & this_pT_bin = this_event.pT_vec[bin];
				const auto & this_pphi_bin = this_event.pphi_vec[bin];
				const auto & this_eta_bin = this_event.eta_vec[bin];
				for (int iParticle = 0; iParticle < n_particles_in_this_bin; iParticle++)
					RNG_output << setprecision(15) << setw(18)
							   /*<< this_pT_bin[iParticle] << "   "
                               << this_eta_bin[iParticle] << "   "*/
                               << this_pphi_bin[iParticle] /*<< "   "
                               << this_event.eventID*/ << endl;
			}
			RNG_output.close();
		}

		return;
	}

}

#endif

ext: h


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <Accelerate/Accelerate.h>

struct _FINCHAIN {
  int length; // Length of state vector
  int currentIter;

  double *currentState, *avgState, *varState, *proposedState, *_M2;
  double logLHDCurrentState;
  double accProb, avgAccProb;

  gsl_rng *r;
};

typedef struct _FINCHAIN FINCHAIN;

/*
  Initialises the finite dimensional chain
*/
void finmcmc_initChain(FINCHAIN *C, int n);

/*
  Frees the finite dimensional chain
*/
void finmcmc_freeChain(FINCHAIN *C);

void finmcmc_printChain(FINCHAIN C);

void finmcmc_writeChain(const FINCHAIN *C, FILE *fp);

void finmcmc_gausRanVec(const FINCHAIN *C, double *x, double stdDev);

/*
  Proposes a standard RWMH move
*/
void finmcmc_proposeRWMH(const FINCHAIN *C, double stdDev, double beta);

/*
  Update chain
*/
void finmcmc_updateRWMH(FINCHAIN *C, double logLHDOfProposal);

//void fin_initChain(FINCHAIN *C, int n, int sizePhi);

//void sampleRMWH(CHAIN *C);
//void sampleIndependenceSampler(CHAIN *C);
//void updateMean(CHAIN *C);
//void updateVar(CHAIN *C);
//void randomPriorDraw(gsl_rng *r, double PRIOR_ALPHA, fftw_complex *randDrawCoeffs);

max_stars_repo_name: dmcdougall/mcmclib

max_issues_repo_name: dmcdougall/mcmclib

max_forks_repo_name: dmcdougall/mcmclib


text:
/*
 *  test.c
 *
 *
 *  Created by ashleyr on 30/06/2010.
 *  Copyright 2010 __MyCompanyName__. All rights reserved.
 *
 */
#include "headera.h"
#include <stdlib.h>
//#include <malloc.h>
//#include <fstream>
//#include <iomanip>
#include <stdio.h>
#include <math.h>
//#include <string>
#include <time.h>
#include <gsl/gsl_math.h>
//#include <gsl/gsl_complex.h>
//#include <gsl/gsl_complex_math.h>
//#include <gsl/gsl_const_num.h>
//#include <gsl/gsl_const_mksa.h>
//#include <gsl/gsl_sf.h>
#include <gsl/gsl_linalg.h>
//#include <gsl/gsl_eigen.h>
//#include <gsl/gsl_blas.h>
//#include <gsl/gsl_permutation.h>
//#include <gsl/gsl_errno.h>
//#include <gsl/gsl_matrix.h>
//#include <gsl/gsl_integration.h>
//#include <gsl/gsl_monte.h>
//#include <gsl/gsl_monte_plain.h>
//#include <gsl/gsl_monte_miser.h>
//#include <gsl/gsl_monte_vegas.h>
//#include <gsl/gsl_rng.h>
//#include <gsl/gsl_randist.h>
//#include <gsl/gsl_odeiv.h>

main(int argc, char *argv[]) {
	int nran = 100000;
	int Nline_max =1000;
	double a;
	double *lcd = dvector(1,1000000);
	double *lcr = dvector(1,1000000);
	double *lsd = dvector(1,1000000);
	double *lsr = dvector(1,1000000);
	double *lzw = dvector(1,1000000);
	double *lpc = dvector(1,1000000);
	int i = 0;
	int k = 0;
	int j = 0;
	int N = 0;
	double ti = time(NULL);
	//int nbin = 99;
	double angm;
	angm = atof(argv[2]);
	//double minang = .2*pi/180.*angm;
	double minang = 0*pi/180.;
	double binsize = 0.15*angm*pi/180.;
	double maxang = angm*12.*pi/180.;
	int nbin = (maxang-minang)/binsize;
	double binl[nbin];
    double angl[nbin];
	double binlpix[nbin];
	double binedges[nbin+2];
	//double g = minang/maxang;
	for(int nb=0;nb<nbin;nb++) {
		binl[nb] = 0;
		binlpix[nb] = 0;
	}
	for(int bn=0;bn<nbin+2;bn++) {
		double num = (bn);
		double angle = maxang - binsize*num;
		if (angle < 0) angle = 0;
		double angb = maxang-binsize*num-binsize/2.;
        binedges[bn] = cos(angle);
        if (bn<nbin) {
            angl[bn] =angb;
        }
		//printf("%lf %lf %lf\n",angle,binedges[bn],cos(angle));
	}
    char fxi_name[200];
    char s1[200];
    sprintf(s1,argv[1]);
    strcpy(fxi_name,s1);
    sprintf(fxi_name,strcat(fxi_name,"odenspczw.dat"));
    FILE *fxi;
    printf("%s\n",s1);
    fxi=fopen(fxi_name,"r");
    if (fxi == NULL) perror ("Error opening file");
    int Nline=0;
    const int bsz=240;
    char buf[bsz];

    while(fgets(buf,bsz,fxi)!= NULL) {
        double cd,cr,sr,sd,zw,pc;
        //if(sscanf(buf,"%*d %f %f %f %f",&cd,&cr,&sd,&sd)!=4)
        //	err_handler("xi input");
        //puts(buf);
        sscanf(buf,"%lf %lf %lf %lf %lf %lf",&sr,&cr,&sd,&cd,&zw,&pc);
        //cd = 1.;
        //cr = 1.;
        //sr = 1.;
        //sd = 1.;
        lcd[Nline] = cd;
        lcr[Nline] = cr;
        lsr[Nline] = sr;
        lsd[Nline] = sd;
        lzw[Nline] = zw;
        lpc[Nline] = pc;
        Nline++;
    }
    printf("%i \n",Nline);
    /*	while (i<nran) {

     l[i] = rand();
     i = i+1;

     }
     */
    float be;
    float r;
    while (k<Nline)  {
        j = k+1;
        while (j<Nline)  {
            a = lcd[k]*lcd[j]*(lcr[k]*lcr[j] + lsr[k]*lsr[j]) + lsd[k]*lsd[j];
            //r = lsr[k];
            be = binedges[0];
            int ba = -1;
            while (a > be) {
                ba++;
                be = binedges[ba+1];
            }
            if (ba > -1 && ba < nbin){
                binlpix[ba] += lpc[k]*lpc[j];
                binl[ba] += lzw[k]*lzw[j]*lpc[k]*lpc[j];
            }
            j++;
            N+=1;

        }
        //printf("%f %f\n",be,r);
        k++;
    }
    double tf = time(NULL)-ti;
    printf("took %g\n",tf);
    char foname[200];
    //char sf1[200];
    strcat(s1,"2ptPixclb.dat");

    sprintf(foname,s1);
    FILE *fo;
    fo =fopen(foname,"w");
    for (int ip=0; ip<nbin; ip++) {
        fprintf(fo,"%g %g \n",angl[ip]*180./pi,binl[ip]/binlpix[ip]);
    }
    //fclose(fxi);
    //exit(0);
}

ext: c

max_stars_repo_name: mehdirezaie/LSSutils

max_issues_repo_name: mehdirezaie/LSSutils

max_forks_repo_name: mehdirezaie/LSSutils


text:
//===--- Sudoku/Location_Utilities.h                                    ---===//
//
// Utilities for class Sudoku::Location
//===----------------------------------------------------------------------===//
#pragma once

#include "Board_Section_traits.h"
#include "Location.h"
#include "Size.h"
#include "traits.h"

#include <gsl/gsl> // index

#include <vector>

#include <algorithm> // minmax_element, is_sorted, all_of
#include <iterator>  // back_inserter
#include <limits>    // numeric_limits

#include "Board.fwd.h" // Forward declarations

namespace Sudoku
{
//===-- function declarations ---------------------------------------------===//
template<int N>
constexpr void valid_dimensions() noexcept;

template<int N>
constexpr bool is_valid(Location<N>) noexcept;
template<int N>
constexpr bool is_valid(const std::vector<Location<N>>& locs) noexcept(true);

template<int N>
constexpr bool is_valid_size(gsl::index elem) noexcept;
template<int N>
constexpr bool is_valid_size(gsl::index row, gsl::index col) noexcept;

template<int N>
constexpr bool is_same_row(Location<N>, Location<N>) noexcept;
template<int N>
constexpr bool is_same_col(Location<N>, Location<N>) noexcept;
template<int N>
constexpr bool is_same_block(Location<N>, Location<N>) noexcept;

template<int N, typename ItrT>
constexpr bool is_same_row(ItrT begin, ItrT end) noexcept;
template<int N, typename ItrT>
constexpr bool is_same_col(ItrT begin, ItrT end) noexcept;
template<int N, typename ItrT>
constexpr bool is_same_block(ItrT begin, ItrT end) noexcept;

template<int N>
std::vector<Location<N>> get_same_row(
	const Location<N>, const std::vector<Location<N>>&) noexcept(true);
template<int N>
std::vector<Location<N>> get_same_col(
	const Location<N>, const std::vector<Location<N>>&) noexcept(true);
template<int N>
std::vector<Location<N>> get_same_block(
	const Location<N>, const std::vector<Location<N>>&) noexcept(true);

template<int N, typename S>
constexpr bool is_same_section(S section, Location<N>) noexcept;

template<typename SectionT, int N>
bool is_same_section(SectionT, std::vector<Location<N>>&) noexcept;

template<
	int N,
	typename S,
	typename = std::enable_if_t<
		Board_Section::traits::is_Row_v<S> ||
		Board_Section::traits::is_Col_v<S>>>
bool intersect_block(S section, Location<N> block_loc) noexcept;

//===----------------------------------------------------------------------===//

// Compile-time only Test
template<int N>
constexpr void valid_dimensions() noexcept
{
	// input check
	static_assert((base_size<N>) > 1, "base_size too small");
	static_assert(
		base_size<N> < elem_size<N> && elem_size<N> <= full_size<N> &&
			base_size<N> < std::numeric_limits<int>::max() && // <limits>
			elem_size<N> < std::numeric_limits<int>::max() &&
			full_size<N> < std::numeric_limits<int>::max(),
		"board size out of bounds");
	// logic check
	static_assert(
		base_size<N> * base_size<N> == elem_size<N> &&
			elem_size<N> * elem_size<N> == full_size<N>,
		"size calculation broken");
}

// Test if Location on Board
template<int N>
[[nodiscard]] inline constexpr bool is_valid(const Location<N> loc) noexcept
{
	return (loc.element() >= 0 && loc.element() < full_size<N>);
}

// Test if Locations on Board and if sorted (ascending)
template<int N>
[[nodiscard]] inline constexpr bool
	is_valid(const std::vector<Location<N>>& locs) noexcept(true)
{ // std::is_sorted can throw std::bad_alloc
	return (
		!locs.empty() && (std::is_sorted(locs.cbegin(), locs.cend()) &&
						  locs.cbegin()->element() >= 0 &&
						  locs.crbegin()->element() < full_size<N>));
}

// Test row/col/block-element
template<int N>
[[nodiscard]] inline constexpr bool
	is_valid_size(const gsl::index elem) noexcept
{
	return (elem >= 0 && elem < elem_size<N>);
}

// Test if location on Board
template<int N>
[[nodiscard]] inline constexpr bool
	is_valid_size(const gsl::index row, const gsl::index col) noexcept
{
	return is_valid_size<N>(row) && is_valid_size<N>(col);
}

//===----------------------------------------------------------------------===//

// check
template<int N>
[[nodiscard]] inline constexpr bool
	is_same_row(const Location<N> left, const Location<N> right) noexcept
{
	return (is_valid<N>(left) && is_valid<N>(right)) &&
		   left.row() == right.row();
}

// check: all in same row
template<int N, typename ItrT>
[[nodiscard]] constexpr bool
	is_same_row(const ItrT begin, const ItrT end) noexcept
{
	{
		static_assert(traits::is_forward<ItrT>);
	}
	if (begin == end)
		return false;

	const auto itr = begin + 1;
	return std::all_of(
		itr, end, [begin](Location<N> i) { return is_same_row<N>(*begin, i); });
}

// return all in same row
template<int N> // NOLINTNEXTLINE(bugprone-exception-escape)
[[nodiscard]] std::vector<Location<N>> get_same_row(
	const Location<N> left,
	const std::vector<Location<N>>& right) noexcept(true)
{ // std::copy_if could throw std::bad_alloc
	std::vector<Location<N>> output{};
	const auto predicate = [&left](Location<N> loc) {
		return is_same_row(left, loc);
	};
	std::copy_if(
		right.cbegin(), right.cend(), std::back_inserter(output), predicate);
	return output;
}

// check
template<int N>
[[nodiscard]] inline constexpr bool
	is_same_col(const Location<N> left, const Location<N> right) noexcept
{
	return (is_valid<N>(left) && is_valid<N>(right)) &&
		   left.col() == right.col();
}

// check: all in same col
template<int N, typename ItrT>
[[nodiscard]] inline constexpr bool
	is_same_col(const ItrT begin, const ItrT end) noexcept
{
	{
		static_assert(traits::is_forward<ItrT>);
	}
	if (begin == end)
		return false;

	const auto itr = begin + 1;
	return std::all_of(
		itr, end, [begin](Location<N> i) { return is_same_col<N>(*begin, i); });
}

// return all in same col
template<int N> // NOLINTNEXTLINE(bugprone-exception-escape)
[[nodiscard]] std::vector<Location<N>> get_same_col(
	const Location<N> left,
	const std::vector<Location<N>>& right) noexcept(true)
{ // std::copy_if could throw std::bad_alloc
	std::vector<Location<N>> output{};
	const auto predicate = [&left](Location<N> loc) {
		return is_same_col(left, loc);
	};
	std::copy_if(
		right.cbegin(), right.cend(), std::back_inserter(output), predicate);
	return output;
}

// check
template<int N>
[[nodiscard]] inline constexpr bool
	is_same_block(const Location<N> left, const Location<N> right) noexcept
{
	return (is_valid<N>(left) && is_valid<N>(right)) &&
		   left.block() == right.block();
}

// check all in same block
template<int N, typename ItrT>
[[nodiscard]] inline constexpr bool
	is_same_block(const ItrT begin, const ItrT end) noexcept
{
	{
		static_assert(traits::is_forward<ItrT>);
	}
	if (begin == end)
		return false;

	const auto itr = begin + 1;
	return std::all_of(itr, end, [begin](Location<N> i) {
		return is_same_block<N>(*begin, i);
	});
}

// return all in same block
template<int N> // NOLINTNEXTLINE(bugprone-exception-escape)
[[nodiscard]] std::vector<Location<N>> get_same_block(
	const Location<N> left,
	const std::vector<Location<N>>& right) noexcept(true)
{ // std::copy_if could throw std::bad_alloc
	std::vector<Location<N>> output{};
	const auto predicate = [&left](Location<N> loc) {
		return is_same_block(left, loc);
	};
	std::copy_if(
		right.cbegin(), right.cend(), std::back_inserter(output), predicate);
	return output;
}

// check: [loc] is in [section]
template<int N, typename S>
[[nodiscard]] inline constexpr bool
	is_same_section(const S section, const Location<N> loc) noexcept
{
	static_assert(Board_Section::traits::is_Section_v<S>);

	if constexpr (Board_Section::traits::is_Row_v<S>)
		return is_same_row(loc, section.cbegin().location());
	else if constexpr (Board_Section::traits::is_Col_v<S>)
		return is_same_col(loc, section.cbegin().location());
	else
	{
		static_assert(Board_Section::traits::is_Block_v<S>);
		return is_same_block(loc, section.cbegin().location());
	}
}

// check: [section] intersects block containing [loc]
template<int N, typename S, typename>
[[nodiscard]] inline bool
	intersect_block(const S section, const Location<N> block_loc) noexcept
{
	for (auto itr = section.cbegin(); itr != section.cend(); ++itr)
	{
		if (is_same_block(block_loc, itr.location()))
		{
			return true;
		}
	}
	return false;
}

// check: at least one [location] inside [section]
template<typename SectionT, int N>
[[nodiscard]] inline bool is_same_section(
	SectionT const section, std::vector<Location<N>> const& locs) noexcept
{
	return std::any_of(locs.cbegin(), locs.cend(), [section](Location<N> L) {
		return is_same_section(section, L);
	});
}

} // namespace Sudoku

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


text:
#pragma once
//=============================================================================
// EXTERNAL DECLARATIONS
//=============================================================================
#include "Core/Components/IComponentInterface.h"
#include "Core/Components/IdType.h"
#include <memory>
#include <gsl/span>

namespace engine
{

//=============================================================================
// FORWARD DECLARATIONS
//=============================================================================
class IGameObject;
using GameObjectRef = std::shared_ptr<IGameObject>;
using GameObjectWeakRef = std::weak_ptr<IGameObject>;

//=============================================================================
// INTERFACE IComponent
//=============================================================================
class IComponent
{
public:
	virtual ~IComponent() {}

public:
	virtual gsl::span<IdType> interfaces() = 0;
	virtual void onAttached(const GameObjectRef &iGameObject) = 0;
	virtual void onDetached(const GameObjectRef &iGameObject) = 0;
};

} // namespace engine


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#ifndef _writeMat_h
#define _writeMat_h
#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>

void bsscr_writeMat(Mat A, char name[], char message[]);
void bsscr_writeVec(Vec V, char name[], char message[]);
#endif