text:
/*
    KGSX: Biomolecular Kino-geometric Sampling and Fitting of Experimental Data
    Yao et al, Proteins. 2012 Jan;80(1):25-43
    e-mail: latombe@cs.stanford.edu, vdbedem@slac.stanford.edu, julie.bernauer@inria.fr

        Copyright (C) 2011-2013 Stanford University

        Permission is hereby granted, free of charge, to any person obtaining a copy of
        this software and associated documentation files (the "Software"), to deal in
        the Software without restriction, including without limitation the rights to
        use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
        of the Software, and to permit persons to whom the Software is furnished to do
        so, subject to the following conditions:

        This entire text, including the above copyright notice and this permission notice
        shall be included in all copies or substantial portions of the Software.

        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
        OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
        FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
        IN THE SOFTWARE.
*/

#ifndef KGS_RELATIVEMSDDIRECTION_H
#define KGS_RELATIVEMSDDIRECTION_H

#include <gsl/gsl_vector.h>

#include <tuple>
#include <vector>

#include "core/graph/KinTree.h"
#include "Direction.h"
#include "Selection.h"

class RelativeMSDDirection: public Direction {
 public:
  RelativeMSDDirection( std::vector< std::tuple<Atom*, Atom*, double> > &relativeDistances);

 protected:
  void computeGradient(Configuration* conf, Configuration* target, gsl_vector* ret);

 private:
  std::vector< std::tuple<Atom*, Atom*, double> > &m_relativeDistances;
};

#endif //KGS_RELATIVEMSDDIRECTION_H


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include <iosfwd>
#include <gsl\gsl>
#include <vector>

#include "Data\Data.h"
#include "Data\TimeUnits.h"

#include <Plat/CameraDevice/CameraSettings.h>

namespace mage
{
    namespace binary_data
    {
        constexpr uint64_t LATEST_BINARY_VERSION = 3;
        constexpr uint64_t LATEST_CALIBRATION_VERSION = 2;

        // portable strings for device information
        struct DeviceInformation
        {
            std::string Id;
            std::string DeviceName;
            std::string DeviceManufacturer;
            std::string DeviceHardwareVersion;
            std::string DeviceFirmwareVersion;
        };

        struct HeaderData
        {
            uint64_t BinaryVersion;
            uint64_t CalibrationVersion;
            calibration::LinearFocalLengthModel Calibration;
            mage::Size ImageSize;
            mage::PixelFormat FormatOfPixels;

            DeviceInformation DeviceInformation;

            HeaderData(const uint64_t& b, const uint64_t& m, const calibration::LinearFocalLengthModel& c, const mage::Size& size, const mage::PixelFormat &fp, const mage::binary_data::DeviceInformation& deviceInfo)
                : BinaryVersion{ b }, CalibrationVersion{ m }, Calibration{ c }, ImageSize{ size }, FormatOfPixels{ fp }, DeviceInformation{ deviceInfo }
            {}

            HeaderData()
                : BinaryVersion{ LATEST_BINARY_VERSION },
                CalibrationVersion{ LATEST_CALIBRATION_VERSION },
                Calibration{},
                ImageSize{},
                FormatOfPixels{ mage::PixelFormat::GRAYSCALE8 },
                DeviceInformation{ "n/a", "n/a", "n/a", "n/a" }
            {}
        };

        struct FileFrameData
        {
            uint64_t CorrelationId;
            hundred_nanoseconds TimeStamp;
            mira::CameraSettings CameraSettings;
            std::vector<uint8_t> Pixels;

            FileFrameData(const uint64_t& correlationId, const hundred_nanoseconds& t, const mira::CameraSettings& cameraSettings, std::vector<uint8_t> p)
                : CorrelationId{ correlationId }, TimeStamp{ t }, CameraSettings{ cameraSettings }, Pixels{ std::move(p) }
            {}
            FileFrameData() = default;
            FileFrameData(FileFrameData&&) = default;
        };

        void SerializeCaptureHeaderData(std::fstream& outputFile, const HeaderData& data);
        void SerializeFrameData(std::fstream& outputFile, const FileFrameData& data);

        bool DeSerializeHeaderData(std::istream& inputFile, HeaderData& data);
        FileFrameData DeSerializeFrameData(std::istream& inputFile, HeaderData& headerData);
    }
}

ext: h

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
// https://stackoverflow.com/questions/66536134/gsl-gsl-linalg-sv-decomp-returning-u-and-v-reversed
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

void pretty_print(const gsl_matrix * M)
{
  // Get the dimension of the matrix.
  int rows = M->size1;
  int cols = M->size2;
  // Now print out the data in a square format.
  for(int i = 0; i < rows; i++){
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_matrix_get(M, i, j));
    }
    printf("\n");
  }
  printf("\n");
}

void pretty_print_vector(const gsl_vector * M)
{
  int cols = M->size;
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_vector_get(M, j));
    }
  printf("\n");
}

int
main()
{
  const size_t M = 4;
  const size_t N = 5;
  double A_data[] = {
    1.0, 0.0, 0.0, 0.0, 2.0,
    0.0, 0.0, 3.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 2.0, 0.0, 0.0, 0.0 };
  gsl_matrix_view A = gsl_matrix_view_array(A_data, 4, 5);
  gsl_matrix * B = gsl_matrix_alloc(N, M);
  gsl_matrix * V = gsl_matrix_alloc(M, M);
  gsl_vector * S = gsl_vector_alloc(M);
  gsl_vector * work = gsl_vector_alloc(M);

  gsl_matrix_transpose_memcpy(B, &A.matrix);

  gsl_linalg_SV_decomp(B, V, S, work);

  printf("U:\n");
  pretty_print(V);

  printf("S:\n");
  pretty_print_vector(S);

  printf("V:\n");
  pretty_print(B);

  gsl_matrix_free(B);
  gsl_matrix_free(V);
  gsl_vector_free(S);
  gsl_vector_free(work);

  return 0;
}

ext: c


text:
#pragma once

// Windows
#if defined(_WIN32)
#if !defined(NOMINMAX)
#define NOMINMAX
#endif
#include <windows.h>
#endif

// Standard
#include <cstdlib>
#include <exception>
#include <stdexcept>
#include <cassert>
#include <string>
#include <vector>
#include <map>
#include <memory>
#include <cstdint>
#include <sstream>
#include <fstream>
#include <chrono>
#include <functional>
#include <algorithm>
#include <iterator>

#if defined(DEBUG) || defined(_DEBUG)
#define _CRTDBG_MAP_ALLOC
#include <crtdbg.h>
#endif

// Guidline Support Library
#include <gsl/gsl>

// OpenGL
#include <GL/gl3w.h>
#define GLFW_INCLUDE_GLCOREARB
#include <GLFW/glfw3.h>

#if defined(_LINUX)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wconversion"
#endif

#define GLM_ENABLE_EXPERIMENTAL
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>

#if defined(_LINUX)
#pragma GCC diagnostic pop
#endif

//#include <SOIL/SOIL.h>

// Assimp
#if defined(_LINUX)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wconversion"
#endif
#include <assimp/scene.h>
#include <assimp/Importer.hpp>
#include <assimp/postprocess.h>
#if defined(_LINUX)
#pragma GCC diagnostic pop
#endif

// Local
#include "ServiceContainer.h"

namespace Library
{
	extern ServiceContainer GlobalServices;
}

max_stars_repo_name: DakkyDaWolf/OpenGL

max_issues_repo_name: DakkyDaWolf/OpenGL

max_forks_repo_name: DakkyDaWolf/OpenGL


text:
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_gegenbauer.h>
#include "scf.h"

/* Given the coefficients stored mat_cos and mat_sin, this code does
 * the summation to compute the potential and accelerations
 *
 * scf_potential() is the full 3D treatment for the host halo
 *
 * scf_potential_spherical() is the spherically symmetric treatment
 * for the subhalos. This spherical code skips the loops over l,m
 * so it's simpler.
 */

float scf_potential(float* a, double* mat_cos, double* mat_sin, int orbit_nmax, int orbit_lmax, float x, float y, float z)
{
    double x2 = (double)x * (double)x;
    double y2 = (double)y * (double)y;
    double r = sqrt(x2+y2+(double)z*(double)z);
    double r_xy = sqrt(x2+y2);
    double phi = atan2((double)y, (double)x);

    double sin_theta = r_xy / r;
    double sin_phi = (double)y / r_xy;
    double cos_theta = (double)z / r;
    double cos_phi = (double)x / r_xy;

    double cos_theta_for_Plm = cos_theta;

#ifdef USE_GSL_GEGENBAUER_HOST
    double xi = (r-1.0)/(r+1.0);
#endif

    double a_r = 0.0;
    double a_t = 0.0;
    double a_p = 0.0;

    int n,l,m;

    float C_lm = 0.0f;
    float D_lm = 0.0f;
    float E_lm = 0.0f;
    float F_lm = 0.0f;

    float this_Phi_nl = 0.f;
    float this_dPhi_nl_dr = 0.f;

#ifdef USE_GSL_GEGENBAUER_HOST
    double* Cna_2l_32_array = (double*)malloc((orbit_nmax+1)*sizeof(double));
    double* Cna_2l_52_array = (double*)malloc(orbit_nmax*sizeof(double));
#endif

    double this_P_lm = 0.0;
    double this_dP_lm = 0.0;
    double* P_l;
    double* dP_l;
    int P_size;

    float cos_mphi;
    float sin_mphi;

    float pot = 0.f;
    a[0] = 0.f;
    a[1] = 0.f;
    a[2] = 0.f;

    for(m=0; m<=orbit_lmax; m++)
    {
        cos_mphi = cos(m*phi);
        sin_mphi = sin(m*phi);

        P_size = orbit_lmax - m + 1;
        P_l = (double*)malloc(P_size*sizeof(double));
        dP_l = (double*)malloc(P_size*sizeof(double));

        gsl_sf_legendre_Plm_deriv_array(orbit_lmax, m, cos_theta_for_Plm, P_l, dP_l);

        for(l=m; l<=orbit_lmax; l++)
        {
            C_lm = 0.f;
            D_lm = 0.f;
            E_lm = 0.f;
            F_lm = 0.f;

#ifdef USE_GSL_GEGENBAUER_HOST
            gsl_sf_gegenpoly_array(orbit_nmax, 2*l+1.5, xi, Cna_2l_32_array);
            if(orbit_nmax > 0)
                gsl_sf_gegenpoly_array(orbit_nmax-1, 2*l+2.5, xi, Cna_2l_52_array);
#endif

            for(n=0; n<=orbit_nmax; n++)
            {
#ifdef USE_GSL_GEGENBAUER_HOST
                this_Phi_nl = Phi_nl(Cna_2l_32_array[n], l, r);
                if(n>0)
                    this_dPhi_nl_dr = dPhi_nl_dr(Cna_2l_52_array[n-1], this_Phi_nl, n, l, r);
                else
                    this_dPhi_nl_dr = dPhi_nl_dr(1e99, this_Phi_nl, n, l, r);
#else
                this_Phi_nl = Phi_nl(n,l,r);
                this_dPhi_nl_dr = dPhi_nl_dr(this_Phi_nl, n, l, r);
#endif

                C_lm += mat_cos[ind(n,l,m)] * this_Phi_nl;
                D_lm += mat_sin[ind(n,l,m)] * this_Phi_nl;

                E_lm += mat_cos[ind(n,l,m)] * this_dPhi_nl_dr;
                F_lm += mat_sin[ind(n,l,m)] * this_dPhi_nl_dr;
            }

            this_P_lm = P_l[l-m];
            this_dP_lm = dP_l[l-m] * -sin_theta;

            pot += this_P_lm * (C_lm*cos_mphi + D_lm*sin_mphi);

            a_r += this_P_lm * (E_lm*cos_mphi + F_lm*sin_mphi);
            a_t += this_dP_lm * (C_lm*cos_mphi + D_lm*sin_mphi);
            a_p += m * this_P_lm * (D_lm*cos_mphi - C_lm*sin_mphi);

        } // l loop

        free(P_l);
        free(dP_l);

    } // m loop

    a_r = -a_r;
    a_t = -a_t/r;
    a_p = -a_p/r_xy;

    a[0] = a_r*sin_theta*cos_phi + a_t*cos_theta*cos_phi - a_p*sin_phi;
    a[1] = a_r*sin_theta*sin_phi + a_t*cos_theta*sin_phi + a_p*cos_phi;
    a[2] = a_r*cos_theta - a_t*sin_theta;

    a[0] *= G;
    a[1] *= G;
    a[2] *= G;
    pot *= G;

#ifdef USE_GSL_GEGENBAUER_HOST
    free(Cna_2l_32_array);
    free(Cna_2l_52_array);
#endif

    return pot;

}

