text:
/**
 *
 * @file qwrapper_zlacpy.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlacpy(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int m, int n, int nb,
                       const PLASMA_Complex64_t *A, int lda,
                       PLASMA_Complex64_t *B, int ldb)
{
    DAG_CORE_LACPY;
    QUARK_Insert_Task(quark, CORE_zlacpy_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &m,     VALUE,
        sizeof(int),                        &n,     VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    A,             INPUT,
        sizeof(int),                        &lda,   VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    B,             OUTPUT,
        sizeof(int),                        &ldb,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlacpy_quark = PCORE_zlacpy_quark
#define CORE_zlacpy_quark PCORE_zlacpy_quark
#endif
void CORE_zlacpy_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    const PLASMA_Complex64_t *A;
    int LDA;
    PLASMA_Complex64_t *B;
    int LDB;

    quark_unpack_args_7(quark, uplo, M, N, A, LDA, B, LDB);
    LAPACKE_zlacpy_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, A, LDA, B, LDB);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlacpy_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum uplo, int m, int n, int nb,
                          const PLASMA_Complex64_t *A, int lda,
                          PLASMA_Complex64_t *B, int ldb,
                          PLASMA_Complex64_t *fake1, int szefake1, int flag1)
{
    DAG_CORE_LACPY;
    if ( fake1 == B ) {
        QUARK_Insert_Task(quark, CORE_zlacpy_quark, task_flags,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &m,     VALUE,
            sizeof(int),                        &n,     VALUE,
            sizeof(PLASMA_Complex64_t)*nb*nb,    A,             INPUT,
            sizeof(int),                        &lda,   VALUE,
            sizeof(PLASMA_Complex64_t)*nb*nb,    B,             OUTPUT | flag1,
            sizeof(int),                        &ldb,   VALUE,
            0);
    }
    else {
        QUARK_Insert_Task(quark, CORE_zlacpy_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &m,     VALUE,
            sizeof(int),                        &n,     VALUE,
            sizeof(PLASMA_Complex64_t)*nb*nb,    A,             INPUT,
            sizeof(int),                        &lda,   VALUE,
            sizeof(PLASMA_Complex64_t)*nb*nb,    B,             OUTPUT,
            sizeof(int),                        &ldb,   VALUE,
            sizeof(PLASMA_Complex64_t)*szefake1, fake1,         flag1,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlacpy_f1_quark = PCORE_zlacpy_f1_quark
#define CORE_zlacpy_f1_quark PCORE_zlacpy_f1_quark
#endif
void CORE_zlacpy_f1_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    const PLASMA_Complex64_t *A;
    int LDA;
    PLASMA_Complex64_t *B;
    int LDB;
    void *fake1;

    quark_unpack_args_8(quark, uplo, M, N, A, LDA, B, LDB, fake1);
    LAPACKE_zlacpy_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, A, LDA, B, LDB);
}


text:
/* Copyright (c) 2011-2012, Jérémy Fix. All rights reserved. */

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

/* * Redistributions of source code must retain the above copyright notice, */
/* this list of conditions and the following disclaimer. */
/* * Redistributions in binary form must reproduce the above copyright notice, */
/* this list of conditions and the following disclaimer in the documentation */
/* and/or other materials provided with the distribution. */
/* * None of the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. */

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

#define DEBUGMODE true

#ifndef UKF_SR_STATE_NDIM_H
#define UKF_SR_STATE_NDIM_H

#include <gsl/gsl_linalg.h> // For the Cholesky decomposition
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>

#include "ukf_types.h"
#include "ukf_math.h"

void print_mat(const gsl_matrix * A)
{
    for(unsigned int i = 0 ; i < A->size1 ; ++i)
    {
        for(unsigned int j = 0 ; j < A ->size2 ; ++j)
            printf("%e ", gsl_matrix_get(A,i,j));

        printf("\n");
    }
}

void print_vec(const gsl_vector * A)
{
    for(unsigned int i = 0 ; i < A->size ; ++i)
    {
        printf("%e ", gsl_vector_get(A,i));
    }
    printf("\n");
}

namespace ukf
{

    /**
      * @short MUST NOT BE USED !!!!!!! Square root UKF for state estimation, additive noise case
      * The notations follow "Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models",p115, PhD, van Der Merwe
      */
    namespace srstate
    {
        /**
          * @short Allocation of the vectors/matrices and initialization
          *
          */
        void ukf_init(ukf_param &p, ukf_state &s)
        {
            // Parameters for the sigma points of the process equation
            p.nbSamples = 2 * p.n + 1;
            p.lambda = p.alpha * p.alpha * (p.n + p.kpa) - p.n;
            p.gamma = sqrt(p.n + p.lambda);

            // Parameters for the sigma points of the observation equation
            p.nbSamplesMeasure = 4 * p.n + 1;
            p.lambda_aug = p.alpha * p.alpha * (2*p.n + p.kpa) - 2*p.n;
            p.gamma_aug = sqrt(2*p.n + p.lambda_aug);

            // Init the matrices used to iterate
            s.xi = gsl_vector_alloc(p.n);
            gsl_vector_set_zero(s.xi);

            // Images of the sigma points
            s.xi_prediction = gsl_matrix_alloc(p.n, p.nbSamples);
            gsl_matrix_set_zero(s.xi_prediction);

            s.xi_mean = gsl_vector_alloc(p.n);
            gsl_vector_set_zero(s.xi_mean);

            // Cholesky factor of the variance/covariance matrix Pxxi
            s.Sxi = gsl_matrix_alloc(p.n,p.n);
            gsl_matrix_set_identity(s.Sxi);
            gsl_matrix_scale(s.Sxi, sqrt(p.prior_x));

            // Cholesky factor of the variance/covariance matrix of the process noise
            s.Svi = gsl_matrix_alloc(p.n,p.n);
            gsl_matrix_set_identity(s.Svi);
            gsl_matrix_scale(s.Svi, sqrt(p.process_noise));

            s.yi_prediction = gsl_matrix_alloc(p.no, p.nbSamplesMeasure);
            gsl_matrix_set_zero(s.yi_prediction);

            s.yi_mean = gsl_vector_alloc(p.no);
            gsl_vector_set_zero(s.yi_mean);

            s.ino_yi = gsl_vector_alloc(p.no);
            gsl_vector_set_zero(s.ino_yi);

            // Cholesky factor of the variance/covariance matrix Pyyi
            s.Syi = gsl_matrix_alloc(p.no, p.no);
            gsl_matrix_set_zero(s.Syi);

            s.Sni = gsl_matrix_alloc(p.no,p.no);
            gsl_matrix_set_identity(s.Sni);
            gsl_matrix_scale(s.Sni, sqrt(p.measurement_noise));

            s.Pxyi = gsl_matrix_alloc(p.n, p.no);
            gsl_matrix_set_zero(s.Pxyi);

            s.sigmaPoint = gsl_vector_alloc(p.n);
            gsl_vector_set_zero(s.sigmaPoint);

            s.sigmaPoints = gsl_matrix_alloc(p.n, p.nbSamples);
            gsl_matrix_set_zero(s.sigmaPoints);

            s.sigmaPointMeasure = gsl_vector_alloc(p.n);
            gsl_vector_set_zero(s.sigmaPoint);

            s.sigmaPointsMeasure = gsl_matrix_alloc(p.n, p.nbSamplesMeasure);
            gsl_matrix_set_zero(s.sigmaPointsMeasure);

            s.U = gsl_matrix_alloc(p.n, p.no);
            gsl_matrix_set_zero(s.U);

            // Weights used to update the statistics
            s.wm_j  = gsl_vector_alloc(p.nbSamples); // Weights used to compute the mean of the sigma points images
            s.wc_j = gsl_vector_alloc(p.nbSamples); // Weights used to update the covariance matrices

            // Set the weights
            gsl_vector_set(s.wm_j, 0, p.lambda / (p.n + p.lambda));
            gsl_vector_set(s.wc_j, 0, p.lambda / (p.n + p.lambda) + (1.0 - p.alpha*p.alpha + p.beta));
            for(int j = 1 ; j < p.nbSamples; j ++)
            {
                gsl_vector_set(s.wm_j, j, 1.0 / (2.0 * (p.n + p.lambda)));
                gsl_vector_set(s.wc_j, j, 1.0 / (2.0 * (p.n + p.lambda)));
            }

            // Set the weights
            s.wm_aug_j  = gsl_vector_alloc(p.nbSamplesMeasure); // Weights used to compute the mean of the sigma points images
            s.wc_aug_j = gsl_vector_alloc(p.nbSamplesMeasure); // Weights used to update the covariance matrices
            gsl_vector_set(s.wm_aug_j, 0, p.lambda_aug / (2*p.n + p.lambda_aug));
            gsl_vector_set(s.wc_aug_j, 0, p.lambda_aug / (2*p.n + p.lambda_aug) + (1.0 - p.alpha*p.alpha + p.beta));
            for(int j = 1 ; j < p.nbSamplesMeasure; j ++)
            {
                gsl_vector_set(s.wm_aug_j, j, 1.0 / (2.0 * (2*p.n + p.lambda_aug)));
                gsl_vector_set(s.wc_aug_j, j, 1.0 / (2.0 * (2*p.n + p.lambda_aug)));
            }

            s.Ki = gsl_matrix_alloc(p.n, p.no);
            s.Ki_T = gsl_matrix_alloc(p.no, p.n);

            // Allocate temporary matrices
            s.temp_n = gsl_vector_alloc(p.n);
            s.temp_no = gsl_vector_alloc(p.no);

            s.temp_n_1 = gsl_matrix_alloc(p.n,1);
            s.temp_1_n = gsl_matrix_alloc(1,p.n);
            s.temp_n_n = gsl_matrix_alloc(p.n, p.n);
            s.temp_n_no = gsl_matrix_alloc(p.n, p.no);
            s.temp_no_1 = gsl_matrix_alloc(p.no,1);
            s.temp_1_no = gsl_matrix_alloc(1,p.no);
            s.temp_no_no = gsl_matrix_alloc(p.no, p.no);

            s.temp_3n_n = gsl_matrix_alloc(3*p.n, p.n);
            s.temp_2nno_no = gsl_matrix_alloc(2*p.n + p.no, p.no);
        }

        /**
          * @short Free of memory allocation
          *
          */
        void ukf_free(ukf_param &p, ukf_state &s)
        {
            gsl_vector_free(s.xi);
            gsl_matrix_free(s.xi_prediction);
            gsl_vector_free(s.xi_mean);
            gsl_matrix_free(s.Sxi);
            gsl_matrix_free(s.Svi);

            gsl_matrix_free(s.yi_prediction);
            gsl_vector_free(s.yi_mean);
            gsl_vector_free(s.ino_yi);
            gsl_matrix_free(s.Syi);
            gsl_matrix_free(s.Sni);

            gsl_matrix_free(s.Pxyi);

            gsl_vector_free(s.sigmaPoint);
            gsl_matrix_free(s.sigmaPoints);

            gsl_vector_free(s.sigmaPointMeasure);
            gsl_matrix_free(s.sigmaPointsMeasure);

            gsl_matrix_free(s.U);

            gsl_vector_free(s.wm_j);
            gsl_vector_free(s.wc_j);

            gsl_vector_free(s.wm_aug_j);
            gsl_vector_free(s.wc_aug_j);

            gsl_matrix_free(s.Ki);
            gsl_matrix_free(s.Ki_T);

            gsl_vector_free(s.temp_n);
            gsl_vector_free(s.temp_no);

            gsl_matrix_free(s.temp_n_1);
            gsl_matrix_free(s.temp_1_n);
            gsl_matrix_free(s.temp_n_n);

            gsl_matrix_free(s.temp_n_no);
            gsl_matrix_free(s.temp_no_1);
            gsl_matrix_free(s.temp_1_no);
            gsl_matrix_free(s.temp_no_no);

            gsl_matrix_free(s.temp_3n_n);
            gsl_matrix_free(s.temp_2nno_no);
        }

        /**
          * @short UKF-additive (zero-mean) noise case, "Kalman Filtering and Neural Networks", p.233
          *
          */
        template <typename FunctProcess, typename FunctObservation>
                void ukf_iterate(ukf_param &p, ukf_state &s, FunctProcess f, FunctObservation h, gsl_vector* yi)
        {
            int i,j;
	    gsl_matrix_view mat_view;
	    gsl_vector_view vec_view;
            // ************************************************** //
            // ************ Compute the sigma points ************ //
            // ************************************************** //

            // 1- Compute the sigma points,
            // Equation (3.209)
            // sigmapoint_j = x_(i-1)
            // sigmapoint_j = x_(i-1) + gamma * Sxi_j for 1 <= j <= n
            // sigmapoint_j = x_(i-1) - gamma * Sxi_(j-(n+1)) for n+1 <= j <= 2n
            gsl_matrix_set_col(s.sigmaPoints, 0, s.xi);
            for(j = 1 ; j < p.n + 1 ; ++j)
                for(i = 0 ; i < p.n ; ++i)
                {
                gsl_matrix_set(s.sigmaPoints,i,j, s.xi->data[i] + p.gamma * gsl_matrix_get(s.Sxi, i, j-1));
                gsl_matrix_set(s.sigmaPoints,i,j+p.n, s.xi->data[i] - p.gamma * gsl_matrix_get(s.Sxi, i, j-1));
            }

            /**********************************/
            /***** Time update equations  *****/
            /**********************************/

            // Time update equations
            // 0 - Compute the image of the sigma points and the mean of these images
            gsl_vector_set_zero(s.xi_mean);
            for(j = 0 ; j < p.nbSamples ; ++j)
            {
                gsl_matrix_get_col(s.sigmaPoint, s.sigmaPoints, j);
                // Compute the image of the sigma points
                // Eq 3.210
		vec_view = gsl_matrix_column(s.xi_prediction,j);
                f(s.params, s.sigmaPoint, &vec_view.vector);

                // Update the mean, Eq (3.211)
                for(i = 0 ; i < p.n ; ++i)
                    s.xi_mean->data[i] += s.wm_j->data[j] * gsl_matrix_get(s.xi_prediction,i,j);
            }

            // 1- Compute the QR decomposition and keeps only the upper triangular part of R
            // Eq. 3.212
            // We first fill the matrix Sxk^-
            // with the transpose of what is in the brackets of qr { ... }
            for(j = 1 ; j < p.nbSamples; ++j)
            {
                for(i = 0 ; i < p.n ; ++i)
                {
                    gsl_matrix_set(s.temp_3n_n, j-1, i, sqrt(s.wc_j->data[1]) * (gsl_matrix_get(s.xi_prediction,i,j) - s.xi_mean->data[i]));
                }
            }
            // The last rows are filled with the cholesky factor of the process noise covariance
            for(j = 0 ; j < p.n ; ++j)
            {
	      mat_view = gsl_matrix_submatrix(s.temp_3n_n, p.nbSamples-1, 0, p.n, p.n);
                gsl_matrix_memcpy(&mat_view.matrix, s.Svi);
            }
            // We now perform the QR decomposition of s.temp_3n_n;
            gsl_linalg_QR_decomp(s.temp_3n_n, s.temp_n);
            // From this QR decomposition, we keep only the upper triangular part of R and copy it in s.Sxi
            for( j = 0 ; j < p.n ; ++j)
            {
                for(i = 0 ; i < j+1 ; ++i)
                    gsl_matrix_set(s.Sxi, i, j, gsl_matrix_get(s.temp_3n_n, i, j));

                for(i = j+1 ; i < p.n ; ++i)
                    gsl_matrix_set(s.Sxi, i, j, 0.0);
            }
            //gsl_matrix_transpose(s.Sxi);
            if(DEBUGMODE) print_mat(s.Sxi);
            if(DEBUGMODE) printf("\n");

            // 2- Perform the cholupdate, which is a rank one update
            // Eq. 3.213
            for(i = 0 ; i < p.n ; ++i)
                s.temp_n->data[i] = gsl_matrix_get(s.xi_prediction,i,0) - s.xi_mean->data[i];

            if(DEBUGMODE) printf("Before cholupdate 1 \n");
            gsl_matrix_transpose(s.Sxi);
            ukf::math::choleskyUpdate(s.Sxi, ukf::math::signof(s.wc_j->data[0]) * sqrt(fabs(s.wc_j->data[0])), s.temp_n);
            // Warning : Here, we keep Sxi as a triangular inferior matrix

            //gsl_matrix_transpose(s.Sxi);
            if(DEBUGMODE) printf("ok\n");
            if(DEBUGMODE) printf("Sxi after cholupdate : \n");
            if(DEBUGMODE) print_mat(s.Sxi);
            // ** Augment the sigma points ** //
            // Eq 3.124

            // 1 - Copy the images of the original sigma points
	    mat_view = gsl_matrix_submatrix(s.sigmaPointsMeasure,0,0, p.n, p.nbSamples);
            gsl_matrix_memcpy(&mat_view.matrix, s.xi_prediction);

            // 2 - Compute the additional sigma points
            for(j = 0 ; j < p.n ; ++j)
            {
                for(i = 0 ; i < p.n ; ++i)
                {
                    gsl_matrix_set(s.sigmaPointsMeasure, i ,j + p.nbSamples , gsl_matrix_get(s.xi_prediction, i, 0) + p.gamma_aug * gsl_matrix_get(s.Svi, i,j));
                    gsl_matrix_set(s.sigmaPointsMeasure, i ,j + p.nbSamples + p.n , gsl_matrix_get(s.xi_prediction, i, 0) - p.gamma_aug * gsl_matrix_get(s.Svi, i,j));
                }
            }

            // Compute the image of the sigma points, i.e. the associated predicted observations
            // and their mean
            gsl_vector_set_zero(s.yi_mean);
            for(j = 0 ; j < p.nbSamplesMeasure ; ++j)
            {
                gsl_matrix_get_col(s.sigmaPointMeasure, s.sigmaPointsMeasure, j);
		vec_view = gsl_matrix_column(s.yi_prediction,j);
                h(s.sigmaPointMeasure, &vec_view.vector);

                // Update the mean , eq (3.176)
                for(i = 0 ; i < p.no ; ++i)
                    s.yi_mean->data[i] += s.wm_aug_j->data[j] * gsl_matrix_get(s.yi_prediction,i,j);
            }

            /*****************************************/
            /***** Measurement update equations  *****/
            /*****************************************/

            // Compute the Cholesky factor of the variance-covariance matrix of the observations
            // Eq. 3.217
            // First put in s.temp_2nno_no the transpose of what is in the brackets of qr{...} ,
            for(j = 0 ; j < 2*p.n ; ++j)
                for(i = 0 ; i < p.no ; ++i)
                    gsl_matrix_set(s.temp_2nno_no, j, i, sqrt(s.wc_aug_j->data[1])*(gsl_matrix_get(s.yi_prediction, i, j+1) - s.yi_mean->data[i]) );

            // Copy the Cholesky factor of the covariance of the observation noise
	    mat_view = gsl_matrix_submatrix(s.temp_2nno_no,2*p.n, 0, p.no, p.no);
            gsl_matrix_memcpy(&mat_view.matrix, s.Sni);

            // We now perform the QR decomposition of s.temp_2nno_no;
            gsl_linalg_QR_decomp(s.temp_2nno_no, s.temp_no);
            // From this QR decomposition, we keep only the upper triangular part of R and copy it in s.Syi
            for( j = 0 ; j < p.no ; ++j)
            {
                for(i = 0 ; i < j+1 ; ++i)
                    gsl_matrix_set(s.Syi, i, j, gsl_matrix_get(s.temp_2nno_no, i, j));

                for(i = j+1 ; i < p.no ; ++i)
                    gsl_matrix_set(s.Syi, i, j, 0.0);
            }
            if(DEBUGMODE) printf("Syi : \n");
            if(DEBUGMODE) print_mat(s.Syi);
            if(DEBUGMODE) printf("\n");
            //gsl_matrix_transpose(s.Syi);

            // 2- Perform the cholupdate, which is a rank one update
            // Eq. 3.218
            // Syi is triangular superior
            for(i = 0 ; i < p.no ; ++i)
                s.temp_no->data[i] = gsl_matrix_get(s.yi_prediction,i,0) - s.yi_mean->data[i];
            if(DEBUGMODE) print_vec(s.temp_no);
            if(DEBUGMODE) printf("Before cholupdate 2 \n");
            gsl_matrix_transpose(s.Syi);
            ukf::math::choleskyUpdate(s.Syi, ukf::math::signof(s.wc_aug_j->data[0])*sqrt(fabs(s.wc_aug_j->data[0])), s.temp_no);
            if(DEBUGMODE) printf("ok \n");
            // Now, Syi is the lower triangular cholesky factor
            //gsl_matrix_transpose(s.Syi);
            if(DEBUGMODE) printf("Syi after cholupdate : \n");
            if(DEBUGMODE) print_mat(s.Syi);

            // Compute the state/observation covariance
            // Eq (3.219)
            gsl_matrix_set_zero(s.Pxyi);
            for(j = 0 ; j < p.nbSamplesMeasure ; ++j)
            {
                for(i = 0 ; i < p.n ; ++i)
                    s.temp_n_1->data[i] = gsl_matrix_get(s.sigmaPointsMeasure,i,j) - s.xi_mean->data[i];

                for(i = 0 ; i < p.no ; ++i)
                    s.temp_1_no->data[i] = gsl_matrix_get(s.yi_prediction,i,j) - s.yi_mean->data[i];

                gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, s.wc_aug_j->data[j] , s.temp_n_1, s.temp_1_no, 1.0, s.Pxyi);
                //gsl_matrix_add(s.Pxyi, s.temp_n_no);
            }

            /*gsl_matrix_set_zero(s.Pxyi);
            for(j = 0 ; j < p.nbSamples ; ++j)
            {
                for(i = 0 ; i < p.n ; ++i)
                    s.temp_n_1->data[i] = gsl_matrix_get(s.xi_prediction,i,j) - s.xi_mean->data[i];

                for(i = 0 ; i < p.no ; ++i)
                    s.temp_1_no->data[i] = gsl_matrix_get(s.yi_prediction,i,j) - s.yi_mean->data[i];

                gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, s.wc_j->data[j] , s.temp_n_1, s.temp_1_no, 1.0, s.Pxyi);
                //gsl_matrix_add(s.Pxyi, s.temp_n_no);
            }*/

            // Compute the Kalman gain,
            // Eq (3.220)

            // We want to solve : Kxy . Sy Sy^T = Pxy
            // Which can be solved using : Sy Sy^T Kxy^T = Pxy^T
            // We have a number of no systems to solve of the form S x = Pxy^T[:,i] ; S^T Kxy[:,i] = x with S being triangular
            // To solve this, we make use of gsl_blas_dtrsv
            for(int i = 0 ; i < p.n ; ++i)
            {
	      vec_view = gsl_matrix_row(s.Pxyi, i);
                gsl_vector_memcpy(s.temp_no, &vec_view.vector);
                gsl_blas_dtrsv(CblasLower,CblasNoTrans,CblasNonUnit,s.Syi,s.temp_no);
                gsl_blas_dtrsv(CblasLower,CblasTrans,CblasNonUnit,s.Syi,s.temp_no);
                // And then copy the result in the ith row of Ki
		vec_view = gsl_matrix_row(s.Ki, i);
                gsl_vector_memcpy(&vec_view.vector, s.temp_no);
            }
            if(DEBUGMODE) printf("Ki : \n");
            if(DEBUGMODE) print_mat(s.Ki);

            // Correction : correct the estimation of the state
            // Eq. 3.221
            // Compute the innovations
            for(i = 0 ; i < p.no ; ++i)
                s.ino_yi->data[i] = yi->data[i] - s.yi_mean->data[i];
            // And correct the current state estimate
            gsl_blas_dgemv(CblasNoTrans, 1.0, s.Ki, s.ino_yi, 1.0, s.xi);

            // Compute the matrix U = Kk Syk
            gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0 , s.Ki, s.Syi, 0, s.U);
            if(DEBUGMODE) printf("Matrix U : \n");
            if(DEBUGMODE) print_mat(s.U);

            // Perform the cholesky updates of Sxi
            // We perform p.no updates rank-1 update with the columns of U
            if(DEBUGMODE) printf("Before cholupdate 3 \n");
            for(i = 0 ; i < p.no ; ++i)
            {
	      vec_view = gsl_matrix_column(s.U, i);
                gsl_vector_memcpy(s.temp_n, &vec_view.vector);
                ukf::math::choleskyUpdate(s.Sxi, -1.0, s.temp_n);
                if(DEBUGMODE) printf("After cholupdate %i\n", i);
                if(DEBUGMODE) print_mat(s.Sxi);
            }

            //printf("After cholupdate %i\n", i);
            //print_mat(s.Sxi);
            //if(DEBUGMODE) printf("\n");

            //
            //            // Update of the process noise
            //            switch(p.process_noise_type)
            //            {
            //            case ukf::UKF_PROCESS_FIXED:
            //                //nothing to do
            //                break;
            //            case ukf::UKF_PROCESS_RLS:
            //                gsl_matrix_memcpy(s.Pvvi, s.Pxxi);
            //                gsl_matrix_scale(s.Pvvi, 1.0/p.process_noise-1.0);
            //                gsl_matrix_memcpy(s.cholPvvi, s.Pvvi);
            //                gsl_linalg_cholesky_decomp(s.cholPvvi);
            //                for(j = 0 ; j < p.n ; j++)
            //                    for(k = j+1 ; k < p.n ; k++)
            //                        gsl_matrix_set(s.cholPvvi,j,k,0.0);
            //                break;
            //            default:
            //                printf("Warning : Unrecognized process noise type\n");
            //            }

        }

        /**
          * @short Evaluation of the output from the sigma points
          *
          */
        template <typename FunctProcess, typename FunctObservation>
                void ukf_evaluate(ukf_param &p, ukf_state &s, FunctProcess f, FunctObservation h, gsl_vector* yi)
        {

            //          int i,j,k;
            //            // ************************************************** //
            //            // ************ Compute the sigma points ************ //
            //            // ************************************************** //
            //
            //            // 0 - Compute the Cholesky decomposition of s.Pxxi
            //            gsl_matrix_memcpy(s.cholPxxi, s.Pxxi);
            //            gsl_linalg_cholesky_decomp(s.cholPxxi);
            //            // Set all the elements of cholPvvi strictly above the diagonal to zero
            //            for(j = 0 ; j < p.n ; j++)
            //                for(k = j+1 ; k < p.n ; k++)
            //                    gsl_matrix_set(s.cholPxxi,j,k,0.0);
            //
            //            // 1- Compute the sigma points,
            //            // Equation (3.170)
            //            // sigmapoint_j = x_(i-1)
            //            // sigmapoint_j = x_(i-1) + gamma * sqrt(P_i-1)_j for 1 <= j <= n
            //            // sigmapoint_j = x_(i-1) - gamma * sqrt(P_i-1)_(j-(n+1)) for n+1 <= j <= 2n
            //            gsl_matrix_set_col(s.sigmaPoints, 0, s.xi);
            //            for(j = 1 ; j < p.n+1 ; ++j)
            //                for(i = 0 ; i < p.n ; ++i)
            //                {
            //                    gsl_matrix_set(s.sigmaPoints,i,j, s.xi->data[i] + p.gamma * gsl_matrix_get(s.cholPxxi, i, j-1));
            //                    gsl_matrix_set(s.sigmaPoints,i,j+p.n, s.xi->data[i] - p.gamma * gsl_matrix_get(s.cholPxxi, i, j-1));
            //                }
            //
            //            /**********************************/
            //            /***** Time update equations  *****/
            //            /**********************************/
            //
            //            // Time update equations
            //            // 0 - Compute the image of the sigma points and the mean of these images
            //            gsl_vector_set_zero(s.xi_mean);
            //            for(j = 0 ; j < p.nbSamples ; ++j)
            //            {
            //                gsl_matrix_get_col(s.sigmaPoint, s.sigmaPoints, j);
            //                f(s.params, s.sigmaPoint, &gsl_matrix_column(s.xi_prediction,j).vector);
            //
            //                // Update the mean, Eq (3.172)
            //                for(i = 0 ; i < p.n ; ++i)
            //                    s.xi_mean->data[i] += s.wm_j->data[j] * gsl_matrix_get(s.xi_prediction,i,j);
            //            }
            //
            //            // 1 - Compute the covariance of the images and add the process noise,
            //            // Equation (3.173)
            //            // Warning, s.Pxxi will now hold P_xk^-
            //            gsl_matrix_set_zero(s.Pxxi);
            //            for(j = 0 ; j < p.nbSamples ; ++j)
            //            {
            //                for(i = 0 ; i < p.n ; ++i)
            //                    s.temp_n_1->data[i] = gsl_matrix_get(s.xi_prediction,i,j) - s.xi_mean->data[i];
            //
            //                gsl_blas_dgemm(CblasNoTrans, CblasTrans, s.wc_j->data[j] , s.temp_n_1, s.temp_n_1, 0, s.temp_n_n);
            //                gsl_matrix_add(s.Pxxi, s.temp_n_n);
            //            }
            //            // Add the covariance of the evolution noise
            //            gsl_matrix_add(s.Pxxi, s.Pvvi);
            //
            //            // Augment sigma points
            //            // Equation 3.174
            //            // First put the images of the initial sigma points
            //            gsl_matrix_memcpy(&gsl_matrix_submatrix(s.sigmaPointsMeasure, 0, 0, p.n, p.nbSamples).matrix, s.xi_prediction);
            //            // And add the additional sigma points eq. (7.56)
            //            for(j = 0 ; j < p.n ; ++j)
            //            {
            //                for(i = 0 ; i < p.n ; ++i)
            //                {
            //                    gsl_matrix_set(s.sigmaPointsMeasure, i, j+p.nbSamples, gsl_matrix_get(s.xi_prediction,i,0)+p.gamma_aug*gsl_matrix_get(s.cholPvvi,i,j));
            //                    gsl_matrix_set(s.sigmaPointsMeasure, i, j+p.nbSamples+p.n, gsl_matrix_get(s.xi_prediction,i,0)-p.gamma_aug*gsl_matrix_get(s.cholPvvi,i,j));
            //                }
            //            }
            //
            //            // Compute the image of the sigma points through the observation equation
            //            // eq (3.175)
            //            gsl_vector_set_zero(yi);
            //            for(j = 0 ; j < p.nbSamplesMeasure ; ++j)
            //            {
            //                gsl_matrix_get_col(s.sigmaPointMeasure, s.sigmaPointsMeasure, j);
            //                h(s.sigmaPointMeasure, &gsl_matrix_column(s.yi_prediction,j).vector);
            //
            //                // Update the mean , eq (3.176)
            //                for(i = 0 ; i < p.no ; ++i)
            //                    yi->data[i] += s.wm_aug_j->data[j] * gsl_matrix_get(s.yi_prediction,i,j);
            //            }
        }
    } // srstate
} // ukf

#endif // UKF_SR_STATE_NDIM_H

ext: h


text:
/** GSLDAPCom.h - <title>GSLDAP: Common</title>

   Copyright (C) 2002-2003 Free Software Foundation, Inc.

   Written by:	Manuel Guesdon <mguesdon@orange-concept.com>
   Date: 	Sep 2002

   $Revision$
   $Date$
   $Id$

   This file is part of the GNUstep LDAP Library.

   <license>
   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Library General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later version.

   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Library General Public License for more details.

   You should have received a copy of the GNU Library General Public
   License along with this library; if not, write to the Free
   Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
   </license>
**/

#ifndef _GSLDAPCom_h__
#define _GSLDAPCom_h__

#include <stdio.h>
#include <string.h>
#include <ldap.h>
#include <lber.h>
#include <sys/time.h>
#include <Foundation/NSObject.h>
#include <Foundation/NSString.h>
#include <Foundation/NSArray.h>
#include <Foundation/NSData.h>
#include <Foundation/NSDictionary.h>
#include <Foundation/NSDebug.h>
#include <Foundation/NSScanner.h>
#include <Foundation/NSCalendarDate.h>
#include <Foundation/NSException.h>
#include <GNUstepBase/GSCategories.h>
#include <gsldap/GSLDAPUtils.h>
#include <gsldap/GSLDAPConnection.h>
#include <gsldap/GSLDAPEntry.h>
#include <gsldap/GSLDAPAttribute.h>
#include <gsldap/GSLDAPObjectClass.h>
#include <gsldap/GSLDAPSyntax.h>
#include <gsldap/GSLDAPMatchingRule.h>
#include <gsldap/GSLDAPSchema.h>

#endif //_GSLDAPCom_h__

max_stars_repo_name: gnustep/gsldap

max_issues_repo_name: gnustep/gsldap

max_forks_repo_name: gnustep/gsldap


text:
#ifndef DIFFEQ_SIM_H
#define DIFFEQ_SIM_H

#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>
#include <iostream>
#include <vector>

using namespace std;

class DiffEq_Sim {
    private:
        double t;      //initial time
        double h;      //time step
        double tmax;   //max time
        double hmin;

    public:
        DiffEq_Sim() {
            t    = 0.0;      //initial time
            h    = 0.1;     //time step
            tmax = 2000;
            hmin = 0.2;
        };

        ~DiffEq_Sim() {};

        int nbins;
        double* y;

        void printY() { for(int i=0; i < nbins; i++) { cout << y[i] << " ";} cout << endl; }

        vector<double> get_state() {
            vector<double> C;
            C.assign(y, y + nbins);
            return C;
        }

        double get_time() { return t; }

        virtual void initialize() {}
        virtual void derivative(const double y[], double dydt[]){}

        static int function(double t, double const y[], double dydt[], void *params) {
            DiffEq_Sim* model = static_cast <DiffEq_Sim*> (params);
            model->derivative(y, dydt);
            return GSL_SUCCESS;
        }

        int run_simulation() {
            gsl_odeiv_evolve*  e = gsl_odeiv_evolve_alloc(nbins);
            gsl_odeiv_control* c = gsl_odeiv_control_y_new(1e-20, 0);
            gsl_odeiv_step*    s = gsl_odeiv_step_alloc(gsl_odeiv_step_rkf45, nbins);
            gsl_odeiv_system sys = {function, NULL, nbins, this };
            while (t < tmax) {  //convergence check here
                int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &t, tmax, &h, y);
                if (status != GSL_SUCCESS) { return status; }
            }
            return 0;
        }

       int step_simulation( double stepsize ) {
            gsl_odeiv_evolve*  e = gsl_odeiv_evolve_alloc(nbins);
            gsl_odeiv_control* c = gsl_odeiv_control_y_new(1e-5, 0);
            gsl_odeiv_step*    s = gsl_odeiv_step_alloc(gsl_odeiv_step_rkf45, nbins);
            gsl_odeiv_system sys = {function, NULL, nbins, this };

            double tstop = t+stepsize;
            while (t < tstop) {
                int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &t, tstop, &h, y);
                if (status != GSL_SUCCESS) { return status; }
            }
            return 0;
        }

        /*
        double* advance_simulation(double I_lim) {
            while (t < 2000) {  //convergence check here
                int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &t, t1, &h, y);
                if (y[1] < I_lim) { return y; }
            }
            printY();
            return y;
        }
*/
};

#endif


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

#ifndef SCILIB_LINALG_LINSOLVE_H
#define SCILIB_LINALG_LINSOLVE_H

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

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

namespace Sci {
namespace Linalg {

// Solve linear system of equations.
template <class Layout>
inline void linsolve(Sci::Matrix_view<double, Layout> a,
                     Sci::Matrix_view<double, Layout> b)
{
    namespace stdex = std::experimental;

    static_assert(a.is_contiguous());
    static_assert(b.is_contiguous());

    assert(a.extent(0) == a.extent(1));
    assert(b.extent(0) == a.extent(1));

    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    const BLAS_INT nrhs = static_cast<BLAS_INT>(b.extent(1));
    const BLAS_INT lda = n;

    Sci::Vector<BLAS_INT, Layout> ipiv(n);

    auto matrix_layout = LAPACK_ROW_MAJOR;
    BLAS_INT ldb = nrhs;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
        ldb = n;
    }
    BLAS_INT info = LAPACKE_dgesv(matrix_layout, n, nrhs, a.data(), lda,
                                  ipiv.data(), b.data(), ldb);
    if (info != 0) {
        throw std::runtime_error("dgesv: factor U is singular");
    }
}

template <class Layout, class Allocator>
inline void linsolve(Sci::Matrix<double, Layout, Allocator>& a,
                     Sci::Matrix<double, Layout, Allocator>& b)
{
    linsolve(a.view(), b.view());
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_LINSOLVE_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_legendre.h>

#include "error.h"

/* Calculate P_m^m(x) from the analytic result:
 *   P_m^m(x) = (-1)^m (2m-1)!! (1-x^2)^(m/2) , m > 0
 *            = 1 , m = 0
 */
static double legendre_Pmm(int m, double x)
{
  if(m == 0)
  {
    return 1.0;
  }
  else
  {
    double p_mm = 1.0;
    double root_factor = sqrt(1.0-x)*sqrt(1.0+x);
    double fact_coeff = 1.0;
    int i;
    for(i=1; i<=m; i++)
    {
      p_mm *= -fact_coeff * root_factor;
      fact_coeff += 2.0;
    }
    return p_mm;
  }
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_legendre_P1_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    result->val = x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
}

int
gsl_sf_legendre_P2_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    result->val = 0.5*(3.0*x*x - 1.0);
    result->err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0);
    return GSL_SUCCESS;
  }
}

int
gsl_sf_legendre_P3_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    result->val = 0.5*x*(5.0*x*x - 3.0);
    result->err = GSL_DBL_EPSILON * (fabs(result->val) + 0.5 * fabs(x) * (fabs(5.0*x*x) + 3.0));
    return GSL_SUCCESS;
  }
}

int
gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(l < 0 || x < -1.0 || x > 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(l == 0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(l == 1) {
    result->val = x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(l == 2) {
    result->val = 0.5 * (3.0*x*x - 1.0);
    result->err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0);
    /*result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
      removed this old bogus estimate [GJ]
      */
    return GSL_SUCCESS;
  }
  else if(x == 1.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x == -1.0) {
    result->val = ( GSL_IS_ODD(l) ? -1.0 : 1.0 );
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(l < 100000) {
    /* upward recurrence: l P_l = (2l-1) z P_{l-1} - (l-1) P_{l-2} */

    double p_ellm2 = 1.0;    /* P_0(x) */
    double p_ellm1 = x;      /* P_1(x) */
    double p_ell = p_ellm1;

    double e_ellm2 = GSL_DBL_EPSILON;
    double e_ellm1 = fabs(x)*GSL_DBL_EPSILON;
    double e_ell = e_ellm1;

    int ell;

    for(ell=2; ell <= l; ell++){
      p_ell = (x*(2*ell-1)*p_ellm1 - (ell-1)*p_ellm2) / ell;
      p_ellm2 = p_ellm1;
      p_ellm1 = p_ell;

      e_ell = 0.5*(fabs(x)*(2*ell-1.0) * e_ellm1 + (ell-1.0)*e_ellm2)/ell;
      e_ellm2 = e_ellm1;
      e_ellm1 = e_ell;
    }

    result->val = p_ell;
    result->err = e_ell + l*fabs(p_ell)*GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else {
    /* Asymptotic expansion.
     * [Olver, p. 473]
     */
    double u  = l + 0.5;
    double th = acos(x);
    gsl_sf_result J0;
    gsl_sf_result Jm1;
    int stat_J0  = gsl_sf_bessel_J0_e(u*th, &J0);
    int stat_Jm1 = gsl_sf_bessel_Jn_e(-1, u*th, &Jm1);
    double pre;
    double B00;
    double c1;

    /* B00 = 1/8 (1 - th cot(th) / th^2
     * pre = sqrt(th/sin(th))
     */
    if(th < GSL_ROOT4_DBL_EPSILON) {
      B00 = (1.0 + th*th/15.0)/24.0;
      pre = 1.0 + th*th/12.0;
    }
    else {
      double sin_th = sqrt(1.0 - x*x);
      double cot_th = x / sin_th;
      B00 = 1.0/8.0 * (1.0 - th * cot_th) / (th*th);
      pre = sqrt(th/sin_th);
    }

    c1 = th/u * B00;

    result->val  = pre * (J0.val + c1 * Jm1.val);
    result->err  = pre * (J0.err + fabs(c1) * Jm1.err);
    result->err += GSL_SQRT_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_J0, stat_Jm1);
  }
}

int
gsl_sf_legendre_Pl_array(const int lmax, const double x, double * result_array)
{
  /* CHECK_POINTER(result_array) */

  if(lmax < 0 || x < -1.0 || x > 1.0) {
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(lmax == 0) {
    result_array[0] = 1.0;
    return GSL_SUCCESS;
  }
  else if(lmax == 1) {
    result_array[0] = 1.0;
    result_array[1] = x;
    return GSL_SUCCESS;
  }
  else {
    /* upward recurrence: l P_l = (2l-1) z P_{l-1} - (l-1) P_{l-2} */

    double p_ellm2 = 1.0;    /* P_0(x) */
    double p_ellm1 = x;    /* P_1(x) */
    double p_ell = p_ellm1;
    int ell;

    result_array[0] = 1.0;
    result_array[1] = x;

    for(ell=2; ell <= lmax; ell++){
      p_ell = (x*(2*ell-1)*p_ellm1 - (ell-1)*p_ellm2) / ell;
      p_ellm2 = p_ellm1;
      p_ellm1 = p_ell;
      result_array[ell] = p_ell;
    }

    return GSL_SUCCESS;
  }
}

int
gsl_sf_legendre_Pl_deriv_array(const int lmax, const double x, double * result_array, double * result_deriv_array)
{
  int stat_array = gsl_sf_legendre_Pl_array(lmax, x, result_array);

  if(lmax >= 0) result_deriv_array[0] = 0.0;
  if(lmax >= 1) result_deriv_array[1] = 1.0;

  if(stat_array == GSL_SUCCESS)
  {
    int ell;

    if(fabs(x - 1.0)*(lmax+1.0)*(lmax+1.0) <  GSL_SQRT_DBL_EPSILON)
    {
      /* x is near 1 */
      for(ell = 2; ell <= lmax; ell++)
      {
        const double pre = 0.5 * ell * (ell+1.0);
        result_deriv_array[ell] = pre * (1.0 - 0.25 * (1.0-x) * (ell+2.0)*(ell-1.0));
      }
    }
    else if(fabs(x + 1.0)*(lmax+1.0)*(lmax+1.0) <  GSL_SQRT_DBL_EPSILON)
    {
      /* x is near -1 */
      for(ell = 2; ell <= lmax; ell++)
      {
        const double sgn = ( GSL_IS_ODD(ell) ? 1.0 : -1.0 ); /* derivative is odd in x for even ell */
        const double pre = sgn * 0.5 * ell * (ell+1.0);
        result_deriv_array[ell] = pre * (1.0 - 0.25 * (1.0+x) * (ell+2.0)*(ell-1.0));
      }
    }
    else
    {
      const double diff_a = 1.0 + x;
      const double diff_b = 1.0 - x;
      for(ell = 2; ell <= lmax; ell++)
      {
        result_deriv_array[ell] = - ell * (x * result_array[ell] - result_array[ell-1]) / (diff_a * diff_b);
      }
    }

    return GSL_SUCCESS;
  }
  else
  {
    return stat_array;
  }
}

int
gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result)
{
  /* If l is large and m is large, then we have to worry
   * about overflow. Calculate an approximate exponent which
   * measures the normalization of this thing.
   */
  const double dif = l-m;
  const double sum = l+m;
  const double t_d = ( dif == 0.0 ? 0.0 : 0.5 * dif * (log(dif)-1.0) );
  const double t_s = ( dif == 0.0 ? 0.0 : 0.5 * sum * (log(sum)-1.0) );
  const double exp_check = 0.5 * log(2.0*l+1.0) + t_d - t_s;

  /* CHECK_POINTER(result) */

  if(m < 0 || l < m || x < -1.0 || x > 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
    /* Bail out. */
    OVERFLOW_ERROR(result);
  }
  else {
    /* Account for the error due to the
     * representation of 1-x.
     */
    const double err_amp = 1.0 / (GSL_DBL_EPSILON + fabs(1.0-fabs(x)));

    /* P_m^m(x) and P_{m+1}^m(x) */
    double p_mm   = legendre_Pmm(m, x);
    double p_mmp1 = x * (2*m + 1) * p_mm;

    if(l == m){
      result->val = p_mm;
      result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mm);
      return GSL_SUCCESS;
    }
    else if(l == m + 1) {
      result->val = p_mmp1;
      result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mmp1);
      return GSL_SUCCESS;
    }
    else{
      /* upward recurrence: (l-m) P(l,m) = (2l-1) z P(l-1,m) - (l+m-1) P(l-2,m)
       * start at P(m,m), P(m+1,m)
       */

      double p_ellm2 = p_mm;
      double p_ellm1 = p_mmp1;
      double p_ell = 0.0;
      int ell;

      for(ell=m+2; ell <= l; ell++){
        p_ell = (x*(2*ell-1)*p_ellm1 - (ell+m-1)*p_ellm2) / (ell-m);
        p_ellm2 = p_ellm1;
        p_ellm1 = p_ell;
      }

      result->val = p_ell;
      result->err = err_amp * (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(p_ell);

      return GSL_SUCCESS;
    }
  }
}

int
gsl_sf_legendre_Plm_array(const int lmax, const int m, const double x, double * result_array)
{
  /* If l is large and m is large, then we have to worry
   * about overflow. Calculate an approximate exponent which
   * measures the normalization of this thing.
   */
  const double dif = lmax-m;
  const double sum = lmax+m;
  const double t_d = ( dif == 0.0 ? 0.0 : 0.5 * dif * (log(dif)-1.0) );
  const double t_s = ( dif == 0.0 ? 0.0 : 0.5 * sum * (log(sum)-1.0) );
  const double exp_check = 0.5 * log(2.0*lmax+1.0) + t_d - t_s;

  /* CHECK_POINTER(result_array) */

  if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(m > 0 && (x == 1.0 || x == -1.0)) {
    int ell;
    for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
    return GSL_SUCCESS;
  }
  else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
    /* Bail out. */
    GSL_ERROR ("overflow", GSL_EOVRFLW);
  }
  else {
    double p_mm   = legendre_Pmm(m, x);
    double p_mmp1 = x * (2.0*m + 1.0) * p_mm;

    if(lmax == m){
      result_array[0] = p_mm;
      return GSL_SUCCESS;
    }
    else if(lmax == m + 1) {
      result_array[0] = p_mm;
      result_array[1] = p_mmp1;
      return GSL_SUCCESS;
    }
    else {
      double p_ellm2 = p_mm;
      double p_ellm1 = p_mmp1;
      double p_ell = 0.0;
      int ell;

      result_array[0] = p_mm;
      result_array[1] = p_mmp1;

      for(ell=m+2; ell <= lmax; ell++){
        p_ell = (x*(2.0*ell-1.0)*p_ellm1 - (ell+m-1)*p_ellm2) / (ell-m);
        p_ellm2 = p_ellm1;
        p_ellm1 = p_ell;
        result_array[ell-m] = p_ell;
      }

      return GSL_SUCCESS;
    }
  }
}

int
gsl_sf_legendre_Plm_deriv_array(
  const int lmax, const int m, const double x,
  double * result_array,
  double * result_deriv_array)
{
  if(m < 0 || m > lmax)
  {
    GSL_ERROR("m < 0 or m > lmax", GSL_EDOM);
  }
  else if(m == 0)
  {
    /* It is better to do m=0 this way, so we can more easily
     * trap the divergent case which can occur when m == 1.
     */
    return gsl_sf_legendre_Pl_deriv_array(lmax, x, result_array, result_deriv_array);
  }
  else
  {
    int stat_array = gsl_sf_legendre_Plm_array(lmax, m, x, result_array);

    if(stat_array == GSL_SUCCESS)
    {
      int ell;

      if(m == 1 && (1.0 - fabs(x) < GSL_DBL_EPSILON))
      {
        /* This divergence is real and comes from the cusp-like
         * behaviour for m = 1. For example, P[1,1] = - Sqrt[1-x^2].
         */
        GSL_ERROR("divergence near |x| = 1.0 since m = 1", GSL_EOVRFLW);
      }
      else if(m == 2 && (1.0 - fabs(x) < GSL_DBL_EPSILON))
      {
        /* m = 2 gives a finite nonzero result for |x| near 1 */
        if(fabs(x - 1.0) < GSL_DBL_EPSILON)
        {
          for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = -0.25 * x * (ell - 1.0)*ell*(ell+1.0)*(ell+2.0);
        }
        else if(fabs(x + 1.0) < GSL_DBL_EPSILON)
        {
          for(ell = m; ell <= lmax; ell++)
          {
            const double sgn = ( GSL_IS_ODD(ell) ? 1.0 : -1.0 );
            result_deriv_array[ell-m] = -0.25 * sgn * x * (ell - 1.0)*ell*(ell+1.0)*(ell+2.0);
          }
        }
        return GSL_SUCCESS;
      }
      else
      {
        /* m > 2 is easier to deal with since the endpoints always vanish */
        if(1.0 - fabs(x) < GSL_DBL_EPSILON)
        {
          for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = 0.0;
          return GSL_SUCCESS;
        }
        else
        {
          const double diff_a = 1.0 + x;
          const double diff_b = 1.0 - x;
          result_deriv_array[0] = - m * x / (diff_a * diff_b) * result_array[0];
          if(lmax-m >= 1) result_deriv_array[1] = (2.0 * m + 1.0) * (x * result_deriv_array[0] + result_array[0]);
          for(ell = m+2; ell <= lmax; ell++)
          {
            result_deriv_array[ell-m] = - (ell * x * result_array[ell-m] - (ell+m) * result_array[ell-1-m]) / (diff_a * diff_b);
          }
          return GSL_SUCCESS;
        }
      }
    }
    else
    {
      return stat_array;
    }
  }
}

int
gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(m < 0 || l < m || x < -1.0 || x > 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(m == 0) {
    gsl_sf_result P;
    int stat_P = gsl_sf_legendre_Pl_e(l, x, &P);
    double pre = sqrt((2.0*l + 1.0)/(4.0*M_PI));
    result->val  = pre * P.val;
    result->err  = pre * P.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_P;
  }
  else if(x == 1.0 || x == -1.0) {
    /* m > 0 here */
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* m > 0 and |x| < 1 here */

    /* Starting value for recursion.
     * Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) (-1)^m (1-x^2)^(m/2) / pi^(1/4)
     */
    gsl_sf_result lncirc;
    gsl_sf_result lnpoch;
    double lnpre_val;
    double lnpre_err;
    gsl_sf_result ex_pre;
    double sr;
    const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
    const double y_mmp1_factor = x * sqrt(2.0*m + 3.0);
    double y_mm, y_mm_err;
    double y_mmp1, y_mmp1_err;
    gsl_sf_log_1plusx_e(-x*x, &lncirc);
    gsl_sf_lnpoch_e(m, 0.5, &lnpoch);  /* Gamma(m+1/2)/Gamma(m) */
    lnpre_val = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
    lnpre_err = 0.25*M_LNPI*GSL_DBL_EPSILON + 0.5 * (lnpoch.err + fabs(m)*lncirc.err);
    /* Compute exp(ln_pre) with error term, avoiding call to gsl_sf_exp_err BJG */
    ex_pre.val = exp(lnpre_val);
    ex_pre.err = 2.0*(sinh(lnpre_err) + GSL_DBL_EPSILON)*ex_pre.val;
    sr     = sqrt((2.0+1.0/m)/(4.0*M_PI));
    y_mm   = sgn * sr * ex_pre.val;
    y_mm_err  = 2.0 * GSL_DBL_EPSILON * fabs(y_mm) + sr * ex_pre.err;
    y_mm_err *= 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-x));
    y_mmp1 = y_mmp1_factor * y_mm;
    y_mmp1_err=fabs(y_mmp1_factor) * y_mm_err;

    if(l == m){
      result->val  = y_mm;
      result->err  = y_mm_err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mm);
      return GSL_SUCCESS;
    }
    else if(l == m + 1) {
      result->val  = y_mmp1;
      result->err  = y_mmp1_err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mmp1);
      return GSL_SUCCESS;
    }
    else{
      double y_ell = 0.0;
      double y_ell_err;
      int ell;

      /* Compute Y_l^m, l > m+1, upward recursion on l. */
      for(ell=m+2; ell <= l; ell++){
        const double rat1 = (double)(ell-m)/(double)(ell+m);
        const double rat2 = (ell-m-1.0)/(ell+m-1.0);
        const double factor1 = sqrt(rat1*(2.0*ell+1.0)*(2.0*ell-1.0));
        const double factor2 = sqrt(rat1*rat2*(2.0*ell+1.0)/(2.0*ell-3.0));
        y_ell = (x*y_mmp1*factor1 - (ell+m-1.0)*y_mm*factor2) / (ell-m);
        y_mm   = y_mmp1;
        y_mmp1 = y_ell;

        y_ell_err = 0.5*(fabs(x*factor1)*y_mmp1_err + fabs((ell+m-1.0)*factor2)*y_mm_err) / fabs(ell-m);
        y_mm_err = y_mmp1_err;
        y_mmp1_err = y_ell_err;
      }

      result->val  = y_ell;
      result->err  = y_ell_err + (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(y_ell);

      return GSL_SUCCESS;
    }
  }
}

int
gsl_sf_legendre_sphPlm_array(const int lmax, int m, const double x, double * result_array)
{
  /* CHECK_POINTER(result_array) */

  if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
    GSL_ERROR ("error", GSL_EDOM);
  }
  else if(m > 0 && (x == 1.0 || x == -1.0)) {
    int ell;
    for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
    return GSL_SUCCESS;
  }
  else {
    double y_mm;
    double y_mmp1;

    if(m == 0) {
      y_mm   = 0.5/M_SQRTPI;          /* Y00 = 1/sqrt(4pi) */
      y_mmp1 = x * M_SQRT3 * y_mm;
    }
    else {
      /* |x| < 1 here */

      gsl_sf_result lncirc;
      gsl_sf_result lnpoch;
      double lnpre;
      const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
      gsl_sf_log_1plusx_e(-x*x, &lncirc);
      gsl_sf_lnpoch_e(m, 0.5, &lnpoch);  /* Gamma(m+1/2)/Gamma(m) */
      lnpre = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
      y_mm   = sqrt((2.0+1.0/m)/(4.0*M_PI)) * sgn * exp(lnpre);
      y_mmp1 = x * sqrt(2.0*m + 3.0) * y_mm;
    }

    if(lmax == m){
      result_array[0] = y_mm;
      return GSL_SUCCESS;
    }
    else if(lmax == m + 1) {
      result_array[0] = y_mm;
      result_array[1] = y_mmp1;
      return GSL_SUCCESS;
    }
    else{
      double y_ell;
      int ell;

      result_array[0] = y_mm;
      result_array[1] = y_mmp1;

      /* Compute Y_l^m, l > m+1, upward recursion on l. */
      for(ell=m+2; ell <= lmax; ell++){
        const double rat1 = (double)(ell-m)/(double)(ell+m);
        const double rat2 = (ell-m-1.0)/(ell+m-1.0);
        const double factor1 = sqrt(rat1*(2*ell+1)*(2*ell-1));
        const double factor2 = sqrt(rat1*rat2*(2*ell+1)/(2*ell-3));
        y_ell = (x*y_mmp1*factor1 - (ell+m-1)*y_mm*factor2) / (ell-m);
        y_mm   = y_mmp1;
        y_mmp1 = y_ell;
        result_array[ell-m] = y_ell;
      }
    }

    return GSL_SUCCESS;
  }
}

int
gsl_sf_legendre_sphPlm_deriv_array(
  const int lmax, const int m, const double x,
  double * result_array,
  double * result_deriv_array)
{
  if(m < 0 || lmax < m || x < -1.0 || x > 1.0)
  {
    GSL_ERROR ("domain", GSL_EDOM);
  }
  else if(m == 0)
  {
    /* m = 0 is easy to trap */
    const int stat_array = gsl_sf_legendre_Pl_deriv_array(lmax, x, result_array, result_deriv_array);
    int ell;
    for(ell = 0; ell <= lmax; ell++)
    {
      const double prefactor = sqrt((2.0 * ell + 1.0)/(4.0*M_PI));
      result_array[ell] *= prefactor;
      result_deriv_array[ell] *= prefactor;
    }
    return stat_array;
  }
  else if(m == 1)
  {
    /* Trapping m = 1 is necessary because of the possible divergence.
     * Recall that this divergence is handled properly in ..._Plm_deriv_array(),
     * and the scaling factor is not large for small m, so we just scale.
     */
    const int stat_array = gsl_sf_legendre_Plm_deriv_array(lmax, m, x, result_array, result_deriv_array);
    int ell;
    for(ell = 1; ell <= lmax; ell++)
    {
      const double prefactor = sqrt((2.0 * ell + 1.0)/(ell + 1.0) / (4.0*M_PI*ell));
      result_array[ell-1] *= prefactor;
      result_deriv_array[ell-1] *= prefactor;
    }
    return stat_array;
  }
  else
  {
    /* as for the derivative of P_lm, everything is regular for m >= 2 */

    int stat_array = gsl_sf_legendre_sphPlm_array(lmax, m, x, result_array);

    if(stat_array == GSL_SUCCESS)
    {
      int ell;

      if(1.0 - fabs(x) < GSL_DBL_EPSILON)
      {
        for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = 0.0;
        return GSL_SUCCESS;
      }
      else
      {
        const double diff_a = 1.0 + x;
        const double diff_b = 1.0 - x;
        result_deriv_array[0] = - m * x / (diff_a * diff_b) * result_array[0];
        if(lmax-m >= 1) result_deriv_array[1] = sqrt(2.0 * m + 3.0) * (x * result_deriv_array[0] + result_array[0]);
        for(ell = m+2; ell <= lmax; ell++)
        {
          const double c1 = sqrt(((2.0*ell+1.0)/(2.0*ell-1.0)) * ((double)(ell-m)/(double)(ell+m)));
          result_deriv_array[ell-m] = - (ell * x * result_array[ell-m] - c1 * (ell+m) * result_array[ell-1-m]) / (diff_a * diff_b);
        }
        return GSL_SUCCESS;
      }
    }
    else
    {
      return stat_array;
    }
  }
}

int
gsl_sf_legendre_array_size(const int lmax, const int m)
{
  return lmax-m+1;
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_legendre_P1(const double x)
{
  EVAL_RESULT(gsl_sf_legendre_P1_e(x, &result));
}

double gsl_sf_legendre_P2(const double x)
{
  EVAL_RESULT(gsl_sf_legendre_P2_e(x, &result));
}

double gsl_sf_legendre_P3(const double x)
{
  EVAL_RESULT(gsl_sf_legendre_P3_e(x, &result));
}

double gsl_sf_legendre_Pl(const int l, const double x)
{
  EVAL_RESULT(gsl_sf_legendre_Pl_e(l, x, &result));
}

double gsl_sf_legendre_Plm(const int l, const int m, const double x)
{
  EVAL_RESULT(gsl_sf_legendre_Plm_e(l, m, x, &result));
}

double gsl_sf_legendre_sphPlm(const int l, const int m, const double x)
{
  EVAL_RESULT(gsl_sf_legendre_sphPlm_e(l, m, x, &result));
}

ext: c

max_stars_repo_name: ruslankuzmin/julia


text:
/*
 * Tema 2 ASC
 * 2020 Spring
 */
#include <cblas.h>
#include "utils.h"

#define IDX2L(i, j, N) (((i) * (N)) + (j))

/*
 * Add your BLAS implementation here
 */

double* my_solver(int N, double *A, double *B) {
	int i, j;
	double alpha = 1.0;
	/* Calculate the final result in C matrix */
	double *C = calloc(N * N, sizeof(double));
	if (C == NULL) {
		return NULL;
	}
	/* Calculate the B*At part in D */
	double *D = calloc(N * N, sizeof(double));
	if (D == NULL) {
		return NULL;
	}
	/* Calculate the A * A in A_square */
	double *A_square = calloc(N * N, sizeof(double));
	if (A_square == NULL) {
		return NULL;
	}

	/* Copy B matrix in D matrix */
	for (i = 0; i != N; i++) {
		for (j = 0; j != N; j++) {
			D[IDX2L(i, j, N)] = B[IDX2L(i, j, N)];
		}
	}
	/* Copy A matrix in A_square matrix */
	for (i = 0; i != N; i++) {
		for (j = 0; j != N; j++) {
			if ( i <= j) {
				A_square[IDX2L(i, j, N)] = A[IDX2L(i, j, N)];
			}
		}
	}
	/*
	 * Multiply matrixes B * At, where At is A transposed
	 * The result is saved in D matrix
	 */
	cblas_dtrmm(CblasRowMajor, CblasRight, CblasUpper, CblasTrans,
				CblasNonUnit, N, N, alpha, A, N, D, N);
	/*
	 * Multiply matrixes A * A, saved result in A_square
	 * Matrix A_square is a triangular matrix
	 */
	cblas_dtrmm(CblasRowMajor, CblasRight, CblasUpper, CblasNoTrans,
			    CblasNonUnit, N, N, alpha, A, N, A_square, N);

	/*
	 * Multiply matrixes A_square * B, where A_square = A * A
	 * The result is saved in B matrix
	 */
	cblas_dtrmm(CblasRowMajor, CblasLeft, CblasUpper, CblasNoTrans,
				CblasNonUnit, N, N, alpha, A_square, N, B, N);

	/* Add matrixes D and B, save the result in C matrix */
	for (i = 0; i != N; i++) {
		for (j = 0; j != N; j++) {
			C[IDX2L(i, j, N)] = D[IDX2L(i, j, N)] + B[IDX2L(i, j, N)];
		}
	}
	/* Free memory from D and A_square */
	free(D);
	free(A_square);
	printf("BLAS SOLVER\n");
	return C;
}

max_stars_repo_name: vladutmargineanu/Code-Optimization

max_issues_repo_name: vladutmargineanu/Code-Optimization

max_forks_repo_name: vladutmargineanu/Code-Optimization


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

/* Based on rk4imp.c by Gerard Jungman */

/* Gaussian implicit 4th order Runge-Kutta method. Error estimation
   is carried out by the step doubling method.
 */

/* References:

   Ascher, U.M., Petzold, L.R., Computer methods for ordinary
   differential and differential-algebraic equations, SIAM,
   Philadelphia, 1998. ISBN 0898714125, 9780898714128
*/

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

#include "odeiv_util.h"
#include "rksubs.c"
#include "modnewton1.c"

/* Stage of method */
#define RK4IMP_STAGE 2

typedef struct
{
  gsl_matrix *A;                /* Runge-Kutta coefficients */
  double *y_onestep;            /* Result with one step */
  double *y_twostep;            /* Result with two half steps */
  double *ytmp;                 /* Temporary work space */
  double *y_save;               /* Backup space */
  double *YZ;                   /* Runge-Kutta points */
  double *fYZ;                  /* Derivatives at YZ */
  gsl_matrix *dfdy;             /* Jacobian matrix */
  double *dfdt;                 /* time derivative of f */
  modnewton1_state_t *esol;     /* nonlinear equation solver */
  double *errlev;               /* desired error level of y */
  const gsl_odeiv2_driver *driver;      /* pointer to driver object */
}
rk4imp_state_t;

static void *
rk4imp_alloc (size_t dim)
{
  rk4imp_state_t *state = (rk4imp_state_t *) malloc (sizeof (rk4imp_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk4imp_state",
                      GSL_ENOMEM);
    }

  state->A = gsl_matrix_alloc (RK4IMP_STAGE, RK4IMP_STAGE);

  if (state->A == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for A", GSL_ENOMEM);
    }

  state->y_onestep = (double *) malloc (dim * sizeof (double));

  if (state->y_onestep == 0)
    {
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  state->y_twostep = (double *) malloc (dim * sizeof (double));

  if (state->y_twostep == 0)
    {
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  state->y_save = (double *) malloc (dim * sizeof (double));

  if (state->y_save == 0)
    {
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_save", GSL_ENOMEM);
    }

  state->YZ = (double *) malloc (dim * RK4IMP_STAGE * sizeof (double));

  if (state->YZ == 0)
    {
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for YZ", GSL_ENOMEM);
    }

  state->fYZ = (double *) malloc (dim * RK4IMP_STAGE * sizeof (double));

  if (state->fYZ == 0)
    {
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for fYZ", GSL_ENOMEM);
    }

  state->dfdt = (double *) malloc (dim * sizeof (double));

  if (state->dfdt == 0)
    {
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdt", GSL_ENOMEM);
    }

  state->dfdy = gsl_matrix_alloc (dim, dim);

  if (state->dfdy == 0)
    {
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdy", GSL_ENOMEM);
    }

  state->esol = modnewton1_alloc (dim, RK4IMP_STAGE);

  if (state->esol == 0)
    {
      gsl_matrix_free (state->dfdy);
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for esol", GSL_ENOMEM);
    }

  state->errlev = (double *) malloc (dim * sizeof (double));

  if (state->errlev == 0)
    {
      modnewton1_free (state->esol);
      gsl_matrix_free (state->dfdy);
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for errlev", GSL_ENOMEM);
    }

  state->driver = NULL;

  return state;
}

static int
rk4imp_apply (void *vstate, size_t dim, double t, double h,
              double y[], double yerr[],
              const double dydt_in[], double dydt_out[],
              const gsl_odeiv2_system * sys)
{
  /* Makes a Gaussian implicit 4th order Runge-Kutta with step size h
     and estimates the local error of the step.
   */

  rk4imp_state_t *state = (rk4imp_state_t *) vstate;

  double *const y_onestep = state->y_onestep;
  double *const y_twostep = state->y_twostep;
  double *const ytmp = state->ytmp;
  double *const y_save = state->y_save;
  double *const YZ = state->YZ; /* Runge-Kutta points */
  double *const fYZ = state->fYZ;
  gsl_matrix *const dfdy = state->dfdy;
  double *const dfdt = state->dfdt;
  double *const errlev = state->errlev;

  const modnewton1_state_t *esol = state->esol;

  /* Runge-Kutta coefficients */

  gsl_matrix *A = state->A;

  const double b[] = { 0.5, 0.5 };
  const double c[] = { (3 - sqrt (3)) / 6, (3 + sqrt (3)) / 6 };

  gsl_matrix_set (A, 0, 0, 1.0 / 4);
  gsl_matrix_set (A, 0, 1, (3 - 2 * sqrt (3)) / 12);
  gsl_matrix_set (A, 1, 0, (3 + 2 * sqrt (3)) / 12);
  gsl_matrix_set (A, 1, 1, 1.0 / 4);

  if (esol == NULL)
    {
      GSL_ERROR ("no non-linear equation solver speficied", GSL_EINVAL);
    }

  /* Get desired error levels via gsl_odeiv2_control object through driver
     object, which is a requirement for this stepper.
   */

  if (state->driver == NULL)
    {
      return GSL_EFAULT;
    }
  else
    {
      size_t i;

      for (i = 0; i < dim; i++)
        {
          if (dydt_in != NULL)
            {
              gsl_odeiv2_control_errlevel (state->driver->c, y[i],
                                           dydt_in[i], h, i, &errlev[i]);
            }
          else
            {
              gsl_odeiv2_control_errlevel (state->driver->c, y[i],
                                           0.0, h, i, &errlev[i]);
            }
        }
    }

  /* Evaluate Jacobian for modnewton1 */

  {
    int s = GSL_ODEIV_JA_EVAL (sys, t, y, dfdy->data, dfdt);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Calculate a single step with size h */

  {
    int s = modnewton1_init ((void *) esol, A, h, dfdy, sys);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = modnewton1_solve ((void *) esol, A, c, t, h, y,
                              sys, YZ, errlev);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    size_t j;

    for (j = 0; j < RK4IMP_STAGE; j++)
      {
        int s =
          GSL_ODEIV_FN_EVAL (sys, t + c[j] * h, &YZ[j * dim], &fYZ[j * dim]);

        if (s != GSL_SUCCESS)
          {
            return s;
          }
      }
  }

  {
    int s = rksubs (y_onestep, h, y, fYZ, b, RK4IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Error estimation by step doubling */

  {
    int s = modnewton1_init ((void *) esol, A, h / 2.0, dfdy, sys);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* 1st half step */

  {
    int s = modnewton1_solve ((void *) esol, A, c, t, h / 2.0, y,
                              sys, YZ, errlev);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    size_t j;

    for (j = 0; j < RK4IMP_STAGE; j++)
      {
        int s = GSL_ODEIV_FN_EVAL (sys, t + c[j] * h / 2.0, &YZ[j * dim],
                                   &fYZ[j * dim]);
        if (s != GSL_SUCCESS)
          {
            return s;
          }
      }
  }

  {
    int s = rksubs (ytmp, h / 2.0, y, fYZ, b, RK4IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      return s;
  }

  /* Save original y values in case of error */

  DBL_MEMCPY (y_save, y, dim);

  /* 2nd half step */

  {
    int s = modnewton1_solve ((void *) esol, A, c, t + h / 2.0, h / 2.0,
                              ytmp, sys, YZ, errlev);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    size_t j;

    for (j = 0; j < RK4IMP_STAGE; j++)
      {
        int s =
          GSL_ODEIV_FN_EVAL (sys, t + h / 2.0 + c[j] * h / 2.0, &YZ[j * dim],
                             &fYZ[j * dim]);
        if (s != GSL_SUCCESS)
          {
            return s;
          }
      }
  }

  {
    int s = rksubs (y_twostep, h / 2.0, ytmp, fYZ, b, RK4IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      {
        DBL_MEMCPY (y, y_save, dim);
        return s;
      }
  }

  /* Note: rk4imp returns y using the results from two half steps
     instead of the single step since the results are freely
     available and more precise.
   */

  DBL_MEMCPY (y, y_twostep, dim);

  /* Error estimation */

  {
    size_t i;
    for (i = 0; i < dim; i++)
      {
        yerr[i] = ODEIV_ERR_SAFETY * 0.5 *
          fabs (y_twostep[i] - y_onestep[i]) / 15.0;
      }
  }

  /* Derivatives at output */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
        {
          /* Restore original values */
          DBL_MEMCPY (y, y_save, dim);

          return s;
        }
    }

  return GSL_SUCCESS;
}

static int
rk4imp_set_driver (void *vstate, const gsl_odeiv2_driver * d)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;

  state->driver = d;

  return GSL_SUCCESS;
}

static int
rk4imp_reset (void *vstate, size_t dim)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;

  DBL_ZERO_MEMSET (state->y_onestep, dim);
  DBL_ZERO_MEMSET (state->y_twostep, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);
  DBL_ZERO_MEMSET (state->y_save, dim);
  DBL_ZERO_MEMSET (state->YZ, dim * RK4IMP_STAGE);
  DBL_ZERO_MEMSET (state->fYZ, dim * RK4IMP_STAGE);

  return GSL_SUCCESS;
}

static unsigned int
rk4imp_order (void *vstate)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;
  state = 0;                    /* prevent warnings about unused parameters */
  return 4;
}

static void
rk4imp_free (void *vstate)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;

  free (state->errlev);
  modnewton1_free (state->esol);
  gsl_matrix_free (state->dfdy);
  free (state->dfdt);
  free (state->fYZ);
  free (state->YZ);
  free (state->y_save);
  free (state->ytmp);
  free (state->y_twostep);
  free (state->y_onestep);
  gsl_matrix_free (state->A);
  free (state);
}

static const gsl_odeiv2_step_type rk4imp_type = {
  "rk4imp",                     /* name */
  1,                            /* can use dydt_in? */
  1,                            /* gives exact dydt_out? */
  &rk4imp_alloc,
  &rk4imp_apply,
  &rk4imp_set_driver,
  &rk4imp_reset,
  &rk4imp_order,
  &rk4imp_free
};

const gsl_odeiv2_step_type *gsl_odeiv2_step_rk4imp = &rk4imp_type;

max_stars_repo_name: WikiGaze/Wikipedia-readers-gaze

max_issues_repo_name: WikiGaze/Wikipedia-readers-gaze

max_forks_repo_name: WikiGaze/Wikipedia-readers-gaze


text:
#ifndef PHI_SVM
#define PHI_SVM
#include <float.h>
#include <vector>
#include <list>
#include <fstream>
#include <iostream>
#include <cmath>
#include <cassert>
#include <cstring>
#include <cstdlib>
#include <sys/time.h>

// BLAS dependency
#ifdef USE_MKL
#include <mkl.h>
#else
typedef int MKL_INT;
#if defined __APPLE__
#include <Accelerate/Accelerate.h>
#else
extern "C" {
#include <cblas.h>
}
#endif
#endif

#include <string>
#include <boost/serialization/serialization.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/back_inserter.hpp>

#include <sstream>

using namespace std;

#define BLOCKSIZE 256
#define REDUCE0 0x00000001
#define REDUCE1 0x00000002
#define MAX_PITCH 262144
#define MAX_POINTS (MAX_PITCH / sizeof(float) - 2)
#define MAXITERATIONS 100000
#define intDivideRoundUp(a, b) (a % b != 0) ? (a / b + 1) : (a / b)

enum KernelType {
  NEWLINEAR,
  NEWGAUSSIAN,
  NEWPRECOMPUTED
};

struct Kernel_params {
  float gamma;
  float coef0;
  int degree;
  float b;
  std::string kernel_type;
};

enum SelectionHeuristic {
  FIRSTORDER,
  SECONDORDER,
  RANDOM,
  ADAPTIVE
};

class PhiSVMModel {
public:
  int nSamples;
  int nDimension;
  float epsilon;
  float bLow;
  float bHigh;
  float* alpha;
  float* f;
  float* kernelDiag;
  float* data;
  float* labels;
  PhiSVMModel() { // default constructor for training from the scratch
  }
  PhiSVMModel(int maxPoints, int maxDimension) {  // constructor for getting an existing model
    alpha = new float[maxPoints];
    f = new float[maxPoints];
    kernelDiag = new float[maxPoints];
    data = new float[maxPoints*maxDimension];
    labels = new float[maxPoints];
  }
  ~PhiSVMModel() {  // leave the data and labels to others to delete
    delete alpha;
    delete f;
    delete kernelDiag;
  }
};

class NewCache {
 public:
  NewCache(int nPointsIn, int cacheSizeIn);
  ~NewCache();
  void findData(const int index, int& offset, bool& compute);
  void search(const int index, int& offset, bool& compute);
  void printCache();
  void printStatistics();

 private:
  int nPoints;
  int cacheSize;
  class DirectoryEntry {
   public:
    enum {
      NEVER,
      EVICTED,
      INCACHE
    };
    DirectoryEntry();
    int status;
    int location;
    std::list<int>::iterator lruListEntry;
  };

  std::vector<DirectoryEntry> directory;
  std::list<int> lruList;
  int occupancy;
  int hits;
  int compulsoryMisses;
  int capacityMisses;
};

class Controller {
 public:
  Controller(float initialGap, SelectionHeuristic currentMethodIn,
             int samplingIntervalIn, int problemSize);
  void addIteration(float gap);
  void print();
  SelectionHeuristic getMethod();

 private:
  bool adaptive;
  int samplingInterval;
  std::vector<float> progress;
  std::vector<int> method;
  SelectionHeuristic currentMethod;
  std::vector<float> rates;
  int timeSinceInspection;
  int inspectionPeriod;
  int beginningOfEpoch;
  int middleOfEpoch;
  int currentInspectionPhase;
  float filter(int begin, int end);
  float findRate(struct timeval* start, struct timeval* finish, int beginning,
                 int end);
  struct timeval start;
  struct timeval mid;
  struct timeval finish;
};

bool hasTwoTypes(float* labels, int n);

float crossValidation(float* data, int nPoints, int nDimension,
                               int nFolds, float* labels, Kernel_params* kp,
                               float cost, SelectionHeuristic heuristicMethod,
                               float epsilon, float tolerance,
                               float* transposedData, bool shuffle);

float incrementalCrossValidation(PhiSVMModel** models, float* data, int nPoints, int nDimension,
                               int nFolds, float* labels, Kernel_params* kp,
                               float cost, SelectionHeuristic heuristicMethod,
                               float epsilon, float tolerance,
                               float* transposedData, int maxNPoints);

void svmGroupClasses(int nPoints, float* labels, int** start_ret,
                     int** count_ret, int* perm);
PhiSVMModel* performTraining(float* data, int nPoints, int nDimension, float* labels,
                     Kernel_params* kp, float cost,
                     SelectionHeuristic heuristicMethod, float epsilon,
                     float tolerance, float* transposedData, PhiSVMModel* curModel);
PhiSVMModel* performOnlineTraining(int nOldPoints, Kernel_params* kp, float cost,
                     SelectionHeuristic heuristicMethod, float epsilon,
                     float tolerance, PhiSVMModel* curModel);

template <int Kernel>
void firstOrder(float* devData, int devDataPitchInFloats,
                float* devTransposedData, int devTransposedDataPitchInFloats,
                float* devLabels, int nPoints, int nDimension, float epsilon,
                float cEpsilon, float* devAlpha, float* devF, float alpha1Diff,
                float alpha2Diff, int iLow, int iHigh, float parameterA,
                float parameterB, float parameterC, float* devCache,
                int devCachePitchInFloats, int iLowCacheIndex,
                int iHighCacheIndex, float* devKernelDiag, void* devResult,
                float cost, bool iLowCompute, bool iHighCompute, int nthreads);
void launchFirstOrder(bool iLowCompute, bool iHighCompute, int kType,
                      int nPoints, int nDimension, float* devData,
                      int devDataPitchInFloats, float* devTransposedData,
                      int devTransposedDataPitchInFloats, float* devLabels,
                      float epsilon, float cEpsilon, float* devAlpha,
                      float* devF, float sAlpha1Diff, float sAlpha2Diff,
                      int iLow, int iHigh, float parameterA, float parameterB,
                      float parameterC, float* devCache,
                      int devCachePitchInFloats, int iLowCacheIndex,
                      int iHighCacheIndex, float* devKernelDiag,
                      void* devResult, float cost, int nthreads);
template <int Kernel>
void secondOrder(float* devData, int devDataPitchInFloats,
                 float* devTransposedData, int devTransposedDataPitchInFloats,
                 float* devLabels, int nPoints, int nDimension, float epsilon,
                 float cEpsilon, float* devAlpha, float* devF, float alpha1Diff,
                 float alpha2Diff, int iLow, int iHigh, float parameterA,
                 float parameterB, float parameterC, float* devCache,
                 int devCachePitchInFloats, int iLowCacheIndex,
                 int iHighCacheIndex, float* devKernelDiag, void* devResult,
                 float cost, float bHigh, bool iHighCompute, int nthreads);
void launchSecondOrder(bool iLowCompute, bool iHighCompute, int kType,
                       int nPoints, int nDimension, float* devData,
                       int devDataPitchInFloats, float* devTransposedData,
                       int devTransposedDataPitchInFloats, float* devLabels,
                       float epsilon, float cEpsilon, float* devAlpha,
                       float* devF, float sAlpha1Diff, float sAlpha2Diff,
                       int iLow, int iHigh, float parameterA, float parameterB,
                       float parameterC, NewCache* kernelCache, float* devCache,
                       int devCachePitchInFloats, int iLowCacheIndex,
                       int iHighCacheIndex, float* devKernelDiag,
                       void* devResult, float cost, int iteration,
                       int nthreads);
template <int Kernel>
void initializeArrays(float* devData, int devDataPitchInFloats, float* devCache,
                      int devCachePitchInFloats, int nPoints, int nDimension,
                      float parameterA, float parameterB, float parameterC,
                      float* devKernelDiag, float* devAlpha, float* devF,
                      float* devLabels, int nthreads);
void launchInitialization(float* devData, int devDataPitchInFloats,
                          float* devCache, int devCachePitchInFloats,
                          int nPoints, int nDimension, int kType,
                          float parameterA, float parameterB, float parameterC,
                          float* devKernelDiag, float* devAlpha, float* devF,
                          float* devLabels, int nthreads);
template <int Kernel>
void takeFirstStep(void* devResult, float* devKernelDiag, float* devData,
                   int devDataPitchInFloats, float* devCache,
                   int devCachePitchInFloats, float* devAlpha, float cost,
                   int nDimension, int iLow, int iHigh, float parameterA,
                   float parameterB, float parameterC);
void launchTakeFirstStep(void* devResult, float* devKernelDiag, float* devData,
                         int devDataPitchInFloats, float* devCache,
                         int devCachePitchInFloats, float* devAlpha, float cost,
                         int nDimension, int iLow, int iHigh, int kType,
                         float parameterA, float parameterB, float parameterC,
                         int nthreads);
void performClassification(float* data, int nData, int nDimension,
                           Kernel_params* kp, float** p_result, PhiSVMModel* model);
void computeKernels(float* devNorms, int devNormsPitchInFloats, float* devAlpha,
                    int nPoints, int nSV, const KernelType kType, int degree,
                    float b, float* devResult);
float kernel(const float v, const int degree, const KernelType kType);
void makeSelfDots(float* devSource, int devSourcePitchInFloats, float* devDest,
                  int sourceCount, int sourceLength);
void makeDots(float* devDots, int devDotsPitchInFloats, float* devSVDots,
              float* devDataDots, int nSV, int nPoints);

class DumpModel {
private:
  friend class boost::serialization::access;
  template<class Archive>
  void serialize(Archive & ar, const unsigned int version) {
    ar & nSamples;
    ar & nDimension;
    for (int i=0; i<nSamples*nDimension; i++) {
      ar & trainingData[i];
    }
    ar & phiSVMModel->nSamples;  // redundant with Model::nSamples
    ar & phiSVMModel->nDimension;  // for precomputed kernel, should be the same as nSamples
    ar & phiSVMModel->epsilon;
    ar & phiSVMModel->bLow;
    ar & phiSVMModel->bHigh;
    for (int i=0; i<phiSVMModel->nSamples; i++) {
      ar & phiSVMModel->alpha[i];
    }
    for (int i=0; i<phiSVMModel->nSamples; i++) {
      ar & phiSVMModel->f[i];
    }
    for (int i=0; i<phiSVMModel->nSamples; i++) {
      ar & phiSVMModel->kernelDiag[i];
    }
    // the data here for precomputed kernel is the kernel matrix
    for (int i=0; i<phiSVMModel->nSamples*phiSVMModel->nDimension; i++) {
      ar & phiSVMModel->data[i];
    }
    for (int i=0; i<phiSVMModel->nSamples; i++) {
      ar & phiSVMModel->labels[i];
    }
  }
public:
  int nSamples;
  int nDimension;
  float* trainingData;
  PhiSVMModel* phiSVMModel;
  DumpModel() {}
  DumpModel(int s, int d) {
    nSamples = s;
    nDimension = d;
    trainingData = new float[s*d];
    phiSVMModel = new PhiSVMModel(MAX_POINTS, MAX_POINTS);
  }
  ~DumpModel() {
    //delete trainingData;
    //delete phiSVMModel->data;
    //delete phiSVMModel->labels;
    //delete phiSVMModel;
  }
};

std::string serialize_DumpModel(DumpModel* model_ptr);
DumpModel* deserialize_DumpModel(std::string, int, int);
void DumpModelToDisk(std::string modelStr);
#endif

max_stars_repo_name: PrincetonUniversity/fcma-toolbox

max_issues_repo_name: PrincetonUniversity/fcma-toolbox

max_forks_repo_name: PrincetonUniversity/fcma-toolbox


text:
# ifndef LANGIL_H
# define LANGIL_H

///////////////////////////////////////////////////////
//////   GILLESPIE CLASS
//////
//////						R. Perez-Carrasco
//////                               Created: 16 Nov '13
///////////////////////////////////////////////////////
/* Class for simulating master equation through langil algorythm. It defines three classes:
		- langil: In charge of the integration
		- species: contains a name a number and methods to change number
		- reaction: contains a name, a stoiciometry associated and a pointer to a propensity computing function
*/

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <string>
#include <fstream>
#include <iostream>
#include <unordered_map>
#include <vector>
#include <math.h>
#include <time.h>
#include <unistd.h>
#include <iomanip>
//#include <omp.h>

#include "species.h"
#include "cell_cycle.h"

#define VERBOSE 1
#define NON_VERBOSE 0

#define TIMESTEPMULT 3
#define TIMESTEP 0
#define ALL 1
#define NOWRITE 2
#define LASTPOINT 4

#define GILLESPIE 0
#define MACROSCOPIC 2
#define CLEEULER 1
#define CLEMILSTEIN 3

#define ADIABATIC true
#define NONADIABATIC false

#define GEOMETRIC true
#define NONGEOMETRIC false

#define DETERMINISTIC_TIME true
#define STOCHASTIC_TIME false

#define EXPONENTIAL_REACTION 0
#define GEOMETRIC_REACTION 1
#define DETERMINISTIC_REACTION 2
#define ADIABATIC_REACTION 3

using namespace std;

struct StructDivisionState{
	double time;
	double species0before;
	double species1before;
	double species0after;
	double species1after;
};

class reaction; // defined after langil
class species; // defined after reaction
class expressionzone; // defined in a separate file

class langil{

	friend class expressionzone;

//	friend class reaction; // reaction have access to the state of the system

	public:

		langil(string outputfilename, double vol=100, int seed=-1);
		void SetState(string aname, float anumber); //Set Initial conditions
		double GetState(string a); // get the state of the selected species
		void SetTime(double time); //Set current time
		void AddSpecies(string aname, float anum); // Add a new molecular species
		void AddSpeciesTimeTracking(string aname, float anum, double lifetime); //  Add a molecular species with a timer
		void AddCellCycleSpecies(); // add a cell cycle to the system
		void addCellPhase(double duration, int type_phase); // add a cell phase to the cell cycle
		void PrintSummarySpecies();
		string StateString();
		void Add_Reaction(string name,vector<int> stoich,double (*prop_f)(v_species&),int typeofrec = EXPONENTIAL_REACTION); // Add a reaction to the system overloaded to create directly
		void UpdateTime(); // Update times for the deterministic events
		void SetRunType(int rt,double dt=0);
		double GillespieStep(double timelimit = -1); // Advance one step in the Gillespie algorithm. Returns time of the step. Timelimit does not react if the reaction would take more than timelimit
		double LangevinEulerStep(double timelimit = -1); // Make an Euler-Maruyama integration step
		double LangevinMilsteinStep(double timelimit = -1); // Make an SDE integration step with the Milstein algorithm (not tested)
		double MacroStep(double timlimit = -1); // Make an Euler (deterministic) step
		double MacroAdiabaticStep(double dtt); //  Make an Euler step if the Adiabaitc species
		double Run(double time, bool verbose=true); // React until reach time time
		double RunTimes(int N,double time, bool verbose=true); // Run Run() N times
		double RunTransition(bool (*prop_f)(v_species&), bool (*prop_g)(v_species&), int loops, double maxT); // Run transition between two points, if loops=1 it also records the transition trajectory
		void WriteState(); // Write state in file
		void WriteTempState(); // Write state in file
		void WritePrevrecState(); // Write previous recorded state in file
		void SetWriteState(int flag, double tl=0);// Switch betweem different states:
		double (langil::*MakeStep)(double); // pointer to the actual integration chosen
			// TIMESTEP: write the state every tl time lapse
			// ALL: write after each reaction takes place
			// NOWRITE: don't write anything in output file
		double SetLangevinTimeStep(double dx); // set dt to a characteristic time that will increase with the size of the system
		void Add_TimeAction(double (*actionfunc)(double));
		void RunTimeAction();
		void setPhaseDuration(int phase, double duration, int type_phase);
		void Set_AdiabaticReaction(string aname);
		void Set_GeometricReaction(string aname);
		void Set_Boundary_Behaviour(void (*b_b)(v_species&, v_species&)); // Boundary behaviour
		void Set_StoreDivisionTimes();// activate the storage of division times in DivisionTimes vector
		void (*Boundary_Behaviour)(v_species&, v_species&); // boundary behaviour function
		vector<StructDivisionState> Get_HistoryDivision();
		void Reset_HistoryDivision();
		void SetDivisionHistory(int b,int a);

	protected:

		v_species x; //value of each species. shared_ptr is used since the vector can contain species and also childrens of that class (i.e. different classes)
		vector<reaction> r;//vector with all the reactions
 		double time; // current time of simulation
		gsl_rng * rng; // allocator for the rng generator
		vector<double> rnd; // rnd numbers
		vector<double> xtemp; // temporal vector for integrations
		v_species xtemp0,xtemp1; // temporal vector for integrations
		double pro,pro0; // auxiliar variables for the program
		int sto;  // auxiliar value for stoichiometry
		double totalprop,cumprop, detprop; // sum of the propensities
		double nexttau; // time of the next step
		vector<reaction>::iterator nextreaction; // reaction selected to react
		vector<int>* nextstoichiometry;
		double Omega; // volume of the system to change between concentrations and absolute numbers
		double signal; // variable that can modulate an external signal its time course
						// may be changed through time actions with Add_TimeAction
		int celldivided; // flag to track possible celldivision
		int geostoichiometry; // dummy variable to set random geometric stoichiometries

		bool record; // Set record of the trajectory
		ofstream trajfile; // File for output traj
		string outputfilename; // File for output name
		int fwrite; // flag for write state when integrating
		double timelapse; // dt for recording if TIMESTEP
		double nextrectime; // t for recording if TIMESTEP
		double totaltime; // total integration time for the trajectory
		double dt; // timestep used in integrations for macroscopic and langevin (if negative, it is used by Langevin as a prediction for dx)

		int cellcyclespeciesidx;  // index of the species vector that tracks cell cycle phases
 		cell_cycle cell; // cell_cycle object to manage cell cycle events

 		bool fStoreDivisionTimes; // flag to store Division Times
 		vector<StructDivisionState> DivisionHistoryVector; // Division history vector. Each component is a vector
 		StructDivisionState divisionstate;
		vector<double(*)(double)> actionfuncvec; // vector to time action functions

};

////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
class reaction{//Class containing the info of each reaction. Each reaction added
// through the method Add_Reaction(reaction) will be stored as an element in a vector
// by langil class. The pointer x, points out to the state of the langil system and
// is managed automatically by langil when a reaction is added to a langil instance
// It has also the property "adiabatic" meaning that it is a fast variable and can always be simulated
// from its deterministic behaviour

	friend class langil;

	public:

	reaction(string aname,vector<int> asto,double (*aprop_f)(v_species&)=NULL,int reacttype= EXPONENTIAL_REACTION,
		    vector<reaction>* x0=NULL){
			name=aname;
			stoichiometry=asto;
			prop_f=aprop_f;
			reactiontype=reacttype;
		}
		double GetPropensity(){
			return prop_f(*x);
		}
		double GetPropensity(v_species& x0){
			return prop_f(x0);
		}
		string GetName(){
			return name;
		}
		void SetPropensity(double (*p_f)(v_species&)){
			prop_f=p_f;
		}
		void SetState(v_species& x0){
			x=&x0;
		}
		vector<int>* GetStoichiometry(){
			return &stoichiometry;
		}
		bool IsAdiabatic(){
			return (reactiontype == ADIABATIC_REACTION);
		}
		bool IsNotAdiabatic(){
			return (reactiontype != ADIABATIC_REACTION);
		}
		bool IsGeometric(){
			return (reactiontype == GEOMETRIC_REACTION);
		}
		bool IsNotGeometric(){
			return (reactiontype != GEOMETRIC_REACTION);
		}
		bool IsNotDetermTime(){
			return (reactiontype != (DETERMINISTIC_REACTION));
		}
		bool IsDetermTime(){
			return (reactiontype == (DETERMINISTIC_REACTION));
		}
		void SetAdiabatic(){
			reactiontype = ADIABATIC_REACTION;
		}
		void SetGeometric(){
			reactiontype = GEOMETRIC_REACTION;
		}
		void SetDeterministicTime(){
			reactiontype = DETERMINISTIC_REACTION;
		}

	protected:

		string name; // name of the reactio
		vector<int> stoichiometry; // stoichiometry
		double (*prop_f)(v_species&); // propensity function
		v_species* x;  //value of each species. shared_ptr is used since the vector can contain species and also childrens of that class (i.e. different classes)
		int reactiontype; // can be any of the macros EXPONENTIAL_REACTION, DETERMINISITIC_REACTION, etc.

 // of the system

};

#endif


text:
/**
 * File: spce_pgp.h
 * Interface to PGplot routines.
 *
 */

#ifndef _SPCE_PGP_H

#define _SPCE_PGP_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_vector.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spce_pathlength.h"
#include "spce_output.h"

extern void
pgp_stamp_image (const ap_pixel * const ap_p, const int n_sub,
		 const object * ob, int beamnum, char filename[],
		 const int negative);
#endif /* !_SPCE_PGP_H */

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe