text:
// Copyright (c) 2013-2017 Anton Kozhevnikov, Thomas Schulthess
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that
// the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the
//    following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
//    and the following disclaimer in the documentation and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

/** \file sbessel.h
 *
 *  \brief Contains implementation of sirius::Spherical_Bessel_functions and sirius::sbessel_approx classes.
 */

#ifndef __SBESSEL_PW_H__
#define __SBESSEL_PW_H__

#include <gsl/gsl_sf_bessel.h>
#include "eigenproblem.h"
#include "Unit_cell/unit_cell.h"

namespace sirius
{

/// Spherical Bessel functions \f$ j_{\ell}(q x) \f$ up to lmax.
class Spherical_Bessel_functions
{
    private:
        int lmax_{-1};
        double q_{0};
        Radial_grid<double> const* rgrid_{nullptr};

        std::vector<Spline<double>> sbessel_;

    public:

        Spherical_Bessel_functions()
        {
        }

        Spherical_Bessel_functions(int lmax__, Radial_grid<double> const& rgrid__, double q__)
            : lmax_(lmax__)
            , q_(q__)
            , rgrid_(&rgrid__)
        {
            assert(q_ >= 0);

            sbessel_ = std::vector<Spline<double>>(lmax__ + 2);
            for (int l = 0; l <= lmax__ + 1; l++) {
                sbessel_[l] = Spline<double>(rgrid__);
            }

            std::vector<double> jl(lmax__ + 2);
            for (int ir = 0; ir < rgrid__.num_points(); ir++) {
                double t = rgrid__[ir] * q__;
                gsl_sf_bessel_jl_array(lmax__ + 1, t, &jl[0]);
                for (int l = 0; l <= lmax__ + 1; l++) {
                    sbessel_[l](ir) = jl[l];
                }
            }

            for (int l = 0; l <= lmax__ + 1; l++) {
                sbessel_[l].interpolate();
            }
        }

        static void sbessel(int lmax__, double t__, double* jl__)
        {
            gsl_sf_bessel_jl_array(lmax__, t__, jl__);
        }

        static void sbessel_deriv_q(int lmax__, double q__, double x__, double* jl_dq__)
        {
            std::vector<double> jl(lmax__ + 2);
            sbessel(lmax__ + 1, x__ * q__, &jl[0]);

            for (int l = 0; l <= lmax__; l++) {
                if (q__ != 0) {
                    jl_dq__[l] = (l / q__) * jl[l] - x__ * jl[l + 1];
                } else {
                    if (l == 1) {
                        jl_dq__[l] = x__ / 3;
                    } else {
                        jl_dq__[l] = 0;
                    }
                }
            }
        }

        Spline<double> const& operator[](int l__) const
        {
            assert(l__ <= lmax_);
            return sbessel_[l__];
        }

        /// Derivative of Bessel function with respect to q.
        /** \f[
         *    \frac{\partial j_{\ell}(q x)}{\partial q} = \frac{\ell}{q} j_{\ell}(q x) - x j_{\ell+1}(q x)
         *  \f]
         */
        Spline<double> deriv_q(int l__)
        {
            assert(l__ <= lmax_);
            assert(q_ >= 0);
            Spline<double> s(*rgrid_);
            if (q_ != 0) {
                for (int ir = 0; ir < rgrid_->num_points(); ir++) {
                    s(ir) = (l__ / q_) * sbessel_[l__](ir) - (*rgrid_)[ir] * sbessel_[l__ + 1](ir);
                }
            } else {
                if (l__ == 1) {
                    for (int ir = 0; ir < rgrid_->num_points(); ir++) {
                        s(ir) = (*rgrid_)[ir] / 3;
                    }
                }
            }
            s.interpolate();
            return std::move(s);
        }
};

class sbessel_approx
{
    private:

        Unit_cell const& unit_cell_;

        int lmax_;

        mdarray<std::vector<double>, 2> qnu_;
        mdarray<double, 4> coeffs_;

        int nqnu_max_;

        static double sbessel_l2norm(double nu, int l, double R)
        {
            if (std::abs(nu) < 1e-10) TERMINATE_NOT_IMPLEMENTED;

            if (l == 0)
            {
                return (nu * R * 2 - std::sin(nu * R * 2)) / 4 / std::pow(nu, 3);
            }
            else
            {
                double jl[l + 2];
                gsl_sf_bessel_jl_array(l + 1, R * nu, &jl[0]);
                return std::pow(R, 3) * (jl[l] * jl[l] - jl[l + 1] * jl[l - 1]) / 2;
            }
        }

    public:

        sbessel_approx(Unit_cell const& unit_cell__,
                       int lmax__,
                       double const qmin__,
                       double const qmax__,
                       double const eps__)
            : unit_cell_(unit_cell__),
              lmax_(lmax__)
        {
            PROFILE("sirius::sbessel_approx");

            qnu_ = mdarray<std::vector<double>, 2>(lmax_ + 1, unit_cell_.num_atom_types());

            for (int l = 0; l <= lmax_; l++) {
                for (int iat = 0; iat < unit_cell_.num_atom_types(); iat++) {
                    qnu_(l, iat) = build_approx_freq(qmin__, qmax__, l, unit_cell_.atom_type(iat).mt_radius(), eps__);
                }
            }

            nqnu_max_ = 0;
            for (int l = 0; l <= lmax_; l++) {
                for (int iat = 0; iat < unit_cell_.num_atom_types(); iat++) {
                    nqnu_max_ = std::max(nqnu_max_, static_cast<int>(qnu_(l, iat).size()));
                }
            }
        }

        void approximate(std::vector<double> const& q__)
        {
            PROFILE("sirius::sbessel_approx::approximate");

            coeffs_ = mdarray<double, 4>(nqnu_max_, q__.size(), lmax_ + 1, unit_cell_.num_atom_types());

            #pragma omp parallel for
            for (int l = 0; l <= lmax_; l++)
            {
                for (int iat = 0; iat < unit_cell_.num_atom_types(); iat++)
                {
                    int n = nqnu(l, iat);

                    mdarray<double, 2> A(n, n);
                    for (int iq = 0; iq < n; iq++)
                    {
                        for (int jq = 0; jq <= iq; jq++)
                        {
                            A(jq, iq) = A(iq, jq) = overlap(qnu_(l, iat)[jq], qnu_(l, iat)[iq], l,
                                                            unit_cell_.atom_type(iat).mt_radius());
                        }

                        for (int j = 0; j < (int)q__.size(); j++)
                        {
                            if (std::abs(q__[j]) < 1e-12)
                            {
                                coeffs_(iq, j, l, iat) = 0;
                            }
                            else
                            {
                                coeffs_(iq, j, l, iat) = overlap(qnu_(l, iat)[iq], q__[j], l,
                                                                 unit_cell_.atom_type(iat).mt_radius());
                            }
                        }
                    }
                    linalg<CPU>::gesv(n, (int)q__.size(), A.at<CPU>(), A.ld(), &coeffs_(0, 0, l, iat), coeffs_.ld());
                }
            }
        }

        inline double qnu(int const iq, int const l, int const iat)
        {
            return qnu_(l, iat)[iq];
        }

        inline int nqnu(int const l, int const iat)
        {
            return static_cast<int>(qnu_(l, iat).size());
        }

        inline int nqnu_max()
        {
            return nqnu_max_;
        }

        inline double coeff(int const iq, int const j, int const l, int const iat)
        {
            return coeffs_(iq, j, l, iat);
        }

        // \int_0^{R} j(nu1 * r) * j(nu2 * r) * r^2 dr
        // this integral can be computed analytically
        static double overlap(double nu1__, double nu2__, int l__, double R__)
        {
            if (std::abs(nu1__) < 1e-10 || std::abs(nu2__) < 1e-10) TERMINATE_NOT_IMPLEMENTED;

            if (std::abs(nu1__ - nu2__) < 1e-12)
            {
                if (l__ == 0)
                {
                    return (nu2__ * R__ * 2 - std::sin(nu2__ * R__ * 2)) / 4 / std::pow(nu2__, 3);
                }
                else
                {
                    double jl[l__ + 2];
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu2__, &jl[0]);
                    return std::pow(R__, 3) * (jl[l__] * jl[l__] - jl[l__ + 1] * jl[l__ - 1]) / 2;
                }
            }
            else
            {
                if (l__ == 0)
                {
                    return (nu2__ * std::cos(nu2__ * R__) * std::sin(nu1__ * R__) - nu1__ * std::cos(nu1__ * R__) * std::sin(nu2__ * R__)) /
                           (std::pow(nu1__, 3) * nu2__ - nu1__ * std::pow(nu2__, 3));
                }
                else
                {
                    double j1[l__ + 2];
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu1__, &j1[0]);

                    double j2[l__ + 2];
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu2__, &j2[0]);

                    return std::pow(R__, 2) * (nu2__ * j2[l__ - 1] * j1[l__] - nu1__ * j1[l__ - 1] * j2[l__]) / (std::pow(nu1__, 2) - std::pow(nu2__, 2));
                }
            }
        }

        std::vector<double> build_approx_freq(double const qmin__,
                                              double const qmax__,
                                              int const l__,
                                              double const R__,
                                              double const eps__)
        {
            std::vector<double> qnu;

            double min_val;
            int n = 2;

            do
            {
                n++;
                qnu.resize(n);
                for (int i = 0; i < n; i++) qnu[i] = qmin__ + (qmax__ - qmin__) * i / (n - 1);

                dmatrix<double_complex> ovlp(n, n);
                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j <= i; j++)
                    {
                        double o = overlap(qnu[j], qnu[i], l__, R__);
                        ovlp(j, i) = o / sbessel_l2norm(qnu[i], l__, R__) / sbessel_l2norm(qnu[j], l__, R__);
                    }
                }

                std::vector<double> eval(n);
                dmatrix<double_complex> z(n, n);

                Eigensolver_lapack<double_complex> solver;
                solver.solve(n, ovlp, &eval[0], z);
                min_val = eval[0];

            } while (min_val > eps__);

            return qnu;
        }
};

class Spherical_Bessel_approximant
{
    private:

        int lmax_;

        double R_;

        /// List of Bessel function scaling factors for each angular momentum.
        std::vector< std::vector<double> > qnu_;

        //mdarray<double, 4> coeffs_;

        int nqnu_max_;

        static double sbessel_l2norm(double nu, int l, double R)
        {
            if (std::abs(nu) < 1e-10)
            {
                if (l == 0) return std::pow(R, 3) / 3.0;
                return 0;
            }

            if (l == 0)
            {
                return (nu * R * 2 - std::sin(nu * R * 2)) / 4 / std::pow(nu, 3);
            }
            else
            {
                double jl[l + 2];
                gsl_sf_bessel_jl_array(l + 1, R * nu, &jl[0]);
                return std::pow(R, 3) * (jl[l] * jl[l] - jl[l + 1] * jl[l - 1]) / 2;
            }
        }

    public:

        Spherical_Bessel_approximant(int lmax__,
                                     double R__,
                                     double const qmin__,
                                     double const qmax__,
                                     double const eps__)
            : lmax_(lmax__),
              R_(R__)
        {
            PROFILE("sirius::Spherical_Bessel_approximant");

            qnu_ = std::vector< std::vector<double> >(lmax_ + 1);

            #pragma omp parallel for
            for (int l = 0; l <= lmax_; l++)
                qnu_[l] = build_approx_freq(qmin__, qmax__, l, R_, eps__);

            nqnu_max_ = 0;
            for (int l = 0; l <= lmax_; l++)
                nqnu_max_ = std::max(nqnu_max_, nqnu(l));
        }

        std::vector<double> approximate(int l__, double nu__)
        {
            int n = nqnu(l__);
            std::vector<double> x(n);
            matrix<double> A(n, n);

            for (int iq = 0; iq < n; iq++)
            {
                for (int jq = 0; jq <= iq; jq++)
                {
                    A(jq, iq) = A(iq, jq) = overlap(qnu(jq, l__), qnu(iq, l__), l__, R_);
                }
                x[iq] = overlap(qnu(iq, l__), nu__, l__, R_);
            }
            linalg<CPU>::gesv(n, 1, A.at<CPU>(), A.ld(), &x[0], n);
            return x;
        }

        void approximate(std::vector<double> const& q__)
        {
            //runtime::Timer t("sirius::sbessel_approx::approximate");

            //coeffs_ = mdarray<double, 4>(nqnu_max_, q__.size(), lmax_ + 1, unit_cell_.num_atom_types());
            //
            //#pragma omp parallel for
            //for (int l = 0; l <= lmax_; l++)
            //{
            //    for (int iat = 0; iat < unit_cell_.num_atom_types(); iat++)
            //    {
            //        int n = nqnu(l, iat);

            //        mdarray<double, 2> A(n, n);
            //        for (int iq = 0; iq < n; iq++)
            //        {
            //            for (int jq = 0; jq <= iq; jq++)
            //            {
            //                A(jq, iq) = A(iq, jq) = overlap(qnu_(l, iat)[jq], qnu_(l, iat)[iq], l,
            //                                                unit_cell_.atom_type(iat).mt_radius());
            //            }

            //            for (int j = 0; j < (int)q__.size(); j++)
            //            {
            //                if (std::abs(q__[j]) < 1e-12)
            //                {
            //                    coeffs_(iq, j, l, iat) = 0;
            //                }
            //                else
            //                {
            //                    coeffs_(iq, j, l, iat) = overlap(qnu_(l, iat)[iq], q__[j], l,
            //                                                     unit_cell_.atom_type(iat).mt_radius());
            //                }
            //            }
            //        }
            //        linalg<CPU>::gesv(n, (int)q__.size(), A.at<CPU>(), A.ld(), &coeffs_(0, 0, l, iat), coeffs_.ld());
            //    }
            //}
        }

        inline double qnu(int const iq, int const l) const
        {
            return qnu_[l][iq];
        }

        inline int nqnu(int const l) const
        {
            return static_cast<int>(qnu_[l].size());
        }

        inline int nqnu_max() const
        {
            return nqnu_max_;
        }

        //inline double coeff(int const iq, int const j, int const l, int const iat)
        //{
        //    return coeffs_(iq, j, l, iat);
        //}
        //

        // \int_0^{R} j(nu1 * r) * j(nu2 * r) * r^2 dr
        // this integral can be computed analytically
        static double overlap(double nu1__, double nu2__, int l__, double R__)
        {
            if (std::abs(nu1__) < 1e-10 && std::abs(nu2__) < 1e-10 && l__ == 0) return std::pow(R__, 3) / 3.0;

            if ((std::abs(nu1__) < 1e-10 || std::abs(nu2__) < 1e-10) && l__ > 0) return 0;

            if ((std::abs(nu1__) < 1e-10 || std::abs(nu2__) < 1e-10) && l__ == 0)
            {
                double nu = std::max(nu1__, nu2__);
                double nuR = nu * R__;
                return (std::sin(nuR) - nuR * std::cos(nuR)) / std::pow(nu, 3);
            }

            if (std::abs(nu1__ - nu2__) < 1e-12)
            {
                if (l__ == 0)
                {
                    return (nu2__ * R__ * 2 - std::sin(nu2__ * R__ * 2)) / 4 / std::pow(nu2__, 3);
                }
                else
                {
                    double jl[l__ + 2];
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu2__, &jl[0]);
                    return std::pow(R__, 3) * (jl[l__] * jl[l__] - jl[l__ + 1] * jl[l__ - 1]) / 2;
                }
            }
            else
            {
                if (l__ == 0)
                {
                    return (nu2__ * std::cos(nu2__ * R__) * std::sin(nu1__ * R__) - nu1__ * std::cos(nu1__ * R__) * std::sin(nu2__ * R__)) /
                           (std::pow(nu1__, 3) * nu2__ - nu1__ * std::pow(nu2__, 3));
                }
                else
                {
                    double j1[l__ + 2];
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu1__, &j1[0]);

                    double j2[l__ + 2];
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu2__, &j2[0]);

                    return std::pow(R__, 2) * (nu2__ * j2[l__ - 1] * j1[l__] - nu1__ * j1[l__ - 1] * j2[l__]) / (std::pow(nu1__, 2) - std::pow(nu2__, 2));
                }
            }
            TERMINATE("this is wrong");
            return -1;
        }

        std::vector<double> build_approx_freq(double const qmin__,
                                              double const qmax__,
                                              int const l__,
                                              double const R__,
                                              double const eps__)
        {
            std::vector<double> qnu;

            double min_val;
            int n = 2;

            do
            {
                n++;
                qnu.resize(n);
                for (int i = 0; i < n; i++) qnu[i] = qmin__ + (qmax__ - qmin__) * i / (n - 1);

                dmatrix<double> ovlp(n, n);
                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j <= i; j++)
                    {
                        double o = overlap(qnu[j], qnu[i], l__, R__);
                        ovlp(j, i) = o / sbessel_l2norm(qnu[i], l__, R__) / sbessel_l2norm(qnu[j], l__, R__);
                    }
                }

                std::vector<double> eval(n);
                dmatrix<double> z(n, n);

                Eigensolver_lapack<double> solver;
                solver.solve(n, ovlp, &eval[0], z);
                min_val = eval[0];

            } while (min_val > eps__);

            return qnu;
        }
};

class Spherical_Bessel_approximant2
{
    private:

        int lmax_;

        double R_;

        /// List of Bessel function scaling factors for each angular momentum.
        std::vector<double> qnu_;

        static double sbessel_l2norm(double nu, int l, double R)
        {
            if (std::abs(nu) < 1e-10)
            {
                if (l == 0) return std::pow(R, 3) / 3.0;
                return 0;
            }

            if (l == 0)
            {
                return (nu * R * 2 - std::sin(nu * R * 2)) / 4 / std::pow(nu, 3);
            }
            else
            {
                double jl[l + 2];
                gsl_sf_bessel_jl_array(l + 1, R * nu, &jl[0]);
                return std::pow(R, 3) * (jl[l] * jl[l] - jl[l + 1] * jl[l - 1]) / 2;
            }
        }

    public:

        Spherical_Bessel_approximant2(int lmax__,
                                     double R__,
                                     double const qmin__,
                                     double const qmax__,
                                     int nq__)
            : lmax_(lmax__),
              R_(R__)
        {
            PROFILE("sirius::Spherical_Bessel_approximant");

            int nq = nfreq(qmin__, qmax__, 0, R__, 1e-12);
            qnu_.resize(nq);
            for (int i = 0; i < nq; i++) qnu_[i] = qmin__ + (qmax__ - qmin__) * i / (nq - 1);

        }

        std::vector<double> approximate(int l__, double nu__)
        {
            int n = nqnu();
            std::vector<double> x(n);
            matrix<double> A(n, n);

            for (int iq = 0; iq < n; iq++)
            {
                for (int jq = 0; jq <= iq; jq++)
                {
                    A(jq, iq) = A(iq, jq) = overlap(qnu(jq), qnu(iq), l__, R_);
                }
                x[iq] = overlap(qnu(iq), nu__, l__, R_);
            }
            linalg<CPU>::gesv(n, 1, A.at<CPU>(), A.ld(), &x[0], n);
            return x;
        }

        inline double qnu(int const iq) const
        {
            return qnu_[iq];
        }

        inline int nqnu() const
        {
            return static_cast<int>(qnu_.size());
        }

        // \int_0^{R} j(nu1 * r) * j(nu2 * r) * r^2 dr
        // this integral can be computed analytically
        static double overlap(double nu1__, double nu2__, int l__, double R__)
        {
            if (std::abs(nu1__) < 1e-10 && std::abs(nu2__) < 1e-10 && l__ == 0) return std::pow(R__, 3) / 3.0;

            if ((std::abs(nu1__) < 1e-10 || std::abs(nu2__) < 1e-10) && l__ > 0) return 0;

            if ((std::abs(nu1__) < 1e-10 || std::abs(nu2__) < 1e-10) && l__ == 0)
            {
                double nu = std::max(nu1__, nu2__);
                double nuR = nu * R__;
                return (std::sin(nuR) - nuR * std::cos(nuR)) / std::pow(nu, 3);
            }

            if (std::abs(nu1__ - nu2__) < 1e-12)
            {
                if (l__ == 0)
                {
                    return (nu2__ * R__ * 2 - std::sin(nu2__ * R__ * 2)) / 4 / std::pow(nu2__, 3);
                }
                else
                {
                    std::vector<double> jl(l__ + 2);
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu2__, &jl[0]);
                    return std::pow(R__, 3) * (jl[l__] * jl[l__] - jl[l__ + 1] * jl[l__ - 1]) / 2;
                }
            }
            else
            {
                if (l__ == 0)
                {
                    return (nu2__ * std::cos(nu2__ * R__) * std::sin(nu1__ * R__) - nu1__ * std::cos(nu1__ * R__) * std::sin(nu2__ * R__)) /
                           (std::pow(nu1__, 3) * nu2__ - nu1__ * std::pow(nu2__, 3));
                }
                else
                {
                    std::vector<double> j1(l__ + 2);
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu1__, &j1[0]);

                    std::vector<double> j2(l__ + 2);
                    gsl_sf_bessel_jl_array(l__ + 1, R__ * nu2__, &j2[0]);

                    return std::pow(R__, 2) * (nu2__ * j2[l__ - 1] * j1[l__] - nu1__ * j1[l__ - 1] * j2[l__]) / (std::pow(nu1__, 2) - std::pow(nu2__, 2));
                }
            }
            TERMINATE("this is wrong");
            return -1;
        }

        int nfreq(double const qmin__,
                  double const qmax__,
                  int const l__,
                  double const R__,
                  double const eps__)
        {
            std::vector<double> qnu;

            double min_val;
            int n = 2;

            do
            {
                n++;
                qnu.resize(n);
                for (int i = 0; i < n; i++) qnu[i] = qmin__ + (qmax__ - qmin__) * i / (n - 1);

                dmatrix<double> ovlp(n, n);
                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j <= i; j++)
                    {
                        double o = overlap(qnu[j], qnu[i], l__, R__);
                        ovlp(j, i) = o / sbessel_l2norm(qnu[i], l__, R__) / sbessel_l2norm(qnu[j], l__, R__);
                    }
                }

                std::vector<double> eval(n);
                dmatrix<double> z(n, n);

                Eigensolver_lapack<double> solver;
                solver.solve(n, ovlp, &eval[0], z);
                min_val = eval[0];

                if (n > 100) return 100;

            } while (min_val > eps__);

            return n;
        }
};

};

#endif


text:
#include <math.h>
#include <stdlib.h>
#if !defined(__APPLE__)
#include <malloc.h>
#endif
#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <string.h>
#include <stdbool.h>

#include <fftw3.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

#include "../cosmolike_core/theory/basics.c"
#include "../cosmolike_core/theory/structs.c"
#include "../cosmolike_core/theory/parameters.c"
#include "../cosmolike_core/emu17/P_cb/emu.c"
#include "../cosmolike_core/theory/recompute.c"
#include "../cosmolike_core/theory/cosmo3D.c"
#include "../cosmolike_core/theory/redshift_spline.c"
#include "../cosmolike_core/theory/halo.c"
#include "../cosmolike_core/theory/HOD.c"
#include "../cosmolike_core/theory/pt.c"
#include "../cosmolike_core/theory/cosmo2D_fourier.c"
#include "../cosmolike_core/theory/IA.c"
#include "../cosmolike_core/theory/BAO.c"
#include "../cosmolike_core/theory/external_prior.c"
#include "../cosmolike_core/theory/covariances_3D.c"
#include "../cosmolike_core/theory/covariances_fourier.c"
#include "../cosmolike_core/theory/CMBxLSS_fourier.c"
#include "../cosmolike_core/theory/covariances_CMBxLSS_fourier.c"

#include "init_LSSxCMB.c"

// Naming convention:
// l = galaxy positions ("l" as in "lens sample")
// k = kappa CMB ("k" as in "kappa")
// s = kappa from source galaxies ("s" as in "source sample")
// And alphabetical order

void run_cov_ls_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_ll_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_ll_ls(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_ll_ll(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_ls_ls(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_ss_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);

void run_cov_ls_kk(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int start);
void run_cov_ls_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int zs2, int start);
void run_cov_lk_lk(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int zl2, int start);
void run_cov_lk_ls(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int n2, int start);
void run_cov_lk_kk(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int start);
void run_cov_lk_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int zs2, int start);
void run_cov_lk_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int n2,int start);
void run_cov_kk_kk(char *OUTFILE, char *PATH, double *ell, double *dell,int start);
void run_cov_kk_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int zs2, int start);
void run_cov_kk_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int n2, int start);
void run_cov_ks_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int zs1, int zs2,int start);
void run_cov_ks_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int zs1, int n2, int start);
void run_cov_ll_kk(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int start);
void run_cov_ll_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int zs2, int start);
void run_cov_ll_lk(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int zl2,int start);

void run_cov_ls_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int zl,zs,z3,z4,nl1,nl2,weight,i,j;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  zl = ZL(n1); zs = ZS(n1);
  printf("\nN_ggl = %d (%d, %d)\n", n1,zl,zs);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n", n2,z3,z4);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 <like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],zl);
      if (weight && ell[nl2] < like.lmax_shear){
        if (test_zoverlap(zl,z3)*test_zoverlap(zl,z4)){
          c_ng = cov_NG_gl_shear_tomo(ell[nl1],ell[nl2],zl,zs,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_gl_shear_tomo(ell[nl1],dell[nl1],zl,zs,z3,z4);
        }
      }
      i=like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1;
      j=like.Ncl*n2+nl2;
      fprintf(F1, "%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell[nl1],ell[nl2],zl,zs,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_ll_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,i,j,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl = %d \n", n1);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n", n2,z3,z4);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 <like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],z1);
      if (weight && ell[nl2] < like.lmax_shear){
        if (test_zoverlap(z1,z3)*test_zoverlap(z1,z4)){
          c_ng = cov_NG_cl_shear_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_cl_shear_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        }
      }
      i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1;
      j=like.Ncl*n2+nl2;
      fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_ll_ls(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,zl,zs,nl1,nl2,i,j,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl_1 = %d \n", n1);
  zl = ZL(n2); zs = ZS(n2);
  printf("N_tomo_2 = %d (%d, %d)\n", n2,zl,zs);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],z1)*test_kmax(ell[nl2],zl);
      if (weight){
        if((z1 == zl)*test_zoverlap(z1,zs)){c_ng = cov_NG_cl_gl_tomo(ell[nl1],ell[nl2],z1,z2,zl,zs);}
        if (nl1 == nl2){
          c_g =  cov_G_cl_gl_tomo(ell[nl1],dell[nl1],z1,z2,zl,zs);
        }
      }
      i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1;
      j=like.Ncl*(tomo.shear_Npowerspectra+n2)+nl2;
      fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],z1,z2,zl,zs,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_ll_ll(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,i,j,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl_1 = %d \n", n1);
  z3 = n2; z4 = n2;
  printf("N_cl_2 = %d\n", n2);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],z1)*test_kmax(ell[nl2],z3);
      if ((z1 == z3) && weight) {
        c_ng = cov_NG_cl_cl_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
      }
      if (nl1 == nl2){
        c_g =  cov_G_cl_cl_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
      }
      i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1;
      j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n2)+nl2;
      fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_ls_ls(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int zl1,zl2,zs1,zs2,nl1,nl2,i,j,weight;
  double c_ng, c_g;
  double fsky = survey.area/41253.0;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");

  zl1 = ZL(n1); zs1 = ZS(n1);
  printf("\nN_tomo_1 = %d (%d, %d)\n", n1,zl1,zs1);
  zl2 = ZL(n2); zs2 = ZS(n2);
  printf("N_tomo_2 = %d (%d, %d)\n", n2,zl2,zs2);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],zl1)*test_kmax(ell[nl2],zl2);
      if (weight && zl1 == zl2) {
        c_ng = cov_NG_gl_gl_tomo(ell[nl1],ell[nl2],zl1,zs1,zl2,zs2);
      }
      if (nl1 == nl2){
        c_g =  cov_G_gl_gl_tomo(ell[nl1],dell[nl1],zl1,zs1,zl2,zs2);
      }
      if (weight ==0 && n2 != n1){
        c_g = 0;
      }
      i=like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1;
      j=like.Ncl*(tomo.shear_Npowerspectra+n2)+nl2;
      fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,zs1,zl2,zs2,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_ss_ss(char *OUTFILE, char *PATH, double *ell, double *dell,int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,i,j,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  z1 = Z1(n1); z2 = Z2(n1);
  printf("N_shear = %d\n", n1);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n",n2,z3,z4);
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      if (ell[nl1] < like.lmax_shear && ell[nl2] < like.lmax_shear){
        c_ng = cov_NG_shear_shear_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
      }
      if (nl1 == nl2){
        c_g =  cov_G_shear_shear_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        if (ell[nl1] > like.lmax_shear && n1!=n2){c_g = 0.;}
      }
      i=like.Ncl*n1+nl1;
      j=like.Ncl*n2+nl2;
      fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
      //printf("%d %d %e %e %d %d %d %d %e %e\n", like.Ncl*n1+nl1,like.Ncl*(n2)+nl2, ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

/*** CMBkappa part ****/
// ls_kk
void run_cov_ls_kk(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int start)
{
   int zl1, zs1,i,j, weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   zl1 = ZL(n1); zs1 = ZS(n1);
   printf("Bin for ls: %d (%d, %d)\n", n1, zl1, zs1);
   for (int nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (int nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1], zl1);
         if (weight && ell[nl2]<like.lmax_kappacmb) {
            c_ng = cov_NG_gs_kk(ell[nl1],ell[nl2], zl1, zs1);
            if (nl1==nl2){
               c_g = cov_G_gs_kk(ell[nl1], dell[nl1], zl1, zs1);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+tomo.shear_Nbin)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,zs1,0,0,c_g,c_ng);
      }
   }
   fclose(F1);
}

// ls_ks
void run_cov_ls_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int zs2, int start)
{
   int zl1, zs1,i,j,weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   zl1 = ZL(n1); zs1 = ZS(n1);
   printf("Bin for ls: %d (%d, %d)\n", n1, zl1, zs1);
   printf("Bin for ks: %d\n", zs2);
   for (int nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (int nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1], zl1);
         if (weight && ell[nl2]<like.lmax_kappacmb) {
            if (test_zoverlap(zl1, zs2)) {
               c_ng = cov_NG_gs_ks(ell[nl1],ell[nl2], zl1, zs1, zs2);
            }
            if (nl1==nl2){
               c_g = cov_G_gs_ks(ell[nl1], dell[nl1], zl1, zs1, zs2);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+zs2)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,zs1,0,zs2,c_g,c_ng);
      }
   }
   fclose(F1);
}
// lk_lk
void run_cov_lk_lk(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int zl2, int start)
{
   int weight,i,j;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   printf("Lens bins: %d, %d\n", zl1, zl2);

   for (int nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (int nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1],zl1)*test_kmax(ell[nl2],zl2);
         bool compute = weight * (ell[nl1]<like.lmax_kappacmb) * (ell[nl2]<like.lmax_kappacmb);
         if (compute && (zl1==zl2)) {
            c_ng = cov_NG_gk_gk(ell[nl1],ell[nl2], zl1, zl2);
         }
         if (nl1==nl2){
            c_g = cov_G_gk_gk(ell[nl1], dell[nl1], zl1, zl2);
         }
         if (!compute && (zl2!=zl1)) {
            c_g = 0;
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+zl1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+zl2)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,0,zl2,0,c_g,c_ng);
      }
   }
   fclose(F1);
}

// lk_ls
void run_cov_lk_ls(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int n2, int start)
{
   int zl2,zs2,nl1,nl2,i,j, weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   // sprintf(filename,"%scov/%s_%s_cov_gkgs_Nell%d_Ns%d_Ng%d_%d", PATH, survey.name, cmb.name, like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   printf("Lens bin for lk: %d\n", zl1);
   zl2 = ZL(n2); zs2 = ZS(n2);
   printf("N_tomo for ls = %d (%d, %d)\n", n2,zl2,zs2);

   for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1],zl1)*test_kmax(ell[nl2],zl2);
         if (weight && ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_kappacmb) {
            if (zl1==zl2) {
               c_ng = cov_NG_gk_gs(ell[nl1],ell[nl2],zl1,zl2,zs2);
            }
            if (nl1==nl2){
               c_g = cov_G_gk_gs(ell[nl1], dell[nl1], zl1,zl2,zs2);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+zl1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+n2)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,0,zl2,zs2,c_g,c_ng);
      }
   }
   fclose(F1);
}

// lk_kk
void run_cov_lk_kk(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int start)
{
   int weight,i,j;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   // sprintf(filename,"%scov/%s_%s_cov_gkkk_Nell%d_Ns%d_Ng%d_%d", PATH, survey.name, cmb.name, like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   printf("Bin for lk: %d\n", zl1);

   for (int nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (int nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1],zl1);
         if (weight && ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_kappacmb) {
            c_ng = cov_NG_gk_kk(ell[nl1],ell[nl2],zl1);
            if (nl1==nl2){
               c_g = cov_G_gk_kk(ell[nl1],dell[nl1],zl1);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+zl1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+tomo.shear_Nbin)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,0,0,0,c_g,c_ng);
      }
   }
   fclose(F1);
}

// lk_ks
void run_cov_lk_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int zs2, int start)
{
   int weight,i,j;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   // sprintf(filename,"%scov/%s_%s_cov_gkks_Nell%d_Ns%d_Ng%d_%d", PATH, survey.name, cmb.name, like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   printf("Bin for lk: %d\n", zl1);
   printf("Bin for ks: %d\n", zs2);

   for (int nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (int nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1], zl1);
         if (weight && ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_kappacmb) {
            if (test_zoverlap(zl1, zs2)) {
               c_ng = cov_NG_gk_ks(ell[nl1],ell[nl2], zl1, zs2);
            }
            if (nl1==nl2){
               c_g =  cov_G_gk_ks(ell[nl1], dell[nl1], zl1, zs2);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+zl1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+zs2)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,0,0,zs2,c_g,c_ng);
      }
   }
   fclose(F1);
}

// lk_ss
void run_cov_lk_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int zl1, int n2,int start)
{
   int zs2,zs3,nl1,nl2,weight,i,j;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   // sprintf(filename,"%scov/%s_%s_cov_gkss_Nell%d_Ns%d_Ng%d_%d",PATH,survey.name, cmb.name, like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   printf("Bin for lk: %d\n", zl1);
   zs2 = Z1(n2); zs3 = Z2(n2);
   printf("Bin for ss: %d (%d, %d)\n", n2, zs2, zs3);
   for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (nl2 = 0; nl2 <like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1],zl1);
         if (weight && ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_shear){
            if (test_zoverlap(zl1, zs2)*test_zoverlap(zl1, zs3)) {
               c_ng = cov_NG_gk_ss(ell[nl1],ell[nl2],zl1,zs2,zs3);
            }
            if (nl1==nl2){
               c_g = cov_G_gk_ss(ell[nl1],dell[nl1],zl1,zs2, zs3);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+zl1)+nl1;
         j=like.Ncl*n2+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],zl1,0,zs2,zs3,c_g,c_ng);
      }
   }
   fclose(F1);
}

// kk_kk
void run_cov_kk_kk(char *OUTFILE, char *PATH, double *ell, double *dell,int start)
{
   int nl1,nl2,i,j,weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   // sprintf(filename,"%scov/%s_%s_cov_kkkk_Nell%d_Ns%d_Ng%d_%d",PATH,survey.name, cmb.name,like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         if (ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_kappacmb){
            c_ng = cov_NG_kk_kk(ell[nl1],ell[nl2]);
         }
         if (nl1==nl2){
            c_g = cov_G_kk_kk(ell[nl1], dell[nl1]);
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+tomo.shear_Nbin)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+tomo.shear_Nbin)+nl2;
         fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],0,0,0,0,c_g,c_ng);
//         printf("l1, l2, Cg, Cng = %le, %le, %le, %le\n", ell[nl1], ell[nl2], c_g, c_ng);
      }
   }
   fclose(F1);
}

// kk_ks
void run_cov_kk_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int zs2, int start)
{
   int nl1,nl2,i,j, weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   // sprintf(filename,"%scov/%s_%s_cov_kkks_Nell%d_Ns%d_Ng%d_%d", PATH, survey.name, cmb.name, like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   printf("Source bin for ks: %d\n", zs2);
   for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         if (ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_kappacmb) {
            c_ng = cov_NG_kk_ks(ell[nl1],ell[nl2],zs2);
            if (nl1==nl2){
               c_g = cov_G_kk_ks(ell[nl1],dell[nl1], zs2);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+tomo.shear_Nbin)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+zs2)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],0,0,0,zs2,c_g,c_ng);
      }
   }
   fclose(F1);
}

// kk_ss
void run_cov_kk_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int n2, int start)
{
   int z2,z3,nl1,nl2,i,j,weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   z2 = Z1(n2); z3 = Z2(n2);
   printf("Bin for ss: %d (%d, %d)\n", n2, z2, z3);
   // sprintf(filename,"%scov/%s_%s_cov_kkss_Nell%d_Ns%d_Ng%d_%d",PATH,survey.name, cmb.name,like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         if (ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_shear){
            c_ng = cov_NG_kk_ss(ell[nl1],ell[nl2],z2,z3);
            if (nl1==nl2){
               c_g = cov_G_kk_ss(ell[nl1],dell[nl1],z2,z3);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+tomo.shear_Nbin)+nl1;
         j=like.Ncl*n2+nl2;
         fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],0,0,z2,z3,c_g,c_ng);
      }
   }
   fclose(F1);
}

// ks_ks
void run_cov_ks_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int zs1, int zs2,int start)
{
   int nl1,nl2,i,j, weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   // sprintf(filename,"%scov/%s_%s_cov_ksks_Nell%d_Ns%d_Ng%d_%d", PATH, survey.name, cmb.name, like.Ncl, tomo.shear_Nbin, tomo.clustering_Nbin, start);
   // printf("Saving to: %s\n",filename);
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   printf("Source bins for ks: %d, %d\n", zs1, zs2);
   for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         if (ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_kappacmb) {
            c_ng = cov_NG_ks_ks(ell[nl1],ell[nl2],zs1,zs2);
         }
         if (nl1==nl2){
            c_g = cov_G_ks_ks(ell[nl1],dell[nl1], zs1, zs2);
         }
         if ((ell[nl1]>like.lmax_kappacmb || ell[nl2]>like.lmax_kappacmb) && zs2!=zs1){
            c_g = 0;
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+zs1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+zs2)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],0,zs1,0,zs2,c_g,c_ng);
      }
   }
   fclose(F1);
}

// ks_ss
void run_cov_ks_ss(char *OUTFILE, char *PATH, double *ell, double *dell, int zs1, int n2, int start)
{
   int z2,z3,nl1,nl2,i,j,weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   z2 = Z1(n2); z3 = Z2(n2);
   printf("Source bin for ks: %d\n", zs1);
   printf("Bin for ss: %d (%d, %d)\n", n2, z2, z3);

   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         if (ell[nl1]<like.lmax_kappacmb && ell[nl2]<like.lmax_shear){
            c_ng = cov_NG_ks_ss(ell[nl1],ell[nl2],zs1,z2,z3);
            if (nl1==nl2){
               c_g = cov_G_ks_ss(ell[nl1],dell[nl1],zs1,z2,z3);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+zs1)+nl1;
         j=like.Ncl*n2+nl2;
         fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],0,zs1,z2,z3,c_g,c_ng);
      }
   }
   fclose(F1);
}

// ll_kk
void run_cov_ll_kk(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int start)
{
   int i,j, weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   // zl1 = ZL(n1); zs1 = ZS(n1);
   printf("Bin for ll: %d \n", n1);
   for (int nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (int nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1], n1);
         if (weight && ell[nl2]<like.lmax_kappacmb) {
            c_ng = cov_NG_gg_kk(ell[nl1],ell[nl2], n1, n1);
            if (nl1==nl2){
               c_g = cov_G_gg_kk(ell[nl1], dell[nl1], n1, n1);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+tomo.shear_Nbin)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],n1,n1,0,0,c_g,c_ng);
      }
   }
   fclose(F1);
}

// ll_ks
void run_cov_ll_ks(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int zs2, int start)
{
   int i,j,weight;
   double c_ng, c_g;
   FILE *F1;
   char filename[300];
   sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
   F1 =fopen(filename,"w");
   // zl1 = ZL(n1); zs1 = ZS(n1);
   printf("Bin for ll: %d\n", n1);
   printf("Bin for ks: %d\n", zs2);
   for (int nl1 = 0; nl1 < like.Ncl; nl1 ++){
      for (int nl2 = 0; nl2 < like.Ncl; nl2 ++){
         c_ng = 0.; c_g = 0.;
         weight = test_kmax(ell[nl1], n1);
         if (weight && ell[nl2]<like.lmax_kappacmb) {
            if (test_zoverlap(n1, zs2)) {
               c_ng = cov_NG_gg_ks(ell[nl1],ell[nl2], n1, n1, zs2);
            }
            if (nl1==nl2){
               c_g = cov_G_gg_ks(ell[nl1], dell[nl1], n1, n1, zs2);
            }
         }
         i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1;
         j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+tomo.clustering_Nbin+zs2)+nl2;
         fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],n1,n1,0,zs2,c_g,c_ng);
      }
   }
   fclose(F1);
}

// ll_lk
void run_cov_ll_lk(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int zl2,int start)
{
  int z1,z2,zl,zs,nl1,nl2,i,j,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl_1 = %d \n", n1);
  // zl = ZL(n2); zs = ZS(n2);
  zl = zl2;
  printf("N_tomo_2 = %d\n", zl2);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],z1)*test_kmax(ell[nl2],zl);
      if (weight){
        if((z1 == zl)) {c_ng = cov_NG_gg_gk(ell[nl1],ell[nl2],z1,z2,zl);}
        if (nl1 == nl2){
          c_g =  cov_G_gg_gk(ell[nl1],dell[nl1],z1,z2,zl);
        }
      }
      i=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1;
      j=like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Nbin+zl)+nl2;
      fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell[nl2],z1,z2,zl,0,c_g,c_ng);
    }
  }
  fclose(F1);
}

int main(int argc, char** argv)
{

  int i,l,m,n,o,s,p,nl1,t,k;
  char OUTFILE[400],filename[400],arg1[400],arg2[400];

  int N_scenarios=2;
  double area_table[2]={18000.0,2000.0};
  double nsource_table[2]={43.0,51.0};
  double nlens_table[2]={50.0,66.0};

  char survey_designation[2][200]={"WFIRSTwidexSO5","WFIRSTstdxSO5"};

  char source_zfile[2][400]={"zdistri_WFIRST_LSST_lensing_fine_bin","zdistri_WFIRST_LSST_lensing_fine_bin"};

#ifdef ONESAMPLE
  char lens_zfile[2][400]={"zdistri_WFIRST_LSST_lensing_fine_bin","zdistri_WFIRST_LSST_lensing_fine_bin"};
  nlens_table[0] = nsource_table[0];
  nlens_table[1] = nsource_table[1];
#else
  char lens_zfile[2][400]={"zdistri_WFIRST_LSST_clustering_fine_bin","zdistri_WFIRST_LSST_clustering_fine_bin"};
#endif

  int hit=atoi(argv[1]);
  Ntable.N_a=100;
  k=1;

  t = atoi(argv[2]);

  //RUN MODE setup
  init_cosmo_runmode("emu");
  // init_binning_fourier(20,30.0,3000.0,3000.0,21.0,10,10);
  init_binning_fourier(15,20.0,3000.0,3000.0,21.0,10,10);
  init_survey(survey_designation[t],nsource_table[t],nlens_table[t],area_table[t]);
  sprintf(arg1,"zdistris/%s",source_zfile[t]);
  sprintf(arg2,"zdistris/%s",lens_zfile[t]);
  init_galaxies(arg1,arg2,"none","none","source_std","WF_SN10");
  init_IA("none","GAMA");
  init_probes("6x2pt");

  init_cmb("so_Y5");

  //set l-bins for shear, ggl, clustering, clusterWL
  double logdl=(log(like.lmax)-log(like.lmin))/like.Ncl;
  double *ell, *dell, *ell_Cluster, *dell_Cluster;
  ell=create_double_vector(0,like.Ncl-1);
  dell=create_double_vector(0,like.Ncl-1);
  ell_Cluster=create_double_vector(0,Cluster.lbin-1);
  dell_Cluster=create_double_vector(0,Cluster.lbin-1);
  int j=0;
  for(i=0;i<like.Ncl;i++){
    ell[i]=exp(log(like.lmin)+(i+0.5)*logdl);
    dell[i]=exp(log(like.lmin)+(i+1)*logdl)-exp(log(like.lmin)+(i*logdl));
    if(ell[i]<like.lmax_shear) printf("%le\n",ell[i]);
  }

  printf("----------------------------------\n");
  sprintf(survey.name,"%s_area%le_ng%le_nl%le",survey_designation[t],survey.area,survey.n_gal,survey.n_lens);
  printf("area: %le n_source: %le n_lens: %le\n",survey.area,survey.n_gal,survey.n_lens);

  // sprintf(covparams.outdir,"/home/u17/timeifler/covparallel/");
#ifdef ONESAMPLE
  sprintf(covparams.outdir,"out_cov_wfirstxso_1sample/");
#else
  sprintf(covparams.outdir,"out_cov_wfirstxso/");
  //sprintf(covparams.outdir,"/halo_nobackup/cosmos/teifler/covparallel/");
#endif

  printf("----------------------------------\n");
  sprintf(OUTFILE,"%s_ssss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.shear_Npowerspectra; l++){
    for (m=l;m<tomo.shear_Npowerspectra; m++){
      if(k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_ss_ss(OUTFILE,covparams.outdir,ell,dell,l,m,k);
      }
      k=k+1;
      //printf("%d\n",k);
    }
  }

  sprintf(OUTFILE,"%s_lsls_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.ggl_Npowerspectra; l++){
    for (m=l;m<tomo.ggl_Npowerspectra; m++){
      if(k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_ls_ls(OUTFILE,covparams.outdir,ell,dell,l,m,k);
      }
     //printf("%d\n",k);
      k=k+1;
    }
  }
  sprintf(OUTFILE,"%s_llll_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Npowerspectra; l++){ //auto bins only for now!
    for (m=l;m<tomo.clustering_Npowerspectra; m++){
      if(k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_ll_ll(OUTFILE,covparams.outdir,ell,dell,l,m,k);
      }
      k=k+1;
      //printf("%d %d %d\n",l,m,k);
    }
  }
  sprintf(OUTFILE,"%s_llss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Npowerspectra; l++){
    for (m=0;m<tomo.shear_Npowerspectra; m++){
      if(k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_ll_ss(OUTFILE,covparams.outdir,ell,dell,l,m,k);
      }
      k=k+1;
      //printf("%d\n",k);
    }
  }
  sprintf(OUTFILE,"%s_llls_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Npowerspectra; l++){
    for (m=0;m<tomo.ggl_Npowerspectra; m++){
      if(k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_ll_ls(OUTFILE,covparams.outdir,ell,dell,l,m,k);
      }
      k=k+1;
      //printf("%d\n",k);
    }
  }
  sprintf(OUTFILE,"%s_lsss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.ggl_Npowerspectra; l++){
    for (m=0;m<tomo.shear_Npowerspectra; m++){
      if(k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_ls_ss(OUTFILE,covparams.outdir,ell,dell,l,m,k);
      }
      k=k+1;
      //printf("%d\n",k);
    }
  }

 // lk_lk
  sprintf(OUTFILE,"%s_lklk_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Nbin; l++){
     for (m=l;m<tomo.clustering_Nbin; m++){
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_lk_lk(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        k=k+1;
     }
  }

  printf("lklk %d\n",k);

  // lk_ls
  sprintf(OUTFILE,"%s_lkls_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Nbin; l++){
     for (m=0;m<tomo.ggl_Npowerspectra; m++){
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_lk_ls(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        k=k+1;
     }
  }
  printf("lkls %d\n",k);

  // lk_kk
  sprintf(OUTFILE,"%s_lkkk_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Nbin; l++){
     if (k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_lk_kk(OUTFILE,covparams.outdir,ell,dell,l,k);
     }
     // printf("%d\n",k);
     k=k+1;
  }
  printf("lkkk %d\n",k);

  // lk_ks
  sprintf(OUTFILE,"%s_lkks_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Nbin; l++){
     for (m=0;m<tomo.shear_Nbin;m++) {
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_lk_ks(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        // printf("%d\n",k);
        k=k+1;
     }
  }
  printf("lkks %d\n",k);

  // lk_ss
  sprintf(OUTFILE,"%s_lkss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Nbin; l++){
     for (m=0;m<tomo.shear_Npowerspectra; m++){
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_lk_ss(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        k=k+1;
        // printf("%d\n",k);
     }
  }
  printf("lkss %d\n",k);

  // ls_kk
  sprintf(OUTFILE,"%s_lskk_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.ggl_Npowerspectra; l++){
     if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ls_kk(OUTFILE,covparams.outdir,ell,dell,l,k);
        }
    // printf("%d\n",k);
     k=k+1;
  }
  printf("lskk%d\n",k);

  // ls_ks
  sprintf(OUTFILE,"%s_lsks_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.ggl_Npowerspectra; l++){
     for (m=0;m<tomo.shear_Nbin;m++) {
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ls_ks(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        // printf("%d\n",k);
        k=k+1;
     }
  }
  printf("lsks %d\n",k);

  // kk_kk
  sprintf(OUTFILE,"%s_kkkk_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  if (k==hit){
    sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
    run_cov_kk_kk(OUTFILE,covparams.outdir,ell,dell,k);
  }
  // printf("%d\n",k);
  k=k+1;
  printf("kkkk %d\n",k);

  // kk_ks
  sprintf(OUTFILE,"%s_kkks_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.shear_Nbin; l++){
     if (k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_kk_ks(OUTFILE,covparams.outdir,ell,dell,l,k);
     }
     // printf("%d\n",k);
     k=k+1;
  }
  printf("kkks %d\n",k);

  // kk_ss
  sprintf(OUTFILE,"%s_kkss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.shear_Npowerspectra; l++){
     if (k==hit){
      sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
      run_cov_kk_ss(OUTFILE,covparams.outdir,ell,dell,l,k);
     }
     k=k+1;
     // printf("%d\n",k);
  }
  printf("kkss %d\n",k);

  // ks_ks
  sprintf(OUTFILE,"%s_ksks_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.shear_Nbin; l++){
     for (m=l;m<tomo.shear_Nbin; m++){
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ks_ks(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        // printf("%d\n",k);
        k=k+1;
     }
  }
  printf("ksks %d\n",k);

 // ks_ss
 sprintf(OUTFILE,"%s_ksss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
 for (l=0;l<tomo.shear_Nbin; l++){
   for (m=0;m<tomo.shear_Npowerspectra; m++){
     if (k==hit){
        sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
        run_cov_ks_ss(OUTFILE,covparams.outdir,ell,dell,l,m,k);
     }
     k=k+1;
     // printf("%d\n",k);
    }
  }
  printf("ksss %d\n",k);

  // ll_kk
  sprintf(OUTFILE,"%s_llkk_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Npowerspectra; l++){
     if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ll_kk(OUTFILE,covparams.outdir,ell,dell,l,k);
        }
    // printf("%d\n",k);
     k=k+1;
  }
  printf("llkk %d\n",k);

  // ll_ks
  sprintf(OUTFILE,"%s_llks_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Npowerspectra; l++){
     for (m=0;m<tomo.shear_Nbin;m++) {
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ll_ks(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        // printf("%d\n",k);
        k=k+1;
     }
  }
  printf("llks %d\n",k);

  // ll_lk
  sprintf(OUTFILE,"%s_lllk_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
  for (l=0;l<tomo.clustering_Npowerspectra; l++){
     for (m=0;m<tomo.clustering_Nbin;m++) {
        if (k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          run_cov_ll_lk(OUTFILE,covparams.outdir,ell,dell,l,m,k);
        }
        // printf("%d\n",k);
        k=k+1;
     }
  }
  printf("lllk %d\n",k);

  printf("number of cov blocks for parallelization: %d\n",k-1);
  printf("-----------------\n");
  printf("PROGRAM EXECUTED\n");
  printf("-----------------\n");
  return 0;
}

ext: c

max_stars_repo_name: CosmoLike/LSSTxSO

max_issues_repo_name: CosmoLike/LSSTxSO

max_forks_repo_name: CosmoLike/LSSTxSO


text:
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <immintrin.h>
#include <omp.h>
#include <Python.h>
#include <numpy/arrayobject.h>

#include <cblas.h>

static PyObject *smmprod_func(PyObject* self, PyObject* args)
{
    PyObject *A_obj = NULL;
    PyObject *A_arr = NULL;
    PyArrayObject *A = NULL;

    PyObject *B_obj = NULL;
    PyObject *B_arr = NULL;
    PyArrayObject *B = NULL;

    PyObject *Omega_obj = NULL;
    PyObject *Omega_row_arr = NULL;
    PyObject *Omega_col_arr = NULL;
    PyObject *Omega_seq = NULL;
    PyArrayObject *Omega_row = NULL;
    PyArrayObject *Omega_col = NULL;

    PyObject *Out_obj = NULL;
    PyObject *Out_arr = NULL;
    PyArrayObject *Out = NULL;

    int result = PyArg_ParseTuple(
        args,
        "OOOO",
        &A_obj,
        &B_obj,
        &Omega_obj,
        &Out_obj
    );

    if (!result) {
        return NULL;
    }

    A_arr = PyArray_FROM_OTF(A_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
    if (A_arr == NULL) {
        return NULL;
    }
    A = (PyArrayObject*) A_arr;

    B_arr = PyArray_FROM_OTF(B_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
    if (B_arr == NULL) {
        return NULL;
    }
    B = (PyArrayObject*) B_arr;

    Omega_seq = PySequence_Fast(Omega_obj, "expected a sequence");

    Omega_row_arr = PyArray_FROM_OTF(PySequence_Fast_GET_ITEM(Omega_seq, 0), NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
    if (Omega_row_arr == NULL) {
        return NULL;
    }
    Omega_row = (PyArrayObject*) Omega_row_arr;

    Omega_col_arr = PyArray_FROM_OTF(PySequence_Fast_GET_ITEM(Omega_seq, 1), NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
    if (Omega_col_arr == NULL) {
        return NULL;
    }
    Omega_col = (PyArrayObject*) Omega_col_arr;

    Out_arr = PyArray_FROM_OTF(Out_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
    if (Out_arr == NULL) {
        return NULL;
    }
    Out = (PyArrayObject*) Out_arr;

    int ndimA = PyArray_NDIM(A);
    if (ndimA != 2) {
        PyErr_SetString(PyExc_ValueError,
            "A should be 2-dimensional");
        return NULL;
    }

    int ndimB = PyArray_NDIM(B);
    if (ndimB != 2) {
        PyErr_SetString(PyExc_ValueError,
            "B should be 2-dimensional");
        return NULL;
    }

//    int ndimOmega = PyArray_NDIM(B);
//    if (ndimOmega != 2) {
//        PyErr_SetString(PyExc_ValueError,
//            "Omega should be 2-dimensional");
//        return NULL;
//    }

    npy_intp *dimsA = PyArray_DIMS(A);
    npy_intp *dimsB = PyArray_DIMS(B);
    npy_intp *dimsOmega_row = PyArray_DIMS(Omega_row);
    npy_intp *dimsOmega_col = PyArray_DIMS(Omega_col);

    if (dimsOmega_row[0] != dimsOmega_col[0]) {
        PyErr_SetString(PyExc_ValueError,
            "index arrays should have equal length");
        return NULL;
    }

    if (dimsA[1] != dimsB[0] ) {
        PyErr_SetString(PyExc_ValueError,
            "Inner dimensions of A and B do not match");
        return NULL;
    }

    double *bufferA = PyArray_DATA(A);
    double *bufferB = PyArray_DATA(B);
    double *bufferOmega_row = PyArray_DATA(Omega_row);
    double *bufferOmega_col = PyArray_DATA(Omega_col);
    double *bufferOut = PyArray_DATA(Out);

    //int m = dimsA[0];
    int k = dimsA[1];
    int n = dimsB[1];
    int j = dimsOmega_row[0];

    // Remove this pragma for single threaded execution
    #pragma omp parallel for
    for (int i=0; i<j; i++) {
        double *iA = bufferA + k * (int)(bufferOmega_row[i]);
        double *iB = bufferB + (int)(bufferOmega_col[i]);

        /*
         * The plain C version, without blas
        double val = 0;
        for (int x=0; x<k; x++) {
            val += *iA * *iB;
            iA += 1;
            iB += n;
        }
        *(bufferOut++) = val;
        */
        bufferOut[i] = cblas_ddot(k, iA, 1, iB, n);
    }

    Py_DECREF(A_arr);
    Py_DECREF(B_arr);
    Py_DECREF(Omega_row_arr);
    Py_DECREF(Omega_col_arr);
    Py_DECREF(Omega_seq);
    Py_DECREF(Out);
    Py_RETURN_NONE;
}

static PyMethodDef smmprod_methods[] =
{
     {"smmprod", smmprod_func, METH_VARARGS, "calculate fast smmprod in pure C"},
     {NULL, NULL, 0, NULL}
};

#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef moduledef = {
    PyModuleDef_HEAD_INIT,
    "_smmprod",
    NULL,
    -1,
    smmprod_methods
};

PyObject *
PyInit__smmprod(void)
{
    import_array();
    PyObject *module = PyModule_Create(&moduledef);
    return module;
}
#else
void init_smmprod(void)
{
    import_array();
    Py_InitModule("_smmprod", smmprod_methods);
    return;
}

#endif

max_stars_repo_name: clemenshage/grslra

max_issues_repo_name: clemenshage/grslra

max_forks_repo_name: clemenshage/grslra


text:
/**
 * @file dereverberation.h
 * @brief Single- and multi-channel dereverberation base on linear prediction in the subband domain.
 * @author John McDonough
 */

#ifndef DEREVERBERATION_H
#define DEREVERBERATION_H

#include <stdio.h>
#include <assert.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_linalg.h>
#include "common/jexception.h"

#include "stream/stream.h"
#include "feature/feature.h"

/*
#include <Eigen/Cholesky>
#include <iostream>
#include <Eigen/Dense>
*/

// ----- definition for class `SingleChannelWPEDereverberationFeature' -----
//
class SingleChannelWPEDereverberationFeature : public VectorComplexFeatureStream {
  typedef vector<gsl_vector_complex*>			Samples_;
  typedef Samples_::iterator				SamplesIterator_;

 public:
  SingleChannelWPEDereverberationFeature(VectorComplexFeatureStreamPtr& samples, unsigned lowerN, unsigned upperN, unsigned iterationsN = 2, double loadDb = -20.0, double bandWidth = 0.0,  double sampleRate = 16000.0, const String& nm = "SingleChannelWPEDereverberationFeature");

  ~SingleChannelWPEDereverberationFeature();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  unsigned estimate_filter(int start_frame_no = 0, int frame_num = -1);
  void reset_filter();
  void next_speaker();
  void print_objective_func(int subbandX){ printing_subbandX_ = subbandX;}

#ifdef ENABLE_LEGACY_BTK_API
  void nextSpeaker(){ next_speaker(); }
#endif

private:
  static const double					subband_floor_;

  void fill_buffer_(int start_frame_no, int frame_num);
  void estimate_Gn_();
  void calc_Rr_(unsigned subbandX);
  void calc_Thetan_();
  void load_R_();
  unsigned set_band_width_(double bandWidth, double sampleRate);

  const gsl_vector_complex* get_lags_(unsigned subbandX, unsigned sampleX);
  VectorComplexFeatureStreamPtr				samples_;

  const unsigned					lowerN_;
  const unsigned					upperN_;
  const unsigned					predictionN_;
  const unsigned					iterationsN_;
  bool							estimated_;
  unsigned						framesN_; // no frames used for filter estimation
  const double						load_factor_;
  const unsigned					lower_bandWidthN_;
  const unsigned					upper_bandWidthN_;

  Samples_						yn_; // buffer to keep observations
  gsl_matrix*						thetan_;
  gsl_vector_complex**					gn_;
  gsl_matrix_complex*					R_;
  gsl_vector_complex*					r_;
  gsl_vector_complex*					lag_samples_;
  int                                                   printing_subbandX_;
};

typedef Inherit<SingleChannelWPEDereverberationFeature, VectorComplexFeatureStreamPtr> SingleChannelWPEDereverberationFeaturePtr;

// ----- definition for class `MultiChannelWPEDereverberation' -----
//
class MultiChannelWPEDereverberation : public Countable {
  typedef vector<VectorComplexFeatureStreamPtr>		SourceList_;
  typedef SourceList_::iterator				SourceListIterator_;

  typedef vector<gsl_vector_complex*>			FrameBrace_;
  typedef FrameBrace_::iterator				FrameBraceIterator_;

  typedef vector<FrameBrace_>				FrameBraceList_;
  typedef FrameBraceList_::iterator			FrameBraceListIterator_;

 public:
  MultiChannelWPEDereverberation(unsigned subbandsN, unsigned channelsN, unsigned lowerN, unsigned upperN, unsigned iterationsN = 2, double loadDb = -20.0, double bandWidth = 0.0, double diagonal_bias = 0.0, double sampleRate = 16000.0);

  ~MultiChannelWPEDereverberation();

  void reset();
  unsigned size() const { return subbandsN_; }
  void set_input(VectorComplexFeatureStreamPtr& samples);
  const gsl_vector_complex* get_output(unsigned channelX);
  gsl_vector_complex** calc_every_channel_output(int frame_no = -5);
  unsigned estimate_filter(int start_frame_no = 0, int frame_num = -1);
  void reset_filter();
  void next_speaker();
  int  frame_no() const { return frame_no_; }
  void print_objective_func(int subbandX){ printing_subbandX_ = subbandX;}

#ifdef ENABLE_LEGACY_BTK_API
  void setInput(VectorComplexFeatureStreamPtr& samples){ set_input(samples); }
  const gsl_vector_complex* getOutput(unsigned channelX, int frame_no = -5){ return get_output(channelX); }
  void nextSpeaker(){ next_speaker(); }
#endif

private:
  static const double					subband_floor_;

  void fill_buffer_(int start_frame_no, int frame_num);
  void estimate_Gn_();
  void calc_Rr_(unsigned subbandX);
  void calc_Thetan_();
  void load_R_();
  unsigned set_band_width_(double bandWidth, double sampleRate);

  void increment_() { frame_no_++; }
  const gsl_vector_complex* get_lags_(unsigned subbandX, unsigned sampleX);

  SourceList_						sources_;
  const unsigned					subbandsN_;
  const unsigned					channelsN_;
  const unsigned					lowerN_;
  const unsigned					upperN_;
  const unsigned					predictionN_;
  const unsigned					iterationsN_;
  const unsigned					totalPredictionN_;

  bool							estimated_;
  unsigned						framesN_;
  const double						load_factor_;
  const unsigned					lower_bandWidthN_;
  const unsigned					upper_bandWidthN_;

  FrameBraceList_					frames_;
  gsl_matrix**						thetan_;
  gsl_vector_complex***					Gn_;
  gsl_matrix_complex**					R_;
  gsl_vector_complex**					r_;
  gsl_vector_complex*					lag_samples_;
  gsl_vector_complex**					output_;

  const int						initial_frame_no_;
  int							frame_no_;

  const double                                          diagonal_bias_;
  int                                                   printing_subbandX_;
};

typedef refcountable_ptr<MultiChannelWPEDereverberation> MultiChannelWPEDereverberationPtr;

// ----- definition for class `MultiChannelWPEDereverberationFeature' -----
//
class MultiChannelWPEDereverberationFeature : public VectorComplexFeatureStream {
public:
  MultiChannelWPEDereverberationFeature(MultiChannelWPEDereverberationPtr& source, unsigned channelX, unsigned primaryChannelX = 0, const String& nm = "MultiChannelWPEDereverberationFeature");

  ~MultiChannelWPEDereverberationFeature();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

private:
  MultiChannelWPEDereverberationPtr			source_;
  const unsigned					channelX_;
  const unsigned					primaryChannelX_; // compute the dereverbed output only when channelX_ == primaryChannelX_. Otherwise copy the precomputed output
};

typedef Inherit<MultiChannelWPEDereverberationFeature, VectorComplexFeatureStreamPtr> MultiChannelWPEDereverberationFeaturePtr;

#endif // DEREVERBERATION_H

ext: h


text:
#pragma once

#include <iosfwd>
#include <map>
#include <memory>
#include <vector>

#include <gsl/gsl>

namespace sframe {

struct openssl_error : std::runtime_error
{
  openssl_error();
};

struct unsupported_ciphersuite_error : std::runtime_error
{
  unsupported_ciphersuite_error();
};

struct authentication_error : std::runtime_error
{
  authentication_error();
};

struct buffer_too_small_error : std::runtime_error
{
  using parent = std::runtime_error;
  using parent::parent;
};

struct invalid_parameter_error : std::runtime_error
{
  using parent = std::runtime_error;
  using parent::parent;
};

enum class CipherSuite : uint16_t
{
  AES_CM_128_HMAC_SHA256_4 = 1,
  AES_CM_128_HMAC_SHA256_8 = 2,
  AES_GCM_128_SHA256 = 3,
  AES_GCM_256_SHA512 = 4,
};

constexpr size_t max_overhead = 17 + 16;

using bytes = std::vector<uint8_t>;
using input_bytes = gsl::span<const uint8_t>;
using output_bytes = gsl::span<uint8_t>;

std::ostream&
operator<<(std::ostream& str, const input_bytes data);

using KeyID = uint64_t;
using Counter = uint64_t;

class SFrame
{
protected:
  CipherSuite suite;

  SFrame(CipherSuite suite_in);
  virtual ~SFrame();

  struct KeyState
  {
    static KeyState from_base_key(CipherSuite suite, const bytes& base_key);

    bytes key;
    bytes salt;
    Counter counter;
  };

  output_bytes _protect(KeyID key_id,
                        output_bytes ciphertext,
                        input_bytes plaintext);
  output_bytes _unprotect(output_bytes ciphertext, input_bytes plaintext);

  virtual KeyState& get_state(KeyID key_id) = 0;
};

class Context : public SFrame
{
public:
  Context(CipherSuite suite);

  void add_key(KeyID kid, const bytes& key);

  output_bytes protect(KeyID key_id,
                       output_bytes ciphertext,
                       input_bytes plaintext);
  output_bytes unprotect(output_bytes plaintext, input_bytes ciphertext);

private:
  std::map<KeyID, KeyState> state;

  KeyState& get_state(KeyID key_id) override;
};

class MLSContext : public SFrame
{
public:
  using EpochID = uint64_t;
  using SenderID = uint32_t;

  MLSContext(CipherSuite suite_in, size_t epoch_bits_in);

  void add_epoch(EpochID epoch_id, const bytes& sframe_epoch_secret);
  void purge_before(EpochID keeper);

  output_bytes protect(EpochID epoch_id,
                       SenderID sender_id,
                       output_bytes ciphertext,
                       input_bytes plaintext);
  output_bytes unprotect(output_bytes plaintext, input_bytes ciphertext);

private:
  const size_t epoch_bits;
  const size_t epoch_mask;

  struct EpochKeys
  {
    const EpochID full_epoch;
    const bytes sframe_epoch_secret;
    std::map<SenderID, KeyState> sender_keys;

    EpochKeys(EpochID full_epoch_in, bytes sframe_epoch_secret_in);
    KeyState& get(CipherSuite suite, SenderID sender_id);
  };

  std::vector<std::unique_ptr<EpochKeys>> epoch_cache;
  KeyState& get_state(KeyID key_id) override;
};

} // namespace sframe

max_stars_repo_name: fanglinliu/sframe

max_issues_repo_name: fanglinliu/sframe

max_forks_repo_name: fanglinliu/sframe