float scf_potential_spherical(float* a, float* mat_cos, int orbit_nmax, float x, float y, float z)
{
    float x2 = x*x;
    float y2 = y*y;

    float r = sqrt(x2+y2+z*z);
    float r_xy = sqrt(x2+y2);

#ifdef USE_GSL_GEGENBAUER_SUBHALO
    float xi = (r-1.0f)/(r+1.0f);
#endif

    float sin_theta = r_xy / r;
    float sin_phi = y / r_xy;
    float cos_theta = z / r;
    float cos_phi = x / r_xy;

    float this_Phi_n0 = 0.f;
    float this_dPhi_n0_dr = 0.f;

#ifdef USE_GSL_GEGENBAUER_SUBHALO
    double* Cna_32_array = (double*)malloc((orbit_nmax+1)*sizeof(double));
    double* Cna_52_array = (double*)malloc(orbit_nmax*sizeof(double));
    gsl_sf_gegenpoly_array(orbit_nmax, 1.5, xi, Cna_32_array);
    gsl_sf_gegenpoly_array(orbit_nmax-1, 2.5, xi, Cna_52_array);
#endif

    float C_lm = 0.0f;
    float E_lm = 0.0f;

    int n;
    for(n=0; n<=orbit_nmax; n++)
    {
#ifdef USE_GSL_GEGENBAUER_SUBHALO
        this_Phi_n0 = Phi_n0(Cna_32_array[n], r);
        if(n>0)
            this_dPhi_n0_dr = dPhi_n0_dr(Cna_52_array[n-1], this_Phi_n0, n, r);
        else
            this_dPhi_n0_dr = dPhi_n0_dr(1e99, this_Phi_n0, n, r);
#else
        this_Phi_n0 = Phi_n0(n,r);
        this_dPhi_n0_dr = dPhi_n0_dr(this_Phi_n0, n, r);
#endif

        C_lm += mat_cos[n] * this_Phi_n0;

        E_lm += mat_cos[n] * this_dPhi_n0_dr;
    }

    float pot = C_lm;

    float a_r = -E_lm;

    a[0] = G * a_r*sin_theta*cos_phi;
    a[1] = G * a_r*sin_theta*sin_phi;
    a[2] = G * a_r*cos_theta;

    pot = G * pot;

#ifdef USE_GSL_GEGENBAUER_SUBHALO
    free(Cna_32_array);
    free(Cna_52_array);
#endif

    return pot;

}


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

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

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

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

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

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

PRIVATE_EXCEPT_IN_TEST:
	/// \brief EolType contains the enum constants that indicate the type of
	/// end-of-line for a text file.
	enum class EolType
	{
		INDETERMINATE, ///< File contains no ends-of-line
		MIXED, ///< Inconsistent end-of-line convention, may be binary
		DOS, ///< End-of-line is "\r\n"
		MACINTOSH, ///< End-of-line is "\r" (now archaic)
		UNIX ///< End-of-line is "\n"
	};

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

	void queryFile(const Path& p) const;
	void translateFile(const Path& p) const;
	static char getIndicatorLetter(EolType eolType);
	static ::std::string toEolString(EolType eolType);
	static ::std::string displayPath(const Path& p);

	/// \brief Scans the input file, counts the end-of-line types, and
	/// optionally translates the ends-of-line into the provided stream.
	///
	/// If pOut is null, then the method simply counts the end-of-line types and
	/// returns the corresponding EolType constant.  In this case, the
	/// targetEolType parameter is ignored.
	///
	/// If pOut is not null, then the method additionally translates the file
	/// into *pOut.  In this case, targetEolType must be one of DOS, MACINTOSH,
	/// or UNIX.
	static EolType scanFile(::std::istream& in,
		/* out */ size_t& numDosEols, /* out */ size_t& numMacEols,
		/* out */ size_t& numUnixEols, /* out */ size_t& totalEols,
		::std::ostream* pOut = nullptr, EolType targetEolType = EolType::INDETERMINATE);

	static void replaceOriginalFileWithTemp(const Path& originalPath,
		const Path& tempPath);

	bool				m_isInQueryMode;
	EolType			m_targetEolType;
	bool				m_forceTranslation;
	FileEnumerator	m_fileEnumerator;
};

#endif // XEOL_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


text:
#pragma once

#include <halley/utils/utils.h>
#include <vector>
#include <gsl/gsl>
#include "halley/data_structures/flat_map.h"
#include "halley/maths/vector2.h"
#include "halley/maths/rect.h"
#include "halley/file/path.h"
#include <map>
#include <unordered_map>
#include <cstdint>
#include <utility>
#include <set>
#include "halley/data_structures/maybe.h"
#include "halley/maths/colour.h"
#include "halley/maths/vector4.h"

namespace Halley
{
    class String;

    class SerializerOptions
    {
    public:
        constexpr static int maxVersion = 1;

        int version = 0;
        bool exhaustiveDictionary = false;
        std::function< std::optional< size_t >( const String& string ) > stringToIndex;
        std::function< const String&( size_t index ) > indexToString;

        SerializerOptions() = default;
        SerializerOptions( int version ) :
            version( version )
        { }
    };

    class SerializerState
    {
    };

    class ByteSerializationBase
    {
    public:
        ByteSerializationBase( SerializerOptions options ) :
            options( std::move( options ) )
        { }

        SerializerState* setState( SerializerState* state );

        template < typename T >
        T* getState() const
        {
            return static_cast< T* >( state );
        }

        int getVersion() const { return version; }
        void setVersion( int v ) { version = v; }

    protected:
        SerializerOptions options;

    private:
        SerializerState* state = nullptr;
        int version = 0;
    };

    class Serializer : public ByteSerializationBase
    {
    public:
        Serializer( SerializerOptions options );
        explicit Serializer( gsl::span< gsl::byte > dst, SerializerOptions options );

        template < typename T, typename std::enable_if< std::is_convertible< T, std::function< void( Serializer& ) > >::value, int >::type = 0 >
        static Bytes toBytes( const T& f, SerializerOptions options = {} )
        {
            auto dry = Serializer( options );
            f( dry );
            Bytes result( dry.getSize() );
            auto s = Serializer( gsl::as_writable_bytes( gsl::span< Halley::Byte >( result ) ), options );
            f( s );
            return result;
        }

        template < typename T, typename std::enable_if< !std::is_convertible< T, std::function< void( Serializer& ) > >::value, int >::type = 0 >
        static Bytes toBytes( const T& value, SerializerOptions options = {} )
        {
            return toBytes( [ &value ]( Serializer& s ) { s << value; }, options );
        }

        size_t getSize() const { return size; }

        Serializer& operator<<( bool val ) { return serializePod( val ); }
        Serializer& operator<<( int8_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( uint8_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( int16_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( uint16_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( int32_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( uint32_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( int64_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( uint64_t val ) { return serializeInteger( val ); }
        Serializer& operator<<( float val ) { return serializePod( val ); }
        Serializer& operator<<( double val ) { return serializePod( val ); }

        Serializer& operator<<( const std::string& str );
        Serializer& operator<<( const String& str );
        Serializer& operator<<( const Path& path );
        Serializer& operator<<( gsl::span< const gsl::byte > span );
        Serializer& operator<<( const Bytes& bytes );

        template < typename T >
        Serializer& operator<<( const std::vector< T >& val )
        {
            unsigned int sz = static_cast< unsigned int >( val.size() );
            *this << sz;
            for ( unsigned int i = 0; i < sz; i++ )
            {
                *this << val[ i ];
            }
            return *this;
        }

        template < typename T, typename U >
        Serializer& operator<<( const FlatMap< T, U >& val )
        {
            *this << static_cast< unsigned int >( val.size() );
            for ( auto& kv : val )
            {
                *this << kv.first << kv.second;
            }
            return *this;
        }

        template < typename T, typename U >
        Serializer& operator<<( const std::map< T, U >& val )
        {
            *this << static_cast< unsigned int >( val.size() );
            for ( auto& kv : val )
            {
                *this << kv.first << kv.second;
            }
            return *this;
        }

        template < typename T, typename U >
        Serializer& operator<<( const std::unordered_map< T, U >& val )
        {
            std::map< T, U > m;
            for ( auto& kv : val )
            {
                m[ kv.first ] = kv.second;
            }
            return ( *this << m );
        }

        template < typename T >
        Serializer& operator<<( const std::set< T >& val )
        {
            unsigned int sz = static_cast< unsigned int >( val.size() );
            *this << sz;
            for ( auto& v : val )
            {
                *this << v;
            }
            return *this;
        }

        template < typename T >
        Serializer& operator<<( const Vector2D< T >& val )
        {
            return *this << val.x << val.y;
        }

        template < typename T >
        Serializer& operator<<( const Vector3D< T >& val )
        {
            return *this << val.x << val.y << val.z;
        }

        template < typename T >
        Serializer& operator<<( const Vector4D< T >& val )
        {
            return *this << val.x << val.y << val.z << val.w;
        }

        template < typename T >
        Serializer& operator<<( const Colour4< T >& val )
        {
            return *this << val.r << val.g << val.b << val.a;
        }

        template < typename T >
        Serializer& operator<<( const Rect2D< T >& val )
        {
            return *this << val.getTopLeft() << val.getBottomRight();
        }

        template < typename T, typename U >
        Serializer& operator<<( const std::pair< T, U >& p )
        {
            return *this << p.first << p.second;
        }

        template < typename T >
        Serializer& operator<<( const std::optional< T >& p )
        {
            if ( p )
            {
                return *this << true << p.value();
            }
            else
            {
                return *this << false;
            }
        }

        template < typename T >
        Serializer& operator<<( const Range< T >& p )
        {
            return *this << p.start << p.end;
        }

        template < typename T, std::enable_if_t< std::is_enum< T >::value == true, int > = 0 >
        Serializer& operator<<( const T& val )
        {
            using B = typename std::underlying_type< T >::type;
            return *this << B( val );
        }

        template < typename T, std::enable_if_t< std::is_enum< T >::value == false, int > = 0 >
        Serializer& operator<<( const T& val )
        {
            val.serialize( *this );
            return *this;
        }

        size_t getPosition() const { return size; }

    private:
        size_t size = 0;
        gsl::span< gsl::byte > dst;
        bool dryRun;

        template < typename T >
        Serializer& serializePod( T val )
        {
            if ( !dryRun )
            {
                memcpy( dst.data() + size, &val, sizeof( T ) );
            }
            size += sizeof( T );
            return *this;
        }

        template < typename T >
        Serializer& serializeInteger( T val )
        {
            if ( options.version >= 1 )
            {
                // Variable-length
                if constexpr ( std::is_signed_v< T > )
                {
                    serializeVariableInteger( static_cast< uint64_t >( val >= 0 ? val : -( val + 1 ) ), val < 0 );
                }
                else
                {
                    serializeVariableInteger( val, {} );
                }
                return *this;
            }
            else
            {
                // Fixed length
                return serializePod( val );
            }
        }

        void serializeVariableInteger( uint64_t val, std::optional< bool > sign );
    };

    class Deserializer : public ByteSerializationBase
    {
    public:
        Deserializer( gsl::span< const gsl::byte > src, SerializerOptions options = {} );
        Deserializer( const Bytes& src, SerializerOptions options = {} );

        template < typename T >
        static T fromBytes( const Bytes& src, SerializerOptions options = {} )
        {
            T result;
            Deserializer s( src, std::move( options ) );
            s >> result;
            return result;
        }

        template < typename T >
        static T fromBytes( gsl::span< const gsl::byte > src, SerializerOptions options = {} )
        {
            T result;
            Deserializer s( src, std::move( options ) );
            s >> result;
            return result;
        }

        template < typename T >
        static void fromBytes( T& target, const Bytes& src, SerializerOptions options = {} )
        {
            Deserializer s( src, std::move( options ) );
            s >> target;
        }

        template < typename T >
        static void fromBytes( T& target, gsl::span< const gsl::byte > src, SerializerOptions options = {} )
        {
            Deserializer s( src, std::move( options ) );
            s >> target;
        }

        Deserializer& operator>>( bool& val ) { return deserializePod( val ); }
        Deserializer& operator>>( int8_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( uint8_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( int16_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( uint16_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( int32_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( uint32_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( int64_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( uint64_t& val ) { return deserializeInteger( val ); }
        Deserializer& operator>>( float& val ) { return deserializePod( val ); }
        Deserializer& operator>>( double& val ) { return deserializePod( val ); }

        Deserializer& operator>>( std::string& str );
        Deserializer& operator>>( String& str );
        Deserializer& operator>>( Path& p );
        Deserializer& operator>>( gsl::span< gsl::byte > span );
        Deserializer& operator>>( Bytes& bytes );

        template < typename T >
        Deserializer& operator>>( std::vector< T >& val )
        {
            unsigned int sz;
            *this >> sz;
            ensureSufficientBytesRemaining( sz ); // Expect at least one byte per vector entry

            val.clear();
            val.reserve( sz );
            for ( unsigned int i = 0; i < sz; i++ )
            {
                val.push_back( T() );
                *this >> val[ i ];
            }
            return *this;
        }

        template < typename T, typename U >
        Deserializer& operator>>( FlatMap< T, U >& val )
        {
            unsigned int sz;
            *this >> sz;
            ensureSufficientBytesRemaining( sz * 2 ); // Expect at least two bytes per map entry

            std::vector< std::pair< T, U > > tmpData( sz );
            for ( unsigned int i = 0; i < sz; i++ )
            {
                *this >> tmpData[ i ].first >> tmpData[ i ].second;
            }
            val = FlatMap< T, U >( boost::container::ordered_unique_range_t(), tmpData.begin(), tmpData.end() );
            return *this;
        }

        template < typename T, typename U >
        Deserializer& operator>>( std::map< T, U >& val )
        {
            unsigned int sz;
            *this >> sz;
            ensureSufficientBytesRemaining( size_t( sz ) * 2 ); // Expect at least two bytes per map entry

            for ( unsigned int i = 0; i < sz; i++ )
            {
                T key;
                U value;
                *this >> key >> value;
                val[ key ] = std::move( value );
            }
            return *this;
        }

        template < typename T, typename U >
        Deserializer& operator>>( std::unordered_map< T, U >& val )
        {
            unsigned int sz;
            *this >> sz;
            ensureSufficientBytesRemaining( sz * 2 ); // Expect at least two bytes per map entry

            for ( unsigned int i = 0; i < sz; i++ )
            {
                T key;
                U value;
                *this >> key >> value;
                val[ key ] = std::move( value );
            }
            return *this;
        }

        template < typename T >
        Deserializer& operator>>( std::set< T >& val )
        {
            unsigned int sz;
            *this >> sz;
            ensureSufficientBytesRemaining( sz ); // Expect at least one byte per set entry

            val.clear();
            for ( unsigned int i = 0; i < sz; i++ )
            {
                T v;
                *this >> v;
                val.insert( std::move( v ) );
            }
            return *this;
        }

        template < typename T >
        Deserializer& operator>>( Vector2D< T >& val )
        {
            *this >> val.x;
            *this >> val.y;
            return *this;
        }

        template < typename T >
        Deserializer& operator>>( Vector3D< T >& val )
        {
            *this >> val.x;
            *this >> val.y;
            *this >> val.z;
            return *this;
        }

        template < typename T >
        Deserializer& operator>>( Vector4D< T >& val )
        {
            *this >> val.x;
            *this >> val.y;
            *this >> val.z;
            *this >> val.w;
            return *this;
        }

        template < typename T >
        Deserializer& operator>>( Colour4< T >& val )
        {
            *this >> val.r;
            *this >> val.g;
            *this >> val.b;
            *this >> val.a;
            return *this;
        }

        template < typename T >
        Deserializer& operator>>( Rect2D< T >& val )
        {
            Vector2D< T > p1, p2;
            *this >> p1;
            *this >> p2;
            val = Rect2D< T >( p1, p2 );
            return *this;
        }

        template < typename T, typename U >
        Deserializer& operator>>( std::pair< T, U >& p )
        {
            return *this >> p.first >> p.second;
        }

        template < typename T >
        Deserializer& operator>>( Range< T >& p )
        {
            return *this >> p.start >> p.end;
        }

        template < typename T >
        Deserializer& operator>>( std::optional< T >& p )
        {
            bool present;
            *this >> present;
            if ( present )
            {
                T tmp;
                *this >> tmp;
                p = tmp;
            }
            else
            {
                p = std::optional< T >();
            }
            return *this;
        }

        template < typename T, std::enable_if_t< std::is_enum< T >::value == true, int > = 0 >
        Deserializer& operator>>( T& val )
        {
            typename std::underlying_type< T >::type tmp;
            *this >> tmp;
            val = T( tmp );
            return *this;
        }

        template < typename T, std::enable_if_t< std::is_enum< T >::value == false, int > = 0 >
        Deserializer& operator>>( T& val )
        {
            val.deserialize( *this );
            return *this;
        }

        template < typename T >
        void peek( T& val )
        {
            const auto oldPos = pos;
            *this >> val;
            pos = oldPos;
        }

        size_t getPosition() const { return pos; }

    private:
        size_t pos = 0;
        gsl::span< const gsl::byte > src;

        template < typename T >
        Deserializer& deserializePod( T& val )
        {
            ensureSufficientBytesRemaining( sizeof( T ) );
            memcpy( &val, src.data() + pos, sizeof( T ) );
            pos += sizeof( T );
            return *this;
        }

        template < typename T >
        Deserializer& deserializeInteger( T& val )
        {
            if ( options.version >= 1 )
            {
                // Variable-length
                bool sign;
                uint64_t temp;
                deserializeVariableInteger( temp, sign, std::is_signed_v< T > );
                if ( sign )
                {
                    int64_t signedTemp = -int64_t( temp ) - 1;
                    val = static_cast< T >( signedTemp );
                }
                else
                {
                    val = static_cast< T >( temp );
                }
                return *this;
            }
            else
            {
                // Fixed length
                return deserializePod( val );
            }
        }

        void deserializeVariableInteger( uint64_t& val, bool& sign, bool isSigned );

        void ensureSufficientBytesRemaining( size_t bytes );
        size_t getBytesRemaining() const;
    };
} // namespace Halley

max_stars_repo_name: chidddy/halley

max_issues_repo_name: chidddy/halley

max_forks_repo_name: chidddy/halley


text:
/*! \file SpecialFunctions.h
 *  \brief header file defining several special functions
    This is a work in progress as generally can always call gsl for most of the functions.
 */

#ifndef SPECIALFUNCTIONS_H
#define SPECIALFUNCTIONS_H

#include <iostream>
#include <complex>
#include <Precision.h>
#include <cmath>
#include <gsl/gsl_sf.h>

using namespace std;
namespace Math
{

    //factorial
    int Factorial(int k);
    //Gamma function not yet idependently implemented (using just gsl but if want arbitrary precision must
	//implement at some point.
	//complete gamma function
	//Double_t GammaFunc(Double_t a);
	//Calculates the normalized (regularized) upper incomplete gamma function (upper integral).
	//Double_t GammaFunc(Double_t a, Double_t x);

	//Again same issue with Gamma Functions, not yet implemented for arbitrary precision
    //calculates spherical harmonics
    //Double_t SphericalHarmonic(Double_t theta, Double_t phi, int l, int m);
    //just real component
    //Double_t SphericalHarmonicR(Double_t theta, Double_t phi, int l, int m);
    //just imaginary component
    //Double_t SphericalHarmonicI(Double_t theta, Double_t phi, int l, int m);

}

#endif

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
#ifndef TNet_Types_h
#define TNet_Types_h

#ifdef HAVE_ATLAS
extern "C"{
  #include <cblas.h>
  #include <clapack.h>
}
#endif

namespace TNet
{
  // TYPEDEFS ..................................................................
#if DOUBLEPRECISION
  typedef double  BaseFloat;
#else
  typedef float   BaseFloat;
#endif

#ifndef UINT_16
  typedef unsigned short  UINT_16   ;
  typedef unsigned        UINT_32   ;
  typedef short           INT_16    ;
  typedef int             INT_32    ;
  typedef float           FLOAT_32  ;
  typedef double          DOUBLE_64 ;
#endif

  // ...........................................................................
  // The following declaration assumes that SSE instructions are enabled
  // and that we are using GNU C/C++ compiler, which defines the __attribute__
  // notation.
  //
  // ENABLE_SSE is defined in <config.h>. Its value depends on options given
  // in the configure phase of builidng the library
#if defined(__GNUC__ )
  // vector of four single floats
  typedef float  v4sf __attribute__((vector_size(16)));
  // vector of two single doubles
  typedef double v2sd __attribute__((vector_size(16)));

  typedef BaseFloat BaseFloat16Aligned __attribute__((aligned(16))) ;

  typedef union
  {
    v4sf    v;
    float   f[4];
  } f4vector;

  typedef union
  {
    v2sd    v;
    double  f[2];
  } d2vector;
#endif // ENABLE_SSE && defined(__GNUC__ )

  typedef enum
  {
#ifdef HAVE_ATLAS
    TRANS    = CblasTrans,
    NO_TRANS = CblasNoTrans
#else
    TRANS    = 'T',
    NO_TRANS = 'N'
#endif
  } MatrixTrasposeType;

} // namespace TNet

#endif // #ifndef TNet_Types_h

max_stars_repo_name: troylee/nnet-asr

max_issues_repo_name: troylee/nnet-asr

max_forks_repo_name: troylee/nnet-asr


text:
/*
Copyright (C) 2019-2020 JingWeiZhangHuai <jingweizhanghuai@163.com>
Licensed under the Apache License, Version 2.0; 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.
*/

#include <cblas.h>
#include "morn_tensor.h"

void _GetRowData0(float *tdata,float *data,int height,int width,int y_locate,int x_locate)
{
    if(y_locate<0)       {memset(data,0,width*sizeof(float)); return;}
    if(y_locate>=height) {memset(data,0,width*sizeof(float)); return;}

    tdata = tdata+y_locate*width;
    if(x_locate<0)
    {
        memset(data,0,(0-x_locate)*sizeof(float));
        memcpy(data-x_locate,tdata,(width+x_locate)*sizeof(float));
        return;
    }
    if(x_locate>0)
    {
        memcpy(data,tdata+x_locate,(width-x_locate)*sizeof(float));
        memset(data+width-x_locate,0,x_locate*sizeof(float));
        return;
    }
    memcpy(data,tdata,width*sizeof(float));
}
void _GetRowData(float *tdata,float *data,int height,int width,int y_locate,int x_locate,int stride)
{
    if(y_locate<0)       {memset(data,0,width/stride*sizeof(float)); return;}
    if(y_locate>=height) {memset(data,0,width/stride*sizeof(float)); return;}

    tdata = tdata+y_locate*width;

    int i=x_locate;int n=0;
    for(;i<0                        ;i+=stride) {data[n]=0;       n++;}
    for(;i<MIN(width,width+x_locate);i+=stride) {data[n]=tdata[i];n++;}
    for(;i<width+x_locate           ;i+=stride) {data[n]=0;       n++;}
}

void _SetRowData(float *tdata,float *data,int height,int width,int y_locate,int x_locate,int stride)
{
    if(y_locate<0)       return;
    if(y_locate>=height) return;

    tdata = tdata+y_locate*width;

    int i=x_locate;int n=0;
    for(;i<0                        ;i+=stride) {                  n++;}
    for(;i<MIN(width,width+x_locate);i+=stride) {tdata[i]+=data[n];n++;}
  //for(;i<width+x_locate           ;i+=stride) {                  n++;}
}

void ConvTensorToMatData0(MTensor *tns,int bc,float *mdata,int knl_height,int knl_width)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    // int mheight= knl_height*knl_width*channel+1;
    int mwidth = width*height;

    int x1=0-(knl_width /2);int x2=x1+knl_width -1;
    int y1=0-(knl_height/2);int y2=y1+knl_height-1;

    float *tdata = tns->data[bc];
    for(int c=0;c<channel;c++)
    {
        int k=0;
        {
            for(int i=x1;i<=x2;i++)
            {
                float *data = mdata + k*mwidth;k++;
                for(int n=y1;n<y1+height;n++)
                {
                    _GetRowData0(tdata,data,height,width,n,i);
                    data+=width;
                }
            }
        }
        for(int j=y1+1+height-1;j<=y2+height-1;j++)
        {
            for(int i=x1;i<=x2;i++)
            {
                float *data = mdata + k*mwidth;k++;
                memcpy(data,data-knl_width*mwidth+width,(mwidth-width)*sizeof(float));
                _GetRowData0(tdata,data+(mwidth-width),height,width,j,i);
            }
        }
        tdata+=width*height;
        mdata+=mwidth*knl_height*knl_width;
    }
    for(int i=0;i<mwidth;i++) mdata[i]=1.0;
}

void ConvTensorToMatData(MTensor *tns,int bc,float *mdata,int knl_height,int knl_width,int y_stride,int x_stride)
{
    if((x_stride==1)&&(y_stride==1)) {ConvTensorToMatData0(tns,bc,mdata,knl_height,knl_width);return;}

    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width = width/x_stride;
    int out_height=height/y_stride;
    // int mheight= knl_height*knl_width*channel+1;
    int mwidth = out_width*out_height;

    int x1=0-(knl_width /2);int x2=x1+knl_width -1;
    int y1=0-(knl_height/2);int y2=y1+knl_height-1;

    float *tdata = tns->data[bc];
    for(int c=0;c<channel;c++)
    {
        int k=0;
        for(int j=y1;j<=y2;j++)for(int i=x1;i<=x2;i++)
        {
            float *data = mdata + k*mwidth;k++;
            for(int n=j;n<j+height;n+=y_stride)
            {
                _GetRowData(tdata,data,height,width,n,i,x_stride);
                data+=out_width;
            }
        }
        tdata+=width*height;
        mdata+=mwidth*knl_height*knl_width;
    }
    for(int i=0;i<mwidth;i++) mdata[i]=1.0;
}

void ConvMatDataToTensor(float *mdata,MTensor *tns,int bc,int knl_height,int knl_width,int y_stride,int x_stride)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width = width/x_stride;
    int out_height=height/y_stride;
    // int mheight= knl_height*knl_width*channel+1;
    int mwidth = out_height*out_width;

    // printf("height=%d,width=%d,channel=%d,out_height=%d,out_width=%d\n",height,width,channel,out_height,out_width);

    int x1=0-(knl_width /2);int x2=x1+knl_width -1;
    int y1=0-(knl_height/2);int y2=y1+knl_height-1;

    float *tdata = tns->data[bc];

    for(int c=0;c<channel;c++)
    {
        int k=0;
        for(int j=y1;j<=y2;j++)for(int i=x1;i<=x2;i++)
        {
            float *data = mdata + k*mwidth;k++;
            for(int n=j;n<j+height;n+=y_stride)
            {
                _SetRowData(tdata,data,height,width,n,i,x_stride);
                data+=out_width;
            }
        }
        tdata+=width*height;
        mdata+=mwidth*knl_height*knl_width;
    }
}

struct TensorConvPara
{
    MLayer *prev;

    int knl_num;
    int knl_height;
    int knl_width;

    int x_stride;
    int y_stride;

    int res_valid;

    float rate;
    float decay;
    float momentum;
};
void *mTensorConvPara(MFile *ini,char *name)
{
    struct TensorConvPara *para = (struct TensorConvPara *)mMalloc(sizeof(struct TensorConvPara));

    para->prev = mNetworkLayer(ini,mINIRead(ini,name,"prev"));
    mException((para->prev == NULL),EXIT,"invalid prev");
    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    para->knl_num   = 1; mINIRead(ini,name,"knl_num"   ,"%d",&(para->knl_num   ));
    para->knl_height= 1; mINIRead(ini,name,"knl_height","%d",&(para->knl_height));
    para->knl_width = 1; mINIRead(ini,name,"knl_width" ,"%d",&(para->knl_width ));
    para->x_stride  = 1; mINIRead(ini,name,"x_stride"  ,"%d",&(para->x_stride  ));
    para->y_stride  = 1; mINIRead(ini,name,"y_stride"  ,"%d",&(para->y_stride  ));

    para->rate    =0.001;if(mINIRead(ini,name,"rate"    ,"%f",&(para->rate    ))==NULL) mINIRead(ini,"para","rate"    ,"%f",&(para->rate    ));
    para->momentum=0.9  ;if(mINIRead(ini,name,"momentum","%f",&(para->momentum))==NULL) mINIRead(ini,"para","momentum","%f",&(para->momentum));
    para->decay   =0.01 ;if(mINIRead(ini,name,"decay"   ,"%f",&(para->decay   ))==NULL) mINIRead(ini,"para","decay"   ,"%f",&(para->decay   ));
    mException((para->decay<0.0f)||(para->decay>=1.0f),EXIT,"invalid para decay");

    return para;
}

struct HandleTensorConv
{
    float *mat;
    float *kernel;
    float *update;
};
void endTensorConv(void *info)
{
    struct HandleTensorConv *handle = (struct HandleTensorConv *)info;
    if(handle->mat   != NULL) mFree(handle->mat);
    if(handle->kernel!= NULL) mFree(handle->kernel);
    if(handle->update!= NULL) mFree(handle->update);
}
#define HASH_TensorConv 0x9087d39c
void TensorConvSet(MLayer *layer)
{
    if(layer->state != DFLT) return;
    // mException(strcmp("Conv",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorConvPara *para = (struct TensorConvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out= layer->tns;

    MHandle *hdl=mHandle(out,TensorConv);
    struct HandleTensorConv *handle = (struct HandleTensorConv *)(hdl->handle);

    int out_height= in->height/para->y_stride;
    int out_width = in->width /para->x_stride;
    int mwidth = (out_height*out_width);
    int mheight= para->knl_height*para->knl_width*in->channel+1;
    int data_size = para->knl_num*mheight;

    mTensorRedefine(out,in->batch,para->knl_num,out_height,out_width,NULL);
    if(morn_network_flag == MORN_TRAIN)
    {
        if(INVALID_TENSOR(res)) mTensorRedefine(res,in->batch,in->channel,in->height,in->width,in->data);
        else
        {
            // printf("llllllllllllllllllllllllll layer->name=%s\n",layer->name);
            mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);
        }

        if(morn_network_flag==MORN_TRAIN)
        {
            if(handle->update != NULL) mFree(handle->update);
            handle->update =(float *)mMalloc(data_size*sizeof(float));
            memset(handle->update,0,data_size*sizeof(float));
        }
    }

    if(handle->kernel !=NULL) mFree(handle->kernel);
    handle->kernel = (float *)mMalloc(data_size*sizeof(float));

    if(morn_network_parafile==NULL)
    {
        float scale = sqrt(2.0f/mheight);
        for(int i=0;i<data_size;i++)
            handle->kernel[i] = mNormalRand(0.0f,1.0f)*scale;
    }
    else
    {
        mNetworkParaRead(layer,"kernel",handle->kernel,data_size*sizeof(float));
    }

    if(handle->mat!=NULL) mFree(handle->mat);
    handle->mat = (float *)mMalloc(mheight*mwidth*sizeof(float));

    hdl->valid = 1;
}

void mTensorConvForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    struct TensorConvPara *para = (struct TensorConvPara *)(layer->para);

    MTensor *in = para->prev->tns;
    MTensor *out=layer->tns;
    // printf("in->data[0][0]=%f\n",in->data[0][0]);

    TensorConvSet(layer);

    MHandle *hdl=mHandle(out,TensorConv);
    struct HandleTensorConv *handle = (struct HandleTensorConv *)(hdl->handle);

    int mwidth = (out->height*out->width);
    int mheight= para->knl_height*para->knl_width*in->channel+1;

    float *kernel_data= handle->kernel;
    float *in_data = handle->mat;

    for(int b=0;b<in->batch;b++)
    {
        // printf("ccccccccccccccc in->width=%d\n",in->width);
        ConvTensorToMatData(in,b,in_data,para->knl_height,para->knl_width,para->y_stride,para->x_stride);
        // printf("ccccccccccccccc in->width=%d\n",in->width);
        float *out_data = out->data[b];

        // printf("\nweight=\n");for(int ii=0;ii<100;ii++) printf("%f,",kernel_data[ii]);
        // printf("\ndata=\n");for(int ii=0;ii<100;ii++) printf("%f,",in_data[ii]*255);
        // printf("\ndata2352=\n");for(int ii=2352;ii<2352+100;ii++) printf("%f,",in_data[ii]*255);
        // printf("\ndata6272=\n");for(int ii=6272;ii<6272+100;ii++) printf("%f,",in_data[ii]*255);
        // printf("\ndata10192=\n");for(int ii=10192;ii<10192+100;ii++) printf("%f,",in_data[ii]*255);
        // printf("\ndata=\n");for(int ii=784*5;ii<784*5+100;ii++) printf("%f,",in_data[ii]*255);
        // printf("\ndata=\n");for(int ii=784*25;ii<784*25+100;ii++) printf("%f,",in_data[ii]*255);
        // printf("\ndata=\n");for(int ii=10000;ii<10000+100;ii++) printf("%f,",in_data[ii]*255);

        // float sum=0;
        // for(int ii=0;ii<75;ii++)
        // {
        //     sum+=(kernel_data[ii]*in_data[mwidth*ii]);
        //     printf("a[ii]=%f,b[n*ii]=%f,sum=%f,ii=%d,ii*mwidth=%d\n",kernel_data[ii],in_data[mwidth*ii]*255,sum,ii,ii*mwidth);
        // }
        // printf("\nsum=%f\n",sum);

        // printf("para->knl_num=%d,mwidth=%d,mheight=%d\n",para->knl_num,mwidth,mheight);
        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
                    para->knl_num,mwidth,mheight,
                    1.0f,
                    kernel_data,mheight,
                        in_data,mwidth,
               0.0f,   out_data,mwidth);

        // printf("\nout=\n");for(int ii=00;ii<200;ii++) printf("%f,",out_data[ii]);
    }

    layer->state = MORN_FORWARD;
}

void mTensorConvBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    struct TensorConvPara *para = (struct TensorConvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->res;

    MHandle *hdl=mHandle(layer->tns,TensorConv);
    struct HandleTensorConv *handle = (struct HandleTensorConv *)(hdl->handle);
    mException((hdl->valid == 0),EXIT,"no forward operate");

    int mwidth = (out->height*out->width);
    int mheight= para->knl_height*para->knl_width*in->channel+1;

    float *kernel_data= handle->kernel;
    float *update_data= handle->update;
    float *    in_data= handle->mat;
    float *   res_data= handle->mat;

    mNetworkParaWrite(layer,"kernel",kernel_data,para->knl_num*mheight*sizeof(float));

    float beta = para->momentum;
    for(int b=0;b<out->batch;b++)
    {
        ConvTensorToMatData(in,b,in_data,para->knl_height,para->knl_width,para->y_stride,para->x_stride);
        float *out_data = out->data[b];

        // printf("m=%d,n=%d,k=%d\n",para->knl_num,mheight,mwidth);
        // printf("\ndelta=\n");for(int ii=00;ii<200;ii++) printf("%f,",0-2*out_data[ii]);
        // printf("\nin=\n");for(int ii=00;ii<200;ii++) printf("%f,",in_data[ii]);
        // printf("\nupdate0=\n");for(int ii=00;ii<75;ii++) printf("%f,",0-update_data[ii]*2*beta);
        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasTrans,
                    para->knl_num,mheight,mwidth,
                    1.0f,//1.0f/mwidth,
                       out_data,mwidth,
                        in_data,mwidth,
                    beta,
                    update_data,mheight);

        // printf("\nupdate1=\n");for(int ii=00;ii<75;ii++) printf("%f,",0-update_data[ii]*2);
        // printf("\nupdate=\n");for(int ii=760;ii<760+75;ii++) printf("%f,",update_data[ii]);
        // printf("\nbias_updates=\n");for(int ii=0;ii<32;ii++) printf("%f,",update_data[ii*76+75]);

        beta = 1.0;
    }
    // for(int i=0;i<28;i++)printf("%f,",in_data[i]);printf("\n");
    // for(int i=0;i<8;i++) printf("%f,",kernel_data[i]);
    // printf("aaaaaaaaaaa,para->knl_num=%d,mwidth=%d\n",para->knl_num,mwidth);

    if(para->res_valid)
    {
        // if(para->prev->state == MORN_FORWARD)
        // {
        //     printf("conv res set 0,conv res set 0,conv res set 0,conv res set 0,conv res set 0res->batch=%d\n",res->batch);
        //     for(int b=0;b<res->batch;b++)

        // }

        for(int b=0;b<in->batch;b++)
        {
            float *out_data = out->data[b];

            if(para->prev->state == MORN_FORWARD)
                memset(res->data[b],0,in->height*in->width*in->channel*sizeof(float));

            // printf("\nconvdelta=\n");for(int ii=200;ii<400;ii++) printf("%f,",out_data[ii]);
            // printf("\nconvkernel=\n");for(int ii=00;ii<200;ii++) printf("%f,",kernel_data[ii]);

            cblas_sgemm(CblasRowMajor,CblasTrans,CblasNoTrans,
                        mheight,mwidth,para->knl_num,
                        1.0f,
                        kernel_data,mheight,
                           out_data,mwidth,
                    0.0,   res_data,mwidth);

            // printf("\nconvdeltaout=\n");for(int ii=00;ii<200;ii++) printf("%f,",0-res_data[ii]*2);

            ConvMatDataToTensor(res_data,res,b,para->knl_height,para->knl_width,para->y_stride,para->x_stride);

            // printf("mheight=%d,mwidth=%d\n",mheight,mwidth);
            // printf("\nconvimd=\n");for(int ii=200;ii<400;ii++) printf("%f,",res->data[0][ii]);
        }

        para->prev->state = MORN_BACKWARD;
    }

    // printf("\nkernel_data0=\n");for(int ii=0;ii<200;ii++) printf("%f,",kernel_data[ii]);
    // printf("\nupdate_data=\n");for(int ii=0;ii<200;ii++) printf("%f,",update_data[ii]);
    // printf("%s:update_data=%f,kernel_data=%f,%f\n",layer->name,update_data[10],kernel_data[10],kernel_data[10]-(para->rate/(float)(in->batch))*update_data[10]);
    cblas_saxpby(para->knl_num*mheight,
                 (0.0f-(para->rate/(float)(in->batch))),update_data,1,
                 (1.0f-(para->decay*para->rate))       ,kernel_data,1);

    // printf("\nkernel_data1=\n");for(int ii=0;ii<200;ii++) printf("%f,",kernel_data[ii]);
}

//////////////////////////////////////////////////////////////////////////////////////////////////

void GroupConvTensorToMatData(MTensor *tns,int bc,float *mdata,int knl_channel,int knl_height,int knl_width,int c_stride,int y_stride,int x_stride)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width  =  width/x_stride;
    int out_height = height/y_stride;
    int out_channel=(channel-knl_channel/2+1)/c_stride;
    int mwidth = knl_height*knl_width*knl_channel+1;
    int mheight= out_height*out_width;

    int x0=(width -(out_width -1)*x_stride)/2;
    int y0=(height-(out_height-1)*y_stride)/2;

    float *tdata = tns->data[bc];
    int tsize = tns->height*tns->width;

    int k,i,j,c;
    for(k=0;k<out_channel;k++)
    {
        for(j=0;j<out_width;j++)
        {
            int n=y0           -knl_height/2;
            int m=x0+j*x_stride-knl_width /2;
            for(i=0;i<mwidth-1;i+=knl_channel)
            {
                int h= i/(knl_width*knl_channel)             +n;if(h<0)h=0;else if(h>=height)h=height-1;
                int w=(i%(knl_width*knl_channel))/knl_channel+m;if(w<0)w=0;else if(w>= width)w= width-1;
                for(c=0;c<knl_channel;c++)
                    mdata[(k*mheight+j)*mwidth+i+c]=tdata[MIN((c+k*c_stride),channel-1)*tsize+h*width+w];
            }
            mdata[(k*mheight+j)*mwidth+mwidth-1]=1.0f;
        }
        for(j=out_width;j<mheight;j++)
        {
            int num=MAX(0,(knl_height-y_stride))*knl_width*knl_channel;
            if(num>0) memcpy(mdata+(k*mheight+j)*mwidth,mdata+((k*mheight+j)-out_width)*mwidth+y_stride*knl_width*knl_channel,num*sizeof(float));

            int n=y0+j/out_width*y_stride-knl_height/2;
            int m=x0+j%out_width*x_stride-knl_width /2;
            for(i=num;i<mwidth-1;i+=knl_channel)
            {
                int h= i/(knl_width*knl_channel)             +n;if(h<0)h=0;else if(h>=height)h=height-1;
                int w=(i%(knl_width*knl_channel))/knl_channel+m;if(w<0)w=0;else if(w>= width)w= width-1;
                for(c=0;c<knl_channel;c++)
                    mdata[(k*mheight+j)*mwidth+i+c]=tdata[MIN((c+k*c_stride),channel-1)*tsize+h*width+w];
            }
            mdata[(k*mheight+j)*mwidth+mwidth-1]=1.0f;
        }
    }
}

void GroupConvMatDataToTensor(float *mdata,MTensor *tns,int bc,int knl_channel,int knl_height,int knl_width,int c_stride,int y_stride,int x_stride)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width = width/x_stride;
    int out_height=height/y_stride;
    int out_channel=(channel-knl_channel/2+1)/c_stride;
    int mwidth = knl_height*knl_width*knl_channel+1;
    int mheight= out_height*out_width;

    int x0=(width -(out_width -1)*x_stride)/2;
    int y0=(height-(out_height-1)*y_stride)/2;

    float *tdata = tns->data[bc];
    int tsize = tns->height*tns->width;

    int k,i,j,c;
    for(k=0;k<out_channel;k++)for(j=0;j<mheight;j++)
    {
        int n=y0+j/out_width*y_stride+knl_height/2+1-knl_height;
        int m=x0+j%out_width*x_stride+knl_width /2+1-knl_width ;
        for(i=0;i<mwidth-1;i+=knl_channel)
        {
            int h= i/(knl_width*knl_channel)             +n;if(h<0)h=0;else if(h>=height)h=height-1;
            int w=(i%(knl_width*knl_channel))/knl_channel+m;if(w<0)w=0;else if(w>= width)w= width-1;
            for(c=0;c<knl_channel;c++)
                tdata[MIN((c+k*c_stride),channel-1)*tsize+h*width+w]+=mdata[(k*mheight+j)*mwidth+i+c];
        }
    }
}

struct TensorGroupConvPara
{
    MLayer *prev;

    int knl_num;
    int knl_height;
    int knl_width;

    int x_stride;
    int y_stride;

    int res_valid;

    float rate;
    float decay;
    float momentum;

    int knl_channel;
    int c_stride;
};
void *mTensorGroupConvPara(MFile *ini,char *name)
{
    struct TensorGroupConvPara *para = (struct TensorGroupConvPara *)mMalloc(sizeof(struct TensorConvPara));

    char *value = mINIRead(ini,name,"prev");
    para->prev = mNetworkLayer(ini,value);
    mException((para->prev == NULL),EXIT,"invalid prev");

    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    value = mINIRead(ini,name,"knl_num");
    if(value != NULL) para->knl_num= atoi(value);else para->knl_num= 1;

    value = mINIRead(ini,name,"knl_channel");
    if(value != NULL) para->knl_channel= atoi(value);else para->knl_channel= DFLT;

    value = mINIRead(ini,name,"knl_height");
    if(value != NULL) para->knl_height= atoi(value);else para->knl_height= 1;

    value = mINIRead(ini,name,"knl_width");
    if(value != NULL) para->knl_width= atoi(value);else para->knl_width= 1;

    value = mINIRead(ini,name,"c_stride");
    if(value != NULL) para->c_stride= atoi(value);else para->c_stride= para->knl_channel;

    value = mINIRead(ini,name,"x_stride");
    if(value != NULL) para->x_stride= atoi(value);else para->x_stride= 1;

    value = mINIRead(ini,name,"y_stride");
    if(value != NULL) para->y_stride= atoi(value);else para->y_stride= 1;

    value = mINIRead(ini,name,"rate");
    if(value != NULL) para->rate = atof(value);
    else
    {
        value = mINIRead(ini,"para","rate");
        if(value != NULL) para->rate = atof(value);
        else              para->rate = 0.001;
    }

    value = mINIRead(ini,name,"decay");
    if(value != NULL) para->decay = atof(value);
    else
    {
        value = mINIRead(ini,"para","decay");
        if(value != NULL) para->decay = atof(value);
        else              para->decay = 0.01;
    }
    mException((para->decay<0.0f)||(para->decay>=1.0f),EXIT,"invalid para decay");

    value = mINIRead(ini,name,"momentum");
    if(value != NULL) para->momentum = atof(value);
    else
    {
        value = mINIRead(ini,"para","momentum");
        if(value != NULL) para->momentum = atof(value);
        else              para->momentum = 0.9;
    }

    return para;
}

struct HandleTensorGroupConv
{
    float *mat;
    float *kernel;
    float *update;
};
void endTensorGroupConv(void *info)
{
    struct HandleTensorGroupConv *handle = (struct HandleTensorGroupConv *)info;
    if(handle->mat   != NULL) mFree(handle->mat);
    if(handle->kernel!= NULL) mFree(handle->kernel);
    if(handle->update!= NULL) mFree(handle->update);
}
#define HASH_TensorGroupConv 0x82866393
void TensorGroupConvSet(MLayer *layer)
{
    if(layer->state != DFLT) return;

    struct TensorGroupConvPara *para = (struct TensorGroupConvPara *)(layer->para);
    if(para->knl_channel<=0)
    {
        layer->type_index=mTensorRegisterIndex("Conv");
        return;
    }

    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out= layer->tns;

    MHandle *hdl=mHandle(out,TensorGroupConv);
    struct HandleTensorGroupConv *handle = (struct HandleTensorGroupConv *)(hdl->handle);

    int out_height= in->height/para->y_stride;
    int out_width = in->width /para->x_stride;
    int out_channel=(in->channel-para->knl_channel/2+1)/para->c_stride;

    int mheight = (out_height*out_width*out_channel);
    int mwidth = para->knl_height*para->knl_width*para->knl_channel+1;
    int data_size = para->knl_num*mwidth;

    mTensorRedefine(out,in->batch,para->knl_num*out_channel,out_height,out_width,NULL);
    if(morn_network_flag == MORN_TRAIN)
    {
        if(INVALID_TENSOR(res)) mTensorRedefine(res,in->batch,in->channel,in->height,in->width,in->data);
        else                    mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);

        if(handle->update != NULL) mFree(handle->update);
        handle->update =(float *)mMalloc(data_size*sizeof(float));
        memset(handle->update,0,data_size*sizeof(float));
    }

    if(handle->kernel !=NULL) mFree(handle->kernel);
    handle->kernel = (float *)mMalloc(data_size*sizeof(float));

    if(morn_network_parafile==NULL)
    {
        float scale = sqrt(2.0f/mwidth);
        for(int i=0;i<data_size;i++)
            handle->kernel[i] = scale*mNormalRand(0.0f,1.0f);
    }
    else
    {
        mNetworkParaRead(layer,"kernel",handle->kernel,data_size*sizeof(float));
    }

    if(handle->mat!=NULL) mFree(handle->mat);
    handle->mat = (float *)mMalloc(mheight*mwidth*sizeof(float));

    hdl->valid = 1;
}

void mTensorGroupConvForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    struct TensorGroupConvPara *para = (struct TensorGroupConvPara *)(layer->para);
    if(para->knl_channel<=0) return mTensorConvForward(layer);

    mException(strcmp("GroupConv",mLayerType(layer)),EXIT,"invalid layer type");

    MTensor *in = para->prev->tns;
    MTensor *out=layer->tns;

    TensorGroupConvSet(layer);

    MHandle *hdl=mHandle(out,TensorGroupConv);
    struct HandleTensorGroupConv *handle = (struct HandleTensorGroupConv *)(hdl->handle);

    int out_channel=(in->channel-para->knl_channel/2+1)/para->c_stride;
    int mheight = (out->height*out->width*out_channel);
    int mwidth = para->knl_height*para->knl_width*para->knl_channel+1;

    float *kernel_data= handle->kernel;
    float *in_data = handle->mat;

    for(int b=0;b<in->batch;b++)
    {
        GroupConvTensorToMatData(in,b,in_data,para->knl_channel,para->knl_height,para->knl_width,para->c_stride,para->y_stride,para->x_stride);
        float *out_data = out->data[b];

        in_data[mwidth-1]=1.0f;

        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasTrans,
                    para->knl_num,mheight,mwidth,
                    1.0f,
                    kernel_data,mwidth,
                        in_data,mwidth,
               0.0f,   out_data,mheight);
    }

    layer->state = MORN_FORWARD;
}

void mTensorGroupConvBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    struct TensorGroupConvPara *para = (struct TensorGroupConvPara *)(layer->para);
    if(para->knl_channel<=0) return mTensorConvBackward(layer);

    mException(strcmp("GroupConv",mLayerType(layer)),EXIT,"invalid layer type");

    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->res;

    MHandle *hdl=mHandle(layer->tns,TensorGroupConv);
    struct HandleTensorGroupConv *handle = (struct HandleTensorGroupConv *)(hdl->handle);
    mException((hdl->valid == 0),EXIT,"no forward operate");

    int out_channel=(in->channel-para->knl_channel/2+1)/para->c_stride;
    int mheight = (out->height*out->width*out_channel);
    int mwidth = para->knl_height*para->knl_width*para->knl_channel+1;

    float *kernel_data= handle->kernel;
    float *update_data= handle->update;
    float *    in_data= handle->mat;
    float *   res_data= handle->mat;

    mNetworkParaWrite(layer,"kernel",kernel_data,para->knl_num*mwidth*sizeof(float));

    for(int b=0;b<in->batch;b++)
    {
        GroupConvTensorToMatData(in,b,in_data,para->knl_channel,para->knl_height,para->knl_width,para->c_stride,para->y_stride,para->x_stride);
        float *out_data = out->data[b];

        in_data[mwidth-1]=1.0f;

        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
                    para->knl_num,mwidth,mheight,
                    1.0f,
                       out_data,mheight,
                        in_data,mwidth,
                    (b==0)?para->momentum:1.0f,
                    update_data,mwidth);
    }
    cblas_saxpby(para->knl_num*mwidth,
                 (0.0f-(para->rate/(float)(in->batch))),update_data,1,
                 (1.0f-(para->rate*para->decay))       ,kernel_data,1);

    if(para->res_valid==0) return;

    if(para->prev->state == MORN_FORWARD)
    {
        for(int b=0;b<res->batch;b++)
            memset(res->data[b],0,in->height*in->width*in->channel*sizeof(float));
        para->prev->state = MORN_BACKWARD;
    }

    for(int b=0;b<in->batch;b++)
    {
        float *out_data = out->data[b];

        cblas_sgemm(CblasRowMajor,CblasTrans,CblasNoTrans,
                    mheight,mwidth,para->knl_num,
                    1.0f,
                       out_data,mheight,
                    kernel_data,mwidth,
                0.0,   res_data,mwidth);

        GroupConvMatDataToTensor(res_data,res,b,para->knl_channel,para->knl_height,para->knl_width,para->c_stride,para->y_stride,para->x_stride);
    }
}


text:
/**
 */
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_multiroots.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spc_trace_functions.h"
#include "disp_conf.h"
#include "trace_conf.h"
#include "crossdisp_utils.h"

/*
 * Function: print_state
 * The function prints a status report
 * on a multi-dimensional root finding
 * solver onto the screen
 *
 * Parameters:
 * @param iter - the iteration number
 * @param s    - the multi-d solver
 */

void print_state (size_t iter, gsl_multiroot_fsolver *s)
{
  printf ("iter = %3zu x = % .3f % .3f "
          "f(x) = % .3e % .3e\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->f, 0),
          gsl_vector_get (s->f, 1));
  // print: current x-value, y-value, x-error, y-error
}

/*
 * Function: drizzle_distort
 * This function can be passed to a gsl multi-d function.
 * It defines the function:
 * x'=f'(x,y) = f(x,y)-x_0
 * y'=g'(x,y) = g(x,y)-y_0
 * with f(x,y) and g(x,y) the equations to undistort
 * the coordinates x and y, respectively and the
 * constants x_0 and y_0.
 * Solving for f'(x,y)=g'(x,y)=0.0
 * means to find (x,y) such that
 * f(x,y) = x_0 and g(x,y) = y_0.
 *
 * Parameters:
 * @param x      - the input x,y
 * @param params - the parameters for the function
 * @param f      - the result values in x,y
 *
 * Returns:
 * @return GSL_SUCCESS - always
 */
int
drizzle_distort(const gsl_vector * x, void *params, gsl_vector * f)
{
  // extract the components from the parameter-structure
  gsl_matrix *coeffs = ((drz_pars *) params)->coeffs;
  px_point  npixels = ((drz_pars *) params)->npixels;
  d_point offset = ((drz_pars *) params)->offset;

  d_point xy_in;
  d_point xy_out;

  // extract the input
  xy_in.x = gsl_vector_get(x,0);
  xy_in.y = gsl_vector_get(x,1);

  // undistort the input
  xy_out = get_drz_position(xy_in, coeffs, npixels);

  // subtract the input
  xy_out.x = xy_out.x-offset.x;
  xy_out.y = xy_out.y-offset.y;

  // set the output
  gsl_vector_set(f, 0, xy_out.x);
  gsl_vector_set(f, 1, xy_out.y);

  // return always the same
  return GSL_SUCCESS;
}

/*
 * Function: distort_point
 * The function finds the distorted coordinates
 * for a given undistorted coordinate and the
 * transformations to get the undistorted coordinates.
 * This function makes the inverse transformation
 * to the drizzle transformation.
 *
 * Parameters:
 * @param coeffs   - the drizzle coefficients
 * @param pixmax   - the image dimensions
 * @param xy_image - the undistorted (x,y)
 *
 * Returns:
 * @return xy_ret - the distorted (x,y)
 */
d_point
distort_point(gsl_matrix *coeffs, const px_point pixmax, d_point xy_image)
{
  const gsl_multiroot_fsolver_type *msolve_type;
  gsl_multiroot_fsolver *msolve;

  int status;
  size_t iter = 0;

  //const size_t n = 2;

  gsl_multiroot_function mult_func;

  drz_pars   *drzpars;
  gsl_vector *xy_in = gsl_vector_alloc(2);
  gsl_vector *xy_out = gsl_vector_alloc(2);
  d_point xy_ret;

  // set up the parameters for the
  // multi-d function
  drzpars = (drz_pars *) malloc(sizeof(drz_pars));
  drzpars->coeffs = coeffs;
  drzpars->offset = xy_image;
  drzpars->npixels = pixmax;

  // set up the multi-d function
  mult_func.f = &drizzle_distort;
  mult_func.n = 2;
  mult_func.params = drzpars;

  // set the starting coordinates
  gsl_vector_set(xy_in, 0, xy_image.x);
  gsl_vector_set(xy_in, 1, xy_image.y);

  // allocate and initialize the multi-d solver
  msolve_type = gsl_multiroot_fsolver_dnewton;
  msolve = gsl_multiroot_fsolver_alloc (msolve_type, 2);
  gsl_multiroot_fsolver_set (msolve, &mult_func, xy_in);

  //  print_state (iter, msolve);

  // iterate
  do
  {
    // count the number of iterations
    iter++;

    // do an iteration
    status = gsl_multiroot_fsolver_iterate (msolve);

    //    print_state (iter, msolve);

    // check if solver is stuck
    if (status)
      break;

    // evaluate the iteration
    status = gsl_multiroot_test_residual (msolve->f, 1e-7);
  }
  while (status == GSL_CONTINUE && iter < 1000);
  // chek for the break conditions

  // transfer the result to the return struct
  xy_ret.x = gsl_vector_get(msolve->x,0);
  xy_ret.y = gsl_vector_get(msolve->x,1);

  // deallocate the different structures
  gsl_multiroot_fsolver_free (msolve);
  gsl_vector_free(xy_in);
  gsl_vector_free(xy_out);

  // return the result
  return xy_ret;
}

/*
 * Function: undistort_point
 *
 * Parameters:
 * @param coeffs   - the drizzle coefficients
 * @param pixmax   - the image dimensions
 * @param xy_image - the undistorted (x,y)
 *
 * Returns:
 * @return xy_new - the undistorted (x,y)
 */
d_point
undistort_point(gsl_matrix *coeffs, const px_point pixmax, d_point xy_image)
{
  d_point xy_new;

  xy_new = get_drz_position(xy_image, coeffs, pixmax);

  return xy_new;
}

/*
 * Function: get_crossdisp_matrix
 * The function extracts the drizzle coefficients stored in the
 * keyword headers of an image extension and creates a matrix
 * to store those coefficients. The matrix is returned.
 * Currently the drizzle coefficients are ALWAYS stored in the
 * primary header of the grism images. For this reason
 * all keywords are read from extension '1' (hardcoded).
 *
 * Parameters:
 * @param  filename   - the image filename
 * @param  sci_numext - the extension number
 *
 * Returns:
 * @return ret - the matrix with the drizzle coefficients
 */
gsl_matrix  *
get_crossdisp_matrix(char * filename, int sci_numext)
{
  int ncoeffs;
  int i;
  gsl_matrix * ret;
  //px_point    pixmax;
  char templt[FLEN_CARD];
  float acoeff, tmp;

  // get the number of coefficients
  tmp = get_float_from_keyword(filename, 1, "DRZCNUM");
  //  tmp = get_float_from_keyword(filename, sci_numext, "DRZCNUM");
  if (isnan(tmp)){
    // in case that the keyword does not exist,
    // find a nice way out of the door
    aXe_message(aXe_M_FATAL, __FILE__, __LINE__,
                "Could not find drizzle keywords in: %s\n", filename);
    ret = gsl_matrix_alloc(1,1);
    gsl_matrix_set(ret, 0,0,GSL_NAN);
    return ret;
  }
  else {
    // cast the value to int
    ncoeffs = (int)tmp;
  }

  // allocate the matrix, set it to the default
  ret = gsl_matrix_alloc(2,ncoeffs);
  gsl_matrix_set_all (ret, 0.0);

  // go over the number of coefficients
  for (i=0; i<ncoeffs; i++){

    // set the name for the x-coefficient, get it and store it
    sprintf (templt, "DRZ%1iX%02i", sci_numext, i+1);
    //    acoeff = get_float_from_keyword(filename, sci_numext, templt);
    acoeff = get_float_from_keyword(filename, 1, templt);
    gsl_matrix_set(ret, 0,i,acoeff);

    // set the name for the y-coefficient, get it and store it
    sprintf (templt, "DRZ%1iY%02i", sci_numext, i+1);
    //    acoeff = get_float_from_keyword(filename, sci_numext, templt);
    acoeff = get_float_from_keyword(filename, 1, templt);
    gsl_matrix_set(ret, 1,i,acoeff);

  }

  // return the matrix
  return ret;
}

/*
 * Function: get_axis_scales
 * The function computes the distance scale along the major and
 * minor half axis of an object. The distance scale is subject
 * to changes due to the image distortion.
 *
 * Parameters:
 * @param actbeam   - the beam to get the drizzle scales for
 * @param drzcoeffs - the drizzle coefficients
 * @param pixmax    - the image dimensions
 *
 * Returns:
 * @return ret - the drizzle scales along major and minor half axis
 */
d_point
get_axis_scales(beam actbeam, gsl_matrix * drzcoeffs, px_point pixmax){

  d_point dispnt;
  d_point drz_ref;
  d_point drz_dis;

  d_point ret;

  // find a dummy point along the major half axis
  dispnt.x = actbeam.refpoint.x + cos(actbeam.aorient);
  dispnt.y = actbeam.refpoint.y + sin(actbeam.aorient);

  // compute the undistorted coos for the refpoint
  drz_ref = get_drz_position(actbeam.refpoint, drzcoeffs, pixmax);

  // compute the undistorted coos for the dummy point
  drz_dis = get_drz_position(dispnt, drzcoeffs, pixmax);

  // compute the scale along the major axis
  ret.x = sqrt((drz_ref.x-drz_dis.x)*(drz_ref.x-drz_dis.x) +
               (drz_ref.y-drz_dis.y)*(drz_ref.y-drz_dis.y));

  // find a dummy point along the minor half axis
  dispnt.x = actbeam.refpoint.x + cos(actbeam.aorient + M_PI/2.0);
  dispnt.y = actbeam.refpoint.y + sin(actbeam.aorient + M_PI/2.0);

  // compute the undistorted coos for the dummy point
  drz_dis = get_drz_position(dispnt, drzcoeffs, pixmax);

  // compute the scale along the minor axis
  ret.y = sqrt((drz_ref.x-drz_dis.x)*(drz_ref.x-drz_dis.x) +
               (drz_ref.y-drz_dis.y)*(drz_ref.y-drz_dis.y));

  // return the result
  return ret;
}

/*
 * Function: get_crossdisp_scale
 * The function computes the drizzle scale in cross-
 * dispersion direction for a given point and the
 * drizzle parameters of an image.
 *
 * Parameters:
 * @param trace     - the object trace
 * @param refpnt    - the reference point
 * @param drzcoeffs - the drizzle coefficients
 * @param pixmax    - the image dimensions
 *
 * Returns:
 * @return crscale - the drizzle scale in crossdispersion
 *                   direction
 */
double
get_crossdisp_scale(trace_func *trace, d_point refpnt, gsl_matrix * drzcoeffs,
                    px_point pixmax)
{
  double crscale;
  double phi_trace;
  //int beamID = 0;
  d_point dispnt;
  d_point drz_ref;
  d_point drz_dis;

  // get the trace angle
  phi_trace = atan2 (trace->deriv (0.0, trace->data), 1.0);

  // compute a dummy point along the cross-dispersion direction
  dispnt.x = refpnt.x + cos(phi_trace + M_PI/2.0);
  dispnt.y = refpnt.y + sin(phi_trace + M_PI/2.0);

  // get the undistorted positions of the refpoint and
  // the dummy point
  drz_ref = get_drz_position(refpnt, drzcoeffs, pixmax);
  drz_dis = get_drz_position(dispnt, drzcoeffs, pixmax);

  // compute the drizzle scale
  crscale = sqrt((drz_ref.x-drz_dis.x)*(drz_ref.x-drz_dis.x) +
                 (drz_ref.y-drz_dis.y)*(drz_ref.y-drz_dis.y));

  // return the drizzle scale
  return crscale;
}

/*
 * Function: evaln
 * The function evaluates the drizzle coefficients to return
 * either the x- or y-value in the undistorted system.
 * This is a translation of a fortran function with
 * the identical name inside the drizzle package.
 *
 * Parameters:
 * @param x         - the x-value in the distorted system
 * @param y         - the y-value in the distorted system
 * @param drzcoeffs - the drizzle coefficients
 * @param row       - which row of the coefficients to evaluate
 *
 * Returns:
 * @return ret - the undistorted coordinate value
 */
double
evaln(double x, double y, gsl_matrix * drzcoeffs, int row){

  double ret=0.0;
  int order;
  int n,m,nc,ncoeffs;
  order = 0;
  ncoeffs = 1;

  // derive the number of 'orders'
  while(ncoeffs+order+2 <= (int)drzcoeffs->size2){
    ncoeffs = ncoeffs+order+2;
    order++;
  }

  // add each individual term
  nc = 0;
  for (n=1; n < order+2; n++){
    for (m=1; m < n+1; m++){
      ret = ret + gsl_matrix_get(drzcoeffs,row,nc) * pow(x,(double)(n-m)) * pow(y,(double)(m-1));
      nc++;
    }
  }

  // give back the result
  return ret;
}

/*
 * Function: devalndx
 * The function numerically differentiates the
 * function 'evaln' with respect to the x-axis.
 *
 * Parameters:
 * @param x         - the x-value in the distorted system
 * @param y         - the y-value in the distorted system
 * @param drzcoeffs - the drizzle coefficients
 * @param row       - which row of the coefficients to evaluate
 *
 * Returns:
 * @return ret - the differential
 */
double
devalndx(double x, double y, gsl_matrix * drzcoeffs, int row){
  double ret=0.0;
  int order;
  int n,m,nc,ncoeffs;

  // derive the number of 'orders'
  order = 0;
  ncoeffs = 1;
  while(ncoeffs+order+2 <= (int)drzcoeffs->size2){
    ncoeffs = ncoeffs+order+2;
    order++;
  }

  // add each individual term
  nc = 0;
  for (n=1; n < order+2; n++){
    for (m=1; m < n; m++){
      ret = ret + gsl_matrix_get(drzcoeffs,row,nc) * (double)(n-m) * pow(x,(double)(n-m-1)) * pow(y,(double)(m-1));
      nc++;
    }
    nc++;
  }

  // return the result
  return ret;
}

/*
 * Function: devalndy
 * The function numerically differentiates the
 * function 'evaln' with respect to the y-axis.
 *
 * Parameters:
 * @param x         - the x-value in the distorted system
 * @param y         - the y-value in the distorted system
 * @param drzcoeffs - the drizzle coefficients
 * @param row       - which row of the coefficients to evaluate
 *
 * Returns:
 * @return ret - the differential
 */
double
devalndy(double x, double y, gsl_matrix * drzcoeffs, int row){
  double ret=0.0;
  int order;
  int n,m,nc,ncoeffs;

  // derive the number of 'orders'
  order = 0;
  ncoeffs = 1;
  while(ncoeffs+order+2 <= (int)drzcoeffs->size2){
    ncoeffs = ncoeffs+order+2;
    order++;
  }

  // add each individual term
  nc = 0;
  for (n=1; n < order+2; n++){
    nc++;
    for (m=2; m < n+1; m++){
      ret = ret + gsl_matrix_get(drzcoeffs,row,nc) * pow(x,(double)(n-m)) * (double)(m-1) * pow(y,(double)(m-2));
      nc++;
    }
  }

  // return the result
  return ret;
}

/*
 * Function: get_jacobian
 *
 * Parameters:
 * @param i         -
 * @param j         -
 * @param drzcoeffs - the drizzle coefficients
 * @param width     -
 * @param height    -
 *
 * Returns:
 * @return ret -
 */
gsl_matrix *
get_jacobian(int i, int j, gsl_matrix * drzcoeffs, int width, int height){

  double xr, yr;
  gsl_matrix * ret;

  ret = gsl_matrix_alloc(2,2);
  gsl_matrix_set_all (ret, 0.0);

  xr = (double)i - ((double)(width/2)+1.0);
  yr = (double)j - ((double)(height/2)+1.0);

  gsl_matrix_set(ret, 0, 0, devalndx( xr, yr, drzcoeffs, 0));
  gsl_matrix_set(ret, 0, 1, devalndy( xr, yr, drzcoeffs, 0));
  gsl_matrix_set(ret, 1, 0, devalndx( xr, yr, drzcoeffs, 1));
  gsl_matrix_set(ret, 1, 1, devalndy( xr, yr, drzcoeffs, 1));

  return ret;
}

/*
 * Function: get_det_jacobian
 *
 * Parameters:
 * @param i         -
 * @param j         -
 * @param drzcoeffs - the drizzle coefficients
 * @param width     -
 * @param height    -
 *
 * Returns:
 * @return ret -
 */
double
get_det_jacobian(int i, int j, gsl_matrix * drzcoeffs, int width, int height){

  gsl_matrix * jacob;
  double ret;

  jacob = get_jacobian(i,j,drzcoeffs,width,height);

  ret = gsl_matrix_get(jacob, 0, 0) * gsl_matrix_get(jacob, 1, 1)
       -gsl_matrix_get(jacob, 0, 1) * gsl_matrix_get(jacob, 1, 0);

  gsl_matrix_free(jacob);

  return ret;
}

/*
 * Function: get_drz_position
 * The function transforms distorted coordinates image to undistorted
 * image coordinates. As an assumption, the undistorted image
 * has the same dimension as the distorted.
 *
 * Parameters:
 * @param in_point  - the coordinates in the distorted frame
 * @param drzcoeffs - the drizzle coefficients
 * @param pixmax    - the image dimension of the distorted frame
 *
 * Returns:
 * @return: ou_point - the coordinates in the undistorted frame
 */
d_point
get_drz_position(d_point in_point, gsl_matrix *drzcoeffs, px_point pixmax)
{
  d_point ou_point;
  double xr, yr;
  double xn, yn;

  // transform the input coos into the drizzle system
  xr = in_point.x - ((double)(pixmax.x/2)+1.0);
  yr = in_point.y - ((double)(pixmax.y/2)+1.0);

  //  ou_point.x = evaln(xr,yr,drzcoeffs,0);
  //  ou_point.y = evaln(xr,yr,drzcoeffs,1);
  // get the undistorted coos in the drizzle system
  xn = evaln(xr,yr,drzcoeffs,0);
  yn = evaln(xr,yr,drzcoeffs,1);

  // transform back into the image system,
  // assuming identical image dimensions
  ou_point.x = xn + ((double)(pixmax.x/2)+1.0);
  ou_point.y = yn + ((double)(pixmax.y/2)+1.0);

  // return the undistorted coos
  return ou_point;
}

/*
 * Function: get_drz_position_free
 * The function transforms distorted coordinates image to undistorted
 * image coordinates. The size of the input and output image are free.
 *
 * Parameters:
 * @param in_point  - the coordinates in the distorted frame
 * @param drzcoeffs - the drizzle coefficients
 * @param pix_in    - the image dimension of the distorted frame
 * @param pix_out   - the image dimension of the distorted frame
 *
 * Returns:
 * @return: ou_point - the coordinates in the undistorted frame
 */
d_point
get_drz_position_free(d_point in_point, gsl_matrix *drzcoeffs,
                      px_point pix_in, px_point pix_out)
{
  d_point ou_point;
  double xr, yr;
  double xn, yn;

  // transform the input coos into the drizzle system
  xr = in_point.x - ((double)(pix_in.x/2)+1.0);
  yr = in_point.y - ((double)(pix_in.y/2)+1.0);

  //  ou_point.x = evaln(xr,yr,drzcoeffs,0);
  //  ou_point.y = evaln(xr,yr,drzcoeffs,1);
  // get the undistorted coos in the drizzle system
  xn = evaln(xr,yr,drzcoeffs,0);
  yn = evaln(xr,yr,drzcoeffs,1);

  // transform back into the image system,
  // assuming identical image dimensions
  ou_point.x = xn + ((double)(pix_out.x/2)+1.0);
  ou_point.y = yn + ((double)(pix_out.y/2)+1.0);

  // return the undistorted coos
  return ou_point;
}

trace_func *
get_tracefunc_at(char *filename, d_point p){

  int beamID = 0;
  trace_func  *trace;
  tracestruct *tstruct;

  tstruct = get_tracestruct_at_pos (filename, beamID, p);
  trace   = vector_to_trace_polyN ( tstruct->pol );

  return trace;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

#include "cosmocalc.h"

#define ODEGROWTH //the integral formulation is onyl correct for LCDM, any wCDM cosmology will not work

#ifdef ODEGROWTH

#define LNAINITGROWTH (log(AEXPN_MIN_GROWTH))
#define ABSERR 1e-8
#define RELERR 1e-8
#define HSTART 1e-4
#define ODE_TYPE gsl_odeiv2_step_rk8pd

static int growth_ode_sys_w0wa(double t, const double y[], double dydt[], void *params);

//ODE for growth from Komatsu et al. 2009, ApJS, 180, 330 (WMAP5 paper)
static int growth_ode_sys_w0wa(double t, const double y[], double dydt[], void *params)
{
  double a = exp(t); //t = log(a)
  double weffa = weff(a);
  double ha = hubble_noscale(a);
  double omegaDEa = cosmoData.OmegaL/ha/ha/pow(a,3.0*(1.0 + weffa));
  double omegaKa = cosmoData.OmegaK/a/a/ha/ha;

  dydt[0] = y[1];
  dydt[1] = (1.5*weffa*omegaDEa - 0.5*omegaKa - 2.5)*y[1] - (2.0*omegaKa + 1.5*(1.0 - weffa)*omegaDEa)*y[0];

  return GSL_SUCCESS;
}

double growth_function_exact(double a)
{
  gsl_odeiv2_system odesys = {growth_ode_sys_w0wa, NULL, 2, NULL};
  gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&odesys,ODE_TYPE,HSTART,ABSERR,RELERR);

  double lna_init;
  double lna_final,y[2];
  double ga,gnorm;
  int status;

  //do g(a) unormalized
  y[0] = 1.0;
  y[1] = 0.0;
  lna_init = LNAINITGROWTH;
  lna_final = log(a);

  status = gsl_odeiv2_driver_apply(d,&lna_init,lna_final,y);
  ga = y[0]*a;
  if(status != GSL_SUCCESS)
    {
      fprintf(stderr,"error in integrating growth function! a = %lf\n",a);
      assert(status == GSL_SUCCESS);
    }

  //do g(a = 1.0) to get normalization
  y[0] = 1.0;
  y[1] = 0.0;
  lna_init = LNAINITGROWTH;
  lna_final = 0.0;

  status = gsl_odeiv2_driver_apply(d,&lna_init,lna_final,y);
  gnorm = y[0];
  if(status != GSL_SUCCESS)
    {
      fprintf(stderr,"error in integrating growth function! a = %lf\n",a);
      assert(status == GSL_SUCCESS);
    }

  gsl_odeiv2_driver_free(d);

  return ga/gnorm;
}

double growth_function_exact_nonorm(double a)
{
  gsl_odeiv2_system odesys = {growth_ode_sys_w0wa, NULL, 2, NULL};
  gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&odesys,ODE_TYPE,HSTART,ABSERR,RELERR);

  double lna_init;
  double lna_final,y[2];
  double ga;
  int status;

  //do g(a) unormalized
  y[0] = 1.0;
  y[1] = 0.0;
  lna_init = LNAINITGROWTH;
  lna_final = log(a);

  status = gsl_odeiv2_driver_apply(d,&lna_init,lna_final,y);
  ga = y[0]*a;
  if(status != GSL_SUCCESS)
    {
      fprintf(stderr,"error in integrating growth function! a = %lf\n",a);
      assert(status == GSL_SUCCESS);
    }

  return ga;
}

double growth_function(double a)
{
  static int initFlag = 1;
  static int currCosmoNum;
  static double growth_function_norm;
  static gsl_spline *cosmocalc_growth_function_spline = NULL;
  static gsl_interp_accel *cosmocalc_growth_function_acc = NULL;

  double growth_function_table[COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH];
  double a_table[COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH];
  long i;

  gsl_odeiv2_system odesys = {growth_ode_sys_w0wa, NULL, 2, NULL};
  gsl_odeiv2_driver *d;
  double afact,lna_init;
  double lna_final,y[2];
  int status;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      //init ODE integrator
      d = gsl_odeiv2_driver_alloc_y_new(&odesys,ODE_TYPE,HSTART,ABSERR,RELERR);

      for(i=0;i<COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH;++i)
        {
	  //get scale factor
	  afact = log((1e-3 + AEXPN_MAX)/AEXPN_MIN_GROWTH)/(COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH-1.0)*((double) i) + log(AEXPN_MIN_GROWTH);
	  a_table[i] = afact;

	  //do growth function integration
	  //init ODE
	  y[0] = 1.0;
	  y[1] = 0.0;
	  lna_init = LNAINITGROWTH;
	  lna_final = afact;
  	  status = gsl_odeiv2_driver_apply(d,&lna_init,lna_final,y);
	  growth_function_table[i] = y[0]*exp(afact);
	  if(status != GSL_SUCCESS)
	    {
	      fprintf(stderr,"error in integrating growth function! a = %lf\n",a);
	      assert(status == GSL_SUCCESS);
	    }
	}

      //get normalization
      y[0] = 1.0;
      y[1] = 0.0;
      lna_init = LNAINITGROWTH;
      lna_final = 0.0;

      status = gsl_odeiv2_driver_apply(d,&lna_init,lna_final,y);
      growth_function_norm = y[0];
      if(status != GSL_SUCCESS)
	{
	  fprintf(stderr,"error in integrating growth function! a = %lf\n",1.0);
	  assert(status == GSL_SUCCESS);
	}

      //free ODE stuff
      gsl_odeiv2_driver_free(d);

      //init the spline and accelerators
      if(cosmocalc_growth_function_spline != NULL)
        gsl_spline_free(cosmocalc_growth_function_spline);
      cosmocalc_growth_function_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_growth_function_spline,a_table,growth_function_table,(size_t) (COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH));
      if(cosmocalc_growth_function_acc != NULL)
        gsl_interp_accel_reset(cosmocalc_growth_function_acc);
      else
        cosmocalc_growth_function_acc = gsl_interp_accel_alloc();
    }

  return gsl_spline_eval(cosmocalc_growth_function_spline,log(a),cosmocalc_growth_function_acc)/growth_function_norm;
}

#else

static double growth_function_integ_funct(double a, void *p);

/* function for integration using gsl integration */
static double growth_function_integ_funct(double a, void *p)
{
  double hubblefac = hubble_noscale(a);
  double adot = a*hubblefac;
  double alim = (*((double*)p));

  return 1.0/adot/adot/adot*hubble_noscale(alim);
}

double growth_function_exact(double a)
{
#define WORKSPACE_NUM 100000
#define ABSERR 0.0
#define RELERR 1e-8

  gsl_integration_workspace *workspace;
  gsl_function F;
  double result,abserr,afact,norm;

  workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);

  afact = a;
  F.function = &growth_function_integ_funct;
  F.params = &(afact);
  gsl_integration_qag(&F,0.0,afact,ABSERR,RELERR,(size_t) WORKSPACE_NUM,GSL_INTEG_GAUSS51,workspace,&result,&abserr);

  afact = 1.0;
  F.params = &(afact);
  gsl_integration_qag(&F,0.0,afact,ABSERR,RELERR,(size_t) WORKSPACE_NUM,GSL_INTEG_GAUSS51,workspace,&norm,&abserr);

  gsl_integration_workspace_free(workspace);

#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM

  return result/norm;
}

double growth_function(double a)
{
#define WORKSPACE_NUM 100000
#define ABSERR 1e-8
#define RELERR 1e-8

  static int initFlag = 1;
  static int currCosmoNum;
  static double growth_function_norm;
  static gsl_spline *cosmocalc_growth_function_spline = NULL;
  static gsl_interp_accel *cosmocalc_growth_function_acc = NULL;

  double growth_function_table[COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH];
  double a_table[COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH];
  long i;
  gsl_integration_workspace *workspace;
  gsl_function F;
  double result,abserr,afact;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);
      F.function = &growth_function_integ_funct;

      for(i=0;i<COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH;++i)
	{
	  afact = (1e-3 + AEXPN_MAX - AEXPN_MIN)/(COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH-1.0)*((double) i) + AEXPN_MIN;
	  a_table[i] = afact;
	  F.params = &(afact);
	  gsl_integration_qag(&F,0.0,afact,ABSERR,RELERR,(size_t) WORKSPACE_NUM,GSL_INTEG_GAUSS51,workspace,&result,&abserr);
	  growth_function_table[i] = result;
	}

      afact = 1.0;
      F.params = &(afact);
      gsl_integration_qag(&F,0.0,1.0,ABSERR,RELERR,(size_t) WORKSPACE_NUM,GSL_INTEG_GAUSS51,workspace,&growth_function_norm,&abserr);

      gsl_integration_workspace_free(workspace);

      //init the spline and accelerators
      if(cosmocalc_growth_function_spline != NULL)
	gsl_spline_free(cosmocalc_growth_function_spline);
      cosmocalc_growth_function_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_growth_function_spline,a_table,growth_function_table,(size_t) (COSMOCALC_GROWTH_FUNCTION_TABLE_LENGTH));
      if(cosmocalc_growth_function_acc != NULL)
	gsl_interp_accel_reset(cosmocalc_growth_function_acc);
      else
	cosmocalc_growth_function_acc = gsl_interp_accel_alloc();
    }

  return gsl_spline_eval(cosmocalc_growth_function_spline,a,cosmocalc_growth_function_acc)/growth_function_norm;

#undef WORKSPACE_NUM
#undef ABSERR
#undef RELERR
}
#endif

#undef ODEGROWTH

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#include <viaio/VImage.h>
#include <string.h>
#include <gsl/gsl_matrix.h>
#include <viaio/Vlib.h>
#include "gsl_utils.h"

void prewhite(gsl_matrix_float** pinvX,gsl_matrix_float** invM,gsl_matrix_float* X,int numlags)
{
  int n = X->size2;

  /* preparation for whitening */

  /* X^T */
  gsl_matrix_float* transX = gsl_matrix_float_alloc(X->size2, X->size1);
  gsl_matrix_float_transpose_memcpy(transX,X);

  /* identity */
  gsl_matrix_float *R = gsl_matrix_float_alloc(n,n);
  gsl_matrix_float_set_identity(R);

  /* pseudo inverse */
  *pinvX = fmat_PseudoInv(transX,NULL);

  /* temporary matrix buffer for binary matrix operations */
  gsl_matrix_float *result;
  gsl_matrix_float *result2;

  result = fmat_x_mat(transX,*pinvX,NULL);
  gsl_matrix_float_sub(R,result);
  gsl_matrix_float_free(result);

  /* result buffer */
  gsl_matrix_float *M;

  if(numlags == 1){
    /* create a n-x-n matrix Diag1 where the k=1 sub- and suberdiagonals
     * are 1 and the remaining entries are 0.
     * Equivalent matlab code:
     *
     * keep = 1:n;
     * indk1=((keep(2:n)-keep(1:n-1))==1);
     * Diag1=diag(indk1,1)+diag(indk1,-1);
     */
    gsl_matrix_float *Diag1 = gsl_matrix_float_alloc(n,n);
    gsl_matrix_float_set_zero(Diag1);
    /* set superdiagonal to 1 */
    gsl_vector_float_view superview = gsl_matrix_float_superdiagonal(Diag1,1);
    gsl_vector_float_set_all(&superview.vector,1);
    /* set subdiagonal to 1 */
    gsl_vector_float_view subview = gsl_matrix_float_subdiagonal(Diag1,1);
    gsl_vector_float_set_all(&subview.vector,1);

    float M11 = trace(R);
    result = fmat_x_mat(R, Diag1, NULL);
    float M12 = trace(result);
    float M21 = M12/2.0;
    fmat_x_mat(R,Diag1,result);
    result2 = fmat_x_mat(result,R,NULL);
    fmat_x_mat(result2,Diag1,result);
    float M22 = trace(result)/2;
    gsl_matrix_float_free(result);
    gsl_matrix_float_free(result2);

    M = gsl_matrix_float_alloc(2,2);

    /* M=[M11 M12; M21 M22] */
    gsl_matrix_float_set(M,0,0,M11);
    gsl_matrix_float_set(M,0,1,M12);
    gsl_matrix_float_set(M,1,0,M21);
    gsl_matrix_float_set(M,1,1,M22);

    gsl_matrix_float_free(Diag1);

  }
  else {

    int i,j;
    M = gsl_matrix_float_alloc(numlags+1,numlags+1);
    gsl_matrix_float_set_zero(M);
    gsl_matrix_float* Di = gsl_matrix_float_alloc(n,n);
    gsl_matrix_float* Dj = gsl_matrix_float_alloc(n,n);

    /* This double loop can be supported by a look-up-table for
     * the diagonal matrixes. */
    for(i=0;i<numlags+1;i++){
      gsl_matrix_float_set_zero(Di);
      /* set maindiagonal to 1 */
      if(i==0){
	gsl_vector_float_view mdiag = gsl_matrix_float_diagonal(Di);
	gsl_vector_float_set_all(&mdiag.vector,1);
      }
      /* set i. super- and subdiagonal to 1 */
      else {
	/* set superdiagonal to 1 */
	gsl_vector_float_view superview = gsl_matrix_float_superdiagonal(Di,i);
	gsl_vector_float_set_all(&superview.vector,1);
	/* set subdiagonal to 1 */
	gsl_vector_float_view subview = gsl_matrix_float_subdiagonal(Di,i);
	gsl_vector_float_set_all(&subview.vector,1);
      }
      for(j=0;j<numlags+1;j++){
	gsl_matrix_float_set_zero(Dj);
	/* set maindiagonal to 1 */
	if(j==0){
	  gsl_vector_float_view mdiag = gsl_matrix_float_diagonal(Dj);
	  gsl_vector_float_set_all(&mdiag.vector,1);
	}
	/* set j. super- and subdiagonal to 1 */
	else {
	  /* set superdiagonal to 1 */
	  gsl_vector_float_view superview = gsl_matrix_float_superdiagonal(Dj,j);
	  gsl_vector_float_set_all(&superview.vector,1);
	  /* set subdiagonal to 1 */
	  gsl_vector_float_view subview = gsl_matrix_float_subdiagonal(Dj,j);
	  gsl_vector_float_set_all(&subview.vector,1);
	}

	result = fmat_x_mat(R,Di,NULL);
	result2 = fmat_x_mat(result,R,NULL);
	fmat_x_mat(result2,Dj,result);
	gsl_matrix_float_free(result2);

	gsl_matrix_float_set(M,i,j,trace(result)/(1+(i>0)));
	gsl_matrix_float_free(result);

      }
    }

    gsl_matrix_float_free(Di);
    gsl_matrix_float_free(Dj);
  }

  /* calculation of the return values */
  *invM = fInv(M,NULL);
  /* pinvX has already been initialized */

  /* free memory */
  gsl_matrix_float_free(transX);
  gsl_matrix_float_free(R);
  gsl_matrix_float_free(M);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia