text:
/**
*\file utils.h
*\author weckyy702 (weckyy702@gmail.com)
*\brief Utility and core Header file
*\date 2021-03-25
*
*MIT License
*Copyright (c) [2021] [Weckyy702 (weckyy702@gmail.com | https://github.com/Weckyy702)]
*Permission is hereby granted, free of charge, to any person obtaining a copy
*of this software and associated documentation files (the "Software"), to deal
*in the Software without restriction, including without limitation the rights
*to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
*copies of the Software, and to permit persons to whom the Software is
*furnished to do so, subject to the following conditions:
*
*The above copyright notice and this permission notice shall be included in all
*copies or substantial portions of the Software.
*
*THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
*IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
*FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
*AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
*LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
*OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
*SOFTWARE.
*
*/
#ifndef RAYCHEL_CORE_H
#define RAYCHEL_CORE_H

#include <cstdint>
#include <limits>
#if !((defined(__cplusplus) && __cplusplus >= 201703L) || (defined(_MSVC_LANG) && _MSVC_LANG >= 201703L))
    #error "C++(17) compilation is required!"
#endif

#include <cassert>
#include <cmath>
#include <cstddef>
#include <gsl/gsl>
#include <type_traits>

//convenience headers
#include <memory>
#include <optional>
#include <utility>
#include <vector>

#include "CMakeSettings.h"
#include "Logger.h"
#include "RaychelMath/constants.h"
#include "RaychelMath/equivalent.h"

#if defined(__clang__) || defined(__GNUC__)
    #define RAYCHEL_FUNC_NAME __PRETTY_FUNCTION__
#elif defined(_MSC_VER)
    #define RAYCHEL_FUNC_NAME __FUNCSIG__
#else
    #error "Unknown compiler detected!"
#endif

#ifdef RAYCHEL_DEBUG
    #define RAYCHEL_LOG(...) Logger::debug(RAYCHEL_FUNC_NAME, ": ", __VA_ARGS__, '\n');
#else
    #define RAYCHEL_LOG(...)
#endif

//terminate the application with the provided message
#ifndef RAYCHEL_NO_LOGGER
    #define RAYCHEL_TERMINATE(...)                                                                                \
        Logger::fatal(RAYCHEL_FUNC_NAME, " at (", __FILE__, ":", __LINE__, "): ", __VA_ARGS__, '\n');                                                \
        std::exit(0x41);
#else
    #define RAYCHEL_TERMINATE(...) std::exit(0x41);
#endif

#if defined(RAYCHEL_DEBUG) || !defined(NDEBUG)
    #define RAYCHEL_ASSERT(exp)                                                                                \
        if (!(exp)) {                                                                                \
            RAYCHEL_TERMINATE("Assertion '", GSL_STRINGIFY(exp), "' failed!");                                                                       \
        }
#else
    #define RAYCHEL_ASSERT(exp) Expects(exp)
#endif

#define RAYCHEL_ASSERT_NOT_REACHED RAYCHEL_TERMINATE("Assertion failed! Expected to not execute ", __FILE__, ":", __LINE__)

//#define RAYCHEL_LOGICALLY_EQUAL //<-- activates logical equivalency for vector-like types

#define RAYCHEL_THROW_EXCEPTION(msg, fatal) throw ::Raychel::exception_context{msg, RAYCHEL_FUNC_NAME, fatal};

#define RAYCHEL_ASSERT_NORMALIZED(vec) RAYCHEL_ASSERT(equivalent(magSq(vec), 1.0f));

namespace Raychel {

    using gsl::byte, gsl::not_null;
    using std::size_t;

    template <typename _num>
    struct vec2Imp;
    template <typename _num>
    struct vec3Imp;
    template <typename _num>
    struct colorImp;
    /*TODO: implement
	template<typename _num>
	class colorAlphaImp;
	*/
    template <typename _num>
    class QuaternionImp;
    template <typename _num>
    struct TransformImp;

    template <typename _number>
    constexpr _number sq(_number x)
    {
        static_assert(std::is_arithmetic_v<_number>, "Raychel::sq<T> requires T to be of arithmetic type!");
        return x * x;
    }

    /**
	*\brief Linearly interpolate between two numbers
	*
	*\tparam _number Type of number to interpolate. Must be arithmetic
	*\param a first number (x=0.0)
	*\param b second number (x=1.0)
	*\param x value of interpolation
	*\return _number the interpolated number
	*/
    template <typename _number>
    constexpr _number lerp(_number a, _number b, long double x)
    {
        return (x * b) + ((1.0 - x) * a);
    }

    template <typename _integral>
    constexpr _integral bit(size_t shift)
    {
        static_assert(std::is_integral_v<_integral>, "Raychel::bit<T> requires T to be of integral type!");
        RAYCHEL_ASSERT(shift < (sizeof(_integral) * 8));
        return static_cast<_integral>(1U << shift);
    }

    /**
	*\brief Get number of digits in a number
	*
	*\param num The number with its digits
	*\return constexpr int
	*/
    constexpr int numDigits(unsigned long long num, unsigned int base = 10U)
    {
        int digits = 0;

        do {
            digits++;
            num /= base;
        } while (num != 0);

        return digits;
    }

} // namespace Raychel

#endif // !RAYCHEL_CORE_H


text:
/**
 *
 * @file testing_spemv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Dulceneia Becker
 * @date 2011-10-06
 * @generated s Tue Jan  7 11:45:19 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

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

#define REAL
#undef COMPLEX

/*--------------------------------------------------------------
 * Check the pemv
 */
static int check_solution(PLASMA_enum trans, PLASMA_enum storev,
                          int M, int N, int L,
                          float alpha, float *A, int LDA,
                                                    float *X, int INCX,
                          float beta,  float *Y0, int INCY0,
                                                    float *Y,  int INCY,
                          float *W, float *Rnorm)
{
    int k;
    float eps = LAPACKE_slamch_work('e');
    float *work;
    float mzone = -1.0;

    /* Copy x to w */
    if ( trans == PlasmaNoTrans ) {
        k = N;
    } else {
        k = M;
    }

    work = (float *)malloc(k * sizeof(float));
    cblas_scopy(k, Y0, INCY0, W, 1);

    /* w = a A x + b w */
    cblas_sgemv(CblasColMajor, (CBLAS_TRANSPOSE)trans,
                M, N,
                (alpha), A,  LDA,
                                    X,  INCX,
                (beta),  W,  1);

    /* y - w */
    cblas_saxpy(k, (mzone), Y, INCY, W, 1);

    /* Max Norm */
    *Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'm', 1, k, W, 1, work);

    if ( (*Rnorm / (M*N)) > eps) {
        return 1;
    } else {
        return 0;
    }
}

/*--------------------------------------------------------------
 * Testing SPEMV
 */
int testing_spemv(int argc, char **argv)
{
    /* Check for number of arguments*/
    if ( argc != 1) {
        USAGE("PEMV", "N",
              "   - N      : number of columns\n");
        return -1;
    }

    /* Args */
    int arg_n = atoi(argv[0]);

    /* Local variables */
    float *A, *X, *Y, *A0, *Y0, *work;
    float alpha, beta, alpha0, beta0;
    int n    = arg_n;
    int lda  = arg_n;

    int info_solution = 0;
    int i, j, k, t;
    int nbtests = 0;
    int nfails = 0;
    int storev;
    int l = 0;
    int m = n;
    int incx = 1;
    int incy = 1;
    char *cstorev;
    float rnorm;
    float eps = LAPACKE_slamch_work('e');

    /* Allocate Data */
    A    = (float *)malloc(lda*n*sizeof(float));
    A0   = (float *)malloc(lda*n*sizeof(float));
    X    = (float *)malloc(lda*n*sizeof(float));
    Y    = (float *)malloc(lda*n*sizeof(float));
    Y0   = (float *)malloc(    n*sizeof(float));
    work = (float *)malloc(  2*n*sizeof(float));

    LAPACKE_slarnv_work(1, ISEED, 1, &alpha0);
    LAPACKE_slarnv_work(1, ISEED, 1, &beta0 );

    /* Check if unable to allocate memory */
    if ( (!A) || (!X) || (!Y0) || (!work) ) {
        printf("Out of Memory \n ");
        exit(0);
    }

     /* Initialize Data */
    PLASMA_splrnt(n, n, A,  lda, 479 );
    PLASMA_splrnt(n, n, X,  lda, 320 );
    PLASMA_splrnt(n, 1, Y0, n,   573 );

    printf("\n");
    printf("------ TESTS FOR PLASMA SPEMV ROUTINE -------  \n");
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf(" The relative machine precision (eps) is %e \n",eps);
    printf(" Computational tests pass if scaled residual is less than eps.\n");
    printf("\n");

    nfails = 0;
    for (i=0; i<6; i++) {

        /* m and n cannot be greater than lda (arg_n) */
        switch (i) {
        case 0: l = 0;       m = arg_n;   n = m;        break;
        case 1: l = 0;       m = arg_n;   n = arg_n/2;  break; /**/
        case 2: l = arg_n;   m = l;       n = l;        break;
        case 3: l = arg_n/2; m = l;       n = arg_n;    break;
        case 4: l = arg_n/2; m = arg_n-l; n = l;        break;
        case 5: l = arg_n/3; m = arg_n-l; n = arg_n/2;  break; /**/
        }

        /* Colwise ConjTrans & Rowwise NoTrans */
#ifdef COMPLEX
        for (t=0; t<3; t++) {
#else
        for (t=0; t<2; t++) {
#endif

            /* Swap m and n for transpose cases */
            if ( t == 1 ) {
                k = m; m = n; n = k;
            }

            LAPACKE_slacpy_work( LAPACK_COL_MAJOR, 'A', m, n,
                                 A, lda, A0, lda);

            if ( trans[t] == PlasmaNoTrans ) {
                storev = PlasmaRowwise;
                cstorev = storevstr[0];

                /* zeroed the upper right triangle */
                int64_t i, j;
                for (j=(n-l); j<n; j++) {
                    for (i=0; i<(j-(n-l)); i++) {
                        A0[i+j*lda] = 0.0;
                    }
                }
            }
            else {
                storev = PlasmaColumnwise;
                cstorev = storevstr[1];

                /* zeroed the lower left triangle */
                int64_t i, j;
                for (j=0; j<(l-1); j++) {
                    for (i=(m-l+1+j); i<m; i++) {
                        A0[i+j*lda] = 0.0;
                    }
                }
            }

            for (j=0; j<3; j++) {

                /* Choose alpha and beta */
                alpha = ( j==1 ) ? 0.0 : alpha0;
                beta  = ( j==2 ) ? 0.0 : beta0;

                /* incx and incy: 1 or lda */
                for (k=0; k<4; k++) {
                    switch (k) {
                    case 0:  incx = 1;    incy = 1;    break;
                    case 1:  incx = 1;    incy = lda;  break;
                    case 2:  incx = lda;  incy = 1;    break;
                    case 3:  incx = lda;  incy = lda;  break;
                    }

                    /* initialize Y with incy */
                    cblas_scopy(n, Y0, 1, Y, incy);

                    /* SPEMV */
                    CORE_spemv( trans[t], storev, m, n, l,
                                alpha, A, lda,
                                       X, incx,
                                beta,  Y, incy,
                                work);

                    /* Check the solution */
                    info_solution = check_solution(trans[t], storev,
                                                   m, n, l,
                                                   alpha, A0,  lda,
                                                          X,   incx,
                                                   beta,  Y0,  1,
                                                          Y,   incy,
                                                   work, &rnorm);

                    if ( info_solution != 0 ) {
                        nfails++;
                        printf("Failed: t=%s, s=%s, M=%3d, N=%3d, L=%3d, alpha=%e, incx=%3d, beta=%e, incy=%3d, rnorm=%e\n",
                               transstr[t], cstorev, m, n, l, (alpha), incx, (beta), incy, rnorm );
                    }
                    nbtests++;
                }
            }
        }
    }

    if ( nfails )
        printf("%d / %d tests failed\n", nfails, nbtests);

    printf("***************************************************\n");
    if (nfails == 0) {
        printf(" ---- TESTING SPEMV ...... PASSED !\n");
    }
    else {
        printf(" ---- TESTING SPEMV ... FAILED !\n");
    }
    printf("***************************************************\n");

    free( A0 );
    free( A );
    free( X );
    free( Y0 );
    free( Y );

    return 0;
}

ext: c


text:
/* im_invfft
 *
 * Copyright: 1990, N. Dessipris.
 *
 * Author: Nicos Dessipris
 * Written on: 12/04/1990
 * Modified on :
 * 28/6/95 JC
 *	- rewritten, based on new im_fwfft() code
 * 10/9/98 JC
 *	- frees memory more quickly
 * 2/4/02 JC
 *	- fftw code added
 * 13/7/02 JC
 *	- Type reset
 * 27/2/03 JC
 *	- tiny speed-up ... save 1 copy on write
 * 22/1/04 JC
 *	- oops, fix for segv on wider than high fftw transforms
 * 3/11/04
 *	- added fftw3 support
 * 7/2/10
 * 	- gtkdoc
 */

/*

    This file is part of VIPS.

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

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

    You should have received a copy of the GNU Lesser General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

 */

/*

    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk

 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /*HAVE_CONFIG_H*/
#include <vips/intl.h>

#include <stdio.h>
#include <math.h>

#ifdef HAVE_FFTW
#include <fftw.h>
#endif /*HAVE_FFTW*/

#ifdef HAVE_FFTW3
#include <fftw3.h>
#endif /*HAVE_FFTW3*/

#include <vips/vips.h>
#include <vips/internal.h>

#ifdef WITH_DMALLOC
#include <dmalloc.h>
#endif /*WITH_DMALLOC*/

#ifdef HAVE_FFTW
/* Call fftw for a 1 band image.
 */
static int
invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
	fftwnd_plan plan;

	IMAGE *cmplx = im_open_local( out, "invfft1:1", "t" );

	/* Make dp complex image.
	 */
	if( !cmplx || im_pincheck( in ) || im_poutcheck( out ) )
		return( -1 );
	if( in->Coding != IM_CODING_NONE || in->Bands != 1 ) {
                im_error( "im_invfft", "%s", _( "one band uncoded only" ) );
                return( -1 );
	}
	if( im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
                return( -1 );

	/* Make the plan for the transform. Yes, they really do use nx for
	 * height and ny for width.
	 */
	if( !(plan = fftw2d_create_plan( in->Ysize, in->Xsize,
		FFTW_BACKWARD,
		FFTW_MEASURE | FFTW_USE_WISDOM | FFTW_IN_PLACE )) ) {
                im_error( "im_invfft",
			"%s", _( "unable to create transform plan" ) );
		return( -1 );
	}

	fftwnd_one( plan, (fftw_complex *) cmplx->data, NULL );

	fftwnd_destroy_plan( plan );

	/* Copy to out.
	 */
        if( im_copy( cmplx, out ) )
                return( -1 );

	return( 0 );
}
#else /*!HAVE_FFTW*/
#ifdef HAVE_FFTW3
/* Complex to complex inverse transform.
 */
static int
invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
	fftw_plan plan;

	IMAGE *cmplx = im_open_local( out, "invfft1:1", "t" );

	/* We have to have a separate buffer for the planner to work on.
	 */
	double *planner_scratch = IM_ARRAY( dummy,
		in->Xsize * in->Ysize * 2, double );

	/* Make dp complex image.
	 */
	if( !cmplx || im_pincheck( in ) || im_poutcheck( out ) )
		return( -1 );
	if( in->Coding != IM_CODING_NONE || in->Bands != 1 ) {
                im_error( "im_invfft",
			"%s", _( "one band uncoded only" ) );
                return( -1 );
	}
	if( im_clip2fmt( in, cmplx, IM_BANDFMT_DPCOMPLEX ) )
                return( -1 );

	/* Make the plan for the transform. Yes, they really do use nx for
	 * height and ny for width.
	 */
	if( !(plan = fftw_plan_dft_2d( in->Ysize, in->Xsize,
		(fftw_complex *) planner_scratch,
		(fftw_complex *) planner_scratch,
		FFTW_BACKWARD,
		0 )) ) {
                im_error( "im_invfft",
			"%s", _( "unable to create transform plan" ) );
		return( -1 );
	}

	fftw_execute_dft( plan,
		(fftw_complex *) cmplx->data, (fftw_complex *) cmplx->data );

	fftw_destroy_plan( plan );

	/* Copy to out.
	 */
        if( im_copy( cmplx, out ) )
                return( -1 );

	return( 0 );
}
#else /*!HAVE_FFTW3*/
/* Fall back to VIPS's built-in fft
 */
static int
invfft1( IMAGE *dummy, IMAGE *in, IMAGE *out )
{
	int bpx = im_ispoweroftwo( in->Xsize );
	int bpy = im_ispoweroftwo( in->Ysize );
	float *buf, *q, *p1, *p2;
	int x, y;

	/* Buffers for real and imaginary parts.
	 */
	IMAGE *real = im_open_local( dummy, "invfft1:1", "t" );
	IMAGE *imag = im_open_local( dummy, "invfft1:2", "t" );

	/* Temps.
	 */
	IMAGE *t1 = im_open_local( dummy, "invfft1:3", "p" );
	IMAGE *t2 = im_open_local( dummy, "invfft1:4", "p" );

	if( !real || !imag || !t1 )
		return( -1 );
        if( im_pincheck( in ) || im_outcheck( out ) )
                return( -1 );
        if( in->Coding != IM_CODING_NONE ||
		in->Bands != 1 || !im_iscomplex( in ) ) {
                im_error( "im_invfft",
			"%s", _( "one band complex uncoded only" ) );
                return( -1 );
	}
	if( !bpx || !bpy ) {
		im_error( "im_invfft",
			"%s", _( "sides must be power of 2" ) );
		return( -1 );
	}

	/* Make sure we have a single-precision complex input image.
	 */
	if( im_clip2fmt( in, t1, IM_BANDFMT_COMPLEX ) )
		return( -1 );

	/* Extract real and imag parts. We have to complement the imaginary.
	 */
	if( im_c2real( t1, real ) )
		return( -1 );
	if( im_c2imag( t1, t2 ) || im_lintra( -1.0, t2, 0.0, imag ) )
		return( -1 );

	/* Transform!
	 */
	if( im__fft_sp( (float *) real->data, (float *) imag->data,
		bpx - 1, bpy - 1 ) ) {
                im_error( "im_invfft",
			"%s", _( "fft_sp failed" ) );
                return( -1 );
	}

	/* WIO to out.
	 */
        if( im_cp_desc( out, in ) )
                return( -1 );
	out->BandFmt = IM_BANDFMT_COMPLEX;
        if( im_setupout( out ) )
                return( -1 );
	if( !(buf = (float *) IM_ARRAY( dummy,
		IM_IMAGE_SIZEOF_LINE( out ), PEL )) )
		return( -1 );

	/* Gather together real and imag parts.
	 */
	for( p1 = (float *) real->data, p2 = (float *) imag->data,
		y = 0; y < out->Ysize; y++ ) {
		q = buf;

		for( x = 0; x < out->Xsize; x++ ) {
			q[0] = *p1++;
			q[1] = *p2++;
			q += 2;
		}

		if( im_writeline( y, out, (PEL *) buf ) )
			return( -1 );
	}

	return( 0 );
}
#endif /*HAVE_FFTW3*/
#endif /*HAVE_FFTW*/

/**
 * im_invfft:
 * @in: input image
 * @out: output image
 *
 * Transform an image from Fourier space to real space. The result is complex.
 * If you are OK with a real result, use im_invfftr() instead, it's quicker.
 *
 * VIPS uses the fftw3 or fftw2 Fourier transform libraries if possible. If
 * they were not available when VIPS was built, it falls back to it's own
 * FFT functions which are slow and only work for square images whose sides
 * are a power of two.
 *
 * See also: im_invfftr(), im_fwfft(), im_disp_ps().
 *
 * Returns: 0 on success, -1 on error.
 */
int
im_invfft( IMAGE *in, IMAGE *out )
{
	IMAGE *dummy = im_open( "im_invfft:1", "p" );

	if( !dummy )
		return( -1 );
	if( im__fftproc( dummy, in, out, invfft1 ) ) {
		im_close( dummy );
		return( -1 );
	}
	im_close( dummy );

	if( out->Bands == 1 )
		out->Type = IM_TYPE_B_W;
	else
		out->Type = IM_TYPE_MULTIBAND;

	return( 0 );
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef __GSL_SORT_VECTOR_H__
#define __GSL_SORT_VECTOR_H__

#include <gsl/gsl_sort_vector_long_double.h>
#include <gsl/gsl_sort_vector_double.h>
#include <gsl/gsl_sort_vector_float.h>

#include <gsl/gsl_sort_vector_ulong.h>
#include <gsl/gsl_sort_vector_long.h>

#include <gsl/gsl_sort_vector_uint.h>
#include <gsl/gsl_sort_vector_int.h>

#include <gsl/gsl_sort_vector_ushort.h>
#include <gsl/gsl_sort_vector_short.h>

#include <gsl/gsl_sort_vector_uchar.h>
#include <gsl/gsl_sort_vector_char.h>

#endif /* __GSL_SORT_VECTOR_H__ */

ext: h

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

// std
#include <optional>

// GSL
#include <gsl/gsl_assert>

// EnTT
#include <entt/entt.hpp>

// Qt
#include <QtWidgets>

max_stars_repo_name: Qt-Widgets/PolyDock

max_issues_repo_name: Qt-Widgets/PolyDock

max_forks_repo_name: Qt-Widgets/PolyDock


text:
/*
 *     Copyright 2021-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */
#pragma once
#define MEMCACHED_ENGINE_H

#include <sys/types.h>
#include <functional>
#include <memory>
#include <string_view>
#include <unordered_set>
#include <utility>

#include <gsl/gsl>
#include <optional>

#include <memcached/visibility.h>

#include "memcached/collections.h"
#include "memcached/engine_common.h"
#include "memcached/thread_pool_config.h"
#include "memcached/types.h"
#include "memcached/vbucket.h"

namespace cb::durability {
class Requirements;
} // namespace cb::durability

namespace cb::mcbp {
class Request;
} // namespace cb::mcbp

namespace cb::prometheus {
enum class Cardinality;
} // namespace cb::prometheus

class BucketStatCollector;
class StatCollector;

/*! \mainpage memcached public API
 *
 * \section intro_sec Introduction
 *
 * The memcached project provides an API for providing engines as well
 * as data definitions for those implementing the protocol in C.  This
 * documentation will explain both to you.
 *
 * \section docs_sec API Documentation
 *
 * Jump right into the modules docs to get started.
 *
 * \example default_engine.cc
 */

/**
 * \defgroup Engine Storage Engine API
 * \defgroup Protex Protocol Extension API
 * \defgroup Protocol Binary Protocol Structures
 *
 * \addtogroup Engine
 * @{
 *
 * Most interesting here is to implement engine_interface_v1 for your
 * engine.
 */

struct DocKey;
struct ServerBucketIface;
struct ServerCoreIface;
struct ServerCallbackIface;
struct ServerLogIface;
struct ServerCookieIface;
struct ServerDocumentIface;
union protocol_binary_request_header;

struct ServerApi {
    ServerCoreIface* core = nullptr;
    ServerCallbackIface* callback = nullptr;
    ServerLogIface* log = nullptr;
    ServerCookieIface* cookie = nullptr;
    ServerDocumentIface* document = nullptr;
    ServerBucketIface* bucket = nullptr;
};

using GET_SERVER_API = ServerApi* (*)();

namespace cb {
using EngineErrorItemPair = std::pair<cb::engine_errc, cb::unique_item_ptr>;

using EngineErrorMetadataPair = std::pair<engine_errc, item_info>;

enum class StoreIfStatus {
    Continue,
    Fail,
    GetItemInfo // please get me the item_info
};

using StoreIfPredicate = std::function<StoreIfStatus(
        const std::optional<item_info>&, cb::vbucket_info)>;

struct EngineErrorCasPair {
    engine_errc status;
    uint64_t cas;
};

/// Result of getVBucketHlcNow()
struct HlcTime {
    enum class Mode { Real, Logical };

    /// Seconds since Unix epoch.
    std::chrono::seconds now;
    Mode mode;
};
} // namespace cb

/**
 * The different compression modes that a bucket supports
 */
enum class BucketCompressionMode : uint8_t {
    Off,     //Data will be stored as uncompressed
    Passive, //Data will be stored as provided by the client
    Active   //Bucket will actively try to compress stored
             //data
};

/* The default minimum compression ratio */
static const float default_min_compression_ratio = 1.2f;

/* The default maximum size for a value */
static const size_t default_max_item_size = 20 * 1024 * 1024;

namespace cb::engine {
/**
 * Definition of the features that an engine can support
 */
enum class Feature : uint16_t {
    Collections = 1,
};

using FeatureSet = std::unordered_set<Feature>;
} // namespace cb::engine

namespace std {
template <>
struct hash<cb::engine::Feature> {
public:
    size_t operator()(const cb::engine::Feature& f) const {
        return static_cast<size_t>(f);
    }
};
} // namespace std

/**
 * Definition of the first version of the engine interface
 */
struct MEMCACHED_PUBLIC_CLASS EngineIface {
    virtual ~EngineIface() = default;

    /**
     * Initialize an engine instance.
     * This is called *after* creation, but before the engine may be used.
     *
     * @param handle the engine handle
     * @param config_str configuration this engine needs to initialize itself.
     */
    virtual cb::engine_errc initialize(const char* config_str) = 0;

    /**
     * Tear down this engine.
     *
     * @param force the flag indicating the force shutdown or not.
     */
    virtual void destroy(bool force) = 0;

    /// Callback that a cookie will be disconnected
    virtual void disconnect(gsl::not_null<const void*> cookie){};

    /**
     * Initiate the bucket shutdown logic (disconnect clients etc)
     */
    virtual void initiate_shutdown() {
        // empty
    }

    // Set the number or reader threads
    virtual void set_num_reader_threads(ThreadPoolConfig::ThreadCount num) {
        // ignored
    }

    // Set the number or writer threads
    virtual void set_num_writer_threads(ThreadPoolConfig::ThreadCount num) {
        // ignored
    }

    // Set the number of storage threads
    virtual void set_num_storage_threads(
            ThreadPoolConfig::StorageThreadCount num) {
        // ignored
    }

    /**
     * Request the engine to cancel all of the ongoing requests which
     * may have cookies in an ewouldblock state.
     *
     * This method is to removed in CC when we'll tighten the logic for
     * bucket deletion to use the following logic:
     *
     * 1) Stop any new requests into the engine.
     * 2) Tell the engine to complete (cancel) anything currently in-flight
     *    (this notification holds a write lock to the engine so it won't
     *    race with any front end threads)
     * 3) wait for in-flight ops (i.e. step B) to finish then delete bucket.
     *
     * In Mad-Hatter initiate_shutdown may race with all of the frontend
     * worker threads, and could lead to operations being added after we've
     * inspected the vbucket. To work around that problem this new method
     * was introduced and will be called multiple times during bucket
     * deletion to work around potential race situations.
     */
    virtual void cancel_all_operations_in_ewb_state() {
        // empty
    }

    /*
     * Item operations.
     */

    /**
     * Allocate an item (extended API)
     *
     * @param cookie The cookie provided by the frontend
     * @param key the item's key
     * @param nbytes the number of bytes that will make up the
     *               value of this item.
     * @param priv_nbytes The number of bytes in nbytes containing
     *                    system data (and may exceed the item limit).
     * @param flags the item's flags
     * @param exptime the maximum lifetime of this item
     * @param vbucket virtual bucket to request allocation from
     * @return pair containing the item and the items information
     * @thows cb::engine_error with:
     *
     *   * `cb::engine_errc::no_bucket` The client is bound to the dummy
     *                                  `no bucket` which don't allow
     *                                  allocations.
     *
     *   * `cb::engine_errc::no_memory` The bucket is full
     *
     *   * `cb::engine_errc::too_big` The requested memory exceeds the
     *                                limit set for items in the bucket.
     *
     *   * `cb::engine_errc::disconnect` The client should be disconnected
     *
     *   * `cb::engine_errc::not_my_vbucket` The requested vbucket belongs
     *                                       to someone else
     *
     *   * `cb::engine_errc::temporary_failure` Temporary failure, the
     *                                          _client_ should try again
     *
     *   * `cb::engine_errc::too_busy` Too busy to serve the request,
     *                                 back off and try again.
     */
    virtual std::pair<cb::unique_item_ptr, item_info> allocateItem(
            gsl::not_null<const void*> cookie,
            const DocKey& key,
            size_t nbytes,
            size_t priv_nbytes,
            int flags,
            rel_time_t exptime,
            uint8_t datatype,
            Vbid vbucket) = 0;

    /**
     * Remove an item.
     *
     * @param cookie The cookie provided by the frontend
     * @param key the key identifying the item to be removed
     * @param cas The cas to value to use for delete (or 0 as wildcard)
     * @param vbucket the virtual bucket id
     * @param durability the durability specifier for the command
     * @param mut_info On a successful remove write the mutation details to
     *                 this address.
     *
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::engine_errc remove(
            gsl::not_null<const void*> cookie,
            const DocKey& key,
            uint64_t& cas,
            Vbid vbucket,
            const std::optional<cb::durability::Requirements>& durability,
            mutation_descr_t& mut_info) = 0;

    /**
     * Indicate that a caller who received an item no longer needs
     * it.
     *
     * @param item the item to be released
     */
    virtual void release(gsl::not_null<ItemIface*> item) = 0;

    /**
     * Retrieve an item.
     *
     * @param cookie The cookie provided by the frontend
     * @param item output variable that will receive the located item
     * @param key the key to look up
     * @param vbucket the virtual bucket id
     * @param documentStateFilter The document to return must be in any of
     *                            of these states. (If `Alive` is set, return
     *                            KEY_ENOENT if the document in the engine
     *                            is in another state)
     *
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::EngineErrorItemPair get(gsl::not_null<const void*> cookie,
                                        const DocKey& key,
                                        Vbid vbucket,
                                        DocStateFilter documentStateFilter) = 0;

    /**
     * Optionally retrieve an item. Only non-deleted items may be fetched
     * through this interface (Documents in deleted state may be evicted
     * from memory and we don't want to go to disk in order to fetch these)
     *
     * @param cookie The cookie provided by the frontend
     * @param key the key to look up
     * @param vbucket the virtual bucket id
     * @param filter callback filter to see if the item should be returned
     *               or not. If filter returns false the item should be
     *               skipped.
     *               Note: the filter is applied only to the *metadata* of the
     *               item_info - i.e. the `value` should not be expected to be
     *               present when filter is invoked.
     * @return A pair of the error code and (optionally) the item
     */
    virtual cb::EngineErrorItemPair get_if(
            gsl::not_null<const void*> cookie,
            const DocKey& key,
            Vbid vbucket,
            std::function<bool(const item_info&)> filter) = 0;

    /**
     * Retrieve metadata for a given item.
     *
     * @param cookie The cookie provided by the frontend
     * @param key the key to look up
     * @param vbucket the virtual bucket id
     *
     * @return  Pair (cb::engine_errc::success, Metadata) if all goes well
     */
    virtual cb::EngineErrorMetadataPair get_meta(
            gsl::not_null<const void*> cookie,
            const DocKey& key,
            Vbid vbucket) = 0;

    /**
     * Lock and Retrieve an item.
     *
     * @param cookie The cookie provided by the frontend
     * @param key the key to look up
     * @param vbucket the virtual bucket id
     * @param lock_timeout the number of seconds to hold the lock
     *                     (0 == use the engines default lock time)
     *
     * @return A pair of the error code and (optionally) the item
     */
    virtual cb::EngineErrorItemPair get_locked(
            gsl::not_null<const void*> cookie,
            const DocKey& key,
            Vbid vbucket,
            uint32_t lock_timeout) = 0;

    /**
     * Unlock an item.
     *
     * @param cookie The cookie provided by the frontend
     * @param key the key to look up
     * @param vbucket the virtual bucket id
     * @param cas the cas value for the locked item
     *
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::engine_errc unlock(gsl::not_null<const void*> cookie,
                                   const DocKey& key,
                                   Vbid vbucket,
                                   uint64_t cas) = 0;

    /**
     * Get and update the expiry time for the document
     *
     * @param cookie The cookie provided by the frontend
     * @param key the key to look up
     * @param vbucket the virtual bucket id
     * @param expirytime the new expiry time for the object
     * @param durability An optional durability requirement
     * @return A pair of the error code and (optionally) the item
     */
    virtual cb::EngineErrorItemPair get_and_touch(
            gsl::not_null<const void*> cookie,
            const DocKey& key,
            Vbid vbucket,
            uint32_t expirytime,
            const std::optional<cb::durability::Requirements>& durability) = 0;

    /**
     * Store an item into the underlying engine with the given
     * state. If the DocumentState is set to DocumentState::Deleted
     * the document shall not be returned unless explicitly asked for
     * documents in that state, and the underlying engine may choose to
     * purge it whenever it please.
     *
     * @param cookie The cookie provided by the frontend
     * @param item the item to store
     * @param cas the CAS value for conditional sets
     * @param semantics the semantics of the store operation
     * @param durability An optional durability requirement
     * @param document_state The state the document should have after
     *                       the update
     * @param preserveTtl if set to true the existing documents TTL should
     *                    be used.
     *
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::engine_errc store(
            gsl::not_null<const void*> cookie,
            gsl::not_null<ItemIface*> item,
            uint64_t& cas,
            StoreSemantics semantics,
            const std::optional<cb::durability::Requirements>& durability,
            DocumentState document_state,
            bool preserveTtl) = 0;

    /**
     * Store an item into the underlying engine with the given
     * state only if the predicate argument returns true when called against an
     * existing item.
     *
     * Optional interface; not supported by all engines.
     *
     * @param cookie The cookie provided by the frontend
     * @param item the item to store
     * @param cas the CAS value for conditional sets
     * @param semantics the semantics of the store operation
     * @param predicate a function that will be called from the engine the
     *                  result of which determines how the store behaves.
     *                  The function is given any existing item's item_info (as
     *                  a std::optional) and a cb::vbucket_info object. In the
     *                  case that the optional item_info is not initialised the
     *                  function can return cb::StoreIfStatus::GetInfo to
     *                  request that the engine tries to get the item_info, the
     *                  engine may ignore this return code if it knows better
     *                  i.e. a memory only engine and no item_info can be
     *                  fetched. The function can also return ::Fail if it
     *                  wishes to fail the store_if (returning predicate_failed)
     *                  or the predicate can return ::Continue and the rest of
     *                  the store_if will execute (and possibly fail for other
     *                  reasons).
     * @param durability An optional durability requirement
     * @param document_state The state the document should have after
     *                       the update
     * @param preserveTtl if set to true the existing documents TTL should
     *                    be used.
     *
     * @return a std::pair containing the engine_error code and new CAS
     */
    virtual cb::EngineErrorCasPair store_if(
            gsl::not_null<const void*> cookie,
            gsl::not_null<ItemIface*> item,
            uint64_t cas,
            StoreSemantics semantics,
            const cb::StoreIfPredicate& predicate,
            const std::optional<cb::durability::Requirements>& durability,
            DocumentState document_state,
            bool preserveTtl) {
        return {cb::engine_errc::not_supported, 0};
    }

    /**
     * Flush the cache.
     *
     * Optional interface; not supported by all engines.
     *
     * @param cookie The cookie provided by the frontend
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::engine_errc flush(gsl::not_null<const void*> cookie) {
        return cb::engine_errc::not_supported;
    }

    /*
     * Statistics
     */

    /**
     * Get statistics from the engine.
     *
     * @param cookie The cookie provided by the frontend
     * @param key optional argument to stats
     * @param value optional value for the given stat group
     * @param add_stat callback to feed results to the output
     *
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::engine_errc get_stats(gsl::not_null<const void*> cookie,
                                      std::string_view key,
                                      std::string_view value,
                                      const AddStatFn& add_stat) = 0;

    /**
     * Get statistics for Prometheus exposition from the engine.
     * Some engines may not support this.
     *
     * @param collector where the bucket may register stats
     *
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::engine_errc get_prometheus_stats(
            const BucketStatCollector& collector,
            cb::prometheus::Cardinality cardinality) {
        return cb::engine_errc::not_supported;
    }

    /**
     * Reset the stats.
     *
     * @param cookie The cookie provided by the frontend
     */
    virtual void reset_stats(gsl::not_null<const void*> cookie) = 0;

    /**
     * Any unknown command will be considered engine specific.
     *
     * @param cookie The cookie provided by the frontend
     * @param request The request from the client
     * @param response function to transmit data
     *
     * @return cb::engine_errc::success if all goes well
     */
    virtual cb::engine_errc unknown_command(const void* cookie,
                                            const cb::mcbp::Request& request,
                                            const AddResponseFn& response) {
        return cb::engine_errc::not_supported;
    }

    /**
     * Set the CAS id on an item.
     */
    virtual void item_set_cas(gsl::not_null<ItemIface*> item,
                              uint64_t cas) = 0;

    /**
     * Set the data type on an item.
     */
    virtual void item_set_datatype(gsl::not_null<ItemIface*> item,
                                   protocol_binary_datatype_t datatype) = 0;

    /**
     * Get information about an item.
     *
     * The loader of the module may need the pointers to the actual data within
     * an item. Instead of having to create multiple functions to get each
     * individual item, this function will get all of them.
     *
     * @param item the item to request information about
     * @param item_info
     * @return true if successful
     */
    virtual bool get_item_info(gsl::not_null<const ItemIface*> item,
                               gsl::not_null<item_info*> item_info) = 0;

    /**
     * The set of collections related functions. May be defined optionally by
     * the engine(s) that support them.
     */
    virtual cb::engine_errc set_collection_manifest(
            gsl::not_null<const void*> cookie, std::string_view json) {
        return cb::engine_errc::not_supported;
    }

    /**
     * Retrieve the last manifest set using set_manifest (a JSON document)
     */
    virtual cb::engine_errc get_collection_manifest(
            gsl::not_null<const void*> cookie, const AddResponseFn& response) {
        return cb::engine_errc::not_supported;
    }

    virtual cb::EngineErrorGetCollectionIDResult get_collection_id(
            gsl::not_null<const void*> cookie, std::string_view path) {
        return cb::EngineErrorGetCollectionIDResult{
                cb::engine_errc::not_supported};
    }

    virtual cb::EngineErrorGetScopeIDResult get_scope_id(
            gsl::not_null<const void*> cookie, std::string_view path) {
        return cb::EngineErrorGetScopeIDResult{cb::engine_errc::not_supported};
    }

    /**
     * Get the scope for the provided key
     *
     * @param cookie cookie object used to identify the request
     * @param key the key to look up
     * @return pair with the manifest UID and if found the scope where the key
     *              belongs.
     */
    virtual cb::EngineErrorGetScopeIDResult get_scope_id(
            gsl::not_null<const void*> cookie,
            const DocKey& key,
            std::optional<Vbid> vbid = std::nullopt) const {
        return cb::EngineErrorGetScopeIDResult(cb::engine_errc::not_supported);
    }

    /**
     * Ask the engine what features it supports.
     */
    virtual cb::engine::FeatureSet getFeatures() = 0;

    /**
     * @returns if XATTRs are enabled for this bucket
     */
    virtual bool isXattrEnabled() {
        return false;
    }

    /**
     * Get the "current" time and mode of the Hybrid Logical Clock for the
     * specified vBucket.
     * @returns seconds since unix epoch of the HLC, along with the current
     * HLC Mode (Real / Logical).
     */
    virtual cb::HlcTime getVBucketHlcNow(Vbid vbucket) = 0;

    /**
     * @returns the compression mode of the bucket
     */
    virtual BucketCompressionMode getCompressionMode() {
        return BucketCompressionMode::Off;
    }

    /**
     * @returns the maximum item size supported by the bucket
     */
    virtual size_t getMaxItemSize() {
        return default_max_item_size;
    };

    /**
     * @returns the minimum compression ratio defined in the bucket
     */
    virtual float getMinCompressionRatio() {
        return default_min_compression_ratio;
    }

    /**
     * Set a configuration parameter in the engine
     *
     * @param cookie The cookie identifying the request
     * @param category The parameter category
     * @param key The name of the parameter
     * @param value The value for the parameter
     * @param vbucket The vbucket specified in the request (only used for
     *                the vbucket sub group)
     * @return The standard engine codes
     */
    virtual cb::engine_errc setParameter(gsl::not_null<const void*> cookie,
                                         EngineParamCategory category,
                                         std::string_view key,
                                         std::string_view value,
                                         Vbid vbucket) {
        return cb::engine_errc::not_supported;
    }

    /**
     * Compact a database
     *
     * @param cookie The cookie identifying the request
     * @param vbid The vbucket to compact
     * @param purge_before_ts The timestamp to purge items before
     * @param purge_before_seq The sequence number to purge items before
     * @param drop_deletes Set to true if deletes should be dropped
     * @return The standard engine codes
     */
    virtual cb::engine_errc compactDatabase(gsl::not_null<const void*> cookie,
                                            Vbid vbid,
                                            uint64_t purge_before_ts,
                                            uint64_t purge_before_seq,
                                            bool drop_deletes) {
        return cb::engine_errc::not_supported;
    }

    /**
     * Get the state of a vbucket
     *
     * @param cookie The cookie identifying the request
     * @param vbid The vbucket to look up
     * @return A pair where the first entry is one of the standard engine codes
     *         and the second is the state when the status is success.
     */
    virtual std::pair<cb::engine_errc, vbucket_state_t> getVBucket(
            gsl::not_null<const void*> cookie, Vbid vbid) {
        return {cb::engine_errc::not_supported, vbucket_state_dead};
    }

    /**
     * Set the state of a vbucket
     *
     * @param cookie The cookie identifying the request
     * @param vbid The vbucket to update
     * @param cas The current value used for CAS swap
     * @param state The new vbucket state
     * @param meta The optional meta information for the state (nullptr if
     *             none is provided)
     * @return The standard engine codes
     */
    virtual cb::engine_errc setVBucket(gsl::not_null<const void*> cookie,
                                       Vbid vbid,
                                       uint64_t cas,
                                       vbucket_state_t state,
                                       nlohmann::json* meta) {
        return cb::engine_errc::not_supported;
    }

    /**
     * Delete a vbucket
     *
     * @param cookie The cookie identifying the request
     * @param vbid The vbucket to delete
     * @param sync If the operation should block (return EWB and be signalled
     *             when the operation is done) or delete in "fire and forget"
     *             mode.
     * @return The standard engine codes
     */
    virtual cb::engine_errc deleteVBucket(gsl::not_null<const void*> cookie,
                                          Vbid vbid,
                                          bool sync) {
        return cb::engine_errc::not_supported;
    }
};

namespace cb {
class ItemDeleter {
public:
    ItemDeleter() : handle(nullptr) {}

    /**
     * Create a new instance of the item deleter.
     *
     * @param handle_ the handle to the the engine who owns the item
     */
    explicit ItemDeleter(EngineIface* handle_) : handle(handle_) {
        if (handle == nullptr) {
            throw std::invalid_argument(
                "cb::ItemDeleter: engine handle cannot be nil");
        }
    }

    /**
     * Create a copy constructor to allow us to use std::move of the item
     */
    ItemDeleter(const ItemDeleter& other) = default;

    void operator()(ItemIface* item) {
        if (handle) {
            handle->release(item);
        } else {
            throw std::invalid_argument("cb::ItemDeleter: item attempted to be "
                                        "freed by null engine handle");
        }
    }

private:
    EngineIface* handle;
};

inline EngineErrorItemPair makeEngineErrorItemPair(cb::engine_errc err) {
    return {err, unique_item_ptr{nullptr, ItemDeleter{}}};
}

inline EngineErrorItemPair makeEngineErrorItemPair(cb::engine_errc err,
                                                   ItemIface* it,
                                                   EngineIface* handle) {
    return {err, unique_item_ptr{it, ItemDeleter{handle}}};
}
}

struct EngineDeletor {
    void operator()(EngineIface* engine) {
        engine->destroy(force);
    }

    bool force = false;
};

using unique_engine_ptr = std::unique_ptr<EngineIface, EngineDeletor>;

/**
 * @}
 */


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include <math.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_voxel_graphics.h"
#include "tz_neurotrace.h"
#include "tz_stack_math.h"
#include "tz_stack_utils.h"
#include "tz_objdetect.h"
#include "tz_voxeltrans.h"
#include "tz_stack_stat.h"
#include "tz_stack_draw.h"
#include "tz_geo3d_vector.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack_bwdist.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
#if 1
  Stack *stack = Read_Stack("../data/lobster_neuron_seed.tif");

  Voxel_List *list = Stack_To_Voxel_List(stack);

  //Print_Voxel_List(list);

  Pixel_Array *pa = Pixel_Array_Read("../data/lobster_neuron_seeds.pa");
  //Print_Pixel_Array(pa);

  Voxel *seed;

  int i;
  double *pa_array = (double *) pa->array;

  printf("mean: %g, std: %g\n", gsl_stats_mean(pa_array, 1, pa->size),
	 sqrt(gsl_stats_variance(pa_array, 1, pa->size)));

  double threshold = gsl_stats_mean(pa_array, 1, pa->size) +
    3.0 * sqrt(gsl_stats_variance(pa_array, 1, pa->size));

  double max_r = gsl_stats_max(pa_array, 1, pa->size);

  printf("%g, %g\n", threshold, max_r);

  max_r = 30.0;
  Set_Neuroseg_Max_Radius(max_r);

  dim_type dim[3];
  dim[0] = stack->width;
  dim[1] = stack->height;
  dim[2] = stack->depth;

  IMatrix *chord = Make_IMatrix(dim, 3);

  Stack *code = Make_Stack(GREY16, stack->width, stack->height, stack->depth);

  Kill_Stack(stack);

  stack = Read_Stack("../data/lobster_neuron_single.tif");

  Rgb_Color color;
  Set_Color(&color, 255, 0, 0);

  Stack *signal = Read_Stack("../data/lobster_neuron.tif");

  Stack *canvas = Translate_Stack(signal, COLOR, 0);

  /************** soma detection *************/
#if 0
  Struct_Element *se = Make_Ball_Se(((int) threshold));
  Stack *stack1 = Stack_Erode_Fast(stack, NULL, se);
  Stack *soma = Stack_Dilate(stack1, NULL, se);
  Kill_Stack(stack1);

  Stack_And(stack, soma, soma);

  Stack_Label_Bwc(canvas, soma, color);

  Kill_Stack(soma);

  Write_Stack("../data/test.tif", canvas);
  return 0;
#endif
  /*******************************************/

  Object_3d *obj = NULL;
  int seed_offset = -1;

  Neurochain *chain = NULL;

  double z_scale = 0.488 / 0.588;

  Set_Zscale(z_scale);

  Stack *traced =
    Make_Stack(GREY, signal->width, signal->height, signal->depth);
  One_Stack(traced);

  for (i = 0; i < pa->size; i++) {
    seed = Voxel_Queue_De(&list);
    printf("------------------------------------------> seed: %d\n", i);
    if (*STACK_PIXEL_8(traced, seed->x, seed->y, seed->z, 0) == 0) {
      continue;
    }

    double width = pa_array[i];

    if (width > max_r) {
      continue;
    }

    chain = New_Neurochain();

    Print_Voxel(seed);
    printf("%g\n", width);

    int max_level = (int) (width + 0.5);
    if (max_level < 6) {
      max_level = 6;
    }

    seed_offset = Stack_Util_Offset(seed->x, seed->y, seed->z,
				    stack->width, stack->height, stack->depth);
    Stack_Level_Code_Constraint(stack, code, chord->array, &seed_offset, 1,
				max_level + 1);

    Voxel_t v;
    Voxel_To_Tvoxel(seed, v);

    Print_Tvoxel(v);

    Stack *tmp_stack = Copy_Stack(stack);
    obj = Stack_Grow_Object_Constraint(tmp_stack, 1, v, chord, code,
				       max_level);
    Free_Stack(tmp_stack);

    Print_Object_3d_Info(obj);

    double vec[3];
    Object_3d_Orientation_Zscale(obj, vec, MAJOR_AXIS, z_scale);

    double theta, psi;
    Geo3d_Vector obj_vec;
    Set_Geo3d_Vector(&obj_vec, vec[0], vec[1], vec[2]);

    Geo3d_Vector_Orientation(&obj_vec, &theta, &psi);

    Set_Neuroseg(&(chain->seg), width, width, 10.0,
		 theta, psi);

    double cpos[3];
    cpos[0] = seed->x;
    cpos[1] = seed->y;
    cpos[2] = seed->z;
    cpos[2] *= z_scale;

    double bpos[3];
    Neuroseg_Pos_Center_To_Bottom(&(chain->seg), cpos, bpos);
    Set_Position(chain->position, bpos[0], bpos[1], bpos[2]);

    if (Initialize_Tracing(signal, chain, NULL, NULL) >= MIN_SCORE) {
      if ((chain->seg.r1 < max_r) && (chain->seg.r2 < max_r)) {
	chain = Trace_Neuron(signal, chain, BOTH, traced);

	Print_Neurochain(chain);

	//Stack_Draw_Object_Bwc(canvas, obj, color);

	Neurochain_Erase(traced, Neurochain_Head(chain));
	Neurochain_Label(canvas, Neurochain_Head(chain));
      }
    }

    Free_Neurochain(chain);

    free(seed);

    Kill_Object_3d(obj);
  }

  Write_Stack("../data/lobster_neuron_traced.tif", canvas);

#endif

  return 0;
}


text:
#include <petscksp.h>
#include "AFS.h"
#include "AFS_ctmc_petsc.h"
#include "im_clam.h"
#include "sparseExp.h"
#include "cs.h"
#include "AFS_ctmc.h"
#include "nlopt.h"
#include "adkGSL.h"
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>
//Functions for jointly computing CTMC transitions and embedded chain

//this function fills two petsc matrices: 1) embedded Discrete MC trans Mat, 2) CTMC trans mat
// also fills the rates vector

void fillPetscTransMats(afsStateSpace *S, double *topol, int *moveType, int *nzCount,
	int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2, Mat *DTMC, Mat *CTMC, gsl_vector *rates){
	int i, j,k, newnz, abcount=0;
	int N = S->nstates,abFlag;
	double  c0r[N],c1r[N],m0r[N],m1r[N],totR[N],rowSums[N], tmp;;

	gsl_vector_set_zero(rates);
	//set rates vector
	for(i=0;i<N;i++){
		c0r[i] = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1)  / theta1 ;
		if(theta2 == 0)
			c1r[i] = 0.0;
		else
			c1r[i] = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1)  / theta2;

		m0r[i] = S->states[i]->aCounts[0] * mig1 ;
		m1r[i] = S->states[i]->aCounts[1] * mig2 ;
		totR[i] = c0r[i] + c1r[i] + m0r[i] + m1r[i];

		rowSums[i] = 0.0;
		if(S->states[i]->nalleles > 1  && c0r[i] >= 0 && c1r[i] >= 0 && totR[i] != 0.0)
			gsl_vector_set(rates,i,1.0/totR[i]);
		else
			gsl_vector_set(rates,i,0);
	}

	newnz=0;
	for(k=0;k<*nzCount;k++){
		i = dim1[k];
		j= dim2[k];
//		printf("i: %d, j:%d, movetype: %d\n",i,j,moveType[k]);
		if(S->states[i]->nalleles ==1 || S->states[j]->nalleles ==1 )abFlag=1;
	//	if(S->states[i]->nalleles ==1 )abFlag=1;
		else abFlag = 0;
		switch(moveType[k]){
			case 0: //coal pop0
			if(totR[i]==0)tmp=0;
			else tmp = c0r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)MatSetValue(*DTMC,i,j,tmp,INSERT_VALUES);
			break;
			case 1: //coal pop1
			if(totR[i]==0)tmp=0;
			else tmp = c1r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)MatSetValue(*DTMC,i,j,tmp,INSERT_VALUES);
			break;
			case 2: //mig pop0
			if(totR[i]==0)tmp=0;
			else tmp = m0r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)MatSetValue(*DTMC,i,j,tmp,INSERT_VALUES);
			break;
			case 3: //mig pop1
			if(totR[i]==0)tmp=0;
			else tmp = m1r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)MatSetValue(*DTMC,i,j,tmp,INSERT_VALUES);
			break;
			case -1:
			abcount+=1;
			break;
		}
	}
	newnz=*nzCount ;
	//now add diagonal elements; topol pre-set to -1
	for(k=0;k<N;k++){
		if((rowSums[k]) != 0.0){
			MatSetValue(*CTMC,k,k,0.0-rowSums[k],INSERT_VALUES);
			newnz++;
		}
		MatSetValue(*DTMC,k,k,0.0,INSERT_VALUES);
	}
	*nzCount = newnz;
	/* Assemble the matrix */
  	MatAssemblyBegin(*CTMC,MAT_FINAL_ASSEMBLY);
  	MatAssemblyEnd(*CTMC,MAT_FINAL_ASSEMBLY);
  	MatAssemblyBegin(*DTMC,MAT_FINAL_ASSEMBLY);
  	MatAssemblyEnd(*DTMC,MAT_FINAL_ASSEMBLY);
}

//fillPetscCsparseTransMats--this function fills one petsc matrix with CTMC transmat and returns pointer to CSparse matrix for DTMC trans mat
// also fills the rates vector
//note that the DTMC matrix needs to be cleared before entry to this function and compressed afterward.

struct cs_di_sparse * fillPetscCsparseTransMats(afsStateSpace *S, double *topol, int *moveType, int *nzCount,
	int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2, Mat *CTMC, gsl_vector *rates){
	int i, j,k, newnz, abcount=0;
	int N = S->nstates,abFlag;
	double  c0r[N],c1r[N],m0r[N],m1r[N],totR[N],rowSums[N], tmp;;
	struct cs_di_sparse *triplet, *tmat;

	//allocate new cs_sparse obj
	triplet = cs_spalloc(N, N, *nzCount + N , 1, 1); //alloc sparse mat with extra space for nonzero identity mats

	gsl_vector_set_zero(rates);
	//set rates vector
	for(i=0;i<N;i++){
		c0r[i] = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1)  / theta1 ;
		if(theta2 == 0)
			c1r[i] = 0.0;
		else
			c1r[i] = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1)  / theta2;
		m0r[i] = S->states[i]->aCounts[0] * mig1 ;
		m1r[i] = S->states[i]->aCounts[1] * mig2 ;
		totR[i] = c0r[i] + c1r[i] + m0r[i] + m1r[i];

		rowSums[i] = 0.0;
		if(S->states[i]->nalleles > 1  && c0r[i] >= 0 && c1r[i] >= 0 && totR[i] != 0.0)
			gsl_vector_set(rates,i,1.0/totR[i]);
		else
			gsl_vector_set(rates,i,0);
	}

	newnz=0;
	for(k=0;k<*nzCount;k++){
		i = dim1[k];
		j= dim2[k];
//		printf("i: %d, j:%d, movetype: %d\n",i,j,moveType[k]);
		if(S->states[i]->nalleles ==1 || S->states[j]->nalleles ==1 )abFlag=1;
	//	if(S->states[i]->nalleles ==1 )abFlag=1;
		else abFlag = 0;
		switch(moveType[k]){
			case 0: //coal pop0
			if(totR[i]==0)tmp=0;
			else tmp = c0r[i] / totR[i] * topol[k];
		//	printf("here--- i:%d j:%d tmp*totR[i]:%f",i,j,tmp*totR[i]);
			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 1: //coal pop1
			if(totR[i]==0)tmp=0;
			else tmp = c1r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 2: //mig pop0
			if(totR[i]==0)tmp=0;
			else tmp = m0r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 3: //mig pop1
			if(totR[i]==0)tmp=0;
			else tmp = m1r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case -1:
			abcount+=1;
			break;
		}
	}
	newnz=*nzCount ;
	//now add diagonal elements; topol pre-set to -1
	for(k=0;k<N;k++){
		if((rowSums[k]) != 0.0){
			MatSetValue(*CTMC,k,k,0.0-rowSums[k],INSERT_VALUES);
			newnz++;
		}
	}
	*nzCount = newnz;
	/* Assemble the matrix */
  	MatAssemblyBegin(*CTMC,MAT_FINAL_ASSEMBLY);
  	MatAssemblyEnd(*CTMC,MAT_FINAL_ASSEMBLY);
	tmat = cs_compress(triplet);
	cs_spfree(triplet);
	return(tmat);
}

//fillPetscCsparseTransMats_prealloc--this function fills one petsc matrix with CTMC transmat and returns pointer to CSparse matrix for DTMC trans mat
// also fills the rates vector
struct cs_di_sparse * fillPetscCsparseTransMats_prealloc(afsStateSpace *S, double *topol, int *moveType, int *nzCount,
	int *dim1, int *dim2, double theta1, double theta2, double mig1, double mig2,
	Mat *CTMC, gsl_vector *rates, struct cs_di_sparse *triplet){

	int i, j,k, newnz, abcount=0;
	int N = S->nstates,abFlag;
	double  *c0r,*c1r,*m0r,*m1r,*totR,*rowSums, tmp;
	struct cs_di_sparse *tmat;

	c0r = malloc(sizeof(double)*N);
	c1r = malloc(sizeof(double)*N);
	m0r = malloc(sizeof(double)*N);
	m1r = malloc(sizeof(double)*N);
	totR = malloc(sizeof(double)*N);
	rowSums = malloc(sizeof(double)*N);

	gsl_vector_set_zero(rates);
	//set rates vector
	for(i=0;i<N;i++){
		c0r[i] = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1)  / theta1 ;
		if(theta2 == 0)
			c1r[i] = 0.0;
		else
			c1r[i] = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1)  / theta2;
		m0r[i] = S->states[i]->aCounts[0] * mig1 ;
		m1r[i] = S->states[i]->aCounts[1] * mig2 ;
		totR[i] = c0r[i] + c1r[i] + m0r[i] + m1r[i];

		rowSums[i] = 0.0;
		if(S->states[i]->nalleles > 1  && c0r[i] >= 0 && c1r[i] >= 0 && totR[i] != 0.0)
			gsl_vector_set(rates,i,1.0/totR[i]);
		else
			gsl_vector_set(rates,i,0);
	}

	newnz=0;
	for(k=0;k<*nzCount;k++){
		i = dim1[k];
		j= dim2[k];
//		printf("i: %d, j:%d, movetype: %d\n",i,j,moveType[k]);
		if(S->states[i]->nalleles ==1 || S->states[j]->nalleles ==1 )abFlag=1;
	//	if(S->states[i]->nalleles ==1 )abFlag=1;
		else abFlag = 0;
		switch(moveType[k]){
			case 0: //coal pop0
			if(totR[i]==0)tmp=0;
			else tmp = c0r[i] / totR[i] * topol[k];
		//	printf("here--- i:%d j:%d tmp*totR[i]:%f",i,j,tmp*totR[i]);
			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 1: //coal pop1
			if(totR[i]==0)tmp=0;
			else tmp = c1r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 2: //mig pop0
			if(totR[i]==0)tmp=0;
			else tmp = m0r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case 3: //mig pop1
			if(totR[i]==0)tmp=0;
			else tmp = m1r[i] / totR[i] * topol[k];

			MatSetValue(*CTMC,i,j,tmp*totR[i],INSERT_VALUES);
			rowSums[i] += tmp* totR[i];
			if(abFlag==0)cs_entry(triplet,i,j,tmp);
			break;
			case -1:
			abcount+=1;
			break;
		}
	}
	newnz=*nzCount ;
	//now add diagonal elements; topol pre-set to -1
	for(k=0;k<N;k++){
		if((rowSums[k]) != 0.0){
			MatSetValue(*CTMC,k,k,0.0-rowSums[k],INSERT_VALUES);
			newnz++;
		}
	}
	*nzCount = newnz;
	/* Assemble the matrix */
  	MatAssemblyBegin(*CTMC,MAT_FINAL_ASSEMBLY);
  	MatAssemblyEnd(*CTMC,MAT_FINAL_ASSEMBLY);
	tmat = cs_compress(triplet);

	free(c0r);
	free(c1r);
	free(m0r);
	free(m1r);
	free(totR);
	free(rowSums);
	return(tmat);
}

//calcLogAFS_IM -- returns the expects log AFS from an IM model
void calcLogAFS_IM(void * p){
	struct clam_lik_params * params = (struct clam_lik_params *) p;
 	PetscInt i,j, N = params->stateSpace->nstates;
  	cs *spMat, *mt, *tmpMat ;
  	double timeV, sum, thetaA, theta2, mig1, mig2;
  	PetscInt Na;
  	gsl_vector *tmpStates;
 	css *S ;
	csn *NN ;
	PetscInt n,Ntmp ;
	double *xx;
 	PetscInt        iStart,iEnd, *idx;
	const PetscInt *idx2;
	PetscScalar    *tmpArray;
	const PetscScalar *tmpArrayC;
	PetscErrorCode ierr;
	MFNConvergedReason reason;

 	//initialize some vectors
 	gsl_vector_set_zero(params->rates);
	MatZeroEntries(params->C);
	MatZeroEntries(params->C2);
 	tmpStates = gsl_vector_alloc(N);
	Ntmp=N;
	Na=params->Na;
 	//For straight MLE the paramVector takes the form [N2,NA,m1,m2,t]
 	theta2 = gsl_vector_get(params->paramVector,0);
 	thetaA = gsl_vector_get(params->paramVector,1);
 	mig1 = gsl_vector_get(params->paramVector,2);
 	mig2 = gsl_vector_get(params->paramVector,3);
 	timeV = gsl_vector_get(params->paramVector,4);
// 	printf("params-> %f %f %f %f %f\n",theta2,thetaA,mig1,mig2,timeV);
// 	printf("nnz %d nnzA%d \n",params->nnz,params->nnzA);

	//fill transMat
 //	tmpMat = fillPetscCsparseTransMats_prealloc(params->stateSpace, params->top, params->move, &params->nnz,
 //		params->dim1, params->dim2, 1, theta2, mig1, mig2,  &params->C, params->rates,params->triplet);

	tmpMat = fillPetscCsparseTransMats(params->stateSpace, params->top, params->move, &params->nnz,
	 		params->dim1, params->dim2, 1, theta2, mig1, mig2,  &params->C, params->rates);
	//using CSparse
	// S = (I-P)^-1
	//add negative ident
//	ident = cs_spalloc(N,N,N,1,1);
//	for(i=0;i<N;i++) cs_entry(ident,i,i,1);
//	eye = cs_compress(ident);
	spMat = cs_add(params->eye,tmpMat,1.0,-1.0);
	//cs_print_adk(spMat);
	mt = cs_transpose(spMat,1);
	//cs_print_adk(mt);

	n = mt->n ;
	S = cs_sqr (0, mt, 0) ;              /* ordering and symbolic analysis */
	NN = cs_lu (mt, S, 1e-12) ;                 /* numeric LU factorization */
	xx = cs_malloc (n, sizeof (double)) ;    /* get workspace */

	MatGetOwnershipRange(params->denseMat1,&iStart,&iEnd);
	PetscMalloc1(N,&idx);
	for(j=0;j<N;j++) idx[j]=j;

	////Compute Entire Inverse Mat
	for(j=iStart;j<iEnd;j++){
	//create unit array for solve
		for(i=0; i<N; i++)params->b[i] = 0.0;
		params->b[j]=1.0;
		//factor outside loop leads to ~40x speedup
		cs_ipvec (NN->pinv, params->b, xx, n) ;       /* x = b(p) */
		cs_lsolve (NN->L, xx) ;               /* x = L\x */
		cs_usolve (NN->U, xx) ;               /* x = U\x */
		cs_ipvec (S->q, xx, params->b, n) ;          /* b(q) = x */
	//	if(!rank)printf("row %d\n",j);
		MatSetValues(params->denseMat1,1,&j,N,idx,params->b,INSERT_VALUES);
	}
	MatAssemblyBegin(params->denseMat1,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(params->denseMat1,MAT_FINAL_ASSEMBLY);
	//temporary clean up
	free(xx);
	cs_spfree(spMat);
	cs_spfree(mt);
//	cs_spfree(ident);
//	cs_spfree(eye);
	cs_spfree(tmpMat);
	cs_nfree(NN);
	cs_sfree(S);
	PetscFree(idx);

//	MatView(params->denseMat1,PETSC_VIEWER_STDOUT_WORLD);

	//broadcast the invMat[0,] as vector to each processor
	if(params->rank==0){
		MatGetRow(params->denseMat1,0,&N,&idx2,&tmpArrayC);
		for (j = 0; j < N; j++) params->b[j] = tmpArrayC[j];
		MatRestoreRow(params->denseMat1,0,&Ntmp,&idx2,&tmpArrayC);	//have to use Ntmp here as MatRestoreRow does something funky
	}
	MPI_Bcast(params->b,N,MPI_DOUBLE,0,PETSC_COMM_WORLD);

	//Get Island Time 0-INF Unnormal
	gsl_matrix_set_zero(params->expAFS);
	fillExpectedAFS_unnorm(params->stateSpace, params->b,params->rates,params->expAFS);
	// printf("////////////////////Island Time 0 - INF unnormalized\n");
	// for(i=0;i< (params->n1+1) ;i++){
	// 	for(j=0;j< (params->n2+1);j++){
	// 		printf("%.5f ",gsl_matrix_get(params->expAFS,i,j));
	// 	}
	// 	printf("\n");
	// }

	//Matrix Exponentiation to get state vector at time t
	//SLEPC Stuff ahead!
	/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	Create the solver and set various options
	- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

	/*
	Set operator matrix, the function to compute, and other options
	*/
	//transpose
	MatTranspose(params->C,MAT_REUSE_MATRIX,&params->C);
	ierr = MFNSetOperator(params->mfn,params->C);CHKERRV(ierr);
	/*
	Set solver parameters at runtime
	*/

	ierr = MFNSetScaleFactor(params->mfn,timeV); CHKERRV(ierr);
	ierr = MFNSetFromOptions(params->mfn);CHKERRV(ierr);

	/* set v = e_1 */
	ierr = MatGetVecs(params->C,PETSC_NULL,&params->y);CHKERRV(ierr);
	ierr = MatGetVecs(params->C,&params->v, PETSC_NULL);CHKERRV(ierr);
	ierr = VecSetValue(params->v,0,1.0,INSERT_VALUES);CHKERRV(ierr);

	ierr = VecAssemblyBegin(params->v);CHKERRV(ierr);
	ierr = VecAssemblyEnd(params->v);CHKERRV(ierr);
	ierr = MFNSolve(params->mfn,params->v,params->y);CHKERRV(ierr);
	ierr = MFNGetConvergedReason(params->mfn,&reason);CHKERRV(ierr);
	if (reason!=MFN_CONVERGED_TOL) ierr = 1; CHKERRV(ierr);

	//retranspose for later use
	MatTranspose(params->C,MAT_REUSE_MATRIX,&params->C);

	//state vector at time t stored in y
	//scatter y and store in expoArray
	VecScatterCreateToAll(params->y,&params->ctx,&params->v_seq);
	VecScatterBegin(params->ctx,params->y,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecScatterEnd(params->ctx,params->y,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecGetArray(params->v_seq,&tmpArray);
	for (j = 0; j < N; j++){
		params->expoArray[j] = tmpArray[j];
		//printf("%.5f\n",expoArray[j]);
	}

	VecRestoreArray(params->v_seq,&tmpArray);
	VecScatterDestroy(&params->ctx);
	VecDestroy(&params->v_seq);
	VecDuplicate(params->y,&params->x);

	//VecView(y, PETSC_VIEWER_STDOUT_WORLD);
	MatMultTranspose(params->denseMat1,params->y,params->x);
	VecScatterCreateToAll(params->x,&params->ctx,&params->v_seq);
	VecScatterBegin(params->ctx,params->x,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecScatterEnd(params->ctx,params->x,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecGetArray(params->v_seq,&tmpArray);
	for (j = 0; j < N; j++){
		params->b[j] = tmpArray[j];
		//printf("rank %d expoArray[%d]:%.5f\n",rank,j,expoArray[j]);
		//printf("rank %d b[%d]:%.5f\n",rank,j,b[j]);
	}

	VecRestoreArray(params->v_seq,&tmpArray);
	VecScatterDestroy(&params->ctx);
	VecDestroy(&params->v_seq);

	fillExpectedAFS_unnorm(params->stateSpace, params->b,params->rates,params->expAFS2);
	//now subtract older AFS (t_t - t_INF) from total AFS (t_0 - t_INF)
	/////////////////
	///////

	gsl_matrix_sub(params->expAFS,params->expAFS2);

	//////////////////////
	//Final phase, collapse popns, reset rates
	//now collapse populations
	//state vector at time t stored in st[] for largerStateSpace; use reverseMap
	for(i=0;i<N;i++)params->st[i]=0.0;
	for(i=0;i<N;i++){params->st[params->map[i]]+=params->expoArray[i];
	//	printf("expoArray[%d]:%f\n",i,expoArray[i]);
	}
	for(i=0; i<params->reducedStateSpace->nstates; i++){
		VecSetValue(params->ancStateVec,i,params->st[params->reverseMap[i]],INSERT_VALUES );
	}

	//Re-Fill trans mats for Ancestral Params

	params->nnzA = coalMarkovChainTopologyMatrix_sparse(params->reducedStateSpace,params->topA,params->moveA, params->dim1A, params->dim2A);

	tmpMat=fillPetscCsparseTransMats(params->reducedStateSpace, params->topA, params->moveA, &params->nnzA,
		params->dim1A, params->dim2A, thetaA, 0, 0, 0, &params->C2, params->rates);
	//////////////////////
	//
	// S = (I-P)^-1
	//
	/////////
	//////
	//add negative ident
//	ident = cs_spalloc(Na,Na,Na,1,1);
//	for(i=0;i<Na;i++) cs_entry(ident,i,i,1);
//	eye = cs_compress(ident);
	spMat = cs_add(params->eyeAnc,tmpMat,1.0,-1.0);

	//cs_print_adk(spMat);
	mt = cs_transpose(spMat,1);
	//VecView(y,PETSC_VIEWER_STDOUT_WORLD);
	n = mt->n ;
	S = cs_sqr (0, mt, 0) ;              /* ordering and symbolic analysis */
	NN = cs_lu (mt, S, 1e-12) ;                 /* numeric LU factorization */
	xx = cs_malloc (n, sizeof (double)) ;    /* get workspace */
	MatGetOwnershipRange(params->denseMat2,&iStart,&iEnd);
//	printf("here\n");
	PetscFree(idx);
	PetscMalloc1(Na,&idx);
	for(j=0;j<Na;j++) idx[j]=j;
	////Compute Entire Inverse Mat
	for(j=iStart;j<iEnd;j++){
	//create unit array for solve
		for(i=0; i<Na; i++)params->b[i] = 0.0;
		params->b[j]=1.0;
		//factor outside loop leads to ~40x speedup
		cs_ipvec (NN->pinv, params->b, xx, n) ;       /* x = b(p) */
		cs_lsolve (NN->L, xx) ;               /* x = L\x */
		cs_usolve (NN->U, xx) ;               /* x = U\x */
		cs_ipvec (S->q, xx, params->b, n) ;          /* b(q) = x */
	//	if(!rank)printf("row %d\n",j);
		MatSetValues(params->denseMat2,1,&j,Na,idx,params->b,INSERT_VALUES);
	}
	MatAssemblyBegin(params->denseMat2,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(params->denseMat2,MAT_FINAL_ASSEMBLY);
	free(xx);

	MatMultTranspose(params->denseMat2,params->ancStateVec,params->ancResVec);
	VecScatterCreateToAll(params->ancResVec,&params->ctx,&params->v_seq);
	VecScatterBegin(params->ctx,params->ancResVec,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecScatterEnd(params->ctx,params->ancResVec,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);

	VecGetArray(params->v_seq,&tmpArray);
	for (j = 0; j < Na; j++){
		params->expoArray[j] = tmpArray[j];
		//printf("%.5f\n",expoArray[j]);
	}

	VecRestoreArray(params->v_seq,&tmpArray);
	VecScatterDestroy(&params->ctx);
	VecDestroy(&params->v_seq);

	//get contribution starting at st
	gsl_matrix_set_zero(params->expAFS2);
	fillExpectedAFS_unnorm(params->reducedStateSpace, params->expoArray,params->rates,params->expAFS2);

	// printf("////////////////////--Ancestral bit\n");
	// for(i=0;i< (n1+1) ;i++){
	// 	for(j=0;j< (n2+1);j++){
	// 		printf("%.5f ",gsl_matrix_get(expAFS2,i,j));
	// 	}
	// 	printf("\n");
	// }

	//if(params->rank==0)printf("////////////////////normalized IM:\n");

	gsl_matrix_add(params->expAFS,params->expAFS2);
	sum = matrixSumDouble(params->expAFS);
	gsl_matrix_scale(params->expAFS, 1.0 / sum);
	params->meanTreeLength = sum;
	// for(i=0;i< (params->n1+1) ;i++){
	// 		for(j=0;j< (params->n2+1);j++){
	// 			if(params->rank==0)printf("%.5f ",gsl_matrix_get(params->expAFS,i,j));
	// 		}
	// 		if(params->rank==0)printf("\n");
	// 	}
	//

	//clean up
	PetscFree(idx);
	VecDestroy(&params->y);
	VecDestroy(&params->x);
	VecDestroy(&params->v);
	PetscFree(tmpArray);
	cs_spfree(spMat);
	cs_spfree(mt);
//	cs_spfree(ident);
//	cs_spfree(eye);
	cs_spfree(tmpMat);
	gsl_vector_free(tmpStates);
	cs_nfree(NN);
	cs_sfree(S);

}

//calcLogAFS_IM_allPETSC -- returns the expects log AFS from an IM model; only uses PETSC
void calcLogAFS_IM_allPETSC(void * p){
	struct clam_lik_params * params = (struct clam_lik_params *) p;
	PetscInt i,j, N = params->stateSpace->nstates;
	cs *spMat, *mt, *ident, *eye, *tmpMat ;
	double timeV, sum, thetaA, theta2, mig1, mig2;
	PetscInt Na;
	gsl_vector *tmpStates;
	css *S ;
	csn *NN ;
	PetscInt n,Ntmp ;
	double *xx;
	PetscInt        iStart,iEnd, *idx;
	const PetscInt *idx2;
	PetscScalar    *tmpArray;
	PetscScalar negOne = -1.0;
//	PetscScalar one = 1.0;

	const PetscScalar *tmpArrayC;
	PetscErrorCode ierr;
	IS perm,iperm;
	MatFactorInfo  info;
	MFNConvergedReason reason;

	//initialize some vectors
	gsl_vector_set_zero(params->rates);
	MatZeroEntries(params->C);
	MatZeroEntries(params->C2);
	MatZeroEntries(params->D);
	tmpStates = gsl_vector_alloc(N);
	Ntmp=N;
	Na=params->Na;

	//For straight MLE the paramVector takes the form [N2,NA,m1,m2,t]
	theta2 = gsl_vector_get(params->paramVector,0);
	thetaA = gsl_vector_get(params->paramVector,1);
	mig1 = gsl_vector_get(params->paramVector,2);
	mig2 = gsl_vector_get(params->paramVector,3);
	timeV = gsl_vector_get(params->paramVector,4);
// 	printf("params-> %f %f %f %f %f\n",theta2,thetaA,mig1,mig2,timeV);
// 	printf("nnz %d nnzA%d \n",params->nnz,params->nnzA);

	//fill transMat
	fillPetscTransMats(params->stateSpace, params->top, params->move, &params->nnz,
		params->dim1, params->dim2, 1, theta2, mig1, mig2, &params->D, &params->C, params->rates);

	//using CSparse
	// S = (I-P)^-1
	//subtract DTMC mat from identity

	MatZeroEntries(params->D_copy);
	MatCopy(params->ident,params->D_copy,DIFFERENT_NONZERO_PATTERN);
	MatAXPY(params->D_copy,negOne,params->D,DIFFERENT_NONZERO_PATTERN);
//	MatTranspose(params->D_copy,MAT_REUSE_MATRIX,&params->D);

	ierr = MatGetOrdering(params->D_copy,  MATORDERINGNATURAL,  &perm,  &iperm);

	ierr = MatFactorInfoInitialize(&info);
	ierr = MatGetFactor(params->D_copy,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&params->F);
	PetscInt icntl_7 = 5;
	//ierr = MatMumpsSetIcntl(params->F,7,icntl_7);
	info.fill = 5.0;
	ierr = MatLUFactorSymbolic(params->F,params->D_copy,perm,iperm,&info);
	ierr = MatLUFactorNumeric(params->F,params->D_copy,&info);

	//ierr = MatLUFactor(params->D_copy, perm, iperm, &info);

	////Compute Entire Inverse Mat
	ierr = MatMatSolve(params->F,params->denseIdent,params->denseMat1);

	//MatView(params->denseMat1,PETSC_VIEWER_STDOUT_SELF);

//	MatGetOwnershipRange(params->denseMat1,&iStart,&iEnd);
//	PetscMalloc1(N,&idx);
//	for(j=0;j<N;j++) idx[j]=j;

// 	for(j=0;j<N;j++){
// 	//create unit array for solve
// 		VecZeroEntries(params->bInv);
// 		VecZeroEntries(params->xInv);
// 		VecSetValue(params->bInv,j,one,INSERT_VALUES);
// 		VecAssemblyBegin(params->bInv);
// 		VecAssemblyEnd(params->bInv);
//
// 		KSPSolve(params->ksp,params->bInv,params->xInv);
// 		VecGetValues(params->x,N,idx,hold);
// 		MatSetValues(params->denseMat1,1,&j,N,idx,hold,INSERT_VALUES);
// 	}
// 	MatAssemblyBegin(params->denseMat1,MAT_FINAL_ASSEMBLY);
// 	MatAssemblyEnd(params->denseMat1,MAT_FINAL_ASSEMBLY);
// 	//temporary clean up
// //	free(xx);
// //	cs_spfree(spMat);
// //	cs_spfree(mt);
// //	cs_spfree(ident);
// //	cs_spfree(eye);
// //	cs_spfree(tmpMat);
// //	cs_nfree(NN);
// //	cs_sfree(S);
// 	PetscFree(idx);
// 	KSPDestroy(&params->ksp);
//	MatView(params->denseMat1,PETSC_VIEWER_STDOUT_WORLD);

	MatDestroy(&params->F);
	//broadcast the invMat[0,] as vector to each processor
	if(params->rank==0){
		MatGetRow(params->denseMat1,0,&N,&idx2,&tmpArrayC);
		for (j = 0; j < N; j++) params->b[j] = tmpArrayC[j];
		MatRestoreRow(params->denseMat1,0,&Ntmp,&idx2,&tmpArrayC);	//have to use Ntmp here as MatRestoreRow does something funky
	}
	MPI_Bcast(params->b,N,MPI_DOUBLE,0,PETSC_COMM_WORLD);

	//Get Island Time 0-INF Unnormal
	gsl_matrix_set_zero(params->expAFS);
	fillExpectedAFS_unnorm(params->stateSpace, params->b,params->rates,params->expAFS);
	// printf("////////////////////Island Time 0 - INF unnormalized\n");
	// for(i=0;i< (params->n1+1) ;i++){
	// 	for(j=0;j< (params->n2+1);j++){
	// 		printf("%.5f ",gsl_matrix_get(params->expAFS,i,j));
	// 	}
	// 	printf("\n");
	// }

	//Matrix Exponentiation to get state vector at time t
	//SLEPC Stuff ahead!
	/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	Create the solver and set various options
	- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

	/*
	Set operator matrix, the function to compute, and other options
	*/
	//transpose
	MatTranspose(params->C,MAT_REUSE_MATRIX,&params->C);
	ierr = MFNSetOperator(params->mfn,params->C);CHKERRV(ierr);
	/*
	Set solver parameters at runtime
	*/

	ierr = MFNSetScaleFactor(params->mfn,timeV); CHKERRV(ierr);
	ierr = MFNSetFromOptions(params->mfn);CHKERRV(ierr);

	/* set v = e_1 */
	ierr = MatGetVecs(params->C,PETSC_NULL,&params->y);CHKERRV(ierr);
	ierr = MatGetVecs(params->C,&params->v, PETSC_NULL);CHKERRV(ierr);
	ierr = VecSetValue(params->v,0,1.0,INSERT_VALUES);CHKERRV(ierr);

	ierr = VecAssemblyBegin(params->v);CHKERRV(ierr);
	ierr = VecAssemblyEnd(params->v);CHKERRV(ierr);
	ierr = MFNSolve(params->mfn,params->v,params->y);CHKERRV(ierr);
	ierr = MFNGetConvergedReason(params->mfn,&reason);CHKERRV(ierr);
	if (reason!=MFN_CONVERGED_TOL) ierr = 1; CHKERRV(ierr);

	//retranspose for later use
	MatTranspose(params->C,MAT_REUSE_MATRIX,&params->C);

	//state vector at time t stored in y
	//scatter y and store in expoArray
	VecScatterCreateToAll(params->y,&params->ctx,&params->v_seq);
	VecScatterBegin(params->ctx,params->y,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecScatterEnd(params->ctx,params->y,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecGetArray(params->v_seq,&tmpArray);
	for (j = 0; j < N; j++){
		params->expoArray[j] = tmpArray[j];
		//printf("%.5f\n",expoArray[j]);
	}

	VecRestoreArray(params->v_seq,&tmpArray);
	VecScatterDestroy(&params->ctx);
	VecDestroy(&params->v_seq);
	VecDuplicate(params->y,&params->x);

	//VecView(y, PETSC_VIEWER_STDOUT_WORLD);
	MatMultTranspose(params->denseMat1,params->y,params->x);
	VecScatterCreateToAll(params->x,&params->ctx,&params->v_seq);
	VecScatterBegin(params->ctx,params->x,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecScatterEnd(params->ctx,params->x,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecGetArray(params->v_seq,&tmpArray);
	for (j = 0; j < N; j++){
		params->b[j] = tmpArray[j];
		//printf("rank %d expoArray[%d]:%.5f\n",rank,j,expoArray[j]);
		//printf("rank %d b[%d]:%.5f\n",rank,j,b[j]);
	}

	VecRestoreArray(params->v_seq,&tmpArray);
	VecScatterDestroy(&params->ctx);
	VecDestroy(&params->v_seq);

	fillExpectedAFS_unnorm(params->stateSpace, params->b,params->rates,params->expAFS2);
	//now subtract older AFS (t_t - t_INF) from total AFS (t_0 - t_INF)
	/////////////////
	///////

	gsl_matrix_sub(params->expAFS,params->expAFS2);

	//////////////////////
	//Final phase, collapse popns, reset rates
	//now collapse populations
	//state vector at time t stored in st[] for largerStateSpace; use reverseMap
	for(i=0;i<N;i++)params->st[i]=0.0;
	for(i=0;i<N;i++){params->st[params->map[i]]+=params->expoArray[i];
	//	printf("expoArray[%d]:%f\n",i,expoArray[i]);
	}
	for(i=0; i<params->reducedStateSpace->nstates; i++){
		VecSetValue(params->ancStateVec,i,params->st[params->reverseMap[i]],INSERT_VALUES );
	}

	//Re-Fill trans mats for Ancestral Params

	params->nnzA = coalMarkovChainTopologyMatrix_sparse(params->reducedStateSpace,params->topA,params->moveA, params->dim1A, params->dim2A);

	tmpMat=fillPetscCsparseTransMats(params->reducedStateSpace, params->topA, params->moveA, &params->nnzA,
		params->dim1A, params->dim2A, thetaA, 0, 0, 0, &params->C2, params->rates);
	//////////////////////
	//
	// S = (I-P)^-1
	//
	/////////
	//////
	//add negative ident
	ident = cs_spalloc(Na,Na,Na,1,1);
	for(i=0;i<Na;i++) cs_entry(ident,i,i,1);
	eye = cs_compress(ident);
	spMat = cs_add(eye,tmpMat,1.0,-1.0);

	//cs_print_adk(spMat);
	mt = cs_transpose(spMat,1);
	//VecView(y,PETSC_VIEWER_STDOUT_WORLD);
	n = mt->n ;
	S = cs_sqr (0, mt, 0) ;              /* ordering and symbolic analysis */
	NN = cs_lu (mt, S, 1e-12) ;                 /* numeric LU factorization */
	xx = cs_malloc (n, sizeof (double)) ;    /* get workspace */
	MatGetOwnershipRange(params->denseMat2,&iStart,&iEnd);
//	printf("here\n");
//	PetscFree(idx);
	PetscMalloc1(Na,&idx);
	for(j=0;j<Na;j++) idx[j]=j;
	////Compute Entire Inverse Mat
	for(j=iStart;j<iEnd;j++){
	//create unit array for solve
		for(i=0; i<Na; i++)params->b[i] = 0.0;
		params->b[j]=1.0;
		//factor outside loop leads to ~40x speedup
		cs_ipvec (NN->pinv, params->b, xx, n) ;       /* x = b(p) */
		cs_lsolve (NN->L, xx) ;               /* x = L\x */
		cs_usolve (NN->U, xx) ;               /* x = U\x */
		cs_ipvec (S->q, xx, params->b, n) ;          /* b(q) = x */
	//	if(!rank)printf("row %d\n",j);
		MatSetValues(params->denseMat2,1,&j,Na,idx,params->b,INSERT_VALUES);
	}
	MatAssemblyBegin(params->denseMat2,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(params->denseMat2,MAT_FINAL_ASSEMBLY);
	free(xx);

	MatMultTranspose(params->denseMat2,params->ancStateVec,params->ancResVec);
	VecScatterCreateToAll(params->ancResVec,&params->ctx,&params->v_seq);
	VecScatterBegin(params->ctx,params->ancResVec,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);
	VecScatterEnd(params->ctx,params->ancResVec,params->v_seq,INSERT_VALUES,SCATTER_FORWARD);

	VecGetArray(params->v_seq,&tmpArray);
	for (j = 0; j < Na; j++){
		params->expoArray[j] = tmpArray[j];
		//printf("%.5f\n",expoArray[j]);
	}

	VecRestoreArray(params->v_seq,&tmpArray);
	VecScatterDestroy(&params->ctx);
	VecDestroy(&params->v_seq);

	//get contribution starting at st
	gsl_matrix_set_zero(params->expAFS2);
	fillExpectedAFS_unnorm(params->reducedStateSpace, params->expoArray,params->rates,params->expAFS2);

	// printf("////////////////////--Ancestral bit\n");
	// for(i=0;i< (n1+1) ;i++){
	// 	for(j=0;j< (n2+1);j++){
	// 		printf("%.5f ",gsl_matrix_get(expAFS2,i,j));
	// 	}
	// 	printf("\n");
	// }

	//if(params->rank==0)printf("////////////////////normalized IM:\n");

	gsl_matrix_add(params->expAFS,params->expAFS2);
	sum = matrixSumDouble(params->expAFS);
	gsl_matrix_scale(params->expAFS, 1.0 / sum);
	// for(i=0;i< (params->n1+1) ;i++){
	// 		for(j=0;j< (params->n2+1);j++){
	// 			if(params->rank==0)printf("%.5f ",gsl_matrix_get(params->expAFS,i,j));
	// 		}
	// 		if(params->rank==0)printf("\n");
	// 	}
	//

	//clean up
	PetscFree(idx);
	VecDestroy(&params->y);
	VecDestroy(&params->x);
	VecDestroy(&params->v);
	ierr = ISDestroy(&perm);
  	ierr = ISDestroy(&iperm);
	PetscFree(tmpArray);
	cs_spfree(spMat);
	cs_spfree(mt);
	cs_spfree(ident);
	cs_spfree(eye);
	cs_spfree(tmpMat);
	gsl_vector_free(tmpStates);
	cs_nfree(NN);
	cs_sfree(S);

}
double calcLikNLOpt(unsigned n, const double *point, double *gradients, void *p){
	int i,j;
	struct clam_lik_params * params = (struct clam_lik_params *) p;
	double localNNZ = params->nnz;
	double output = 666.0;
	double lik=0.0;
	double x[5];

	for(i = 0;i<5;i++) x[i] = point[i];
	MPI_Bcast(x,5,MPI_DOUBLE,0,PETSC_COMM_WORLD);

	for(i = 0;i<5;i++) gsl_vector_set(params->paramVector, i, x[i]);

	//fill in the expAFS table
	calcLogAFS_IM(p);
	params->nnz = localNNZ;
	//compute lik
	for(i=0;i<params->obsData->size1;i++){
		for(j=0;j<params->obsData->size2;j++){
			if (gsl_matrix_get(params->expAFS,i,j) > 0.0) //corners of AFS are zero prob
				lik += gsl_matrix_get(params->obsData,i,j) * (log(gsl_matrix_get(params->expAFS,i,j)));
		}
	}
	output = -1.0* lik;
	MPI_Bcast(&output,1,MPI_DOUBLE,0,PETSC_COMM_WORLD);
	if(params->rank == 0 && vbse){
		for(i = 0;i<5;i++) printf("x[%d]: %lf\t",i,x[i]);
		printf("lik: %lf\n",lik);
	}
	params->fEvals += 1;
//	gsl_matrix_prettyPrint(params->expAFS);
	return(output);
}

double calcLikNLOpt_gradients(unsigned n, const double *point, double *gradients, void *p){
	int i,j;
	//struct clam_lik_params * params = (struct clam_lik_params *) p;
	//double localNNZ = params->nnz;
	double output = 666.0;
	double lik=0.0;
	double newF;
	double x[5];
	double temp[5];
	double epsilon=1e-5;

	for(i = 0;i<5;i++) x[i] = point[i];
	MPI_Bcast(x,5,MPI_DOUBLE,0,PETSC_COMM_WORLD);

	lik = calcLikNLOpt(n, point, NULL, p);
	if (gradients != NULL){
		for(i=0; i<n ;i++){
			for( j=0; j< n; j++)
				temp[j] = x[j];

			temp[i] +=  epsilon;
			newF = calcLikNLOpt(n,temp,NULL,p);
			gradients[i] = (newF -  lik)/epsilon;
		}
	}

	output = lik;
	MPI_Bcast(&output,1,MPI_DOUBLE,0,PETSC_COMM_WORLD);
//	gsl_matrix_prettyPrint(params->expAFS);
	return(output);
}

//maximizeLik-- maximizes the logLikFunction function using constraints via
// nlOpt library

void maximizeLikNLOpt(double *lik, void *p, double *mle){
	//struct clam_lik_params * params = (struct clam_lik_params *) p;
	int np = 5;
	int err;
	//double x[5] = {1.0,1.0,1.0,1.0,1.5};
	double minimum=666.6;
	nlopt_opt opt;

	opt = nlopt_create(NLOPT_LD_LBFGS, np);
	//define constraints

	nlopt_set_lower_bounds(opt, lowerBounds);
	nlopt_set_upper_bounds(opt, upperBounds);
	nlopt_set_min_objective(opt, calcLikNLOpt_gradients, (void *) p);

	nlopt_set_xtol_rel(opt, 1e-6);
	//set the seed for deterministic sequence in COBYLA alg
	//nlopt_srand(1232333);
	err = nlopt_optimize(opt, mle, &minimum);
	if (err < 0) {
	    fprintf(stderr,"nlopt failed with code %d!\n",err);
	}
//	MPI_Barrier(MPI_COMM_WORLD);
	*lik=minimum;
	nlopt_destroy(opt);

}

//maximizeLikNLOpt_twoStage-- begins with coarse global, then goes local

void maximizeLikNLOpt_twoStage(double *lik, void *p, double *mle){
	//struct clam_lik_params * params = (struct clam_lik_params *) p;
	int np = 5;
	int err;
	//double x[5] = {1.0,1.0,1.0,1.0,1.5};
	double minimum=666.6;
	nlopt_opt opt;

	opt = nlopt_create(NLOPT_GN_DIRECT_L_RAND, np);

	nlopt_set_lower_bounds(opt, lowerBounds);
	nlopt_set_upper_bounds(opt, upperBounds);
	nlopt_set_min_objective(opt, calcLikNLOpt_gradients, (void *) p);

	//nlopt_set_xtol_rel(opt, 0.25);
	//set the seed for deterministic sequence in COBYLA alg
	//nlopt_srand(1232333);
	err = nlopt_optimize(opt, mle, &minimum);
	if (err < 0) {
	    fprintf(stderr,"nlopt failed with code %d!\n",err);
	}
//	MPI_Barrier(MPI_COMM_WORLD);
	*lik=minimum;
	nlopt_destroy(opt);

}

//maximizeLikNLOpt_MLSL-- multi-level restarts

void maximizeLikNLOpt_MLSL(double *lik, void *p, double *mle){
	int np = 5;
	int err;
	//double x[5] = {1.0,1.0,1.0,1.0,1.5};
	double minimum=666.6;
	nlopt_opt opt, l_opt;

	opt = nlopt_create(NLOPT_G_MLSL, np);

	nlopt_set_lower_bounds(opt, lowerBounds);
	nlopt_set_upper_bounds(opt, upperBounds);
	nlopt_set_min_objective(opt, calcLikNLOpt_gradients, (void *) p);

	l_opt = nlopt_create(NLOPT_LD_LBFGS, np);
	nlopt_set_lower_bounds(l_opt, lowerBounds);
	nlopt_set_upper_bounds(l_opt, upperBounds);
	nlopt_set_min_objective(l_opt, calcLikNLOpt_gradients, (void *) p);
	//define constraints

	nlopt_set_local_optimizer(opt, l_opt);

	//nlopt_set_xtol_rel(opt, 0.25);
	//set the seed for deterministic sequence in COBYLA alg
	//nlopt_srand(1232333);
	err = nlopt_optimize(opt, mle, &minimum);
	if (err < 0) {
	    fprintf(stderr,"nlopt failed with code %d!\n",err);
	}
//	MPI_Barrier(MPI_COMM_WORLD);
	*lik=minimum;
	nlopt_destroy(l_opt);
	nlopt_destroy(opt);

}

///////////////// Hessian stuff ///////////////////

double hessianMatrix_element(double lik, double *mle, int i, int j, double hi, double hj, void *p){
	//struct clam_lik_params * params = (struct clam_lik_params *) p;
	double Pi0, Pj0,output,fpp,fpm,fmp,fmm;

	Pi0 = mle[i]; Pj0 = mle[j];

	//printf("mle: %f\n",lik);
	if(i==j){
		mle[i] = Pi0 + hi;
		fpp = calcLikNLOpt(5,mle,NULL,p);
		mle[i] = Pj0 - hi;
		fmm = calcLikNLOpt(5,mle,NULL,p);
		output = ((fpp - lik) + (fmm - lik))/(hi*hi);
	}
	else{
		// f(xi + hi, xj + h)
		mle[i] = Pi0 + hi;
		mle[j] = Pj0 + hj;
		fpp = calcLikNLOpt(5,mle,NULL,p);

	        // f(xi + hi, xj - hj)
		mle[i] = Pi0 + hi;
		mle[j] = Pj0 - hj;
		fpm = calcLikNLOpt(5,mle,NULL,p);

	        // f(xi - hi, xj + hj)
		mle[i] = Pi0 - hi;
		mle[j] = Pj0 + hj;
		fmp = calcLikNLOpt(5,mle,NULL,p);

	        // f(xi - hi, xj - hj)
		mle[i] = Pi0 - hi;
		mle[j] = Pj0 - hj;
		fmm = calcLikNLOpt(5,mle,NULL,p);

		output = (fpp - fpm - fmp + fmm)/(4 * hi * hj);
	}
	mle[i]=Pi0; mle[j]=Pj0;
	return output;
}

gsl_matrix *hessian(double *mle, double lik, void *p){
	int i, j;
	double eps;
	gsl_matrix *H;

	H = gsl_matrix_alloc(5,5);
	eps = 0.01;
	for(i=0;i<5;i++){
		for (j=i;j<5;j++){
			gsl_matrix_set(H,i,j,hessianMatrix_element(lik, mle, i, j, eps, eps, p));
			gsl_matrix_set(H,j,i,gsl_matrix_get(H,i,j));
		}
	}
	return(H);
}

gsl_matrix *getFisherInfoMatrix(double *mle, double lik, void *p){
	gsl_matrix *H;
	gsl_matrix *fi;
	gsl_permutation *perm = gsl_permutation_alloc(5);
	int s;

	H = hessian(mle,lik,p);
	//gsl_matrix_prettyPrint(H);printf("\n");
	fi = gsl_matrix_alloc(5,5);
	gsl_matrix_scale(H,-1.0);
	gsl_linalg_LU_decomp (H, perm, &s);
	gsl_linalg_LU_invert (H, perm, fi);
	//gsl_matrix_prettyPrint(fi);printf("\n");
	return(fi);
}

gsl_matrix *getGodambeInfoMatrix(double *mle, double lik, void *p){
	gsl_matrix *H;
	gsl_matrix *J,*Jtemp, *Jinv, *origData, *boot;
	gsl_vector *cU, *grad_temp;
	struct clam_lik_params * params = (struct clam_lik_params *) p;
	gsl_permutation *perm = gsl_permutation_alloc(5);
	int s,i, nBoot;

	origData = params->obsData;
	nBoot = 100;

	H = hessian(mle,lik,p);
	gsl_matrix_scale(H,-1.0);
	//gsl_matrix_prettyPrint(H);printf("\n");
	J = gsl_matrix_alloc(5,5);
	Jinv = gsl_matrix_alloc(5,5);
	Jtemp = gsl_matrix_alloc(5,5);
	boot = gsl_matrix_alloc(origData->size1,origData->size2);

	gsl_matrix_set_zero(J);
	gsl_matrix_set_zero(Jinv);
	gsl_matrix_set_zero(Jtemp);
	cU = gsl_vector_alloc(5);
	gsl_vector_set_zero(cU);

	//do bootstraps, estimate J
	for(i=0;i<nBoot;i++){
		gsl_matrix_bootstrap(origData, boot, params->rng);
		params->obsData = boot;
		grad_temp = getGradient(mle,p);
		gsl_vector_outer_product(grad_temp,grad_temp,Jtemp);
		gsl_matrix_add(J,Jtemp);
		gsl_vector_add(cU,grad_temp);
		//gsl_matrix_prettyPrint(Jtemp);
		gsl_vector_free(grad_temp);
	}
	gsl_vector_scale(cU,1.0/(float)nBoot);
	gsl_matrix_scale(J,1.0/(float)nBoot);
	//printf("cU:\n");
	//gsl_vector_fprintf(stdout,cU,"%f");
	//printf("J:\n");
	//gsl_matrix_prettyPrint(J);printf("\n");
	//G = H*J^-1*H
	gsl_linalg_LU_decomp (J, perm, &s);
	gsl_linalg_LU_invert (J, perm, Jinv);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, H, Jinv, 0.0, Jtemp);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Jtemp, H, 0.0, J);
	//printf("godambe:\n");
	//gsl_matrix_prettyPrint(J);printf("\n");
	//need to invert this to get the Godambe information matrix
	gsl_linalg_LU_decomp (J, perm, &s);
	gsl_linalg_LU_invert (J, perm, Jinv);
	return(Jinv);
}

gsl_vector *getGradient(double *mle, void *p){
	int i, j;
	double eps, epsVals[5], paramTemp[5], fp, fm;
	gsl_vector *gradient;

	gradient = gsl_vector_alloc(5);
	eps = 0.01;
	for(i=0;i<5;i++){
		if( mle[i] != 0){
			epsVals[i] = mle[i] * eps;
			//printf("eps[%d]=%f\n",i,epsVals[i]);
		}
		else{
			epsVals[i] = eps;
		}
		paramTemp[i] = mle[i];
	}
	for(i=0;i<5;i++){
		for(j=0;j<5;j++) paramTemp[j] = mle[j];
		if( mle[i] != 0){
			paramTemp[i] = mle[i] + epsVals[i];
			fp = calcLikNLOpt(5,paramTemp,NULL,p);
			paramTemp[i] = mle[i] - epsVals[i];
			fm = calcLikNLOpt(5,paramTemp,NULL,p);
			gsl_vector_set(gradient,i,(fp-fm)/(2*epsVals[i]));
		}
		else{
			paramTemp[i] = mle[i] + epsVals[i];
			fp = calcLikNLOpt(5,paramTemp,NULL,p);
			paramTemp[i] = mle[i];
			fm = calcLikNLOpt(5,paramTemp,NULL,p);
			gsl_vector_set(gradient,i,(fp-fm)/epsVals[i]);
		}
	}
	return(gradient);
}

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA