text:
//
//  ger.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 12/23/11.
//  Copyright (c) 2011 University of Colorado Denver. All rights reserved.
//

#ifndef _ger_h
#define _ger_h

/// @file ger.h Performs vector outer product.

#include "latl.h"

namespace LATL
{
   /// @brief Performs a vector outer product of two real vectors.
   ///
   /// For a real matrix A, real vectors x and y, and real scalar alpha,
   ///
   ///        A := alpha * x * y' + A
   ///
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param m Specifies the number of rows of the matrix A.  m>=0
   /// @param n Specifies the number of coumns of the matrix A.  n>=0
   /// @param alpha Real scalar.
   /// @param x Pointer to real vector x.
   /// @param incx Increment of the vector x.  x!=0
   /// @param y Pointer to real vector y.
   /// @param incy Increment of the vector y.  y!=0
   /// @param A Pointer to real m-by-n matrix A.
   /// @param ldA Column length of matrix A.  ldA>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int GER(int_t m, int_t n, real_t alpha, real_t *x, int_t incx, real_t *y, int_t incy, real_t *A, int_t ldA)
   {
      const real_t zero(0.0);
      int_t i,j,kx,jy,ix;
      if(m<0)
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(incy==0)
         return -7;
      else if(ldA<m)
         return -9;
      else if((m==0)||(n==0)||(alpha==zero))
         return 0;

      jy=(incy>0)?0:(1-n)*incy;
      kx=(incx>0)?0:(1-m)*incx;

      if(incx==1)
      {
         for(j=0;j<n;j++)
         {
            for(i=0;i<m;i++)
               A[i]+=x[i]*alpha*y[jy];
            jy+=incy;
            A+=ldA;
         }
      }
      else
      {
         for(j=0;j<n;j++)
         {
            ix=kx;
            for(i=0;i<m;i++)
            {
               A[i]+=x[ix]*alpha*y[jy];
               ix+=incx;
            }
            A+=ldA;
            jy+=incy;
         }
      }
      return 0;
   }

#ifdef __latl_cblas
#include <cblas.h>

   template <> int GER<float>(int_t m, int_t n, float alpha, float *x, int_t incx, float *y, int_t incy, float *A, int_t ldA)
   {
      if(m<0)
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(incy==0)
         return -7;
      else if(ldA<m)
         return -9;

      cblas_sger(CblasColMajor,m,n,alpha,x,incx,y,incy,A,ldA);

      return 0;
   }

   template <> int GER<double>(int_t m, int_t n, double alpha, double *x, int_t incx, double *y, int_t incy, double *A, int_t ldA)
   {
      if(m<0)
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(incy==0)
         return -7;
      else if(ldA<m)
         return -9;

      cblas_dger(CblasColMajor,m,n,alpha,x,incx,y,incy,A,ldA);

      return 0;
   }

#endif
}
#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
//#include </usr/local/include/gsl/gsl_matrix.h>
#include <gsl/gsl_matrix.h> // suprisingly, both work
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

/*
 *  * Program that reads reference spectra and measure spectra data
 *  and linearly unmixes the data using method of least squares fit.
 *  Assumes readings is set to 400 and file path is preset to:
 *  "/mnt/c/Users/jakei/'OneDrive - McGill University'/'Year 2'/1BIEN210/'Group Project'/".
 *
 * syntax: ./linearUnmixer MUfilename ref1filename ref2filename ref3filename ...
 * returns: completion status, prints composition of ref1, ref2, ... , refn
 *
 *  * ****************************************************************
 *  * Author       Dept.           Date            Notes
 *  * ****************************************************************
 *  * Jake Z       Bio. Eng.       Mar 12 2020     Initial version
 *  * Jake Z		"	   Apr 03 2020     It works now! TODO: investigate the discrepancy between what you set $readings and the amount of space you allocate to the matrices as you forgot to account for the initial data clearing of the title headers in the data.
 *  * */
const int READINGS = 1300;		//note, if the readings specified here is more then the number of readings in file, funky stuff happens, TODO: implent check for this
const int SKIP = 50;
int normalizeVector( gsl_vector *vector )
{
	double min = gsl_vector_min( vector );
	double norm = gsl_vector_max( vector ) - min;
	gsl_vector_add_constant( vector, -min );
	gsl_vector_scale( vector, 1/norm );
	return 0;
}

int readFileVector( gsl_vector *vector, char *filePath, char *fileName )		//read file vals into vector
{
	char tmpFilepath[120] = "\0";
	char tmpReading[20];
	char *tmp2;
	double reading;					//WARNING this is similar to global var $readings, consider renaming

	strcat( tmpFilepath, filePath );		//TODO replace this and the first instance of strcat into strcpy, that should work
	strcat( tmpFilepath, fileName ); 	//combines the filepath and name into one searchable /path/file
	FILE *file = fopen( tmpFilepath, "rt");
	printf("%s\n", tmpFilepath);			//TODO DELETE THIS
	if ( file == NULL ) {			//check if file exists
		fprintf( stderr, "Error, unable to locate the Mu data file\n");
		exit(100);
	}

	for ( int i = 0; i < SKIP; i++)		//move past header of read files
	{
		fgets( tmpReading, 19, file );
	}
	for ( int i = 0; i < READINGS; i++ )
	{
		fgets( tmpReading, 19, file );

		tmp2 = strrchr( tmpReading, '\t' );
		tmp2 = tmp2 + 1;
		reading = atof(tmp2);
		gsl_vector_set( vector, i, reading);
	}
	fclose(file);
	printf("successful vector \n");			//TODO delete this
	return 0;
}

int readFileMatrix( gsl_matrix *matrix, char *filePath, int argc, char *argv[] )
{
	char *tmp2;
	double reading;
	for ( int j = 2; j < argc; j++ )
	{
		gsl_vector_view column;
		char tmpFilepath[120] = "\0";				//TODO replace this and the first instance of strcat into strcpy, that should work
		char tmpReading[20];
		strcat( tmpFilepath, filePath );
		strcat( tmpFilepath, argv[j] ); 	//combines the filepath and name into one searchable /path/file
		FILE *file = fopen( tmpFilepath, "rt");
		printf("%s\n", tmpFilepath);				//TODO DELETE THIS
		if ( file == NULL ) {			//check if file exists
			fprintf( stderr, "Error, unable to locate a reference data file\n");
			exit(100);
		}
		for ( int i = 0; i < SKIP; i++)		//move past header of read files
		{
			fgets( tmpReading, 20, file);
		}
		for ( int i = 0; i < READINGS; i++ ) //fill the ith element of the jth column
		{
			fgets( tmpReading, 19, file );

			tmp2 = strrchr( tmpReading, '\t' );
			tmp2 = tmp2 + 1;
			reading = atof(tmp2);
			gsl_matrix_set( matrix, i, j-2, reading );
		}
		fclose(file);
		column = gsl_matrix_column( matrix, j-2 );				//NORMALISE EACH REFERENCE VECOTR.
		normalizeVector( &column.vector );
	}
	printf("success, matrix filled \n");			//TODO DELETE THIS
	return 0;
}

int main(int argc, char *argv[])
{
	if ( argc <= 2 )
	{
		fprintf( stderr, "Error, incorrect usage.\n./linearUnmixer MU_SPECTRA REF_SPECTRA\n");
		return 3;
	}
	//spectrum files labelled  /mnt/c/Users/jakei/'OneDrive - McGill University'/'Year 2'/1BIEN210/'Group Project'/SpectrumFile_00N.txt
	//intialize and allocate memory for variabels
	int numRefs = argc - 2;
	//char *filePath = "";
	char *filePath = "/mnt/c/Users/jakei/OneDrive - McGill University/Year 2/1BIEN210/Group Project/";
	gsl_vector *muVector = gsl_vector_alloc(READINGS);
	gsl_vector *fVector = gsl_vector_alloc(numRefs);
	gsl_vector *x = gsl_vector_alloc(numRefs);			//TODO RENAME x AND PROBECOMPVECTOR
	gsl_matrix *refMatrix = gsl_matrix_alloc(READINGS, numRefs);
	gsl_matrix *C = gsl_matrix_alloc( numRefs, numRefs );
	gsl_matrix *Cinverse = gsl_matrix_alloc( numRefs, numRefs );

	int s;
	gsl_permutation *p = gsl_permutation_alloc(numRefs);			//CHANGED ALLOC TO CALLOC

	//initailize the variables to 0
	gsl_vector_set_zero(muVector);
	gsl_vector_set_zero(fVector);
	gsl_vector_set_zero(x);
	gsl_matrix_set_zero(refMatrix);
	gsl_matrix_set_zero(C);

	gsl_permutation_init(p);			//INITIALISED THE PERMUTATION

	readFileVector( muVector, filePath, argv[1] ); //call function which reads in values of first arguement file into the vector mu
	readFileMatrix( refMatrix, filePath, argc, argv ); //obtain matrix of reference vals, (400, 2) in 400 READINGS w/ 2 probes
	//normalize values, vector is already normalised
	normalizeVector( muVector );

	//print vector to csv file for plotting
	FILE *csvfile = fopen( "/mnt/c/Users/jakei/OneDrive - McGill University/Year 2/1BIEN210/Group Project/muVector.csv", "wt" );
	gsl_vector_fprintf( csvfile, muVector, "%f" );
	fclose(csvfile);
	FILE *csvfile2 = fopen( "/mnt/c/Users/jakei/OneDrive - McGill University/Year 2/1BIEN210/Group Project/refVectors.txt", "wt" );
	gsl_matrix_fprintf( csvfile2, refMatrix, "%f" );
	fclose(csvfile2);

	//now we have variable muVector with vals and refMatrix with vals.
	gsl_blas_dgemm( CblasTrans, CblasNoTrans, 1, refMatrix, refMatrix, 0, C);	//matrix C is now A^T A
	gsl_blas_dgemv( CblasTrans, 1, refMatrix, muVector, 0, fVector);	//fVector is 'b'
	gsl_linalg_LU_decomp( C, p, &s);
//	gsl_linalg_LU_solve( C, p, fVector, x );	//can replace with gsl_linalg_LU_svx and delete x vector from program

	gsl_linalg_LU_invert( C, p, Cinverse );		//TODO maybe implement the other way you can solve for basis using inverse on both sides, (A^T A )^-1 A^T A x = (a^t a)^-1 (a^t b) = x
	gsl_blas_dgemv( CblasNoTrans, 1, Cinverse, fVector, 0, x);
/*
	gsl_vector_fprintf( stdout, muVector, "%f" );
	printf("***********************************************\n");
*/
	gsl_vector_fprintf( stdout, x, "%f" );

	gsl_matrix_free(Cinverse);
	gsl_matrix_free(C);
	gsl_matrix_free(refMatrix);
	gsl_vector_free(muVector);
	gsl_vector_free(fVector);
	gsl_vector_free(x);
	return 0;

}

ext: c

max_stars_repo_name: jakeinater/spectralILU

max_issues_repo_name: jakeinater/spectralILU

max_forks_repo_name: jakeinater/spectralILU


text:
#ifndef RL_MDL_NLOPTINVERSEKINEMATICS_H
#define RL_MDL_NLOPTINVERSEKINEMATICS_H

#include <chrono>
#include <nlopt.h>
#include <random>
#include <utility>

#include "InverseKinematics.h"

namespace rl
{
	namespace mdl
	{
		class NloptInverseKinematics : public InverseKinematics
		{
		public:
			NloptInverseKinematics(Kinematic* kinematic);

			virtual ~NloptInverseKinematics();

			bool solve();

