text:
/**
 * @file stream.h
 * @brief Representation of feature streams.
 * @author John McDonough.
 */

#ifndef STREAM_H
#define STREAM_H

#include <gsl/gsl_vector.h>
#include <gsl/gsl_complex.h>
#include "common/refcount.h"

// ----- interface class for 'FeatureStream' -----
//
template <typename Type, typename item_type>
class FeatureStream : public Countable {
 public:
  virtual ~FeatureStream();

  const String& name() const { return name_; }
  unsigned      size() const { return size_; }

  virtual const Type* next(int frame_no = -5) = 0;
  const Type* current() {
    if (frame_no_ < 0)
      throw jconsistency_error("Frame index (%d) < 0.", frame_no_);
    return next(frame_no_);
  }

  bool is_end(){ return is_end_; };

  virtual void reset() {
    frame_no_ = frame_reset_no_;
    is_end_ = false;
  }

  virtual int frame_no() const { return frame_no_; }
  size_t itemsize() { return sizeof(item_type); };

 protected:
  FeatureStream(unsigned sz, const String& nm);
  void gsl_vector_set_(Type *vector, int index, item_type value);
  void increment_() { frame_no_++; }

  const int					frame_reset_no_;
  const unsigned				size_;
  int						frame_no_; /*!< lapse time after reset() */
  Type*						vector_;
  bool                                          is_end_;

 private:
  const String					name_;
};

// ----- declare partial specializations of 'FeatureStream' -----
//
template<> FeatureStream<gsl_vector_char, char>::FeatureStream(unsigned sz, const String& nm);
template<> FeatureStream<gsl_vector_short, short>::FeatureStream(unsigned sz, const String& nm);
template<> FeatureStream<gsl_vector_float, float>::FeatureStream(unsigned sz, const String& nm);
template<> FeatureStream<gsl_vector, double>::FeatureStream(unsigned sz, const String& nm);
template<> FeatureStream<gsl_vector_complex, gsl_complex>::FeatureStream(unsigned sz, const String& nm);
template<> FeatureStream<gsl_vector_char, char>::~FeatureStream();
template<> FeatureStream<gsl_vector_short, short>::~FeatureStream();
template<> FeatureStream<gsl_vector_float, float>::~FeatureStream();
template<> FeatureStream<gsl_vector, double>::~FeatureStream();
template<> FeatureStream<gsl_vector_complex, gsl_complex>::~FeatureStream();
template<> void FeatureStream<gsl_vector_char, char>::gsl_vector_set_(gsl_vector_char *vector, int index, char value);
template<> void FeatureStream<gsl_vector_short, short>::gsl_vector_set_(gsl_vector_short *vector, int index, short value);
template<> void FeatureStream<gsl_vector_float, float>::gsl_vector_set_(gsl_vector_float *vector, int index, float value);
template<> void FeatureStream<gsl_vector, double>::gsl_vector_set_(gsl_vector *vector, int index, double value);
template<> void FeatureStream<gsl_vector_complex, gsl_complex>::gsl_vector_set_(gsl_vector_complex *vector, int index, gsl_complex value);

typedef FeatureStream<gsl_vector_char, char>		VectorCharFeatureStream;
typedef FeatureStream<gsl_vector_short, short>		VectorShortFeatureStream;
typedef FeatureStream<gsl_vector_float, float>		VectorFloatFeatureStream;
typedef FeatureStream<gsl_vector, double>		VectorFeatureStream;
typedef FeatureStream<gsl_vector_complex, gsl_complex>	VectorComplexFeatureStream;

typedef refcountable_ptr<VectorCharFeatureStream>	VectorCharFeatureStreamPtr;
typedef refcountable_ptr<VectorShortFeatureStream>	VectorShortFeatureStreamPtr;
typedef refcountable_ptr<VectorFloatFeatureStream>	VectorFloatFeatureStreamPtr;
typedef refcountable_ptr<VectorFeatureStream>		VectorFeatureStreamPtr;
typedef refcountable_ptr<VectorComplexFeatureStream>	VectorComplexFeatureStreamPtr;

#endif


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

#pragma once

#include <gsl\gsl>

#include "Data\Data.h"
#include "Data\Pose.h"
#include "opencv\cv.h"
#include "MageSettings.h"

namespace mage
{
    struct VOIKeyframe
    {
        const VolumeOfInterestSettings& VoiSettings;
        cv::Vec3f WorldPosition;
        cv::Vec3f Forward;
        cv::Vec3f Right;
        cv::Vec3f Up;
        float NearDepth;
        float FarDepth;

        // Cached values for performance.
        cv::Vec3f Centroid;
        float DistanceAlphaToXi;
        float ModifiedDistanceAlphaToOmega;

        VOIKeyframe(const Pose& pose, const Depth& depth, const VolumeOfInterestSettings& voiSettings);
        float TeardropScore(const cv::Vec3f&) const;
    };

    bool CalculateVolumeOfInterest(gsl::span<const VOIKeyframe> keyframes, const VolumeOfInterestSettings& voiSettings, AxisAlignedVolume& voi);
}

ext: h

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#ifndef Single_Pin_Conduction_h
#define Single_Pin_Conduction_h

#include <vector>

// Trilinos includes
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_SerialDenseSolver.hpp"
#include "Teuchos_SerialDenseVector.hpp"

// vendored includes
#include <gsl/gsl>

namespace enrico {

//===========================================================================//
/*!
 * \class Single_Pin_Conduction
 * \brief Solve thermal conduction in single fuel pin.
 *
 * This class implements a simple 1-D radial thermal conduction equation
 * for each axial level of a fuel pin.  Axial conduction is ignored and there
 * is currently no gap conductance model.  This model should only be used for
 * very approximate qualitative behavior.
 */
//===========================================================================//

class Single_Pin_Conduction {
public:
  //@{
  //! Typedefs
  using RCP_PL = Teuchos::RCP<Teuchos::ParameterList>;
  using Matrix = Teuchos::SerialDenseMatrix<int, double>;
  using Vector = Teuchos::SerialDenseVector<int, double>;
  using Solver = Teuchos::SerialDenseSolver<int, double>;
  //@}

private:
  // >>> DATA

  // Geometry
  double d_fuel_radius;
  double d_clad_radius;
  double d_delta_r_fuel;
  double d_delta_r_clad;
  std::vector<double> d_delta_z;

  // Thermal conductivities
  double d_k_fuel;
  double d_k_clad;

public:
  // Constructor
  Single_Pin_Conduction(RCP_PL& parameters, const std::vector<double>& delta_z);

  // Set fuel radius (cm)
  void set_fuel_radius(double r)
  {
    Expects(r > 0.0);
    Expects(r < 2.0);
    d_fuel_radius = r;
  }

  // Set clad radius (cm)
  void set_clad_radius(double r)
  {
    Expects(r > 0.0);
    Expects(r < 2.0);
    d_clad_radius = r;
  }

  // Solve conduction equation at each axial level
  void solve(const std::vector<double>& power,
             const std::vector<double>& channel_temperature,
             std::vector<double>& fuel_temperature);
};

//---------------------------------------------------------------------------//
} // end namespace enrico

//---------------------------------------------------------------------------//
#endif // Single_Pin_Conduction_h

//---------------------------------------------------------------------------//
// end of Single_Pin_Conduction.h
//---------------------------------------------------------------------------//

max_stars_repo_name: pshriwise/enrico

max_issues_repo_name: pshriwise/enrico

max_forks_repo_name: pshriwise/enrico


text:
#ifndef XYCO_RUNTIME_UTILS_SELECT_H
#define XYCO_RUNTIME_UTILS_SELECT_H

#include <gsl/pointers>
#include <utility>

#include "runtime/runtime.h"
#include "type_wrapper.h"

namespace xyco::runtime {
template <typename T1, typename T2>
class SelectFuture
    : public Future<std::variant<TypeWrapper<T1>, TypeWrapper<T2>>> {
  using CoOutput = std::variant<TypeWrapper<T1>, TypeWrapper<T2>>;

 public:
  auto poll(Handle<void> self) -> Poll<CoOutput> override {
    if (!ready_) {
      ready_ = true;
      auto [wrapper1, wrapper2] = std::pair<Future<void>, Future<void>>(
          future_wrapper<T1, 0>(std::move(futures_.first)),
          future_wrapper<T2, 1>(std::move(futures_.second)));
      wrappers_ = {wrapper1.get_handle(), wrapper2.get_handle()};
      RuntimeCtx::get_ctx()->spawn(std::move(wrapper1));
      RuntimeCtx::get_ctx()->spawn(std::move(wrapper2));
      return Pending();
    }

    if (result_.index() == 0) {
      RuntimeCtx::get_ctx()->cancel_future(wrappers_.second);
      return Ready<CoOutput>{
          CoOutput(std::in_place_index<0>, std::get<0>(result_))};
    }

    RuntimeCtx::get_ctx()->cancel_future(wrappers_.first);
    return Ready<CoOutput>{
        CoOutput(std::in_place_index<1>, std::get<1>(result_))};
  }

  SelectFuture(Future<T1> &&future1, Future<T2> &&future2)
      : Future<CoOutput>(nullptr),
        ready_(false),
        result_(std::in_place_index<2>, true),
        futures_(std::move(future1), std::move(future2)) {}

 private:
  template <typename T, int Index>
  auto future_wrapper(Future<T> future) -> Future<void>
  requires(!std::is_same_v<T, void>) {
    auto result = co_await future;
    if (result_.index() == 2) {
      result_ = std::variant<TypeWrapper<T1>, TypeWrapper<T2>, bool>(
          std::in_place_index<Index>, TypeWrapper<T>{result});
      RuntimeCtx::get_ctx()->register_future(this);
    }
  }

  template <typename T, int Index>
  auto future_wrapper(Future<T> future) -> Future<void>
  requires(std::is_same_v<T, void>) {
    co_await future;
    if (result_.index() == 2) {
      result_ = std::variant<TypeWrapper<T1>, TypeWrapper<T2>, bool>(
          std::in_place_index<Index>, TypeWrapper<T>());
      RuntimeCtx::get_ctx()->register_future(this);
    }
  }

  bool ready_;
  std::variant<TypeWrapper<T1>, TypeWrapper<T2>, bool> result_;
  std::pair<Future<T1>, Future<T2>> futures_;
  std::pair<Handle<PromiseBase>, Handle<PromiseBase>> wrappers_;
};

template <typename T1, typename T2>
auto select(Future<T1> future1, Future<T2> future2)
    -> Future<std::variant<TypeWrapper<T1>, TypeWrapper<T2>>> {
  co_return co_await SelectFuture<T1, T2>(std::move(future1),
                                          std::move(future2));
}
}  // namespace xyco::runtime

#endif  // XYCO_RUNTIME_UTILS_SELECT_H


text:
/* randist/dirichlet.c
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>

/* The Dirichlet probability distribution of order K-1 is

     p(\theta_1,...,\theta_K) d\theta_1 ... d\theta_K =
        (1/Z) \prod_i=1,K \theta_i^{alpha_i - 1} \delta(1 -\sum_i=1,K \theta_i)

   The normalization factor Z can be expressed in terms of gamma functions:

      Z = {\prod_i=1,K \Gamma(\alpha_i)} / {\Gamma( \sum_i=1,K \alpha_i)}

   The K constants, \alpha_1,...,\alpha_K, must be positive. The K parameters,
   \theta_1,...,\theta_K are nonnegative and sum to 1.

   The random variates are generated by sampling K values from gamma
   distributions with parameters a=\alpha_i, b=1, and renormalizing.
   See A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991).

   Gavin E. Crooks <gec@compbio.berkeley.edu> (2002)
*/

void
gsl_ran_dirichlet (const gsl_rng * r, const size_t K,
                   const double alpha[], double theta[])
{
  size_t i;
  double norm = 0.0;

  for (i = 0; i < K; i++)
    {
      theta[i] = gsl_ran_gamma (r, alpha[i], 1.0);
    }

  for (i = 0; i < K; i++)
    {
      norm += theta[i];
    }

  for (i = 0; i < K; i++)
    {
      theta[i] /= norm;
    }
}

double
gsl_ran_dirichlet_pdf (const size_t K,
                       const double alpha[], const double theta[])
{
  return exp (gsl_ran_dirichlet_lnpdf (K, alpha, theta));
}

double
gsl_ran_dirichlet_lnpdf (const size_t K,
                         const double alpha[], const double theta[])
{
  /*We calculate the log of the pdf to minimize the possibility of overflow */
  size_t i;
  double log_p = 0.0;
  double sum_alpha = 0.0;

  for (i = 0; i < K; i++)
    {
      log_p += (alpha[i] - 1.0) * log (theta[i]);
    }

  for (i = 0; i < K; i++)
    {
      sum_alpha += alpha[i];
    }

  log_p += gsl_sf_lngamma (sum_alpha);

  for (i = 0; i < K; i++)
    {
      log_p -= gsl_sf_lngamma (alpha[i]);
    }

  return log_p;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef _UI_HELPERS_H_
#define _UI_HELPERS_H_

#define OEMRESOURCE

#include <algorithm>
#include <functional>
#include <memory>
#include <optional>
#include <span>
#include <string>
#include <type_traits>
#include <unordered_map>
#include <vector>

#include <ppl.h>

#include <gsl/gsl>

// Included before windows.h, because pfc.h includes winsock2.h
#include "../pfc/pfc.h"

#include <windows.h>
#include <WindowsX.h>
#include <SHLWAPI.H>
#include <vssym32.h>
#include <uxtheme.h>
#include <shldisp.h>
#include <ShlObj.h>
#include <gdiplus.h>
#include <Usp10.h>
#include <CommonControls.h>

#include <wil/resource.h>

// Windows SDK headers define min and max macros, which are used by gdiplus.h.
//
// This undefines them temporarily to avoid conflicts with the C++ standard
// library. They are redefined at the end of this header.
#ifdef max
#define _UI_HELPERS_OLD_MAX_
#undef max
#endif

#ifdef min
#define _UI_HELPERS_OLD_MIN_
#undef min
#endif

#include "../mmh/stdafx.h"

#ifndef RECT_CX
#define RECT_CX(rc) (rc.right - rc.left)
#endif

#ifndef RECT_CY
#define RECT_CY(rc) (rc.bottom - rc.top)
#endif

#include "literals.h"

#include "handle.h"
#include "win32_helpers.h"
#include "ole.h"

#include "dialog.h"
#include "dpi.h"
#include "container_window.h"
#include "message_hook.h"
#include "trackbar.h"
#include "solid_fill.h"
#include "theming.h"

#include "gdi.h"
#include "text_drawing.h"
#include "list_view/list_view.h"
#include "drag_image.h"
#include "info_box.h"
#include "menu.h"

#include "ole/data_object.h"
#include "ole/enum_format_etc.h"

#ifdef _UI_HELPERS_OLD_MAX_
#define max(a, b) (((a) > (b)) ? (a) : (b))
#undef _UI_HELPERS_OLD_MAX_
#endif

#ifdef _UI_HELPERS_OLD_MIN_
#define min(a, b) (((a) < (b)) ? (a) : (b))
#undef _UI_HELPERS_OLD_MIN_
#endif

#endif


text:
#ifndef __UTILS_H_INCLUDED__
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <set>
#include <tuple>
#include <math.h>
#include <Eigen/Dense>
#include <gsl/gsl>
#define __UTILS_H_INCLUDED__
#include <cosan/base/CosanBO.h>
namespace Cosan{
    /**
     * @brief General string to number conversion function
     * @details As we allow for user-determined numeric type `NumericType`, the detailed implementation or functions needed maybe slightly
     * different among each other. For instance, when we are trying to read data from `csv` file. Different data type requires different
     * string-to-numeric function. For `double`, one is required `std::stod` while 'std::stof' is the candidate function if
     * `float` is chosen. To take care of this variant before running time, we conside the static-i. The feature
     * allows us to discard branches of an if statement at compile-time based on a constant expression condition. In the following
     * code as an example, we define a template function with input requiring `NumericType` as numeric and then the implementation
     * is decided via `if constexpr` statement.
     **/
    template<Numeric NumericType>
    NumericType StringToNum(const std::string& arg, std::size_t* pos = 0) {
        static_assert(std::is_arithmetic<NumericType>::value, "NumericType must be numeric");
        if constexpr (std::is_same_v<NumericType, unsigned long>) {
            return std::stoul(arg,pos);
        }
        else if constexpr (std::is_same_v<NumericType, unsigned long long>){
            return std::stoull(arg,pos);
        }
        else if constexpr (std::is_same_v<NumericType, int>){
            return std::stoi(arg,pos);
        }
        else if constexpr (std::is_same_v<NumericType, long>){
            return std::stol(arg,pos);
        }
        else if constexpr (std::is_same_v<NumericType, long long>){
            return std::stoll(arg,pos);
        }
        else if constexpr (std::is_same_v<NumericType, float>){
            return std::stof(arg,pos);
        }
        else if constexpr (std::is_same_v<NumericType, double>){
            return std::stod(arg,pos);
        }
        else{
            return std::stold(arg,pos);
        }
    }
}

template<typename Matrix>
Matrix load_csv1 (const std::string & path) {
    std::ifstream indata;

    indata.open(path);
    std::string line;
    std::vector<double> values;
    gsl::index rows = 0;
    while (std::getline(indata, line)) {
        std::stringstream lineStream(line);
        std::string cell;
        while (getline(lineStream, cell, ',')) {
            values.push_back(stod(cell));
        }
        ++rows;
    }

    return Eigen::Map<const Eigen::Matrix<typename Matrix::Scalar, Matrix::RowsAtCompileTime, Matrix::ColsAtCompileTime, Eigen::RowMajor> >(values.data(), rows, values.size()/rows);
}

template<typename Matrix>
void save_csv(const std::string & path, const Matrix & matrix){
//    if (path.substr(path.size()-4)!=".csv"){
//        path.append(".csv");
//    }
    std::ofstream file(path,std::ios::out);
    if (file.is_open()){
        Eigen::IOFormat csvFmt(Eigen::FullPrecision,0,",");
        file<<matrix.format(csvFmt);
        file.close();
    }
}

#endif

max_stars_repo_name: zhxinyu/cosan

max_issues_repo_name: zhxinyu/cosan

max_forks_repo_name: zhxinyu/cosan


text:
#include <complex.h>
#include <math.h>
#include <petsc.h>
#include <stdint.h>
#define LOOPY_CALL_WITH_INTEGER_TYPES(MACRO_NAME) \
    MACRO_NAME(int8, char) \
    MACRO_NAME(int16, short) \
    MACRO_NAME(int32, int) \
    MACRO_NAME(int64, long)
#define LOOPY_DEFINE_FLOOR_DIV_POS_B(SUFFIX, TYPE) \
    static inline TYPE loopy_floor_div_pos_b_##SUFFIX(TYPE a, TYPE b) \
    { \
        if (a<0) \
            a = a - (b-1); \
        return a/b; \
    }
LOOPY_CALL_WITH_INTEGER_TYPES(LOOPY_DEFINE_FLOOR_DIV_POS_B)
#undef LOOPY_DEFINE_FLOOR_DIV_POS_B
#undef LOOPY_CALL_WITH_INTEGER_TYPES
#define BYTES8 (8*8)
typedef double double8 __attribute__ ((vector_size (BYTES8)));
#define BYTES4 (8*4)
typedef int int8 __attribute__ ((vector_size (BYTES4)));

void wrap_form0_cell_integral_otherwise(int const start, int const end, double *__restrict__ dat2, double const *__restrict__ dat1, double const *__restrict__ glob0, double const *__restrict__ dat0, int const *__restrict__ map0, int const *__restrict__ map1)
{
  double8 form_form0_cell_integral_otherwise __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_0 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_1 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_10 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_11 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_12 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_13 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_14 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_15 __attribute__ ((aligned (64)));
  double const form_form0_cell_integral_otherwise_16[4] __attribute__ ((aligned (64))) = { 0.041666666666666664, 0.041666666666666664, 0.041666666666666664, 0.041666666666666664 };
  double const form_form0_cell_integral_otherwise_17[4 * 4] __attribute__ ((aligned (64))) = { 0.13819660112500914, 0.585410196624969, 0.138196601125011, 0.13819660112501098, 0.13819660112500912, 0.13819660112501092, 0.585410196624969, 0.13819660112501098, 0.13819660112500914, 0.13819660112501095, 0.138196601125011, 0.585410196624969, 0.5854101966249672, 0.13819660112501092, 0.138196601125011, 0.13819660112501098 };
  double8 form_form0_cell_integral_otherwise_18 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_19 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_2 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_20 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_21 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_22[4] __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_23 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_24 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_25 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_26 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_27 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_28 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_29 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_3 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_30 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_31 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_32 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_33 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_34 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_35 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_36 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_37 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_38 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_39 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_4 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_40 __attribute__ ((aligned (64)));
  double const form_form0_cell_integral_otherwise_41[4] __attribute__ ((aligned (64))) = { -1.0, 0.0, 0.0, 1.0 };
  double8 form_form0_cell_integral_otherwise_42 __attribute__ ((aligned (64)));
  double const form_form0_cell_integral_otherwise_43[4] __attribute__ ((aligned (64))) = { -1.0, 0.0, 1.0, 0.0 };
  double8 form_form0_cell_integral_otherwise_44 __attribute__ ((aligned (64)));
  double const form_form0_cell_integral_otherwise_45[4] __attribute__ ((aligned (64))) = { -1.0, 1.0, 0.0, 0.0 };
  double8 form_form0_cell_integral_otherwise_5 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_6 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_7 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_8 __attribute__ ((aligned (64)));
  double8 form_form0_cell_integral_otherwise_9 __attribute__ ((aligned (64)));
  double8 t0[4] __attribute__ ((aligned (64)));
  double8 t1[4 * 3] __attribute__ ((aligned (64)));
  double8 t2[4] __attribute__ ((aligned (64)));
  double8 const zero_vec_float64 = { 0.0 };

  {
    /* initial slab for 'n_outer' */
    /*   */
  }
  {
    int const n_outer = start + -1 * ((7 + 7 * start) / 8);

    if (-1 + end + -1 * start >= 0)
    {
      for (int i2 = 0; i2 <= 3; ++i2)
        t2[i2] = zero_vec_float64;
      for (int i0 = 0; i0 <= 3; ++i0)
      {
        for (int i1 = 0; i1 <= 2; ++i1)
          for (int n_batch = start + -8 * n_outer; n_batch <= (-8 + end + -8 * n_outer >= 0 ? 7 : -1 + end + -8 * n_outer); ++n_batch)
            t1[(3 * i0 + i1)][n_batch] = dat1[3 * map1[32 * n_outer + 4 * n_batch + i0] + i1];
        for (int n_batch = start + -8 * n_outer; n_batch <= (-8 + end + -8 * n_outer >= 0 ? 7 : -1 + end + -8 * n_outer); ++n_batch)
          t0[i0][n_batch] = dat0[map0[32 * n_outer + 4 * n_batch + i0]];
      }
      for (int n_batch = start + -8 * n_outer; n_batch <= (-8 + end + -8 * n_outer >= 0 ? 7 : -1 + end + -8 * n_outer); ++n_batch)
      {
        /* no-op (insn=form__start) */
        /*   */
      }
      form_form0_cell_integral_otherwise_20 = zero_vec_float64;
      form_form0_cell_integral_otherwise = -1.0 * t1[0];
      form_form0_cell_integral_otherwise_0 = form_form0_cell_integral_otherwise + t1[3];
      form_form0_cell_integral_otherwise_1 = -1.0 * t1[1];
      form_form0_cell_integral_otherwise_2 = form_form0_cell_integral_otherwise_1 + t1[7];
      form_form0_cell_integral_otherwise_3 = -1.0 * t1[2];
      form_form0_cell_integral_otherwise_4 = form_form0_cell_integral_otherwise_3 + t1[11];
      form_form0_cell_integral_otherwise_5 = form_form0_cell_integral_otherwise_1 + t1[10];
      form_form0_cell_integral_otherwise_6 = form_form0_cell_integral_otherwise_3 + t1[8];
      form_form0_cell_integral_otherwise_7 = form_form0_cell_integral_otherwise_2 * form_form0_cell_integral_otherwise_4 + -1.0 * form_form0_cell_integral_otherwise_5 * form_form0_cell_integral_otherwise_6;
      form_form0_cell_integral_otherwise_8 = form_form0_cell_integral_otherwise + t1[6];
      form_form0_cell_integral_otherwise_9 = form_form0_cell_integral_otherwise_3 + t1[5];
      form_form0_cell_integral_otherwise_10 = form_form0_cell_integral_otherwise_5 * form_form0_cell_integral_otherwise_9;
      form_form0_cell_integral_otherwise_11 = form_form0_cell_integral_otherwise_1 + t1[4];
      form_form0_cell_integral_otherwise_12 = form_form0_cell_integral_otherwise + t1[9];
      form_form0_cell_integral_otherwise_13 = form_form0_cell_integral_otherwise_11 * form_form0_cell_integral_otherwise_6 + -1.0 * form_form0_cell_integral_otherwise_2 * form_form0_cell_integral_otherwise_9;
      form_form0_cell_integral_otherwise_14 = form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_13 + form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_7 + form_form0_cell_integral_otherwise_8 * (form_form0_cell_integral_otherwise_10 + -1.0 * form_form0_cell_integral_otherwise_11 * form_form0_cell_integral_otherwise_4);
      #pragma omp simd
      for (int n_batch_simd = start + -8 * n_outer; n_batch_simd <= (-8 + end + -8 * n_outer >= 0 ? 7 : -1 + end + -8 * n_outer); ++n_batch_simd)
        form_form0_cell_integral_otherwise_15[n_batch_simd] = fabs(form_form0_cell_integral_otherwise_14[n_batch_simd]);
      for (int form_j = 0; form_j <= 3; ++form_j)
        form_form0_cell_integral_otherwise_22[form_j] = zero_vec_float64;
      for (int form_ip = 0; form_ip <= 3; ++form_ip)
      {
        form_form0_cell_integral_otherwise_18 = zero_vec_float64;
        for (int form_i = 0; form_i <= 3; ++form_i)
          form_form0_cell_integral_otherwise_18 = form_form0_cell_integral_otherwise_18 + form_form0_cell_integral_otherwise_17[4 * form_ip + form_i] * t0[form_i];
        form_form0_cell_integral_otherwise_19 = form_form0_cell_integral_otherwise_16[form_ip] * form_form0_cell_integral_otherwise_15;
        form_form0_cell_integral_otherwise_20 = form_form0_cell_integral_otherwise_20 + form_form0_cell_integral_otherwise_19;
        form_form0_cell_integral_otherwise_21 = form_form0_cell_integral_otherwise_19 * glob0[0] * form_form0_cell_integral_otherwise_18;
        for (int form_j_0 = 0; form_j_0 <= 3; ++form_j_0)
          form_form0_cell_integral_otherwise_22[form_j_0] = form_form0_cell_integral_otherwise_22[form_j_0] + form_form0_cell_integral_otherwise_17[4 * form_ip + form_j_0] * form_form0_cell_integral_otherwise_21;
      }
      form_form0_cell_integral_otherwise_23 = 1.0 / form_form0_cell_integral_otherwise_14;
      form_form0_cell_integral_otherwise_24 = form_form0_cell_integral_otherwise_7 * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_25 = -1.0 * t0[0];
      form_form0_cell_integral_otherwise_26 = form_form0_cell_integral_otherwise_25 + t0[1];
      form_form0_cell_integral_otherwise_27 = (form_form0_cell_integral_otherwise_10 + form_form0_cell_integral_otherwise_11 * -1.0 * form_form0_cell_integral_otherwise_4) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_28 = form_form0_cell_integral_otherwise_25 + t0[2];
      form_form0_cell_integral_otherwise_29 = form_form0_cell_integral_otherwise_13 * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_30 = form_form0_cell_integral_otherwise_25 + t0[3];
      form_form0_cell_integral_otherwise_31 = form_form0_cell_integral_otherwise_29 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_24 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_27 * form_form0_cell_integral_otherwise_28;
      form_form0_cell_integral_otherwise_32 = (form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_6 + form_form0_cell_integral_otherwise_4 * -1.0 * form_form0_cell_integral_otherwise_8) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_33 = (form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_4 + form_form0_cell_integral_otherwise_9 * -1.0 * form_form0_cell_integral_otherwise_12) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_34 = (form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_9 + -1.0 * form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_6) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_35 = form_form0_cell_integral_otherwise_34 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_32 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_33 * form_form0_cell_integral_otherwise_28;
      form_form0_cell_integral_otherwise_36 = (form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_5 + -1.0 * form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_2) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_37 = (form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_11 + -1.0 * form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_5) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_38 = (form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_2 + -1.0 * form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_11) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_39 = form_form0_cell_integral_otherwise_38 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_36 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_37 * form_form0_cell_integral_otherwise_28;
      form_form0_cell_integral_otherwise_40 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_38 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_29 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_34) * form_form0_cell_integral_otherwise_20;
      form_form0_cell_integral_otherwise_42 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_37 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_27 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_33) * form_form0_cell_integral_otherwise_20;
      form_form0_cell_integral_otherwise_44 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_36 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_24 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_32) * form_form0_cell_integral_otherwise_20;
      for (int form_j_1 = 0; form_j_1 <= 3; ++form_j_1)
        t2[form_j_1] = t2[form_j_1] + form_form0_cell_integral_otherwise_41[form_j_1] * form_form0_cell_integral_otherwise_40 + form_form0_cell_integral_otherwise_43[form_j_1] * form_form0_cell_integral_otherwise_42 + form_form0_cell_integral_otherwise_22[form_j_1] + form_form0_cell_integral_otherwise_45[form_j_1] * form_form0_cell_integral_otherwise_44;
      for (int n_batch = start + -8 * n_outer; n_batch <= (-8 + end + -8 * n_outer >= 0 ? 7 : -1 + end + -8 * n_outer); ++n_batch)
      {
        /* no-op (insn=statement4) */
        /*   */
      }
      for (int i3 = 0; i3 <= 3; ++i3)
        for (int n_batch = start + -8 * n_outer; n_batch <= (-8 + end + -8 * n_outer >= 0 ? 7 : -1 + end + -8 * n_outer); ++n_batch)
          dat2[map0[32 * n_outer + 4 * n_batch + i3]] = dat2[map0[32 * n_outer + 4 * n_batch + i3]] + t2[i3][n_batch];
    }
  }
  {
    /* bulk slab for 'n_outer' */
    /*   */
  }
  for (int n_outer = 1 + start + -1 * ((7 + 7 * start) / 8); n_outer <= -2 + (7 + end) / 8; ++n_outer)
  {
    for (int i2 = 0; i2 <= 3; ++i2)
      t2[i2] = zero_vec_float64;
    for (int i0 = 0; i0 <= 3; ++i0)
    {
      for (int i1 = 0; i1 <= 2; ++i1)
        for (int n_batch = 0; n_batch <= 7; ++n_batch)
          t1[(3 * i0 + i1)][n_batch] = dat1[3 * map1[32 * n_outer + 4 * n_batch + i0] + i1];
      for (int n_batch = 0; n_batch <= 7; ++n_batch)
        t0[i0][n_batch] = dat0[map0[32 * n_outer + 4 * n_batch + i0]];
    }
    for (int n_batch = 0; n_batch <= 7; ++n_batch)
    {
      /* no-op (insn=form__start) */
      /*   */
    }
    form_form0_cell_integral_otherwise_20 = zero_vec_float64;
    form_form0_cell_integral_otherwise = -1.0 * t1[0];
    form_form0_cell_integral_otherwise_0 = form_form0_cell_integral_otherwise + t1[3];
    form_form0_cell_integral_otherwise_1 = -1.0 * t1[1];
    form_form0_cell_integral_otherwise_2 = form_form0_cell_integral_otherwise_1 + t1[7];
    form_form0_cell_integral_otherwise_3 = -1.0 * t1[2];
    form_form0_cell_integral_otherwise_4 = form_form0_cell_integral_otherwise_3 + t1[11];
    form_form0_cell_integral_otherwise_5 = form_form0_cell_integral_otherwise_1 + t1[10];
    form_form0_cell_integral_otherwise_6 = form_form0_cell_integral_otherwise_3 + t1[8];
    form_form0_cell_integral_otherwise_7 = form_form0_cell_integral_otherwise_2 * form_form0_cell_integral_otherwise_4 + -1.0 * form_form0_cell_integral_otherwise_5 * form_form0_cell_integral_otherwise_6;
    form_form0_cell_integral_otherwise_8 = form_form0_cell_integral_otherwise + t1[6];
    form_form0_cell_integral_otherwise_9 = form_form0_cell_integral_otherwise_3 + t1[5];
    form_form0_cell_integral_otherwise_10 = form_form0_cell_integral_otherwise_5 * form_form0_cell_integral_otherwise_9;
    form_form0_cell_integral_otherwise_11 = form_form0_cell_integral_otherwise_1 + t1[4];
    form_form0_cell_integral_otherwise_12 = form_form0_cell_integral_otherwise + t1[9];
    form_form0_cell_integral_otherwise_13 = form_form0_cell_integral_otherwise_11 * form_form0_cell_integral_otherwise_6 + -1.0 * form_form0_cell_integral_otherwise_2 * form_form0_cell_integral_otherwise_9;
    form_form0_cell_integral_otherwise_14 = form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_13 + form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_7 + form_form0_cell_integral_otherwise_8 * (form_form0_cell_integral_otherwise_10 + -1.0 * form_form0_cell_integral_otherwise_11 * form_form0_cell_integral_otherwise_4);
    #pragma omp simd
    for (int n_batch_simd = 0; n_batch_simd <= 7; ++n_batch_simd)
      form_form0_cell_integral_otherwise_15[n_batch_simd] = fabs(form_form0_cell_integral_otherwise_14[n_batch_simd]);
    for (int form_j = 0; form_j <= 3; ++form_j)
      form_form0_cell_integral_otherwise_22[form_j] = zero_vec_float64;
    for (int form_ip = 0; form_ip <= 3; ++form_ip)
    {
      form_form0_cell_integral_otherwise_18 = zero_vec_float64;
      for (int form_i = 0; form_i <= 3; ++form_i)
        form_form0_cell_integral_otherwise_18 = form_form0_cell_integral_otherwise_18 + form_form0_cell_integral_otherwise_17[4 * form_ip + form_i] * t0[form_i];
      form_form0_cell_integral_otherwise_19 = form_form0_cell_integral_otherwise_16[form_ip] * form_form0_cell_integral_otherwise_15;
      form_form0_cell_integral_otherwise_20 = form_form0_cell_integral_otherwise_20 + form_form0_cell_integral_otherwise_19;
      form_form0_cell_integral_otherwise_21 = form_form0_cell_integral_otherwise_19 * glob0[0] * form_form0_cell_integral_otherwise_18;
      for (int form_j_0 = 0; form_j_0 <= 3; ++form_j_0)
        form_form0_cell_integral_otherwise_22[form_j_0] = form_form0_cell_integral_otherwise_22[form_j_0] + form_form0_cell_integral_otherwise_17[4 * form_ip + form_j_0] * form_form0_cell_integral_otherwise_21;
    }
    form_form0_cell_integral_otherwise_23 = 1.0 / form_form0_cell_integral_otherwise_14;
    form_form0_cell_integral_otherwise_24 = form_form0_cell_integral_otherwise_7 * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_25 = -1.0 * t0[0];
    form_form0_cell_integral_otherwise_26 = form_form0_cell_integral_otherwise_25 + t0[1];
    form_form0_cell_integral_otherwise_27 = (form_form0_cell_integral_otherwise_10 + form_form0_cell_integral_otherwise_11 * -1.0 * form_form0_cell_integral_otherwise_4) * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_28 = form_form0_cell_integral_otherwise_25 + t0[2];
    form_form0_cell_integral_otherwise_29 = form_form0_cell_integral_otherwise_13 * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_30 = form_form0_cell_integral_otherwise_25 + t0[3];
    form_form0_cell_integral_otherwise_31 = form_form0_cell_integral_otherwise_29 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_24 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_27 * form_form0_cell_integral_otherwise_28;
    form_form0_cell_integral_otherwise_32 = (form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_6 + form_form0_cell_integral_otherwise_4 * -1.0 * form_form0_cell_integral_otherwise_8) * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_33 = (form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_4 + form_form0_cell_integral_otherwise_9 * -1.0 * form_form0_cell_integral_otherwise_12) * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_34 = (form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_9 + -1.0 * form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_6) * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_35 = form_form0_cell_integral_otherwise_34 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_32 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_33 * form_form0_cell_integral_otherwise_28;
    form_form0_cell_integral_otherwise_36 = (form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_5 + -1.0 * form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_2) * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_37 = (form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_11 + -1.0 * form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_5) * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_38 = (form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_2 + -1.0 * form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_11) * form_form0_cell_integral_otherwise_23;
    form_form0_cell_integral_otherwise_39 = form_form0_cell_integral_otherwise_38 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_36 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_37 * form_form0_cell_integral_otherwise_28;
    form_form0_cell_integral_otherwise_40 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_38 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_29 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_34) * form_form0_cell_integral_otherwise_20;
    form_form0_cell_integral_otherwise_42 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_37 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_27 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_33) * form_form0_cell_integral_otherwise_20;
    form_form0_cell_integral_otherwise_44 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_36 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_24 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_32) * form_form0_cell_integral_otherwise_20;
    for (int form_j_1 = 0; form_j_1 <= 3; ++form_j_1)
      t2[form_j_1] = t2[form_j_1] + form_form0_cell_integral_otherwise_41[form_j_1] * form_form0_cell_integral_otherwise_40 + form_form0_cell_integral_otherwise_43[form_j_1] * form_form0_cell_integral_otherwise_42 + form_form0_cell_integral_otherwise_22[form_j_1] + form_form0_cell_integral_otherwise_45[form_j_1] * form_form0_cell_integral_otherwise_44;
    for (int n_batch = 0; n_batch <= 7; ++n_batch)
    {
      /* no-op (insn=statement4) */
      /*   */
    }
    for (int i3 = 0; i3 <= 3; ++i3)
      for (int n_batch = 0; n_batch <= 7; ++n_batch)
        dat2[map0[32 * n_outer + 4 * n_batch + i3]] = dat2[map0[32 * n_outer + 4 * n_batch + i3]] + t2[i3][n_batch];
  }
  {
    /* final slab for 'n_outer' */
    /*   */
  }
  {
    int const n_outer = loopy_floor_div_pos_b_int32(-1 + end, 8);

    if (-1 + 8 * n_outer + -1 * start >= 0)
    {
      for (int i2 = 0; i2 <= 3; ++i2)
        t2[i2] = zero_vec_float64;
      for (int i0 = 0; i0 <= 3; ++i0)
      {
        for (int i1 = 0; i1 <= 2; ++i1)
          for (int n_batch = 0; n_batch <= -1 + end + -8 * n_outer; ++n_batch)
            t1[(3 * i0 + i1)][n_batch] = dat1[3 * map1[32 * n_outer + 4 * n_batch + i0] + i1];
        for (int n_batch = 0; n_batch <= -1 + end + -8 * n_outer; ++n_batch)
          t0[i0][n_batch] = dat0[map0[32 * n_outer + 4 * n_batch + i0]];
      }
      for (int n_batch = 0; n_batch <= -1 + end + -8 * n_outer; ++n_batch)
      {
        /* no-op (insn=form__start) */
        /*   */
      }
      form_form0_cell_integral_otherwise_20 = zero_vec_float64;
      form_form0_cell_integral_otherwise = -1.0 * t1[0];
      form_form0_cell_integral_otherwise_0 = form_form0_cell_integral_otherwise + t1[3];
      form_form0_cell_integral_otherwise_1 = -1.0 * t1[1];
      form_form0_cell_integral_otherwise_2 = form_form0_cell_integral_otherwise_1 + t1[7];
      form_form0_cell_integral_otherwise_3 = -1.0 * t1[2];
      form_form0_cell_integral_otherwise_4 = form_form0_cell_integral_otherwise_3 + t1[11];
      form_form0_cell_integral_otherwise_5 = form_form0_cell_integral_otherwise_1 + t1[10];
      form_form0_cell_integral_otherwise_6 = form_form0_cell_integral_otherwise_3 + t1[8];
      form_form0_cell_integral_otherwise_7 = form_form0_cell_integral_otherwise_2 * form_form0_cell_integral_otherwise_4 + -1.0 * form_form0_cell_integral_otherwise_5 * form_form0_cell_integral_otherwise_6;
      form_form0_cell_integral_otherwise_8 = form_form0_cell_integral_otherwise + t1[6];
      form_form0_cell_integral_otherwise_9 = form_form0_cell_integral_otherwise_3 + t1[5];
      form_form0_cell_integral_otherwise_10 = form_form0_cell_integral_otherwise_5 * form_form0_cell_integral_otherwise_9;
      form_form0_cell_integral_otherwise_11 = form_form0_cell_integral_otherwise_1 + t1[4];
      form_form0_cell_integral_otherwise_12 = form_form0_cell_integral_otherwise + t1[9];
      form_form0_cell_integral_otherwise_13 = form_form0_cell_integral_otherwise_11 * form_form0_cell_integral_otherwise_6 + -1.0 * form_form0_cell_integral_otherwise_2 * form_form0_cell_integral_otherwise_9;
      form_form0_cell_integral_otherwise_14 = form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_13 + form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_7 + form_form0_cell_integral_otherwise_8 * (form_form0_cell_integral_otherwise_10 + -1.0 * form_form0_cell_integral_otherwise_11 * form_form0_cell_integral_otherwise_4);
      #pragma omp simd
      for (int n_batch_simd = 0; n_batch_simd <= -1 + end + -8 * n_outer; ++n_batch_simd)
        form_form0_cell_integral_otherwise_15[n_batch_simd] = fabs(form_form0_cell_integral_otherwise_14[n_batch_simd]);
      for (int form_j = 0; form_j <= 3; ++form_j)
        form_form0_cell_integral_otherwise_22[form_j] = zero_vec_float64;
      for (int form_ip = 0; form_ip <= 3; ++form_ip)
      {
        form_form0_cell_integral_otherwise_18 = zero_vec_float64;
        for (int form_i = 0; form_i <= 3; ++form_i)
          form_form0_cell_integral_otherwise_18 = form_form0_cell_integral_otherwise_18 + form_form0_cell_integral_otherwise_17[4 * form_ip + form_i] * t0[form_i];
        form_form0_cell_integral_otherwise_19 = form_form0_cell_integral_otherwise_16[form_ip] * form_form0_cell_integral_otherwise_15;
        form_form0_cell_integral_otherwise_20 = form_form0_cell_integral_otherwise_20 + form_form0_cell_integral_otherwise_19;
        form_form0_cell_integral_otherwise_21 = form_form0_cell_integral_otherwise_19 * glob0[0] * form_form0_cell_integral_otherwise_18;
        for (int form_j_0 = 0; form_j_0 <= 3; ++form_j_0)
          form_form0_cell_integral_otherwise_22[form_j_0] = form_form0_cell_integral_otherwise_22[form_j_0] + form_form0_cell_integral_otherwise_17[4 * form_ip + form_j_0] * form_form0_cell_integral_otherwise_21;
      }
      form_form0_cell_integral_otherwise_23 = 1.0 / form_form0_cell_integral_otherwise_14;
      form_form0_cell_integral_otherwise_24 = form_form0_cell_integral_otherwise_7 * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_25 = -1.0 * t0[0];
      form_form0_cell_integral_otherwise_26 = form_form0_cell_integral_otherwise_25 + t0[1];
      form_form0_cell_integral_otherwise_27 = (form_form0_cell_integral_otherwise_10 + form_form0_cell_integral_otherwise_11 * -1.0 * form_form0_cell_integral_otherwise_4) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_28 = form_form0_cell_integral_otherwise_25 + t0[2];
      form_form0_cell_integral_otherwise_29 = form_form0_cell_integral_otherwise_13 * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_30 = form_form0_cell_integral_otherwise_25 + t0[3];
      form_form0_cell_integral_otherwise_31 = form_form0_cell_integral_otherwise_29 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_24 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_27 * form_form0_cell_integral_otherwise_28;
      form_form0_cell_integral_otherwise_32 = (form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_6 + form_form0_cell_integral_otherwise_4 * -1.0 * form_form0_cell_integral_otherwise_8) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_33 = (form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_4 + form_form0_cell_integral_otherwise_9 * -1.0 * form_form0_cell_integral_otherwise_12) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_34 = (form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_9 + -1.0 * form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_6) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_35 = form_form0_cell_integral_otherwise_34 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_32 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_33 * form_form0_cell_integral_otherwise_28;
      form_form0_cell_integral_otherwise_36 = (form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_5 + -1.0 * form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_2) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_37 = (form_form0_cell_integral_otherwise_12 * form_form0_cell_integral_otherwise_11 + -1.0 * form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_5) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_38 = (form_form0_cell_integral_otherwise_0 * form_form0_cell_integral_otherwise_2 + -1.0 * form_form0_cell_integral_otherwise_8 * form_form0_cell_integral_otherwise_11) * form_form0_cell_integral_otherwise_23;
      form_form0_cell_integral_otherwise_39 = form_form0_cell_integral_otherwise_38 * form_form0_cell_integral_otherwise_30 + form_form0_cell_integral_otherwise_36 * form_form0_cell_integral_otherwise_26 + form_form0_cell_integral_otherwise_37 * form_form0_cell_integral_otherwise_28;
      form_form0_cell_integral_otherwise_40 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_38 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_29 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_34) * form_form0_cell_integral_otherwise_20;
      form_form0_cell_integral_otherwise_42 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_37 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_27 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_33) * form_form0_cell_integral_otherwise_20;
      form_form0_cell_integral_otherwise_44 = (form_form0_cell_integral_otherwise_39 * form_form0_cell_integral_otherwise_36 + form_form0_cell_integral_otherwise_31 * form_form0_cell_integral_otherwise_24 + form_form0_cell_integral_otherwise_35 * form_form0_cell_integral_otherwise_32) * form_form0_cell_integral_otherwise_20;
      for (int form_j_1 = 0; form_j_1 <= 3; ++form_j_1)
        t2[form_j_1] = t2[form_j_1] + form_form0_cell_integral_otherwise_41[form_j_1] * form_form0_cell_integral_otherwise_40 + form_form0_cell_integral_otherwise_43[form_j_1] * form_form0_cell_integral_otherwise_42 + form_form0_cell_integral_otherwise_22[form_j_1] + form_form0_cell_integral_otherwise_45[form_j_1] * form_form0_cell_integral_otherwise_44;
      for (int n_batch = 0; n_batch <= -1 + end + -8 * n_outer; ++n_batch)
      {
        /* no-op (insn=statement4) */
        /*   */
      }
      for (int i3 = 0; i3 <= 3; ++i3)
        for (int n_batch = 0; n_batch <= -1 + end + -8 * n_outer; ++n_batch)
          dat2[map0[32 * n_outer + 4 * n_batch + i3]] = dat2[map0[32 * n_outer + 4 * n_batch + i3]] + t2[i3][n_batch];
    }
  }
}


text:
/* multifit/gsl_multifit.h
 *
 * Copyright (C) 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MULTIFIT_H__
#define __GSL_MULTIFIT_H__

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

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t n; /* number of observations */
  size_t p; /* number of parameters */
  gsl_matrix * A;
  gsl_matrix * Q;
  gsl_matrix * QSI;
  gsl_vector * S;
  gsl_vector * t;
  gsl_vector * xt;
  gsl_vector * D;
}
gsl_multifit_linear_workspace;

gsl_multifit_linear_workspace *
gsl_multifit_linear_alloc (size_t n, size_t p);

void
gsl_multifit_linear_free (gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear (const gsl_matrix * X,
                     const gsl_vector * y,
                     gsl_vector * c,
                     gsl_matrix * cov,
                     double * chisq,
                     gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_svd (const gsl_matrix * X,
                         const gsl_vector * y,
                         double tol,
                         size_t * rank,
                         gsl_vector * c,
                         gsl_matrix * cov,
                         double *chisq,
                         gsl_multifit_linear_workspace * work);

int
gsl_multifit_wlinear (const gsl_matrix * X,
                      const gsl_vector * w,
                      const gsl_vector * y,
                      gsl_vector * c,
                      gsl_matrix * cov,
                      double * chisq,
                      gsl_multifit_linear_workspace * work);

int
gsl_multifit_wlinear_svd (const gsl_matrix * X,
                          const gsl_vector * w,
                          const gsl_vector * y,
                          double tol,
                          size_t * rank,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double *chisq,
                          gsl_multifit_linear_workspace * work);

int
gsl_multifit_linear_est (const gsl_vector * x,
                         const gsl_vector * c,
                         const gsl_matrix * cov, double *y, double *y_err);

__END_DECLS

#endif /* __GSL_MULTIFIT_H__ */


text:
#ifndef _DRAGONCELLO_H
#define _DRAGONCELLO_H

#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>

#include "constants.h"

#include <gsl/gsl_integration.h>
#include <gsl/gsl_matrix.h>

using namespace std;
using namespace DRAGON;

struct diagEntry {
  double upper;
  double central;
  double lower;
};

struct norms {
    double max;
    double mean;
    double norm2;
};

class DRAGONCELLO {

public:

  DRAGONCELLO(int,int);
  ~DRAGONCELLO() { cout << "I am the destructor " << endl; }

  void buildXGrid();
  void buildYGrid();
  void buildZGrid();
  void buildPGrid();
  void initEmptyGrids();
  void initGrids();
  void initSource();
  void initBfield();
  void initAnisoSpatialDiffusion();
  void initDriftLikeVelocities();
  void initEloss();

  void dumpSpectra();
  void dumpSpectra_2(int counter);
  void dumpProfiles();
  void dumpProfiles2D();

  void compute_epsilon(int counter);

  void computeResidualLosses(int);
  void computeResidualDiffAniso(int);

  void computeSourceEnergyIntegral(int, int, int);
  void computeSolutionEnergyIntegral(int, int, int);

  void initCN2Daniso();

  void Run2D();

  void propagateInZaniso(double dt, double dtbar, int counter);
  void propagateInRaniso(double dt, double dtbar, int counter);
  void propagateMixedDerivatives(double dt, double dtbar, int counter);
  void propagateEloss(double dt, double dtbar);

  void printLocalFlux();
  void dumpSolution(int counter);

  inline long int index(int ix, int iy, int iz, int ip) {
    return  (((ix * ny + iy) * nz + iz) * np + ip);
  }

  inline long int indexSpatial(int ix, int iy, int iz) {
    return  (((ix * ny + iy) * nz + iz));
  }

  inline double computeVelocity(double momentum) {
	  return c_light*sqrt( 1. - mass*mass*c_light*c_light/(momentum*momentum)  );
  }

  inline vector<double> getP() {
	  return p_vec;
  }

  void solveTridiagonalSystem(vector<double>& a, vector<double>& b, vector<double>& c, vector<double>& r, vector<double>& u, int n);
  int gsl_linalg_solve_tridiag(const vector<double> & diag, const vector<double> & abovediag, const vector<double> & belowdiag,
                                const vector<double> & rhs, vector<double> & x);
  int solve_tridiag_nonsym(const vector<double> & diag, const vector<double> & abovediag, const vector<double> & belowdiag,
                            const vector<double> & rhs, vector<double> & x, size_t N);

private:

  int Z, A;
  int species_id;

  double deltax, deltay, deltaz, deltaLogp;
  vector<double> x_vec, y_vec, z_vec, p_vec;
  vector<double> deltaxUp, deltayUp, deltazUp;
  vector<double> deltaxDown, deltayDown, deltazDown;
  vector<double> deltaxCentral, deltayCentral, deltazCentral;

  unsigned int ixsun, iysun, izsun, ipLocal;

  vector<double> N;          //(r,z,p) or (x,y,z,p)
  vector<double> N_previous; //(r,z,p) or (x,y,z,p)

  double sourceEnergyIntegral;
  double totalInjectedEnergy;
  double solutionEnergyIntegral;
  double reaccelerationPower;
  double totalReaccelerationEnergy;

  vector<double> analyticalSolution;

  vector<double> sourceTerm;

  vector<double> b_r;  // (r,z) or (r,phi,z); b_r = B_r/|B| magnetic field unit vector projected along r
  vector<double> b_phi;  // (r,z) or (r,phi,z);
  vector<double> b_x;  // b_x = B_x/|B| magnetic field unit vector projected along x
  vector<double> b_y;  // b_y = B_y/|B|
  vector<double> b_z;  // same as before, but along z

  vector<double> diffusionCoefficient;      //(r,z,p) or (x,y,z,p)
  vector<double> diffusionCoefficient_xx;      //(r,z,p) or (x,y,z,p)
  vector<double> diffusionCoefficient_zz;      //(r,z,p) or (x,y,z,p)

  vector<double> Dpara;  //(r,z,p) or (x,y,z,p)
  vector<double> Dperp;  //(r,z,p) or (x,y,z,p)

  vector<double> D_rr;  //(r,z,p) or (r,phi,z,p)
  vector<double> D_rz;  //(r,z,p) or (r,phi,z,p)
  vector<double> D_xx;  //(x,z,p) or (x,y,z,p)
  vector<double> D_xz;  //(x,z,p) or (x,y,z,p)
  vector<double> D_zz;  //(r,z,p) or (x,z,p) or ...
  vector<double> u_r;  //(r,z,p) or (x,y,z,p)
  vector<double> u_x;  //(x,z,p) or (x,y,z,p)
  vector<double> u_z;  //(r,z,p) or (x,z,p) or ...

  vector<double> reaccelerationCoefficient;      //(r,z,p) or (x,y,z,p)
  vector<double> energyLossTerm;      //(r,z,p) or (x,y,z,p)

  vector<double> upperDiagonalZ;
  vector<double> lowerDiagonalZ;
  vector<double> centralDiagonalZ;

  vector<double> upperDiagonalX;
  vector<double> lowerDiagonalX;
  vector<double> centralDiagonalX;

  vector<double> upperDiagonalR;
  vector<double> lowerDiagonalR;
  vector<double> centralDiagonalR;

  vector<double> upperDiagonalDpp;
  vector<double> lowerDiagonalDpp;
  vector<double> centralDiagonalDpp;

  vector<double> upperDiagonalEloss;
  vector<double> lowerDiagonalEloss;
  vector<double> centralDiagonalEloss;
};

#endif

max_stars_repo_name: sscerr/DRAGONCELLO

max_issues_repo_name: sscerr/DRAGONCELLO

max_forks_repo_name: sscerr/DRAGONCELLO


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_movstat.h>

/* compute moving variance by explicitely constructing window and computing variance */
int
slow_movvar(const gsl_movstat_end_t etype, const gsl_vector * x, gsl_vector * y,
            const int H, const int J)
{
  const size_t n = x->size;
  const int K = H + J + 1;
  double *window = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
      double variance = (wsize > 1) ? gsl_stats_variance(window, 1, wsize) : 0.0;

      gsl_vector_set(y, i, variance);
    }

  free(window);

  return GSL_SUCCESS;
}

/* compute moving variance by explicitely constructing window and computing variance */
int
slow_movsd(const gsl_movstat_end_t etype, const gsl_vector * x, gsl_vector * y,
           const int H, const int J)
{
  const size_t n = x->size;
  const int K = H + J + 1;
  double *window = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
      double sd = (wsize > 1) ? gsl_stats_sd(window, 1, wsize) : 0.0;

      gsl_vector_set(y, i, sd);
    }

  free(window);

  return GSL_SUCCESS;
}

static double
func_var(const size_t n, double x[], void * params)
{
  (void) params;
  if (n > 1)
    return gsl_stats_variance(x, 1, n);
  else
    return 0.0;
}

static double
func_sd(const size_t n, double x[], void * params)
{
  (void) params;
  if (n > 1)
    return gsl_stats_sd(x, 1, n);
  else
    return 0.0;
}

static void
test_variance_proc(const double tol, const size_t n, const size_t H, const size_t J,
                   const gsl_movstat_end_t etype, gsl_rng * rng_p)
{
  gsl_movstat_workspace * w = gsl_movstat_alloc2(H, J);
  gsl_vector * x = gsl_vector_alloc(n);
  gsl_vector * y = gsl_vector_alloc(n);
  gsl_vector * z = gsl_vector_alloc(n);
  gsl_movstat_function F1, F2;
  char buf[2048];

  F1.function = func_var;
  F1.params = NULL;

  F2.function = func_sd;
  F2.params = NULL;

  random_vector(x, rng_p);

  /* test variance */

  /* y = variance(x) with slow brute force method */
  slow_movvar(etype, x, y, H, J);

  /* y = variance(x) with fast method */
  gsl_movstat_variance(etype, x, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u variance random", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = variance(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_movstat_variance(etype, z, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u variance random in-place", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = variance(x) with user-defined function */
  gsl_movstat_apply(etype, &F1, x, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u variance user", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* test standard deviation */

  /* y = stddev(x) with slow brute force method */
  slow_movsd(etype, x, y, H, J);

  /* y = stddev(x) with fast method */
  gsl_movstat_sd(etype, x, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u stddev random", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = stddev(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_movstat_sd(etype, z, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u stddev random in-place", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = stddev(x) with user-defined function */
  gsl_movstat_apply(etype, &F2, x, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u stddev user", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  gsl_movstat_free(w);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
}

static void
test_variance(gsl_rng * rng_p)
{
  const double eps = 1.0e-10;

  test_variance_proc(eps, 1000, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 1000, 3, 3, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 1000, 0, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 1000, 5, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 2000, 10, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 2000, 5, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_variance_proc(eps, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_variance_proc(eps, 1000, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 1000, 3, 3, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 1000, 1, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 1000, 5, 1, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 2000, 10, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 2000, 5, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_variance_proc(eps, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_variance_proc(eps, 1000, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 1000, 3, 3, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 1000, 0, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 1000, 5, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 2000, 10, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 2000, 5, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_variance_proc(eps, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "cosmocalc.h"
#include "haloprofs.h"
#include "weaklens.h"

double eps_fun(double x, void *p)
{
  double *v = (double*)p;
  return NFWprof_menc(x,v[0],v[1],v[2]) - v[3];
}

double solve_scale_nfw(gsl_root_fsolver *s, gsl_function F, double *p)
{
  int status;
  int iter,max_iter = 100;
  double x_lo = 0.0,x_hi = 5.0;
  double r;

  iter = 0;
  gsl_root_fsolver_set(s,&F,x_lo,x_hi);
  do
    {
      iter++;
      status = gsl_root_fsolver_iterate(s);
      r = gsl_root_fsolver_root(s);
      x_lo = gsl_root_fsolver_x_lower(s);
      x_hi = gsl_root_fsolver_x_upper(s);
      status = gsl_root_test_interval(x_lo, x_hi,0,0.001);
    }
  while (status == GSL_CONTINUE && iter < max_iter);

  return r;
}

int main(int argc, char **argv)
{
  //init
  cosmoData.cosmoNum = 1;
  cosmoData.OmegaM = 0.27;
  cosmoData.OmegaL = 0.73;
  cosmoData.OmegaB = 0.045;
  cosmoData.OmegaNu = 0.0;
  cosmoData.OmegaK = 0.0;
  cosmoData.h = 0.7;
  cosmoData.w0 = -1.0;
  cosmoData.wa = 0.0;
  cosmoData.SpectralIndex = 0.96;
  cosmoData.Sigma8 = 0.8;
  cosmoData.useSmoothTransFunc = 0;

  //prints mass function in bins to stdout
  double a = atof(argv[1]);
  double mp = atof(argv[2]);
  double np;
  double mmin = 1e1;
  double mmax = 1e18;
  long Nm = 2000;
  double m,dlnm = log(mmax/mmin)/Nm;
  long i;
  double rvir,rs,c;

  cosmoData.delta = 200.0*RHO_CRIT*hubble_noscale(a)*hubble_noscale(a)/(cosmoData.OmegaM*RHO_CRIT/a/a/a);
  //fprintf(stderr,"delta = %f\n",cosmoData.delta);

  gsl_root_fsolver *s;
  gsl_function F;
  double p[4],r1,rN,rN2;

  F.function = &eps_fun;
  F.params = p;
  s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);

  fprintf(stdout,"# m rs rvir r1 rN rN2\n");
  for(i=0;i<Nm;++i)
    {
      m = exp(dlnm*i)*mmin;
      rvir = pow(m/(cosmoData.delta*RHO_CRIT*cosmoData.OmegaM*4.0/3.0*M_PI),1.0/3.0);
      c = duffy2008_concNFW(m,a);
      //if(c < 4.0)
      //c = 4.0;
      rs = rvir/c;

      np = m/mp;

      if(np >= 1.0)
	{
	  p[0] = m;
	  p[1] = rvir;
	  p[2] = c;
	  p[3] = mp;
	  r1 = solve_scale_nfw(s,F,p);
	}
      else
	r1 = -1.0;

      if(np > 100)
	{
	  p[0] = m;
	  p[1] = rvir;
	  p[2] = c;
	  p[3] = mp*100.0;
	  rN = solve_scale_nfw(s,F,p);
	}
      else
	rN = -1.0;

      if(np > 10)
	{
	  p[0] = m;
	  p[1] = rvir;
	  p[2] = c;
	  p[3] = mp*10.0;
	  rN2 = solve_scale_nfw(s,F,p);
	}
      else
	rN2 = -1.0;

      fprintf(stdout,"%e %e %e %e %e %e\n",m,rs,rvir,r1,rN,rN2);
    }

  gsl_root_fsolver_free(s);

  return 0;
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_wavelet.h>

int
main (int argc, char **argv)
{
  (void)(argc); /* avoid unused parameter warning */
  int i, n = 256, nc = 20;
  double *orig_data = malloc (n * sizeof (double));
  double *data = malloc (n * sizeof (double));
  double *abscoeff = malloc (n * sizeof (double));
  size_t *p = malloc (n * sizeof (size_t));

  FILE * f;
  gsl_wavelet *w;
  gsl_wavelet_workspace *work;

  w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
  work = gsl_wavelet_workspace_alloc (n);

  f = fopen (argv[1], "r");
  for (i = 0; i < n; i++)
    {
      fscanf (f, "%lg", &orig_data[i]);
      data[i] = orig_data[i];
    }
  fclose (f);

  gsl_wavelet_transform_forward (w, data, 1, n, work);

  for (i = 0; i < n; i++)
    {
      abscoeff[i] = fabs (data[i]);
    }

  gsl_sort_index (p, abscoeff, 1, n);

  for (i = 0; (i + nc) < n; i++)
    data[p[i]] = 0;

  gsl_wavelet_transform_inverse (w, data, 1, n, work);

  for (i = 0; i < n; i++)
    {
      printf ("%g %g\n", orig_data[i], data[i]);
    }

  gsl_wavelet_free (w);
  gsl_wavelet_workspace_free (work);

  free (data);
  free (orig_data);
  free (abscoeff);
  free (p);
  return 0;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 * \brief Function Definition Library to implement Sensor Fusion Algorithm
 * \author jaspalsingh@cmail.carleton.ca
 * \version 1.0
 * \date 2019-12-01
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include "sensor_fusion.h"

/**
 * \brief Structure pointer for computed support degree matrix
 * \details Structure pointer for calculating support degree matrix
 * \param[in] sensor_readings array of sensor readings
 * \return Success: pointer to the structure support_degree_matrix
 * \return Failure: NULL
 */
struct support_degree_matrix *compute_support_degree_matrix(double *sensor_readings) {
    struct support_degree_matrix *spd;
    spd = (struct support_degree_matrix *)malloc(sizeof(struct support_degree_matrix));
    // count the number of sensor readings
    int length = sizeof(sensor_readings) / sizeof(sensor_readings[0]);
    spd->sensor_count = length;
    double **arrptr = (double **)malloc(length * sizeof(double *));
    int count = 0;
    for (int i = 0; i < length; i++) {
        arrptr[i] = (double *)malloc(length * sizeof(double));
    }
    // allocate memory space for support degree matrix structure
    spd->sd_matrix = (double *)malloc(sizeof(double) * length * length);
    if (arrptr == NULL || spd == NULL || spd->sd_matrix == NULL) {
        printf("ERROR: Failed to allocate memory at %s\n", __func__);
        return NULL;
    }
    // compute support degree matrix
    for (int i = 0; i < length; i++) {
        for (int j = 0; j < length; j++) {
            arrptr[i][j] = exp(-1 * fabs(sensor_readings[i] - sensor_readings[j]));
            spd->sd_matrix[count] = arrptr[i][j];
            count++;
        }
    }
    printf("INFO : Computed support degree matrix\n");
    for (int i = 0; i < length * length; i++) {
        printf("DEBUG: Value => %f\n", spd->sd_matrix[i]);
    }
    free(arrptr);
    return spd;
}

/**
 * \brief Structure pointer calculate eigen value and eigen vector
 * \details Structure pointer for calculating eigen value and eigen vector
 * \param[in] support_degree_matrix pointer to structure support degree matrix
 * \return Success: pointer to structure eigen_value_vector
 * \return Failure: NULL
 */
struct eigen_value_vector *compute_eigen(struct support_degree_matrix *spd) {
    if ((spd->sd_matrix == NULL) || (spd->sensor_count == 0)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return NULL;
    }
    int length = spd->sensor_count;
    struct eigen_value_vector *eigen;
    // allocation memory space for eigen value vector structure
    eigen = (struct eigen_value_vector *)malloc(sizeof(struct eigen_value_vector));
    eigen->eigen_value = (double *)malloc(sizeof(double) * length);
    eigen->eigen_vector = (double **)malloc(length * sizeof(double *));
    for (int i = 0; i < length; i++) {
        eigen->eigen_vector[i] = (double *)malloc(length * sizeof(double));
    }
    if (eigen == NULL || eigen->eigen_vector == NULL || eigen->eigen_value == NULL) {
        printf("ERROR: Failed to allocate memory at %s\n", __func__);
        return NULL;
    }
    //Refer: https://www.gnu.org/software/gsl/doc/html/eigen.html
    // compute eigen values
    gsl_matrix_view m = gsl_matrix_view_array(spd->sd_matrix, length, length);
    gsl_vector *eval = gsl_vector_alloc(length);
    gsl_matrix *evec = gsl_matrix_alloc(length, length);
    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(length);
    gsl_eigen_symmv(&m.matrix, eval, evec, w);
    gsl_eigen_symmv_free(w);
    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);
    {
        int i;
        for (i = 0; i < length; i++) {
            eigen->eigen_value[i] = gsl_vector_get(eval, i);
            gsl_vector_view evec_i = gsl_matrix_column(evec, i);
            for (int j = 0; j < length; j++) {
                eigen->eigen_vector[i][j] = *(*(&evec_i.vector.data) + j * length);
            }
        }
    }
    for (int i = 0; i < length; i++) {
        printf("DEBUG: Eigen Value => %g\n", eigen->eigen_value[i]);
        printf("DEBUG: Eigen Vector => \n");
        for (int j = 0; j < length; j++) {
            printf("%g\n", eigen->eigen_vector[i][j]);
        }
    }
    gsl_vector_free(eval);
    gsl_matrix_free(evec);
    return eigen;
}

/**
 * \brief Function to calculate the contribution rate of principal component
 * \details Compute the contribution rate of kth and mth principal components
 * \param[in] eigen_value_vector pointer to structure of eigen_value_vector
 * \param[in] sensor_count total count of sensors
 * \return Success: pointer to contribution_rate of type double
 * \return Failure: NULL
 */
double *compute_contribution_rate(struct eigen_value_vector *eigen,
                                  int sensor_count) {
    if ((sensor_count == 0) || (eigen->eigen_value == NULL) || (eigen->eigen_vector == NULL)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return NULL;
    }
    double *contribution_rate = (double *)malloc(sensor_count * sizeof(double));
    if (contribution_rate == NULL) {
        printf("ERROR: Failed to allocate memory at %s\n", __func__);
        return NULL;
    }
    double sum = 0;
    for (int i = 0; i < sensor_count; i++) {
        sum += eigen->eigen_value[i];
    }
    printf("INFO : Compute contribution rate\n");
    for (int j = 0; j < sensor_count; j++) {
        contribution_rate[j] = eigen->eigen_value[j] / sum;
        printf("DEBUG: Rate => %g\n", contribution_rate[j]);
    }
    return contribution_rate;
}

/**
 * \brief Function to compute contribution rate
 * \details Compute the count of contribution rates to be used for sensors
 * \param[in] contribution_rate
 * \param[in] threshold
 * \param[in] sensor_count number of sensors to be used
 * \return Success: contribution rate count,
 * \return Failure: -1
 */
int compute_contrib_rate_count(double *contribution_rate,
                               float threshold,
                               int sensor_count) {
    if ((contribution_rate == NULL) || (threshold == 0) || (sensor_count == 0)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return -1;
    }
    double sum = 0;
    int m;
    for (int k = 0; k < sensor_count; k++) {
        for (int j = 0; j <= k; j++) {
            sum += contribution_rate[j];
        }
        if (sum <= threshold) {
            m = k - 1;
            printf("DEBUG: Contribution rate to use %d\n", m);
            return m;
        }
    }
    printf("DEBUG: Count of Contribution rates => %d\n", sensor_count);
    return sensor_count;
}

/**
 * \brief Function to compute principal component
 * \details Computed principal component for each of eigen vector
 * \param[in] support_degree_matrix pointer to structure support degree matrix
 * \param[in] eigen_vector pointer to pointer of eigen vector of type double
 * \param[in] contrib_rate_count number of contribution rate to be used
 * \return Success: principal of support degree matrix,
 * \return Failure: NULL
 */
double **compute_principal_component(struct support_degree_matrix *spd,
                                     double **eigen_vector,
                                     int contrib_rate_count) {
    if ((spd->sensor_count == 0) || (contrib_rate_count == 0) || (spd->sd_matrix == NULL) || (eigen_vector == NULL)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return NULL;
    }
    int n = spd->sensor_count;
    int m = contrib_rate_count;
    int count = 0;
    double **arrptr = (double **)malloc(n * sizeof(double *));
    for (int i = 0; i < n; i++) {
        arrptr[i] = (double *)malloc(n * sizeof(double));
    }
    if (arrptr == NULL) {
        printf("ERROR: Failed to allocate memory at %s\n", __func__);
        return NULL;
    }
    printf("INFO : Compute support degree matrix\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            arrptr[i][j] = spd->sd_matrix[count];
            printf("%f\n", arrptr[i][j]);
            count++;
        }
    }
    double **principal_components_matrix = (double **)malloc(m * sizeof(double *));
    for (int i = 0; i < n; i++) {
        principal_components_matrix[i] = (double *)malloc(n * sizeof(double));
    }
    if (principal_components_matrix == NULL) {
        printf("ERROR: Failed to allocate memory at %s\n", __func__);
        return NULL;
    }
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            principal_components_matrix[i][j] = 0;
            for (int k = 0; k < n; k++) {
                principal_components_matrix[i][j] += eigen_vector[i][k] * arrptr[k][j];
            }
        }
    }
    free(arrptr);
    return principal_components_matrix;
}

/**
 * \brief Function to compute integrated support degree
 * \details Computed support degree matrix for each of the sensors
 * \param[in] principle_components pointer to pointer of principle component of type double
 * \param[in] contribution_rate pointer to contribution rate of type double
 * \param[in] contrib_rate_count number of contribution rate to be used
 * \param[in] sensor_count number of sensors
 * \return Success: array of support degree matrix,
 * \return Failure: NULL
 */
double *compute_integrated_support_degree(double **principle_components,
                                          double *contribution_rate,
                                          int contrib_rate_count,
                                          int sensor_count) {
    if ((sensor_count == 0) || (contrib_rate_count == 0) || (principle_components == NULL) || (contribution_rate == NULL)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return NULL;
    }
    double *arr = (double *)malloc(sensor_count * sizeof(double));
    if (arr == NULL) {
        printf("ERROR: Failed to allocate memory at %s\n", __func__);
        return NULL;
    }
    // Calculate integrated support degree matrix
    for (int i = 0; i < sensor_count; i++) {
        for (int j = 0; j < contrib_rate_count; j++) {
            printf("INFO: Principal Component => %d-%f, Contribution Rate => %d-%f\n", i, principle_components[i][j], j, contribution_rate[i]);
            arr[i] += principle_components[i][j] * contribution_rate[i];
        }
    }
    return arr;
}

/**
 * \brief Function to eliminate invalid sensor readings
 * \details Discard the faulty sensor readings from the support degree matrix
 * \param[in,out] integrated_support_degree_matrix pointer to integrated support degree matrix
 * \param[in] fault_tolerance_threshold threshold value to cut off invalid readings
 * \param[in] sensor_count total count of sensors
 * \return Success: 0,
 * \return Failure: -1
 */
int eliminate_incorrect_data(double *integrated_support_degree_matrix,
                             double fault_tolerance_threshold,
                             int sensor_count) {
    int i;
    double mean = 0;
    double sum = 0;
    double *arr = integrated_support_degree_matrix;
    if ((sensor_count == 0) || (arr == NULL)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return -1;
    }
    for (i = 0; i < sensor_count; i++) {
        sum += arr[i];
    }
    mean = sum / (i + 1);
    for (i = 0; i < sensor_count; i++) {
        if (fabs(arr[i]) < fabs(fault_tolerance_threshold * mean)) {
            arr[i] = 0;
        }
    }
    return 0;
}

/**
 * \brief Function to compute weight coefficient for each sensor
 * \details Compute weighted coefficients for each sensors from the support degree matrix
 * \param[in] integrated_support_degree_matrix pointer to integrated support degree matrix
 * \param[in] sensor_count total count of sensors
 * \return Success: array of weighted coefficients,
 * \return Failure: NULL
 */
double *compute_weight_coefficient(double *integrated_support_degree_matrix,
                                   int sensor_count) {
    double sum_weight_coefficient = 0;
    if ((sensor_count == 0) || (integrated_support_degree_matrix == NULL)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return NULL;
    }
    double *arr = (double *)malloc(sensor_count * sizeof(double));
    if (arr == NULL) {
        printf("ERROR: Failed to allocate memory at %s\n", __func__);
        return NULL;
    }
    memset(arr, 0, sensor_count * sizeof(double));
    for (int i = 0; i < sensor_count; i++) {
        sum_weight_coefficient += integrated_support_degree_matrix[i];
    }
    for (int i = 0; i < sensor_count; i++) {
        arr[i] = integrated_support_degree_matrix[i] / sum_weight_coefficient;
    }
    return arr;
}

/**
 * \brief Function to compute fused output
 * \details Compute the aggregated value of sensors from the weighted coefficients
 * \param[in] weight_coefficient pointer to weighted coefficient of sensors
 * \param[in] sensor_data pointer to data of sensors
 * \param[in] sensor_count total count of sensors
 * \return Success: fused value of sensor data
 * \return Failure: -1
 */
double compute_fused_output(double *weight_coefficient,
                            double *sensor_data,
                            int sensor_count) {
    double fused_output = 0;
    if ((sensor_count == 0) || (weight_coefficient == NULL) || (sensor_data == NULL)) {
        printf("ERROR: Invalid param passed at %s\n", __func__);
        return -1;
    }
    for (int i = 0; i < sensor_count; i++) {
        fused_output += weight_coefficient[i] * sensor_data[i];
    }
    return fused_output;
}

max_stars_repo_name: jaspal-carleton/SensorFusionAlgorithm

max_issues_repo_name: jaspal-carleton/SensorFusionAlgorithm

max_forks_repo_name: jaspal-carleton/SensorFusionAlgorithm


text:
#ifndef ENTITY_H_
#define ENTITY_H_

#include <vector>
#include <memory>
#include <utility>
#include <gsl.h>
#include "Component.h"

using std::vector;
using std::unique_ptr;
using std::make_unique;
using std::move;
using std::forward;

class Entity
{
public:

    static void test();

    template <typename C, typename... Args>
    void addComponent(Args &&... args);

    template <typename C>
    C * getComponent() const;

private:

    Component * getComponent(Component::Id _id, Component::Id _mask) const;

    vector<unique_ptr<Component>> m_components;
    vector<Component::Id> m_componentIds;      // invariant : m_components.size() == m_componentIds.size()
};

template<typename C, typename ...Args>
inline void Entity::addComponent(Args && ...args)
{
    Expects(IdInfoOf<C>::value.isRegistered() != 0);
    m_components.emplace_back(make_unique<C>(forward<Args>(args)...));
    m_componentIds.emplace_back(IdInfoOf<C>::value.getId());
}

template<typename C>
inline C * Entity::getComponent() const
{
    auto & idInfo = IdInfoOf<C>::value;
    auto id = idInfo.getId();
    auto mask = idInfo.getMaskId();
    return static_cast<C*>(getComponent(id, mask));
}

#endif // ENTITY_H_

max_stars_repo_name: Onduril/GetComponent

max_issues_repo_name: Onduril/GetComponent

max_forks_repo_name: Onduril/GetComponent