			::rl::math::Real delta;

			::std::chrono::nanoseconds duration;

			::rl::math::Real epsilonRotation;

			::rl::math::Real epsilonTranslation;

			double tolerance;

		protected:

		private:
			static void check(const nlopt_result& ret);

			::rl::math::Real error(const ::rl::math::Vector& q);

			static ::rl::math::Real f(unsigned n, const double* x, double* grad, void* data);

			::std::uniform_real_distribution< ::rl::math::Real> randDistribution;

			::std::mt19937 randEngine;
		};
	}
}

#endif // RL_MDL_NLOPTINVERSEKINEMATICS_H

ext: h

max_stars_repo_name: Roboy/rl

max_issues_repo_name: Roboy/rl

max_forks_repo_name: Roboy/rl


text:
/**
 * This file is part of the "libterminal" project
 *   Copyright (c) 2019-2020 Christian Parpart <christian@parpart.family>
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once

#include <terminal/Color.h>
#include <terminal/RenderBuffer.h>
#include <terminal/Screen.h>

#include <terminal_renderer/Atlas.h>
#include <terminal_renderer/BoxDrawingRenderer.h>
#include <terminal_renderer/RenderTarget.h>

#include <text_shaper/font.h>
#include <text_shaper/shaper.h>

#include <crispy/FNV.h>
#include <crispy/LRUCache.h>
#include <crispy/point.h>
#include <crispy/size.h>

#include <unicode/convert.h>
#include <unicode/run_segmenter.h>

#include <gsl/span>
#include <gsl/span_ext>

#include <functional>
#include <list>
#include <memory>
#include <unordered_map>
#include <vector>

namespace terminal::renderer
{

enum class TextStyle
{
    Invalid = 0x00,
    Regular = 0x10,
    Bold = 0x11,
    Italic = 0x12,
    BoldItalic = 0x13,
};

constexpr TextStyle operator|(TextStyle a, TextStyle b) noexcept
{
    return static_cast<TextStyle>(static_cast<unsigned>(a) | static_cast<unsigned>(b));
}

constexpr bool operator<(TextStyle a, TextStyle b) noexcept
{
    return static_cast<unsigned>(a) < static_cast<unsigned>(b);
}

struct TextCacheKey
{
    std::u32string_view text;
    TextStyle style;

    constexpr bool operator<(TextCacheKey const& _rhs) const noexcept
    {
        if (text < _rhs.text)
            return true;

        return text < _rhs.text || style < _rhs.style;
    }

    constexpr bool operator==(TextCacheKey const& _rhs) const noexcept
    {
        return text == _rhs.text && style == _rhs.style;
    }

    constexpr bool operator!=(TextCacheKey const& _rhs) const noexcept { return !(*this == _rhs); }
};

} // end namespace terminal::renderer

namespace std
{
template <>
struct hash<terminal::renderer::TextCacheKey>
{
    size_t operator()(terminal::renderer::TextCacheKey const& _key) const noexcept
    {
        auto fnv = crispy::FNV<char32_t> {};
        return static_cast<size_t>(fnv(fnv.basis(), _key.text, static_cast<char32_t>(_key.style)));
    }
};
} // namespace std

namespace terminal::renderer
{

struct GridMetrics;

enum class TextShapingEngine
{
    OpenShaper, //!< Uses open-source implementation: harfbuzz/freetype/fontconfig
    DWrite,     //!< native platform support: Windows
    CoreText,   //!< native platform support: OS/X
};

enum class FontLocatorEngine
{
    Mock,       //!< mock font locator API
    FontConfig, //!< platform independant font locator API
    DWrite,     //!< native platform support: Windows
    CoreText,   //!< native font locator on OS/X
};

std::unique_ptr<text::font_locator> createFontLocator(FontLocatorEngine _engine);

struct FontDescriptions
{
    double dpiScale = 1.0;
    crispy::Point dpi = { 0, 0 }; // 0 => auto-fill with defaults
    text::font_size size;
    text::font_description regular;
    text::font_description bold;
    text::font_description italic;
    text::font_description boldItalic;
    text::font_description emoji;
    text::render_mode renderMode;
    TextShapingEngine textShapingEngine = TextShapingEngine::OpenShaper;
    FontLocatorEngine fontLocator = FontLocatorEngine::FontConfig;
    bool builtinBoxDrawing = true;
};

inline bool operator==(FontDescriptions const& a, FontDescriptions const& b) noexcept
{
    return a.size.pt == b.size.pt && a.regular == b.regular && a.bold == b.bold && a.italic == b.italic
           && a.boldItalic == b.boldItalic && a.emoji == b.emoji && a.renderMode == b.renderMode;
}

inline bool operator!=(FontDescriptions const& a, FontDescriptions const& b) noexcept
{
    return !(a == b);
}

struct FontKeys
{
    text::font_key regular;
    text::font_key bold;
    text::font_key italic;
    text::font_key boldItalic;
    text::font_key emoji;
};

/// Text Rendering Pipeline
class TextRenderer: public Renderable
{
  public:
    TextRenderer(GridMetrics const& _gridMetrics,
                 text::shaper& _textShaper,
                 FontDescriptions& _fontDescriptions,
                 FontKeys const& _fontKeys);

    void setRenderTarget(RenderTarget& _renderTarget) override;

    void debugCache(std::ostream& _textOutput) const;
    void clearCache() override;

    void updateFontMetrics();

    void setPressure(bool _pressure) noexcept { pressure_ = _pressure; }

    /// Must be invoked before a new terminal frame is rendered.
    void beginFrame();

    /// Renders a given terminal's grid cell that has been
    /// transformed into a RenderCell.
    void renderCell(RenderCell const& _cell);

    /// Must be invoked when rendering the terminal's text has finished for this frame.
    void endFrame();

  private:
    /// Puts a sequence of codepoints that belong to the same grid cell at @p _pos
    /// at the end of the currently filled line.
    void appendCell(gsl::span<char32_t const> _codepoints, TextStyle _style, RGBColor _color);
    text::shape_result const& cachedGlyphPositions();
    text::shape_result requestGlyphPositions();
    text::shape_result shapeRun(unicode::run_segmenter::range const& _run);
    void endSequence();

    void renderRun(crispy::Point _startPos,
                   gsl::span<text::glyph_position const> _glyphPositions,
                   RGBColor _color);

    /// Renders an arbitrary texture.
    void renderTexture(crispy::Point const& _pos,
                       RGBAColor const& _color,
                       atlas::TextureInfo const& _textureInfo);

    // rendering
    //
    struct GlyphMetrics
    {
        ImageSize bitmapSize;  // glyph size in pixels
        crispy::Point bearing; // offset baseline and left to top and left of the glyph's bitmap
    };

    using TextureAtlas = atlas::MetadataTextureAtlas<text::glyph_key, GlyphMetrics>;
    using DataRef = TextureAtlas::DataRef;

    std::optional<DataRef> getTextureInfo(text::glyph_key const& _id,
                                          unicode::PresentationStyle _presentation);

    void renderTexture(crispy::Point const& _pos,
                       RGBAColor const& _color,
                       atlas::TextureInfo const& _textureInfo,
                       GlyphMetrics const& _glyphMetrics,
                       text::glyph_position const& _gpos);

    TextureAtlas* atlasForBitmapFormat(text::bitmap_format _format) noexcept
    {
        switch (_format)
        {
        case text::bitmap_format::alpha_mask: return monochromeAtlas_.get();
        case text::bitmap_format::rgba: return colorAtlas_.get();
        case text::bitmap_format::rgb: return lcdAtlas_.get();
        default: return nullptr; // Should NEVER EVER happen.
        }
    }

    // general properties
    //
    GridMetrics const& gridMetrics_;
    FontDescriptions& fontDescriptions_;
    FontKeys const& fonts_;

    // performance optimizations
    //
    bool pressure_ = false;
    std::unordered_map<text::glyph_key, text::bitmap_format> glyphToTextureMapping_;
    std::list<std::u32string> cacheKeyStorage_;
    crispy::LRUCache<TextCacheKey, text::shape_result> cache_;

    // target surface rendering
    //
    text::shaper& textShaper_;
    std::unique_ptr<TextureAtlas> monochromeAtlas_;
    std::unique_ptr<TextureAtlas> colorAtlas_;
    std::unique_ptr<TextureAtlas> lcdAtlas_;

    // sub-renderer
    //
    BoxDrawingRenderer boxDrawingRenderer_;

    // render states
    TextStyle style_ = TextStyle::Invalid;
    RGBColor color_ {};

    crispy::Point textPosition_;
    std::vector<char32_t> codepoints_;
    std::vector<unsigned> clusters_;
    unsigned cellCount_ = 0;
    bool textStartFound_ = false;
    bool forceCellGroupSplit_ = false;

    // output fields
    //
    std::vector<text::shape_result> shapedLines_;
};

} // namespace terminal::renderer

namespace fmt
{ // {{{
template <>
struct formatter<terminal::renderer::FontDescriptions>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(terminal::renderer::FontDescriptions const& fd, FormatContext& ctx)
    {
        return format_to(ctx.out(),
                         "({}, {}, {}, {}, {}, {})",
                         fd.size,
                         fd.regular,
                         fd.bold,
                         fd.italic,
                         fd.boldItalic,
                         fd.emoji,
                         fd.renderMode);
    }
};

template <>
struct formatter<terminal::renderer::FontLocatorEngine>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(terminal::renderer::FontLocatorEngine value, FormatContext& ctx)
    {
        using terminal::renderer::FontLocatorEngine;
        switch (value)
        {
        case FontLocatorEngine::Mock: return format_to(ctx.out(), "Mock");
        case FontLocatorEngine::FontConfig: return format_to(ctx.out(), "FontConfig");
        case FontLocatorEngine::DWrite: return format_to(ctx.out(), "DirectWrite");
        case FontLocatorEngine::CoreText: return format_to(ctx.out(), "CoreText");
        }
        return format_to(ctx.out(), "UNKNOWN");
    }
};

template <>
struct formatter<terminal::renderer::TextShapingEngine>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(terminal::renderer::TextShapingEngine value, FormatContext& ctx)
    {
        using terminal::renderer::TextShapingEngine;
        switch (value)
        {
        case TextShapingEngine::OpenShaper: return format_to(ctx.out(), "OpenShaper");
        case TextShapingEngine::DWrite: return format_to(ctx.out(), "DirectWrite");
        case TextShapingEngine::CoreText: return format_to(ctx.out(), "CoreText");
        }
        return format_to(ctx.out(), "UNKNOWN");
    }
};

template <>
struct formatter<terminal::renderer::TextCacheKey>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(terminal::renderer::TextCacheKey value, FormatContext& ctx)
    {
        return format_to(ctx.out(), "({}, \"{}\")", value.style, unicode::convert_to<char>(value.text));
    }
};

} // namespace fmt

max_stars_repo_name: contour-terminal/contour

max_issues_repo_name: dankamongmen/contour

max_forks_repo_name: dankamongmen/contour


text:
#ifndef GBNET_GRAPHORNOR
#define GBNET_GRAPHORNOR

#include <iostream>
#include <string>
#include <vector>
#include <utility>
#include <map>
#include <set>
#include <stdio.h>

#include <gsl/gsl_rng.h>

#include "RVNode.h"
#include "XNode.h"
#include "HNodeORNOR.h"
#include "YDataNode.h"
#include "SNode.h"
#include "TNode.h"
#include "ZNode.h"
#include "NodeDictionary.h"
#include "GraphBase.h"
#include "YNoiseNode.h"

namespace gbn
{
    const double SPRIOR[3 * 3] = {0.40, 0.40, 0.20,
                                  0.30, 0.40, 0.30,
                                  0.20, 0.40, 0.40};

    class GraphORNOR: public GraphBase
    {
        private:
        protected:
        public:
            const double XPROB[2] = {0.99, 0.01};
            const double XPROB_ACTIVE[2] = {0.10, 0.90};
            double SPROB[3][3] = {{0.900, 0.090, 0.010},
                                        {0.050, 0.900, 0.050},
                                        {0.010, 0.090, 0.900}};
            double YPROB[3] = {0.05, 0.90, 0.05};
            double YNOISE[3][3] = {{0.945, 0.050, 0.005},
                                   {0.050, 0.900, 0.050},
                                   {0.005, 0.050, 0.945}};
            const double YNOISE_ALPHA[3][3] = {{2000,  100,   10},
                                               { 100, 2000,  100},
                                               {  10,  100, 2000}};

            GraphORNOR (unsigned int);
            virtual ~GraphORNOR ();

            void build_structure (network_t, evidence_dict_t, prior_active_tf_set_t, bool = true,
                                  const double [3 * 3] = SPRIOR, double z_alpha = 25., double z_beta = 25., double z0_alpha = 25., double z0_beta = 25., double t_alpha = 25., double t_beta = 25., bool comp_yprob = true, bool const_params = false,
                                  double t_focus = 2., double t_lmargin = 2., double t_hmargin = 2., double zn_focus = 2., double zn_lmargin = 8., double zn_hmargin = 2.);
    };
}

#endif

max_stars_repo_name: umbibio/gbnet

max_issues_repo_name: umbibio/gbnet

max_forks_repo_name: umbibio/gbnet


text:
//
// beta.h
//
// code using beta distribution
//
// copyright 2020 Peter Andrews
//

#ifndef PAA_BETA_H_
#define PAA_BETA_H_

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

#include <cmath>
#include <mutex>
#include <utility>

#include "error.h"
#include "utility.h"

namespace paa {

class GSLRNG {
 public:
  explicit GSLRNG(const uint64_t seed) : gen{gsl_rng_alloc(gsl_rng_mt19937)} {
    if (!gen) throw Error("Problem creating gsl generator");
    gsl_rng_set(gen, seed);
  }
  ~GSLRNG() { gsl_rng_free(gen); }
  operator gsl_rng * () const { return gen; }
  GSLRNG(const GSLRNG &) = delete;
  GSLRNG & operator=(const GSLRNG &) = delete;
  void lock() { mutex.lock(); }
  void unlock() { mutex.unlock(); }

 private:
  gsl_rng * gen;
  std::mutex mutex{};
};

template<class Val>
class BetaT {
 public:
  using Params = std::pair<Val, Val>;
  BetaT(const Val alpha__, const Val beta__) :
      alpha_{alpha__}, beta_{beta__} {
        if (alpha_ <= 0) throw Error("Nonpositive alpha in Beta distribution");
        if (beta_ <= 0) throw Error("Nonpositive beta in Beta distribution");
      }
  explicit BetaT(const Params params) :
      BetaT{params.first, params.second} {}
  BetaT(const Val mean_, const Val stdev_, const bool) :
      alpha_{get_alpha(mean_, stdev_)},
      beta_{get_beta(mean_, stdev_)} {
        if (stdev_ >= max_stdev(mean_))
          throw Error("Maximum stdev for basic Beta distribution is")
              << max_stdev(mean_);
      }
  Val alpha() const { return alpha_; }
  Val beta() const { return beta_; }
  Val low() const { return 0.0; }
  Val high() const { return 1.0; }
  Val range() const { return 1.0; }
  Val middle() const { return 0.5; }
  static Val get_alpha(const Val mean_, const Val stdev_) {
    return mean_ * (mean_ * (1 - mean_) / sqr(stdev_) - 1);
  }
  static Val get_beta(const Val mean_, const Val stdev_) {
    return (1 - mean_) * (mean_ * (1 - mean_) / sqr(stdev_) - 1);
  }
  static Val max_stdev(const Val mean_) {
    return sqrt(mean_ * (1 - mean_));
  }
  Val mean() const { return alpha_ / (alpha_ + beta_); }
  Val stdev() const {
    return sqrt((alpha_ * beta_) /
                (sqr(alpha_ + beta_) * (alpha_ + beta_ + 1)));
  }
  Val operator()(const Val val) const {
    return gsl_ran_beta_pdf(val, alpha_, beta_);
  }
  Val operator()(gsl_rng * gen) const {
    return gsl_ran_beta(gen, alpha_, beta_);
  }

 private:
  Val alpha_;
  Val beta_;
};
using Beta = BetaT<double>;
using fBeta = BetaT<float>;

template<class Val>
class Beta4T {
 public:
  Beta4T(const Val alpha__, const Val beta__,
        const Val low__, const Val high__) :
      beta_dist{alpha__, beta__}, low_{low__}, high_{high__} {
        if (mean() < low())
          throw Error("Mean less than low in Beta distribution");
        if (mean() > high())
          throw Error("Mean greater than high in Beta distribution");
      }
  Beta4T(const Val mean_, const Val stdev_,
        const Val low__, const Val high__,
        const bool) :
      beta_dist{get_alpha_beta(mean_, stdev_, low__, high__)},
      low_{low__}, high_{high__} {}

  Val alpha() const { return beta_dist.alpha(); }
  Val beta() const { return beta_dist.beta(); }
  Val low() const { return low_; }
  Val high() const { return high_; }
  Val range() const { return (high_ - low_); }
  Val middle() const { return (high_ - low_) / 2; }
  Val mean() const { return low_ + scale * beta_dist.mean(); }
  Val stdev() const { return scale * beta_dist.stdev(); }
  static Val max_stdev(const Val mean_, const Val low__, const Val high__) {
    const Val sf{high__ - low__};
    const Val m{(mean_ - low__) / sf};
    return sf * BetaT<Val>::max_stdev(m);
  }
  Val operator()(const Val val) const {
    return static_cast<Val>(beta_dist((val - low_) / scale)) / scale;
  }
  Val operator()(gsl_rng * gen) const {
    return low_ + scale * static_cast<Val>(beta_dist(gen));
  }

 private:
  static typename BetaT<Val>::Params get_alpha_beta(
      const Val mean_, const Val stdev_,
      const Val low__, const Val high__) {
    const Val sf{high__ - low__};
    const Val m{(mean_ - low__) / sf};
    const Val sd{stdev_ / sf};
    return typename BetaT<Val>::Params{
      BetaT<Val>::get_alpha(m, sd), BetaT<Val>::get_beta(m, sd)};
  }

  BetaT<Val> beta_dist;
  Val low_;
  Val high_;
  Val scale{high_ - low_};
};
using Beta4 = Beta4T<double>;
using fBeta4 = Beta4T<float>;

}  // namespace paa

#endif  // PAA_BETA_H_

max_stars_repo_name: docpaa/mumdex

max_stars_repo_licenses: Zlib MIT

max_issues_repo_name: docpaa/mumdex

max_issues_repo_licenses: Zlib MIT

max_forks_repo_name: docpaa/mumdex

max_forks_repo_licenses: Zlib MIT


text:
/**
 *
 * @file example_cposv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example of solving a linear system with a Cholesky factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:45:21 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

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

int check_solution(int, int, PLASMA_Complex32_t*, int, PLASMA_Complex32_t*, PLASMA_Complex32_t*, int);

int IONE=1;
int ISEED[4] = {0,0,0,1};   /* initial seed for clarnv() */

int main ()
{

    int cores = 2;
    int N     = 10;
    int LDA   = 10;
    int NRHS  = 5;
    int LDB   = 10;
    int info;
    int info_solution;

    PLASMA_Complex32_t *A1   = (PLASMA_Complex32_t *)malloc(LDA*N*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *A2   = (PLASMA_Complex32_t *)malloc(LDA*N*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *B1   = (PLASMA_Complex32_t *)malloc(LDB*NRHS*sizeof(PLASMA_Complex32_t));
    PLASMA_Complex32_t *B2   = (PLASMA_Complex32_t *)malloc(LDB*NRHS*sizeof(PLASMA_Complex32_t));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)){
        printf("Out of Memory \n ");
        return EXIT_SUCCESS;
    }

    /* Plasma Initialize */
    PLASMA_Init(cores);
    printf("-- PLASMA is initialized to run on %d cores. \n",cores);

    /*-------------------------------------------------------------
    *  TESTING CPOSV
    */

    /* Initialize A1 and A2 for Symmetric Positif Matrix (Hessenberg in the complex case) */
    PLASMA_cplghe( (float)N, N, A1, LDA, 51 );
    PLASMA_clacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

    /* Initialize B1 and B2 */
    PLASMA_cplrnt( N, NRHS, B1, LDB, 371 );
    PLASMA_clacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

    /* PLASMA CPOSV */
    info = PLASMA_cposv(PlasmaUpper, N, NRHS, A2, LDA, B2, LDB);

    /* Check the solution */
    info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB);

    if ((info_solution != 0)|(info != 0))
       printf("-- Error in CPOSV example ! \n");
    else
       printf("-- Run of CPOSV example successful ! \n");

    free(A1); free(A2); free(B1); free(B2);

    PLASMA_Finalize();

    return EXIT_SUCCESS;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the solution of the linear system
 */

int check_solution(int N, int NRHS, PLASMA_Complex32_t *A1, int LDA, PLASMA_Complex32_t *B1, PLASMA_Complex32_t *B2, int LDB)
{
    int info_solution;
    float Rnorm, Anorm, Xnorm, Bnorm;
    PLASMA_Complex32_t alpha, beta;
    float *work = (float *)malloc(N*sizeof(float));
    float eps;

    eps = LAPACKE_slamch_work('e');

    alpha = 1.0;
    beta  = -1.0;

    Xnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B2, LDB, work);
    Anorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A1, LDA, work);
    Bnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B1, LDB, work);

    cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, CBLAS_SADDR(alpha), A1, LDA, B2, LDB, CBLAS_SADDR(beta), B1, LDB);
    Rnorm = LAPACKE_clange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B1, LDB, work);

    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n",Rnorm/((Anorm*Xnorm+Bnorm)*N*eps));

    if (Rnorm/((Anorm*Xnorm+Bnorm)*N*eps) > 10.0){
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
     }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }

    free(work);

    return info_solution;
}


text:
#pragma once

#include <gsl\gsl>
#include <winrt\base.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "FullScreenQuad.h"

namespace Rendering
{
	class ComputeShaderMaterial;

	class ComputeShaderDemo final : public Library::DrawableGameComponent
	{
	public:
		ComputeShaderDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		ComputeShaderDemo(const ComputeShaderDemo&) = delete;
		ComputeShaderDemo(ComputeShaderDemo&&) = default;
		ComputeShaderDemo& operator=(const ComputeShaderDemo&) = default;
		ComputeShaderDemo& operator=(ComputeShaderDemo&&) = default;
		~ComputeShaderDemo();

		bool AnimationEnabled() const;
		void SetAnimationEnabled(bool enabled);

		float BlueColor() const;
		void SetBlueColor(float blueColor);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		std::unique_ptr<ComputeShaderMaterial> mMaterial;
		Library::FullScreenQuad mFullScreenQuad;
		winrt::com_ptr<ID3D11ShaderResourceView> mColorTexture;
		bool mAnimationEnabled{ true };
	};
}


text:
/*****************************************************************************
 *                                                                           *
 *  Copyright 2018 Rice University                                           *
 *                                                                           *
 *  Licensed under the Apache License, Version 2.0 (the "License");          *
 *  you may not use this file except in compliance with the License.         *
 *  You may obtain a copy of the License at                                  *
 *                                                                           *
 *      http://www.apache.org/licenses/LICENSE-2.0                           *
 *                                                                           *
 *  Unless required by applicable law or agreed to in writing, software      *
 *  distributed under the License is distributed on an "AS IS" BASIS,        *
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
 *  See the License for the specific language governing permissions and      *
 *  limitations under the License.                                           *
 *                                                                           *
 *****************************************************************************/

#ifndef LDA_DOC_TOPIC_PROB_SELECT_H
#define LDA_DOC_TOPIC_PROB_SELECT_H

#include "Lambda.h"
#include "LambdaCreationFunctions.h"
#include "SelectionComp.h"
#include "PDBVector.h"
#include "IntIntVectorPair.h"
#include "IntDoubleVectorPair.h"
#include "DocAssignment.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <random>
#include <gsl/gsl_vector.h>

/* This class samples the topic probability for the document */
using namespace pdb;
class LDADocTopicProbSelection : public SelectionComp<IntDoubleVectorPair, DocAssignment> {

private:
    Vector<double> prior;
    Handle<Vector<char>> myMem;
    int topicNum;

public:
    ENABLE_DEEP_COPY

    LDADocTopicProbSelection() {}

    LDADocTopicProbSelection(Vector<double>& fromPrior) {

        this->prior = fromPrior;
        this->topicNum = fromPrior.size();

	/* Set up the random number generator */
        gsl_rng* src = gsl_rng_alloc(gsl_rng_mt19937);
        std::random_device rd;
        std::mt19937 gen(rd());
        gsl_rng_set(src, gen());

        int spaceNeeded = sizeof(gsl_rng) + src->type->size;
        myMem = makeObject<Vector<char>>(spaceNeeded, spaceNeeded);

        memcpy(myMem->c_ptr(), src, sizeof(gsl_rng));
        memcpy(myMem->c_ptr() + sizeof(gsl_rng), src->state, src->type->size);

        gsl_rng_free(src);
    }

    Lambda<bool> getSelection(Handle<DocAssignment> checkMe) override {
        return makeLambda(checkMe, [](Handle<DocAssignment>& checkMe) { return true; });
    }

    Lambda<Handle<IntDoubleVectorPair>> getProjection(Handle<DocAssignment> checkMe) override {
        return makeLambda(checkMe, [&](Handle<DocAssignment>& checkMe) {

            gsl_rng* rng = getRng();
            Vector<unsigned>& topicCount = checkMe->getVector();
            Handle<Vector<double>> mySamples = makeObject<Vector<double>>(topicNum, topicNum);
            double* totalProb = new double[topicNum];
            for (int i = 0; i < topicNum; i++) {
                totalProb[i] = (this->prior)[i] + topicCount[i];
            }

	    /* Sample the topic probability */
            gsl_ran_dirichlet(rng, topicNum, totalProb, mySamples->c_ptr());

            Handle<IntDoubleVectorPair> result =
                makeObject<IntDoubleVectorPair>(checkMe->getKey(), mySamples);

            delete[] totalProb;

            return result;
        });
    }

    gsl_rng* getRng() {
        gsl_rng* dst = (gsl_rng*)myMem->c_ptr();
        dst->state = (void*)(myMem->c_ptr() + sizeof(gsl_rng));
        dst->type = gsl_rng_mt19937;
        return dst;
    }
};

#endif

max_stars_repo_name: dimitrijejankov/pdb

max_issues_repo_name: dcbdan/plinycompute

max_forks_repo_name: dcbdan/plinycompute


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_deriv.h>

static void
central_deriv (const gsl_function * f, double x, double h,
               double *result, double *abserr_round, double *abserr_trunc)
{
  /* Compute the derivative using the 5-point rule (x-h, x-h/2, x,
     x+h/2, x+h). Note that the central point is not used.

     Compute the error using the difference between the 5-point and
     the 3-point rule (x-h,x,x+h). Again the central point is not
     used. */

  double fm1 = GSL_FN_EVAL (f, x - h);
  double fp1 = GSL_FN_EVAL (f, x + h);

  double fmh = GSL_FN_EVAL (f, x - h / 2);
  double fph = GSL_FN_EVAL (f, x + h / 2);

  double r3 = 0.5 * (fp1 - fm1);
  double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3;

  double e3 = (fabs (fp1) + fabs (fm1)) * GSL_DBL_EPSILON;
  double e5 = 2.0 * (fabs (fph) + fabs (fmh)) * GSL_DBL_EPSILON + e3;

  /* The next term is due to finite precision in x+h = O (eps * x) */

  double dy = GSL_MAX (fabs (r3 / h), fabs (r5 / h)) *(fabs (x) / h) * GSL_DBL_EPSILON;

  /* The truncation error in the r5 approximation itself is O(h^4).
     However, for safety, we estimate the error from r5-r3, which is
     O(h^2).  By scaling h we will minimise this estimated error, not
     the actual truncation error in r5. */

  *result = r5 / h;
  *abserr_trunc = fabs ((r5 - r3) / h); /* Estimated truncation error O(h^2) */
  *abserr_round = fabs (e5 / h) + dy;   /* Rounding error (cancellations) */
}

int
gsl_deriv_central (const gsl_function * f, double x, double h,
                   double *result, double *abserr)
{
  double r_0, round, trunc, error;
  central_deriv (f, x, h, &r_0, &round, &trunc);
  error = round + trunc;

  if (round < trunc && (round > 0 && trunc > 0))
    {
      double r_opt, round_opt, trunc_opt, error_opt;

      /* Compute an optimised stepsize to minimize the total error,
         using the scaling of the truncation error (O(h^2)) and
         rounding error (O(1/h)). */

      double h_opt = h * pow (round / (2.0 * trunc), 1.0 / 3.0);
      central_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt);
      error_opt = round_opt + trunc_opt;

      /* Check that the new error is smaller, and that the new derivative
         is consistent with the error bounds of the original estimate. */

      if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error)
        {
          r_0 = r_opt;
          error = error_opt;
        }
    }

  *result = r_0;
  *abserr = error;

  return GSL_SUCCESS;
}

static void
forward_deriv (const gsl_function * f, double x, double h,
               double *result, double *abserr_round, double *abserr_trunc)
{
  /* Compute the derivative using the 4-point rule (x+h/4, x+h/2,
     x+3h/4, x+h).

     Compute the error using the difference between the 4-point and
     the 2-point rule (x+h/2,x+h).  */

  double f1 = GSL_FN_EVAL (f, x + h / 4.0);
  double f2 = GSL_FN_EVAL (f, x + h / 2.0);
  double f3 = GSL_FN_EVAL (f, x + (3.0 / 4.0) * h);
  double f4 = GSL_FN_EVAL (f, x + h);

  double r2 = 2.0*(f4 - f2);
  double r4 = (22.0 / 3.0) * (f4 - f3) - (62.0 / 3.0) * (f3 - f2) +
    (52.0 / 3.0) * (f2 - f1);

  /* Estimate the rounding error for r4 */

  double e4 = 2 * 20.67 * (fabs (f4) + fabs (f3) + fabs (f2) + fabs (f1)) * GSL_DBL_EPSILON;

  /* The next term is due to finite precision in x+h = O (eps * x) */

  double dy = GSL_MAX (fabs (r2 / h), fabs (r4 / h)) * fabs (x / h) * GSL_DBL_EPSILON;

  /* The truncation error in the r4 approximation itself is O(h^3).
     However, for safety, we estimate the error from r4-r2, which is
     O(h).  By scaling h we will minimise this estimated error, not
     the actual truncation error in r4. */

  *result = r4 / h;
  *abserr_trunc = fabs ((r4 - r2) / h); /* Estimated truncation error O(h) */
  *abserr_round = fabs (e4 / h) + dy;
}

int
gsl_deriv_forward (const gsl_function * f, double x, double h,
                   double *result, double *abserr)
{
  double r_0, round, trunc, error;
  forward_deriv (f, x, h, &r_0, &round, &trunc);
  error = round + trunc;

  if (round < trunc && (round > 0 && trunc > 0))
    {
      double r_opt, round_opt, trunc_opt, error_opt;

      /* Compute an optimised stepsize to minimize the total error,
         using the scaling of the estimated truncation error (O(h)) and
         rounding error (O(1/h)). */

      double h_opt = h * pow (round / (trunc), 1.0 / 2.0);
      forward_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt);
      error_opt = round_opt + trunc_opt;

      /* Check that the new error is smaller, and that the new derivative
         is consistent with the error bounds of the original estimate. */

      if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error)
        {
          r_0 = r_opt;
          error = error_opt;
        }
    }

  *result = r_0;
  *abserr = error;

  return GSL_SUCCESS;
}

int
gsl_deriv_backward (const gsl_function * f, double x, double h,
                    double *result, double *abserr)
{
  return gsl_deriv_forward (f, x, -h, result, abserr);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <gsl/gsl_chebyshev.h>
const struct _GSLMethods
chebyshev_f   = { (void_m_t) gsl_cheb_free,
		  /* gsl_multimin_fminimizer_restart */  (void_m_t) NULL,
		  NULL, NULL};

static int
PyGSL_cheb_init(PyGSL_Solver *self, PyObject *args, PyObject *kw)
{
     ;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/*
Copyright (c) 2015, Patrick Weltevrede
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#define _FILE_OFFSET_BITS 64
#define _USE_LARGEFILE 1
#define _LARGEFILE_SOURCE 1
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_sort.h>
#include "psrsalsa.h"
#define MaxNrOnpulseRegions 10
int NrOnpulseRegions, OnPulseRegion[MaxNrOnpulseRegions][2];
void PlotProfile(int NrBins, float *Ipulse, char *xlabel, char *ylabel, char *title, int Highlight, int color, int clearPage);
void ShiftProfile(int shift, int NrBins, float *Iprofile, float *outputProfile);
int main(int argc, char **argv)
{
  char output_fname[1000], PlotDevice[100], singlechar, *inputname;
  int circularShift, noinput, onlyI, memsave, currentfilenumber, dummy_int;
  int shift, bin1, bin2, subintWritten, poladd, freqadd;
  long i, j, pol, fchan, nsub, binnr, nrinputfiles, subintslost, currentOutputSubint, sumNsub, curNrInsubint;
  float *Iprofile, *Iprofile_firstfile, *shiftedProfile, *subint, x, y, *float_ptr, *float_ptr2;
  datafile_definition **fin;
  datafile_definition fout, clone;
  psrsalsaApplication application;
  initApplication(&application, "padd", "[options] inputfiles");
  application.switch_blocksize = 1;
  application.switch_verbose = 1;
  application.switch_debug = 1;
  application.switch_filelist = 1;
  application.switch_iformat = 1;
  application.switch_oformat = 1;
  application.switch_formatlist = 1;
  application.switch_fscr = 1;
  application.switch_FSCR = 1;
  application.switch_nocounters = 1;
  application.switch_changeRefFreq = 1;
  application.oformat = FITS_format;
  strcpy(PlotDevice, "?");
  onlyI = 0;
  strcpy(output_fname, "addedfile.gg");
  NrOnpulseRegions = 0;
  circularShift = 0;
  noinput = 0;
  shift = 0;
  memsave = 0;
  sumNsub = 1;
  poladd = 0;
  freqadd = 0;
  x = y = singlechar = 0;
  if(argc < 2) {
    printf("Program to add data files together. Usage:\n\n");
    printApplicationHelp(&application);
    printf("Optional options:\n");
    printf("-w                      Output name. Default is \"%s\"\n", output_fname);
    printf("-I                      Only process the first polarization channel\n");
    printf("-c                      Turn on circular shifting (so that no subint is lost).\n");
    printf("                        The first and last subintegrations are spilling over\n");
    printf("                        into each other.\n");
    printf("-n                      No graphical input, circular shifting by this number of\n");
    printf("                        bins, so this option implies -c. So -n 0 results in a\n");
    printf("                        simple concatenation of the input files.\n");
    printf("-nsub                   This number of subintegrations are summed before being\n");
    printf("                        written out\n");
    printf("-poladd                 The input files are to be interpretted as separate\n");
    printf("                        polarization channels (option implies -n 0).\n");
    printf("-freqadd                The input files are to be interpretted as separate\n");
    printf("                        frequency bands (option implies -n 0).\n");
    printf("-memsave                Only one full input data-set exists in memory at a time,\n");
    printf("                        but every input file will be opened twice.\n");
    printf("\n");
    printCitationInfo();
    terminateApplication(&application);
    return 0;
  }else {
    for(i = 1; i < argc; i++) {
      dummy_int = i;
      if(processCommandLine(&application, argc, argv, &dummy_int)) {
 i = dummy_int;
      }else if(strcmp(argv[i], "-w") == 0 || strcmp(argv[i], "-W") == 0) {
 strcpy(output_fname,argv[i+1]);
        i++;
      }else if(strcmp(argv[i], "-memsave") == 0) {
 memsave = 1;
      }else if(strcmp(argv[i], "-c") == 0 || strcmp(argv[i], "-C") == 0) {
 circularShift = 1;
      }else if(strcmp(argv[i], "-I") == 0) {
 onlyI = 1;
      }else if(strcmp(argv[i], "-poladd") == 0) {
 poladd = 1;
 circularShift = 1;
 noinput = 1;
 shift = 0;
      }else if(strcmp(argv[i], "-freqadd") == 0) {
 freqadd = 1;
 circularShift = 1;
 noinput = 1;
 shift = 0;
      }else if(strcmp(argv[i], "-n") == 0) {
 circularShift = 1;
 noinput = 1;
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%d", &shift, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR padd: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 i++;
      }else if(strcmp(argv[i], "-nsub") == 0) {
 if(parse_command_string(application.verbose_state, argc, argv, i+1, 0, -1, "%ld", &sumNsub, NULL) == 0) {
   printerror(application.verbose_state.debug, "ERROR padd: Cannot parse '%s' option.", argv[i]);
   return 0;
 }
 if(sumNsub < 1) {
   fflush(stdout);
   printerror(application.verbose_state.debug, "ERROR padd: Cannot parse option %s, expected one integer number > 1", argv[i]);
   return 0;
 }
 i++;
      }else {
 if(argv[i][0] == '-') {
   printerror(application.verbose_state.debug, "ERROR padd: Unknown option %s. Run padd without command-line options to get help.", argv[i]);
   terminateApplication(&application);
   return 0;
 }else {
   if(applicationAddFilename(i, application.verbose_state) == 0)
     return 0;
 }
      }
    }
  }
  if(applicationFilenameList_checkConsecutive(argv, application.verbose_state) == 0) {
    return 0;
  }
  nrinputfiles = numberInApplicationFilenameList(&application, argv, application.verbose_state);
  if(nrinputfiles < 2) {
    printerror(application.verbose_state.debug, "ERROR padd: Need at least two input files");
    return 0;
  }
  if(freqadd && poladd) {
    printerror(application.verbose_state.debug, "ERROR padd: Cannot use the -freqadd and -poladd flag simultaneously.");
    return 0;
  }
  if(sumNsub != 1 && poladd) {
    printerror(application.verbose_state.debug, "ERROR padd: Cannot use the -nsub and -poladd flag simultaneously.");
    return 0;
  }
  if(sumNsub != 1 && freqadd) {
    printerror(application.verbose_state.debug, "ERROR padd: Cannot use the -nsub and -freqadd flag simultaneously.");
    return 0;
  }
  if(circularShift != 1 && poladd) {
    printerror(application.verbose_state.debug, "ERROR padd: You must use the -c option together with -poladd.");
    return 0;
  }
  if(circularShift != 1 && freqadd) {
    printerror(application.verbose_state.debug, "ERROR padd: You must use the -c option together with -freqadd.");
    return 0;
  }
  if((noinput != 1 || shift != 0) && poladd) {
    printerror(application.verbose_state.debug, "ERROR padd: You must use the -n 0 option together with -poladd.");
    return 0;
  }
  if((noinput != 1 || shift != 0) && freqadd) {
    printerror(application.verbose_state.debug, "ERROR padd: You must use the -n 0 option together with -poladd.");
    return 0;
  }
  fin = malloc(nrinputfiles*sizeof(datafile_definition *));
  if(fin == NULL) {
    printerror(application.verbose_state.debug, "ERROR padd: Memory allocation error");
    return 0;
  }
  for(i = 0; i < nrinputfiles; i++) {
    fin[i] = malloc(sizeof(datafile_definition));
    if(fin[i] == NULL) {
      printerror(application.verbose_state.debug, "ERROR padd: Memory allocation error");
      return 0;
    }
  }
  currentfilenumber = 0;
  while((inputname = getNextFilenameFromList(&application, argv, application.verbose_state)) != NULL) {
    verbose_definition verbose2;
    copyVerboseState(application.verbose_state, &verbose2);
    verbose2.indent = application.verbose_state.indent + 2;
    if(memsave == 0 || (currentfilenumber == 0 && noinput == 0)) {
      if(currentfilenumber == 0) {
 printf("Read in input files:\n");
      }
      if(openPSRData(fin[currentfilenumber], inputname, application.iformat, 0, 1, 0, verbose2) == 0) {
 printerror(application.verbose_state.debug, "ERROR padd: Cannot open %s\n", inputname);
 return 0;
      }
      if(currentfilenumber == 0) {
 for(i = 1; i < argc; i++) {
   if(strcmp(argv[i], "-header") == 0) {
     fflush(stdout);
     printwarning(application.verbose_state.debug, "WARNING: If using the -header option, be aware it applied BEFORE the preprocessing.");
   }
 }
      }
      if(preprocessApplication(&application, fin[currentfilenumber]) == 0) {
 printerror(application.verbose_state.debug, "ERROR padd: preprocess option failed on file %s\n", inputname);
 return 0;
      }
    }else {
      if(currentfilenumber == 0)
 printf("Read in headers of input files:\n");
      if(openPSRData(fin[currentfilenumber], inputname, application.iformat, 0, 0, 0, verbose2) == 0) {
 printerror(application.verbose_state.debug, "ERROR padd: Cannot open %s\n", inputname);
 return 0;
      }
      if(readHeaderPSRData(fin[currentfilenumber], 1, 0, verbose2) == 0) {
 printerror(application.verbose_state.debug, "ERROR padd: Cannot read header of file %s\n", inputname);
 return 0;
      }
    }
    if(currentfilenumber == 0 && noinput == 0) {
      Iprofile_firstfile = (float *)malloc(fin[0]->NrPols*fin[0]->NrBins*sizeof(float));
      shiftedProfile = (float *)malloc(fin[0]->NrPols*fin[0]->NrBins*sizeof(float));
      if(Iprofile_firstfile == NULL || shiftedProfile == NULL) {
 printerror(application.verbose_state.debug, "ERROR padd: Cannot allocate memory.");
 return 0;
      }
      if(read_profilePSRData(*fin[0], Iprofile_firstfile, NULL, 0, application.verbose_state) != 1) {
 printerror(application.verbose_state.debug, "ERROR padd: Reading pulse profile of first input file failed.");
 return 0;
      }
    }
    if(memsave) {
      if(closePSRData(fin[currentfilenumber], 1, application.verbose_state) != 0) {
 printerror(application.verbose_state.debug, "ERROR padd: Closing file %s failed\n", inputname);
 return 0;
      }
    }
    if(fin[currentfilenumber]->NrBins != fin[0]->NrBins) {
      printerror(application.verbose_state.debug, "ERROR padd: Number of pulse longitude bins not equal in input files.");
      return 0;
    }
    if(fin[0]->NrPols != fin[currentfilenumber]->NrPols && onlyI == 0) {
      printerror(application.verbose_state.debug, "ERROR padd: Number of polarization channels are different. Maybe you want to use the -I option?");
      return 0;
    }
    if(freqadd == 0) {
      if(fin[currentfilenumber]->NrFreqChan != fin[0]->NrFreqChan) {
 printerror(application.verbose_state.debug, "ERROR padd: Number of frequency channels not equal in input files.");
 return 0;
      }
    }
    if(poladd || freqadd) {
      if(fin[0]->NrSubints != fin[currentfilenumber]->NrSubints) {
 printerror(application.verbose_state.debug, "ERROR padd: Number of subints are different. This is not allowed with the -poladd or -freqadd options.");
 return 0;
      }
    }
    if(poladd) {
      if(fin[currentfilenumber]->NrPols != 1) {
 printerror(application.verbose_state.debug, "ERROR padd: Number of polarization channels in input files should be 1 if using the -poladd option.");
 return 0;
      }
    }
    if(fin[currentfilenumber]->isDeDisp != fin[0]->isDeDisp) {
      printerror(application.verbose_state.debug, "ERROR padd: Dedispersion state is not equal in input files.");
      return 0;
    }
    if(fin[currentfilenumber]->isDeFarad != fin[0]->isDeFarad) {
      printerror(application.verbose_state.debug, "ERROR padd: De-Faraday rotation state is not equal in input files.");
      return 0;
    }
    if(fin[currentfilenumber]->isDePar != fin[0]->isDePar) {
      printerror(application.verbose_state.debug, "ERROR padd: Parallactic angle state is not equal in input files.");
      return 0;
    }
    currentfilenumber++;
  }
  printf("Reading of input files done\n");
  unsigned long *sort_indx;
  sort_indx = (unsigned long *)malloc(nrinputfiles*sizeof(unsigned long));
  if(sort_indx == NULL) {
    printerror(application.verbose_state.debug, "ERROR padd: Memory allocation error\n");
    return 0;
  }
    for(currentfilenumber = 0; currentfilenumber < nrinputfiles; currentfilenumber++) {
      sort_indx[currentfilenumber] = currentfilenumber;
    }
  cleanPSRData(&fout, application.verbose_state);
  copy_params_PSRData(*(fin[0]), &fout, application.verbose_state);
  if(onlyI) {
    fout.NrPols = 1;
  }
  fout.format = application.oformat;
  fout.NrSubints = 0;
  subintslost = 0;
  if(poladd) {
    fout.NrPols = nrinputfiles;
    fout.NrSubints = fin[0]->NrSubints;
  }else if(freqadd) {
    fout.NrSubints = fin[0]->NrSubints;
    fout.NrFreqChan = 0;
    for(currentfilenumber = 0; currentfilenumber < nrinputfiles; currentfilenumber++) {
      fout.NrFreqChan += fin[currentfilenumber]->NrFreqChan;
    }
    fout.freqMode = FREQMODE_FREQTABLE;
    if(fout.freqlabel_list != NULL) {
      free(fout.freqlabel_list);
    }
    fout.freqlabel_list = (double *)malloc(fout.NrFreqChan*fout.NrSubints*sizeof(double));
    if(fout.freqlabel_list == NULL) {
      printerror(application.verbose_state.debug, "ERROR padd: Memory allocation error");
      return 0;
    }
    long currentOutputChan = 0;
    for(i = 0; i < nrinputfiles; i++) {
      currentfilenumber = sort_indx[i];
      for(fchan = 0; fchan < fin[currentfilenumber]->NrFreqChan; fchan++) {
 for(nsub = 0; nsub < fin[currentfilenumber]->NrSubints; nsub++) {
   double freq;
   freq = get_weighted_channel_freq(*(fin[currentfilenumber]), nsub, fchan, application.verbose_state);
   if(set_weighted_channel_freq(&fout, nsub, currentOutputChan, freq, application.verbose_state) == 0) {
     printerror(application.verbose_state.debug, "ERROR padd: Constructing frequency table failed");
     return 0;
   }
 }
 currentOutputChan++;
      }
    }
  }else {
    for(i = 0; i < nrinputfiles; i++) {
      fout.NrSubints += fin[i]->NrSubints;
      if(circularShift == 0) {
 fout.NrSubints -= 1;
 subintslost += 1;
      }
    }
  }
  if(freqadd || poladd) {
    printf("\nOutput data will be %ld subints, %ld phase bins %ld frequency channels and %ld polarizations.\n", fout.NrSubints, fout.NrBins, fout.NrFreqChan, fout.NrPols);
  }else {
    printf("\nInput data contains %ld subints, %ld phase bins %ld frequency channels and %ld polarizations.\n", fout.NrSubints+subintslost, fout.NrBins, fout.NrFreqChan, fout.NrPols);
    printf("%ld subints are lost because of the alignment of input data", subintslost);
    if(subintslost > 0)
      printf(" (consider using -c option)");
  }
  fout.tsubMode = TSUBMODE_TSUBLIST;
  if(fout.tsub_list != NULL) {
    free(fout.tsub_list);
  }
  fout.tsub_list = (double *)malloc(fout.NrSubints*sizeof(double));
  if(fout.tsub_list == NULL) {
    printerror(application.verbose_state.debug, "ERROR padd: Memory allocation error");
    return 0;
  }
  currentOutputSubint = 0;
  fout.tsub_list[0] = 0;
  subintWritten = 0;
  curNrInsubint = 0;
  if(poladd || freqadd) {
    for(nsub = 0; nsub < fin[0]->NrSubints; nsub++) {
      fout.tsub_list[nsub] = get_tsub(*(fin[0]), nsub, application.verbose_state);
    }
  }else {
    for(i = 0; i < nrinputfiles; i++) {
      currentfilenumber = sort_indx[i];
      for(nsub = 0; nsub < fin[currentfilenumber]->NrSubints; nsub++) {
 if(currentOutputSubint < fout.NrSubints)
   fout.tsub_list[currentOutputSubint] += get_tsub(*(fin[currentfilenumber]), nsub, application.verbose_state);
 if(sumNsub == 1) {
   subintWritten = 1;
 }else if(curNrInsubint == sumNsub - 1) {
   subintWritten = 1;
 }
 curNrInsubint++;
 if(subintWritten) {
   subintWritten = 0;
   curNrInsubint = 0;
   currentOutputSubint++;
   if(currentOutputSubint < fout.NrSubints)
     fout.tsub_list[currentOutputSubint] = 0;
 }
      }
    }
  }
  if(freqadd == 0 && poladd == 0) {
    dummy_int = fout.NrSubints % sumNsub;
    fout.NrSubints = fout.NrSubints/sumNsub;
    printf("\nOutput data will contain %ld subints", fout.NrSubints);
    if(sumNsub > 1)
      printf(" after summing every %ld input subints", sumNsub);
    if(dummy_int)
      printf(" (%d input subints lost because of incomplete last subint)", dummy_int);
    printf("\n\n");
  }
  if(fout.gentype == GENTYPE_PULSESTACK && sumNsub != 1) {
    if(fout.NrSubints != 1) {
      fout.gentype = GENTYPE_SUBINTEGRATIONS;
    }else {
      fout.gentype = GENTYPE_PROFILE;
    }
  }
  if(openPSRData(&fout, output_fname, fout.format, 1, 0, 0, application.verbose_state) == 0) {
    printerror(application.verbose_state.debug, "ERROR padd: Cannot open %s", output_fname);
    return 0;
  }
  if(writeHeaderPSRData(&fout, argc, argv, application.history_cmd_only, application.verbose_state) != 1) {
    printerror(application.verbose_state.debug, "ERROR padd: Cannot write header to %s", output_fname);
    return 0;
  }
  Iprofile = (float *)malloc(fout.NrPols*fout.NrBins*sizeof(float));
  if(Iprofile == NULL) {
    printerror(application.verbose_state.debug, "ERROR padd: Cannot allocate memory.");
    return 0;
  }
  if(sumNsub > 1) {
    subint = (float *)malloc(fout.NrPols*fout.NrBins*fout.NrFreqChan*sizeof(float));
    if(subint == NULL) {
      printerror(application.verbose_state.debug, "ERROR padd: Cannot allocate memory.");
      return 0;
    }
  }
  if(noinput == 0) {
    ppgopen(PlotDevice);
    ppgask(0);
    ppgslw(1);
  }
  currentfilenumber = 0;
  currentOutputSubint = 0;
  curNrInsubint = 0;
  subintWritten = 0;
  rewindFilenameList(&application);
  long currentfilenumber_index;
  long fchan_offset = 0;
  while((inputname = getNextFilenameFromList(&application, argv, application.verbose_state)) != NULL) {
    currentfilenumber_index = sort_indx[currentfilenumber];
    if(memsave) {
      closePSRData(fin[currentfilenumber_index], 0, application.verbose_state);
      if(openPSRData(fin[currentfilenumber_index], inputname, application.iformat, 0, 1, 0, application.verbose_state) == 0) {
 printerror(application.verbose_state.debug, "ERROR padd: Cannot open %s\n", inputname);
 return 0;
      }
      if(currentfilenumber == 0) {
 for(i = 1; i < argc; i++) {
   if(strcmp(argv[i], "-header") == 0) {
     fflush(stdout);
     printwarning(application.verbose_state.debug, "WARNING: If using the -header option, be aware it applied BEFORE the preprocessing.");
   }
 }
      }
      if(preprocessApplication(&application, fin[currentfilenumber_index]) == 0) {
 printerror(application.verbose_state.debug, "ERROR padd: preprocess option failed on file %s\n", inputname);
 return 0;
      }
    }
    if(shift >= fout.NrBins)
      shift -= fout.NrBins;
    if(shift < 0)
      shift += fout.NrBins;
    if(noinput == 0 && currentfilenumber != 0) {
      if(read_profilePSRData(*fin[currentfilenumber_index], Iprofile, NULL, 0, application.verbose_state) != 1) {
 printerror(application.verbose_state.debug, "ERROR padd: Reading pulse profile failed.");
 return 0;
      }
      do {
 ShiftProfile(shift, fout.NrBins, Iprofile, shiftedProfile);
 PlotProfile(fout.NrBins, Iprofile_firstfile, "Bin number", "Intensity", "Click to shift profile, press s to stop", 0, 1, 1);
 PlotProfile(fout.NrBins, shiftedProfile, "", "", "", 0, 2, 0);
 j = 0;
 do {
   if(j == 0)
     ppgband(0, 0, 0.0, 0.0, &x, &y, &singlechar);
   else
     ppgband(4, 0, bin1, 0.0, &x, &y, &singlechar);
   if(singlechar == 65) {
     if(j == 0)
       bin1 = x;
     else
       bin2 = x;
     j++;
   }else if(singlechar == 115 || singlechar ==83 ) {
     j = 10;
   }
 }while(j < 2);
 if(j < 10) {
   shift += bin2 - bin1;
   if(shift >= fout.NrBins)
     shift -= fout.NrBins;
   if(shift < 0)
     shift += fout.NrBins;
 }
      }while(j < 10);
    }
    if(shift != 0) {
      if(continuous_shift(*fin[currentfilenumber_index], &clone, shift, circularShift, "padd", MEMORY_format, 0, NULL, application.verbose_state, application.verbose_state.debug) != 1) {
 printerror(application.verbose_state.debug, "ERROR padd: circular shift failed.");
      }
      swap_orig_clone(fin[currentfilenumber_index], &clone, application.verbose_state);
    }
    long nrPulsesInCurFile;
    nrPulsesInCurFile = fin[currentfilenumber_index]->NrSubints;
    if(shift == 0 && circularShift == 0)
      nrPulsesInCurFile -= 1;
    for(nsub = 0; nsub < nrPulsesInCurFile; nsub++) {
      int nrpolsinloop;
      nrpolsinloop = fout.NrPols;
      if(poladd) {
 nrpolsinloop = 1;
      }
      for(pol = 0; pol < nrpolsinloop; pol++) {
 for(fchan = 0; fchan < fin[currentfilenumber_index]->NrFreqChan; fchan++) {
   if(readPulsePSRData(fin[currentfilenumber_index], nsub, pol, fchan, 0, fin[currentfilenumber_index]->NrBins, Iprofile, application.verbose_state) != 1) {
     printerror(application.verbose_state.debug, "ERROR padd: Read error");
     return 0;
   }
   if(sumNsub == 1) {
     if(poladd == 0 && freqadd == 0) {
  if(writePulsePSRData(&fout, currentOutputSubint, pol, fchan, 0, fout.NrBins, Iprofile, application.verbose_state) != 1) {
    printerror(application.verbose_state.debug, "ERROR padd: Write error");
    return 0;
  }
     }else {
       if(poladd) {
  if(writePulsePSRData(&fout, nsub, currentfilenumber, fchan, 0, fout.NrBins, Iprofile, application.verbose_state) != 1) {
    printerror(application.verbose_state.debug, "ERROR padd: Write error");
    return 0;
  }
       }else if(freqadd) {
  if(writePulsePSRData(&fout, nsub, pol, fchan+fchan_offset, 0, fout.NrBins, Iprofile, application.verbose_state) != 1) {
    printerror(application.verbose_state.debug, "ERROR padd: Write error");
    return 0;
  }
       }
     }
     subintWritten = 1;
   }else {
     float_ptr = &subint[fout.NrBins*(pol+fout.NrPols*fchan)];
     float_ptr2 = Iprofile;
     if(curNrInsubint == 0) {
       for(binnr = 0; binnr < fout.NrBins; binnr++) {
  *float_ptr = *float_ptr2;
  float_ptr++;
  float_ptr2++;
       }
     }else {
       for(binnr = 0; binnr < fout.NrBins; binnr++) {
  *float_ptr += *float_ptr2;
  float_ptr++;
  float_ptr2++;
       }
     }
     if(curNrInsubint == sumNsub - 1) {
       if(writePulsePSRData(&fout, currentOutputSubint, pol, fchan, 0, fout.NrBins, float_ptr-fout.NrBins, application.verbose_state) != 1) {
  printerror(application.verbose_state.debug, "ERROR padd: Write error");
  return 0;
       }
       subintWritten = 1;
     }
   }
 }
      }
      curNrInsubint++;
      if(subintWritten) {
 subintWritten = 0;
 curNrInsubint = 0;
 currentOutputSubint++;
      }
      if(application.verbose_state.nocounters == 0) {
 printf("Processing input file %d: %.1f%%     \r", currentfilenumber+1, (100.0*(nsub+1))/(float)(fin[currentfilenumber_index]->NrSubints));
 fflush(stdout);
      }
    }
    if(freqadd) {
      fchan_offset += fin[currentfilenumber_index]->NrFreqChan;
    }
    if(application.verbose_state.nocounters == 0) {
      printf("Processing file %d is done (%s).                                \n", currentfilenumber+1, fin[currentfilenumber_index]->filename);
    }
    closePSRData(fin[currentfilenumber_index], 0, application.verbose_state);
    currentfilenumber++;
  }
  closePSRData(&fout, 0, application.verbose_state);
  if(noinput == 0)
    ppgend();
  free(Iprofile);
  if(noinput == 0) {
    free(shiftedProfile);
    free(Iprofile_firstfile);
  }
  if(sumNsub > 1)
    free(subint);
  for(i = 0; i < nrinputfiles; i++) {
    free(fin[i]);
  }
  free(fin);
  free(sort_indx);
  terminateApplication(&application);
  return 0;
}
int CheckOnPulse(int bin, int NrRegions, int Regions[MaxNrOnpulseRegions][2])
{
  int i;
  for(i = 0; i < NrRegions; i++) {
    if(bin >= Regions[i][0] && bin <= Regions[i][1])
      return i+1;
  }
  return 0;
}
void PlotProfile(int NrBins, float *Ipulse, char *xlabel, char *ylabel, char *title, int Highlight, int color, int clearPage)
{
  long j;
  float ymin, ymax;
  ymin = ymax = Ipulse[0];
  for(j = 1; j < NrBins; j++) {
    if(Ipulse[j] > ymax)
      ymax = Ipulse[j];
    if(Ipulse[j] < ymin)
      ymin = Ipulse[j];
  }
  if(clearPage) {
    ppgpage();
    ppgsci(1);
    ppgsvp(0.1, 0.9, 0.1, 0.9);
    ppgswin(0,NrBins-1,-0.1,1.1);
    ppgbox("bcnsti",0.0,0,"bcnti",0.0,0);
    ppglab(xlabel, ylabel, title);
  }
  ppgsci(color);
  ppgmove(0, Ipulse[0]/ymax);
  for(j = 1; j < NrBins; j++) {
    if(Highlight != 0)
      ppgsci(color+CheckOnPulse(j,NrOnpulseRegions,OnPulseRegion));
    else
      ppgsci(color);
    ppgdraw(j, Ipulse[j]/ymax);
  }
  ppgsci(1);
}
void ShiftProfile(int shift, int NrBins, float *Iprofile, float *outputProfile)
{
  int b, b2;
  for(b = 0; b < NrBins; b++) {
    b2 = b+shift;
    if(b2 < 0)
      b2 += NrBins;
    if(b2 >= NrBins)
      b2 -= NrBins;
    outputProfile[b2] = Iprofile[b];
  }
}

max_stars_repo_name: weltevrede/psrsalsa

max_issues_repo_name: weltevrede/psrsalsa

max_forks_repo_name: weltevrede/psrsalsa


text:
// =============================================================================
// == niw_sampled.h
// == --------------------------------------------------------------------------
// == A class for a Normal Inverse-Wishart distribution
// == --------------------------------------------------------------------------
// == Copyright 2013. MIT. All Rights Reserved.
// == Written by Jason Chang 11-03-2013
// == --------------------------------------------------------------------------
// == If this code is used, the following should be cited:
// ==
// == [1] J. Chang and J. W. Fisher II, "Parallel Sampling of DP Mixtures
// ==     Models using Sub-Cluster Splits". Neural Information Processing
// ==     Systems (NIPS 2013), Lake Tahoe, NV, USA, Dec 2013.
// =============================================================================

#ifndef _NIW_SAMPLED_H_INCLUDED_
#define _NIW_SAMPLED_H_INCLUDED_

//#include "matrix.h"
//#include "mex.h"
#include <math.h>
//#include "array.h"

//#include "helperMEX.h"
//#include "debugMEX.h"
#include "dpmmSubclusters/normal.h"

#include "dpmmSubclusters/linear_algebra.h"
#include "dpmmSubclusters/myfuncs.h"

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

#include <vector>
using std::vector;

#ifndef pi
#define pi 3.14159265
#endif

#ifndef logpi
#define logpi 1.144729885849
#endif

class niw_sampled
{
   gsl_rng *r;

public:
   // prior hyperparameters
   double kappah;
   double nuh;
   arr(double) thetah;
   arr(double) Deltah;
   int D;
   int D2;

   // sufficient statistics of the observed data
   arr(double) t;
   arr(double) T;
   int N;

   // posterior hyperparameters
   double kappa;
   double nu;
   arr(double) theta;
   arr(double) Delta;

   arr(double) tempVec;
   arr(double) tempMtx;

   // instantiated gaussian parameters
   normal param;
   double logpmu_prior, logpmu_posterior;

public:
   // --------------------------------------------------------------------------
   // -- niw_sampled
   // --   constructor; initializes to empty
   // --------------------------------------------------------------------------
   niw_sampled();
   // --------------------------------------------------------------------------
   // -- niw_sampled
   // --   copy constructor;
   // --------------------------------------------------------------------------
   niw_sampled(const niw_sampled& that);
   // --------------------------------------------------------------------------
   // -- operator=
   // --   assignment operator
   // --------------------------------------------------------------------------
   niw_sampled& operator=(const niw_sampled& that);
   // --------------------------------------------------------------------------
   // -- copy
   // --   returns a copy of this
   // --------------------------------------------------------------------------
   void copy(const niw_sampled& that);
   // --------------------------------------------------------------------------
   // -- niw_sampled
   // --   constructor; intializes to all the values given
   // --------------------------------------------------------------------------
   niw_sampled(int _D, double _kappah, double _nuh, arr(double) _thetah, arr(double) _Deltah);

   // --------------------------------------------------------------------------
   // -- ~niw_sampled
   // --   destructor
   // --------------------------------------------------------------------------
   virtual ~niw_sampled();

   // --------------------------------------------------------------------------
   // -- ~cleanup
   // --   deletes all the memory allocated by this
   // --------------------------------------------------------------------------
   virtual void cleanup();

public:
   // --------------------------------------------------------------------------
   // -- empty
   // --   Empties out the statistics of the niw_sampled (i.e. no data).
   // --------------------------------------------------------------------------
   void clear();
   void empty();

   bool isempty() const;
   int getN() const;
   int getD() const;
   arr(double) get_mean() const;
   arr(double) get_cov() const;
   arr(double) get_prec() const;
   gsl_rng* get_r();

   void set_normal(normal &other);
   void set_normal(arr(double) _mean, arr(double) _cov);
   normal* get_normal();

   // --------------------------------------------------------------------------
   // -- update_posteriors
   // --   Updates the posterior hyperparameters
   // --------------------------------------------------------------------------
   void update_posteriors();
   void update_posteriors_sample();

   // --------------------------------------------------------------------------
   // -- add_data
   // --   functions to add an observation to the niw_sampled. Updates the sufficient
   // -- statistics, posterior hyperparameters, and predictive parameters
   // --
   // --   parameters:
   // --     - data : the new observed data point of size [1 D]
   // --------------------------------------------------------------------------
   void add_data_init(arr(double) data);
   void add_data(arr(double) data);
   void merge_with(niw_sampled* &other, bool doSample);
   void merge_with(niw_sampled* &other1, niw_sampled* &other2, bool doSample);
   void merge_with(niw_sampled &other, bool doSample);
   void merge_with(niw_sampled &other1, niw_sampled &other2, bool doSample);
   void set_stats(int _N, arr(double) _t, arr(double) _T);

   double Jdivergence(const niw_sampled &other);

   double predictive_loglikelihood(arr(double) data) const;
   double predictive_loglikelihood_mode(arr(double) data) const;
   double data_loglikelihood() const;
   double data_loglikelihood_marginalized() const;

   double data_loglikelihood_marginalized_testmerge(niw_sampled *other) const;
   void sample(normal &_param);
   void sample_scale(normal &_param);
   void sample();
   void find_mode();
   double logmu_posterior(const normal &_param) const;
   double logmu_posterior() const;
   double logmu_prior(const normal &_param) const;
   double logmu_prior() const;
};

inline double niw_sampled::predictive_loglikelihood(arr(double) data) const
{
   return param.predictive_loglikelihood(data);
}

#endif

max_stars_repo_name: jstraub/dpMM

max_issues_repo_name: jstraub/dpMM

max_forks_repo_name: jstraub/dpMM


text:
#include <limits.h>

#include <asf.h>
#include <asf_export.h>
#include <asf_raster.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>

/* Get sample size in bytes of the data types represented by the
   meta_parameters_t.  */
size_t
get_sample_size (meta_parameters *metadata)
{
  size_t sample_size;
  switch ( metadata->general->data_type ) {
  case BYTE:
    sample_size = sizeof (int8_t);
    break;
  case INTEGER16:
    sample_size = sizeof (int16_t);
    break;
  case INTEGER32:
    sample_size = sizeof (int32_t);
    break;
  case REAL32:
    sample_size = sizeof (float);
    break;
  case REAL64:
    sample_size = sizeof (double);
    break;
  default:
    /* Other types aren't handled.  */
    asfPrintError ("This program can handle byte, int16, int32, real32, and real64 data types.\n");
    sample_size = 0; // not reached
  }

  return sample_size;
}

/* Get the image data in data file image_data_file_name, using
   metadata metadata.  A pointer to new memory containing the image
   data is returned.  */
void *
get_image_data (meta_parameters *metadata, const char *image_data_file)
{
  size_t sample_size = get_sample_size (metadata);
  size_t pixel_count;
  void *data;
  size_t read_count;
  int return_code;

  /* Read the image data itself.  */
  FILE *ifp = fopen (image_data_file, "r");
  if ( ifp == NULL )
    asfPrintError("Failed to open %s: %s", image_data_file, strerror(errno));

  /* Total number of samples in image.  */
  pixel_count = metadata->general->line_count * metadata->general->sample_count;
  data = MALLOC (pixel_count * sample_size);
  read_count = fread (data, sample_size, pixel_count, ifp);
  if ( read_count != pixel_count ) {
    if ( feof (ifp) ) {
      asfPrintError("Read wrong amount of data from %s", image_data_file);
    }
    else if ( ferror (ifp) ) {
      asfPrintError("Read of file %s failed: %s", image_data_file,
                    strerror(errno));
    }
    else {
      /* Shouldn't get here.  */
      asfPrintError("Unknown error reading %s\n", image_data_file);
    }
    exit (EXIT_FAILURE);
  }

  return_code = fclose (ifp);
  asfRequire ((return_code==0), "Error closing file %s\n", image_data_file);

  return data;
}

/* Return the average of all the kernel_size * kernel_size elements
   centered around element i, j.  The i and j arguments must be within
   the bounds of img, and the kernel_size must be odd.  The image is
   reflected at the edges and corners for the purpose of determining
   this average.  */
unsigned char
averaging_kernel (gsl_matrix_uchar *img, int kernel_size, size_t i, size_t j)
{
  int i_idx, j_idx;
  int i_min = i - kernel_size / 2;
  int i_max = i + kernel_size / 2;
  int j_min = j - kernel_size / 2;
  int j_max = j + kernel_size / 2;
  int sum = 0;
  int average; /* Truncated average.  */

  asfRequire(kernel_size%2 != 0, "Odd-sized kernels only.\n");

  for ( i_idx = i_min ; i_idx < i_max ; i_idx++ ) {
    /* The i index to use, adjusted in case we are off the edge of the
       image.  This choice (and the corresponding choice for j)
       implement a kernel that pretends that a mirror image of the
       image exists at the image edges (and corners). */
    int itu = i_idx;
    if ( itu < 0 )
      itu = -itu - 1;
    else if ( itu >= img->size1 )
      itu = img->size1 - (itu - img->size1) - 1;
    for ( j_idx = j_min ; j_idx < j_max ; j_idx++ ) {
      /* See the comment for variable itu above.  */
      int jtu = j_idx;
      if ( jtu < 0 )
        jtu = -jtu - 1;
      else if ( jtu >= img->size2 )
        jtu = img->size2 - (jtu - img->size2) - 1;
      sum += gsl_matrix_uchar_get (img, itu, jtu);
    }
  }

  average = sum / pow (kernel_size, 2); /* Truncated average.  */
  /* Since we are averaging unsigned char values, this should always
     be true.  */
  asfRequire(average<=UCHAR_MAX,
             "The average value of the image is above its maximum value.\n");

  return average;
}

/* Average together tiles of kernel_size pixels in the *width by
   *height image at pixels, reducing its size by a factor of about
   kernel size in each dimension.  The new image replaces the old in
   pixels (extra memory is freed) and a pointer to the possibly
   relocated pixel data is returned.  The new image width and height
   replace the input width and height arguments.  */
unsigned char *average_unsigned_char_pixels (unsigned char *pixels,
                                             unsigned long *width,
                                             unsigned long *height,
					     int kernel_size)
{
  /* Input width and height.  */
  size_t iwidth = *width, iheight = *height;
  gsl_matrix_uchar *iimg;
  size_t ii, jj;
  size_t owidth;
  size_t oheight;
  gsl_matrix_uchar *oimg;
  unsigned char *reallocated_pixels;

  asfRequire (kernel_size%2 != 0, "Odd-sized kernels only.\n");

  /* Make a matrix form of the input image for easy indexing.  */
  iimg = gsl_matrix_uchar_alloc (iheight, iwidth);
  for ( ii = 0 ; ii < iheight ; ii++ ) {
    for ( jj = 0 ; jj < iwidth ; jj++ ) {
      gsl_matrix_uchar_set (iimg, ii, jj, pixels[ii * iwidth + jj]);
    }
  }

  /* Dimensions of the averaged image.  */
  owidth = ceil (iwidth / kernel_size);
  oheight = ceil (iheight / kernel_size);

  /* Form the output image.  */
  oimg = gsl_matrix_uchar_alloc (oheight, owidth);
  for ( ii = 0 ; ii < oheight ; ii++ ) {
    for ( jj = 0 ; jj < owidth ; jj++ ) {
      gsl_matrix_uchar_set
        (oimg, ii, jj, averaging_kernel (iimg, kernel_size,
                                         ii * kernel_size + kernel_size / 2,
                                         jj * kernel_size + kernel_size / 2));
    }
  }

  /* Write output image back into pixel memory.  */
  for ( ii = 0 ; ii < oheight ; ii++ ) {
    for ( jj = 0 ; jj < owidth ; jj++ ) {
      pixels[ii * owidth + jj] = gsl_matrix_uchar_get (oimg, ii, jj);
    }
  }

  /* Done with matrix forms of the image.  */
  gsl_matrix_uchar_free (oimg);
  gsl_matrix_uchar_free (iimg);

  /* Modify the input/output arguments and resize the allocated space
     as promised.  */
  *width = owidth;
  *height = oheight;
  reallocated_pixels = realloc (pixels, owidth * oheight);

  return reallocated_pixels;
}

/* Given a block of pixel_count floats, map them linearly from range
   [max (minsample, mean - 3 * sigma), min (maxsample, mean + 3 * sigma)]
   into the unsigned char range, with input samples outside the above
   range clamped.  Return the data in new malloc()ed memory.  */
unsigned char *map_floats_to_unsigned_bytes (float *daf, size_t pixel_count)
{
  unsigned char *pixels = malloc (pixel_count * sizeof (unsigned char));

  /* Minimum and maximum values in the input data.  */
  float imin = gsl_stats_float_min (daf, 1, pixel_count);
  float imax = gsl_stats_float_max (daf, 1, pixel_count);
  /* Mean value of input data.  */
  float imean = gsl_stats_float_mean (daf, 1, pixel_count);
  /* Standard deviation of input data.  */
  float isdev = gsl_stats_float_sd (daf, 1, pixel_count);
  /* Minimum and maximum after clamping.  */
  float omin = GSL_MAX (imean - 3 * isdev, imin);
  float omax = GSL_MIN (imean + 3 * isdev, imax);

  /* Shift we need to apply to the data to get it to fall in the
     range of the unsigned chars.  */
  float bias = -omin + 0.25;

  /* Compute all the output pixels.  */
  size_t ii;
  for ( ii = 0 ; ii < pixel_count ; ii++ ) {
    if ( daf[ii] < omin ) {
      pixels[ii] = 0;           /* Clamp low.  */
    }
    else if ( daf[ii] > omax ) {
      pixels[ii] = UCHAR_MAX;   /* Clamp high.  */
    }
    else {
      pixels[ii] = (daf[ii] + bias) * (UCHAR_MAX / (omax - omin));
    }
  }

  return pixels;
}

/* Scale the *width x *height image at pixels st its large dimension
   is less than or equal to max_large_dimension.  The memory pointed
   to by pixels is resized and possible relocated, with the new
   location being returned.  The new image width and height are
   returned in *width and *height.  */
unsigned char *scale_unsigned_char_image_dimensions (unsigned char *pixels,
                                      unsigned long max_large_dimension,
                                      unsigned long *width,
                                      unsigned long *height)
{
  /* This assertion is pretty obvious, but since the algorithm needs
     it to work correctly, its included.  */
  asfRequire(max_large_dimension>1,
             "Image needs to be at least one pixel by one pixel.\n");
  if ( GSL_MAX (*width, *height) > max_large_dimension ) {
    int kernel_size = GSL_MAX (*width, *height) / max_large_dimension + 1;
    if ( kernel_size % 2 != 1 ) {
      kernel_size++;
    }
    pixels = average_unsigned_char_pixels (pixels, width, height, kernel_size);
  }

  return pixels;
}

/* Figure basic statistics for the image as needed for the different data
   scaling techniques available to the user */
void
get_statistics (FloatImage *si, scale_t sample_mapping, int sampling_stride,
                float *mean, float *standard_deviation,
                float *min_sample, float *max_sample,
                float *omin, float *omax, gsl_histogram **hist,
                float no_data)
{
  asfRequire(*max_sample >= *min_sample,
              "get_statistics() error ...max_sample parameter < min_sample parameter\n");

  if ( sample_mapping == SIGMA ) {
    float_image_approximate_statistics (si, sampling_stride, mean,
                                        standard_deviation, no_data);
    *omin = *mean - 2 * (*standard_deviation);
    *omax = *mean + 2 * (*standard_deviation);
  }

  else if ( sample_mapping == MINMAX ) {
    float_image_statistics (si, min_sample, max_sample, mean,
                            standard_deviation, no_data);
    *omin = *min_sample;
    *omax = *max_sample;
  }

  else if ( sample_mapping == HISTOGRAM_EQUALIZE ) {
    float_image_statistics (si, min_sample, max_sample, mean,
                            standard_deviation, no_data);
    // Add a little bit of tail padding on the histgram to avoid problems
    // when searching for image min or max in the histogram (gsl histogram
    // limitation)
    //*omin = *min_sample * 0.025 * (*min_sample);
    //*omax = *max_sample * 0.025 * (*max_sample);
    //
    // NOTE:  The above doesn't work properly for data with negative
    // values and also, I think the first multiply is a typo and that
    // a '+' was probably intended.  If so, then the affect of the
    // above code is for the GSL increment (bin counting) function to
    // ignore values outside the adjusted min/max range, e.g. ignore
    // boundary case values (black and white noise tends to sum in the
    // limits, so ignoring the boundaries produces a higher quality
    // image.)

    // NOTE: The following assumes that the intent of the above code was
    // to ignore boundary case values, but the new histogram limits are
    // now calculated based on the full data range.
    float bin_range;

    bin_range = (*max_sample - *min_sample) / NUM_HIST_BINS;
    *omin = *min_sample + (0.025 * bin_range);
    *omax = *max_sample - (0.025 * bin_range);
    *hist = float_image_gsl_histogram (si, *omin, *omax, NUM_HIST_BINS);
  }
}

#ifndef linux
#ifndef darwin
#ifndef win32
static double
round (double arg)
{
    return floor (arg + 0.5);
}
#endif // #ifndef win32
#endif // #ifndef darwin
#endif // #ifndef linux

/* Create a byte value based on the floating point value and whatever scaling
   method the user chose */
unsigned char
pixel_float2byte(float paf, scale_t sample_mapping,
                 float omin, float omax, gsl_histogram *hist,
                 gsl_histogram_pdf *hist_pdf, float no_data_value)
{
  const unsigned char min_brightness=0;          // Minimum brightess via byte
  const unsigned char max_brightness=UCHAR_MAX;  // Maximum brightess via byte
  unsigned char pab = 0;                         // Pixel As Byte.
  size_t hist_bin;                       // histogram bin for given float value
  double slope, offset;

  // This case is easy -- always map the "no data" value to 0.
  if (paf == no_data_value)
      return 0;

  switch ( sample_mapping ) {
  case TRUNCATE:
    if ( paf <= (float) min_brightness ) {
      pab = min_brightness;
    }
    else if ( paf >= (float) max_brightness ) {
      pab = max_brightness;
    }
    else {
      pab = (unsigned char) paf;
    }
    break;
  case MINMAX:
  case MINMAX_MEDIAN:
  case SIGMA:
    slope = max_brightness / (omax-omin);
    offset = -slope * omin;
    if ((slope * paf + offset) < (float)min_brightness)
      pab = min_brightness;
    else if ((slope * paf + offset) > (float)max_brightness)
      pab = max_brightness;
    else
      pab = slope * paf + offset;
    break;
  case HISTOGRAM_EQUALIZE:
    // Since we used a mask when we called float_image_statistics() earlier,
    // we've got to account for it when writing the image out. The only way we
    // could think to do that was casting the mask value to byte :(.
    //
    // FIXME: The following cast of FLOAT_IMAGE_DEFAULT_MASK (currently equal
    // to 0.0) is incorrect for images such as DEMs and db's that contain
    // negative data.  I'm voting for something similar, but force data values
    // sufficiently close to omin TO omin ...if we do any forcing at all.
    // SIDE EFFECTS: As-is, the following will result in black shot-noise and
    // 'topo line' type artifacts in images that contain negative data (in regions
    // where zero-crossings occur or in noise near 0.0)
    //
    if (0==gsl_fcmp(paf, FLOAT_IMAGE_DEFAULT_MASK, 0.00000000001)) {
      pab = (unsigned char)FLOAT_IMAGE_DEFAULT_MASK;
    }
    else {
      // Determine which bin the pixel-as-float is in
      if (paf <= gsl_histogram_min(hist)) {
	hist_bin = 0;
      }
      else if (paf >= gsl_histogram_max(hist)) {
        hist_bin = gsl_histogram_bins(hist) - 1;
      }
      else {
        // Pixel is within a valid range, so the following should not complain
	gsl_histogram_find (hist, paf, &hist_bin);
      }

      double pdf_at_index = gsl_histogram_get ((gsl_histogram*)hist_pdf,
					       hist_bin);
      pab = (unsigned char)(max_brightness * pdf_at_index);
    }
    break;
  default:
    asfPrintError("Impossible: sample_mapping=%d\n", sample_mapping);
    break;
  }
  return pab;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady