text:
//**************************************************************************\
//* This file is property of and copyright by the ALICE Project            *\
//* ALICE Experiment at CERN, All rights reserved.                         *\
//*                                                                        *\
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *\
//*                  for The ALICE HLT Project.                            *\
//*                                                                        *\
//* Permission to use, copy, modify and distribute this software and its   *\
//* documentation strictly for non-commercial purposes is hereby granted   *\
//* without fee, provided that the above copyright notice appears in all   *\
//* copies and that both the copyright notice and this permission notice   *\
//* appear in the supporting documentation. The authors make no claims     *\
//* about the suitability of this software for any purpose. It is          *\
//* provided "as is" without express or implied warranty.                  *\
//**************************************************************************

/// \file GPUQA.h
/// \author David Rohr

#ifndef GPUQA_H
#define GPUQA_H

#include "GPUSettings.h"
struct AliHLTTPCClusterMCWeight;
class TH1F;
class TH2F;
class TCanvas;
class TPad;
class TLegend;
class TPad;
class TH1;
class TFile;
class TH1D;
class TObjArray;
class TColor;
typedef short int Color_t;

#if !defined(GPUCA_BUILD_QA) || defined(GPUCA_GPUCODE)

namespace GPUCA_NAMESPACE
{
namespace gpu
{
class GPUChainTracking;

class GPUQA
{
 public:
  GPUQA(GPUChainTracking* chain) {}
  ~GPUQA() = default;

  int InitQA(int tasks = 0) { return 1; }
  void RunQA(bool matchOnly = false) {}
  int DrawQAHistograms() { return 1; }
  void SetMCTrackRange(int min, int max) {}
  bool SuppressTrack(int iTrack) const { return false; }
  bool SuppressHit(int iHit) const { return false; }
  int HitAttachStatus(int iHit) const { return false; }
  int GetMCTrackLabel(unsigned int trackId) const { return -1; }
  bool clusterRemovable(int attach, bool prot) const { return false; }
  static bool QAAvailable() { return false; }
  static bool IsInitialized() { return false; }
};
} // namespace gpu
} // namespace GPUCA_NAMESPACE

#else

#include "GPUTPCDef.h"
#include <cmath>
#include <vector>
#include <memory>
#ifdef GPUCA_TPC_GEOMETRY_O2
#include <gsl/span>
#endif

namespace o2
{
class MCCompLabel;
namespace tpc
{
class TrackTPC;
struct ClusterNativeAccess;
} // namespace tpc
} // namespace o2

struct AliHLTTPCClusterMCLabel;

namespace GPUCA_NAMESPACE::gpu
{
class GPUChainTracking;
class GPUParam;
struct GPUTPCMCInfo;
struct GPUQAGarbageCollection;

class GPUQA
{
 public:
  GPUQA();
  GPUQA(GPUChainTracking* chain, const GPUSettingsQA* config = nullptr, const GPUParam* param = nullptr);
  ~GPUQA();

  int InitQA(int tasks = -1);
  void RunQA(bool matchOnly = false, const std::vector<o2::tpc::TrackTPC>* tracksExternal = nullptr, const std::vector<o2::MCCompLabel>* tracksExtMC = nullptr, const o2::tpc::ClusterNativeAccess* clNative = nullptr);
  int DrawQAHistograms(TObjArray* qcout = nullptr);
  void DrawQAHistogramsCleanup(); // Needed after call to DrawQAHistograms with qcout != nullptr when GPUSettingsQA.shipToQCAsCanvas = true to clean up the Canvases etc.
  void SetMCTrackRange(int min, int max);
  bool SuppressTrack(int iTrack) const;
  bool SuppressHit(int iHit) const;
  int HitAttachStatus(int iHit) const;
  int GetMCTrackLabel(unsigned int trackId) const;
  bool clusterRemovable(int attach, bool prot) const;
  static bool QAAvailable() { return true; }
  bool IsInitialized() { return mQAInitialized; }

  const std::vector<TH1F>& getHistograms1D() const { return *mHist1D; }
  const std::vector<TH2F>& getHistograms2D() const { return *mHist2D; }
  const std::vector<TH1D>& getHistograms1Dd() const { return *mHist1Dd; }
  void resetHists();
  int loadHistograms(std::vector<TH1F>& i1, std::vector<TH2F>& i2, std::vector<TH1D>& i3, int tasks = -1);

  static constexpr int N_CLS_HIST = 8;
  static constexpr int N_CLS_TYPE = 3;

  static constexpr int MC_LABEL_INVALID = -1e9;

  enum QA_TASKS {
    taskTrackingEff = 1,
    taskTrackingRes = 2,
    taskTrackingResPull = 4,
    taskClusterAttach = 8,
    taskTrackStatistics = 16,
    taskClusterCounts = 32,
    taskDefault = 63,
    taskDefaultPostprocess = 31
  };

 private:
  struct additionalMCParameters {
    float pt, phi, theta, eta, nWeightCls;
  };

  struct additionalClusterParameters {
    int attached, fakeAttached, adjacent, fakeAdjacent;
    float pt;
  };

  int InitQACreateHistograms();
  int DoClusterCounts(unsigned long long int* attachClusterCounts, int mode = 0);
  void PrintClusterCount(int mode, int& num, const char* name, unsigned long long int n, unsigned long long int normalization);

  void SetAxisSize(TH1F* e);
  void SetLegend(TLegend* l);
  double* CreateLogAxis(int nbins, float xmin, float xmax);
  void ChangePadTitleSize(TPad* p, float size);
  void DrawHisto(TH1* histo, char* filename, char* options);
  void doPerfFigure(float x, float y, float size);
  void GetName(char* fname, int k);
  template <class T>
  T* GetHist(T*& ee, std::vector<std::unique_ptr<TFile>>& tin, int k, int nNewInput);

#ifdef GPUCA_TPC_GEOMETRY_O2
  using mcLabels_t = gsl::span<const o2::MCCompLabel>;
  using mcLabel_t = o2::MCCompLabel;
  using mcLabelI_t = mcLabel_t;
  using mcInfo_t = GPUTPCMCInfo;
  mcLabels_t GetMCLabel(unsigned int i);
  mcLabel_t GetMCLabel(unsigned int i, unsigned int j);
#else
  using mcLabels_t = AliHLTTPCClusterMCLabel;
  using mcLabel_t = AliHLTTPCClusterMCWeight;
  struct mcLabelI_t {
    int getTrackID() const { return AbsLabelID(track); }
    int getEventID() const { return 0; }
    long int getTrackEventSourceID() const { return getTrackID(); }
    bool isFake() const { return track < 0; }
    bool isValid() const { return track != MC_LABEL_INVALID; }
    void invalidate() { track = MC_LABEL_INVALID; }
    void setFakeFlag(bool v = true) { track = v ? FakeLabelID(track) : AbsLabelID(track); }
    void setNoise() { track = MC_LABEL_INVALID; }
    bool operator==(const mcLabel_t& l);
    bool operator!=(const mcLabel_t& l) { return !(*this == l); }
    mcLabelI_t() = default;
    mcLabelI_t(const mcLabel_t& l);
    int track = MC_LABEL_INVALID;
  };
  using mcInfo_t = GPUTPCMCInfo;
  const mcLabels_t& GetMCLabel(unsigned int i);
  const mcLabel_t& GetMCLabel(unsigned int i, unsigned int j);
  const mcInfo_t& GetMCTrack(const mcLabelI_t& label);
  static int FakeLabelID(const int id);
  static int AbsLabelID(const int id);
#endif
  template <class T>
  static auto& GetMCTrackObj(T& obj, const mcLabelI_t& l);

  unsigned int GetNMCCollissions();
  unsigned int GetNMCTracks(int iCol);
  unsigned int GetNMCLabels();
  const mcInfo_t& GetMCTrack(unsigned int iTrk, unsigned int iCol);
  const mcInfo_t& GetMCTrack(const mcLabel_t& label);
  int GetMCLabelNID(const mcLabels_t& label);
  int GetMCLabelNID(unsigned int i);
  int GetMCLabelID(unsigned int i, unsigned int j);
  int GetMCLabelCol(unsigned int i, unsigned int j);
  static int GetMCLabelID(const mcLabels_t& label, unsigned int j);
  static int GetMCLabelID(const mcLabel_t& label);
  float GetMCLabelWeight(unsigned int i, unsigned int j);
  float GetMCLabelWeight(const mcLabels_t& label, unsigned int j);
  float GetMCLabelWeight(const mcLabel_t& label);
  const auto& GetClusterLabels();
  bool mcPresent();

  static bool MCComp(const mcLabel_t& a, const mcLabel_t& b);

  GPUChainTracking* mTracking;
  const GPUSettingsQA& mConfig;
  const GPUParam& mParam;

  const char* str_perf_figure_1 = "ALICE Performance 2018/03/20";
  // const char* str_perf_figure_2 = "2015, MC pp, #sqrt{s} = 5.02 TeV";
  const char* str_perf_figure_2 = "2015, MC Pb-Pb, #sqrt{s_{NN}} = 5.02 TeV";
  //-------------------------

  std::vector<mcLabelI_t> mTrackMCLabels;
#ifdef GPUCA_TPC_GEOMETRY_O2
  std::vector<std::vector<int>> mTrackMCLabelsReverse;
  std::vector<std::vector<int>> mRecTracks;
  std::vector<std::vector<int>> mFakeTracks;
  std::vector<std::vector<additionalMCParameters>> mMCParam;
  std::vector<std::vector<mcInfo_t>> mMCInfos;
  std::vector<int> mNColTracks;
#else
  std::vector<int> mTrackMCLabelsReverse[1];
  std::vector<int> mRecTracks[1];
  std::vector<int> mFakeTracks[1];
  std::vector<additionalMCParameters> mMCParam[1];
#endif
  std::vector<additionalClusterParameters> mClusterParam;
  int mNTotalFakes = 0;

  TH1F* mEff[4][2][2][5][2]; // eff,clone,fake,all - findable - secondaries - y,z,phi,eta,pt - work,result
  TCanvas* mCEff[6];
  TPad* mPEff[6][4];
  TLegend* mLEff[6];

  TH1F* mRes[5][5][2]; // y,z,phi,lambda,pt,ptlog res - param - res,mean
  TH2F* mRes2[5][5];
  TCanvas* mCRes[7];
  TPad* mPRes[7][5];
  TLegend* mLRes[6];

  TH1F* mPull[5][5][2]; // y,z,phi,lambda,pt,ptlog res - param - res,mean
  TH2F* mPull2[5][5];
  TCanvas* mCPull[7];
  TPad* mPPull[7][5];
  TLegend* mLPull[6];

  enum CL_types { CL_attached = 0,
                  CL_fake = 1,
                  CL_att_adj = 2,
                  CL_fakeAdj = 3,
                  CL_tracks = 4,
                  CL_physics = 5,
                  CL_prot = 6,
                  CL_all = 7 };
  TH1D* mClusters[N_CLS_TYPE * N_CLS_HIST - 1]; // attached, fakeAttached, attach+adjacent, fakeAdjacent, physics, protected, tracks, all / count, rel, integral
  TCanvas* mCClust[N_CLS_TYPE];
  TPad* mPClust[N_CLS_TYPE];
  TLegend* mLClust[N_CLS_TYPE];

  struct counts_t {
    long long int nRejected = 0, nTube = 0, nTube200 = 0, nLoopers = 0, nLowPt = 0, n200MeV = 0, nPhysics = 0, nProt = 0, nUnattached = 0, nTotal = 0, nHighIncl = 0, nAbove400 = 0, nFakeRemove400 = 0, nFullFakeRemove400 = 0, nBelow40 = 0, nFakeProtect40 = 0, nMergedLooper = 0;
    double nUnaccessible = 0;
  } mClusterCounts;

  TH1F* mTracks;
  TCanvas* mCTracks;
  TPad* mPTracks;
  TLegend* mLTracks;

  TH1F* mNCl;
  TCanvas* mCNCl;
  TPad* mPNCl;
  TLegend* mLNCl;

  std::vector<TH2F*> mHistClusterCount;

  std::vector<TH1F>* mHist1D = nullptr;
  std::vector<TH2F>* mHist2D = nullptr;
  std::vector<TH1D>* mHist1Dd = nullptr;
  bool mHaveExternalHists = false;
  std::vector<TH1F**> mHist1D_pos{};
  std::vector<TH2F**> mHist2D_pos{};
  std::vector<TH1D**> mHist1Dd_pos{};
  template <class T>
  auto getHistArray();
  template <class T, typename... Args>
  void createHist(T*& h, const char* name, Args... args);

  std::unique_ptr<GPUQAGarbageCollection> mGarbageCollector;
  template <class T, typename... Args>
  T* createGarbageCollected(Args... args);
  void clearGarbagageCollector();

  int mNEvents = 0;
  bool mQAInitialized = false;
  int mQATasks = 0;
  std::vector<std::vector<int>> mcEffBuffer;
  std::vector<std::vector<int>> mcLabelBuffer;
  std::vector<std::vector<bool>> mGoodTracks;
  std::vector<std::vector<bool>> mGoodHits;

  static std::vector<TColor*> mColors;
  static int initColors();

  int mMCTrackMin = -1, mMCTrackMax = -1;

  const o2::tpc::ClusterNativeAccess* mClNative = nullptr;
};

inline bool GPUQA::SuppressTrack(int iTrack) const { return (mConfig.matchMCLabels.size() && !mGoodTracks[mNEvents][iTrack]); }
inline bool GPUQA::SuppressHit(int iHit) const { return (mConfig.matchMCLabels.size() && !mGoodHits[mNEvents - 1][iHit]); }
inline int GPUQA::HitAttachStatus(int iHit) const { return (mClusterParam.size() && mClusterParam[iHit].fakeAttached ? (mClusterParam[iHit].attached ? 1 : 2) : 0); }

} // namespace GPUCA_NAMESPACE::gpu

#endif
#endif

ext: h

max_stars_repo_name: iarsene/AliRoot

max_issues_repo_name: iarsene/AliRoot

max_forks_repo_name: iarsene/AliRoot


text:
/* rng/slatec.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/**

* ======================================================================
* NIST Guide to Available Math Software.
* Source for module RAND from package CMLIB.
* Retrieved from TIBER on Fri Oct 11 11:43:42 1996.
* ======================================================================
      FUNCTION RAND(R)
C***BEGIN PROLOGUE  RAND
C***DATE WRITTEN   770401   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  L6A21
C***KEYWORDS  RANDOM NUMBER,SPECIAL FUNCTION,UNIFORM
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Generates a uniformly distributed random number.
C***DESCRIPTION
C
C      This pseudo-random number generator is portable among a wide
C variety of computers.  RAND(R) undoubtedly is not as good as many
C readily available installation dependent versions, and so this
C routine is not recommended for widespread usage.  Its redeeming
C feature is that the exact same random numbers (to within final round-
C off error) can be generated from machine to machine.  Thus, programs
C that make use of random numbers can be easily transported to and
C checked in a new environment.
C      The random numbers are generated by the linear congruential
C method described, e.g., by Knuth in Seminumerical Methods (p.9),
C Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
C sequence, the I+1 -st number is generated from
C             X(I+1) = (A*X(I) + C) MOD M,
C where here M = 2**22 = 4194304, C = 1731 and several suitable values
C of the multiplier A are discussed below.  Both the multiplier A and
C random number X are represented in double precision as two 11-bit
C words.  The constants are chosen so that the period is the maximum
C possible, 4194304.
C      In order that the same numbers be generated from machine to
C machine, it is necessary that 23-bit integers be reducible modulo
C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
C integers be multiplied exactly.  Furthermore, if the restart option
C is used (where R is between 0 and 1), then the product R*2**22 =
C R*4194304 must be correct to the nearest integer.
C      The first four random numbers should be .0004127026,
C .6750836372, .1614754200, and .9086198807.  The tenth random number
C is .5527787209, and the hundredth is .3600893021 .  The thousandth
C number should be .2176990509 .
C      In order to generate several effectively independent sequences
C with the same generator, it is necessary to know the random number
C for several widely spaced calls.  The I-th random number times 2**22,
C where I=K*P/8 and P is the period of the sequence (P = 2**22), is
C still of the form L*P/8.  In particular we find the I-th random
C number multiplied by 2**22 is given by
C I   =  0  1*P/8  2*P/8  3*P/8  4*P/8  5*P/8  6*P/8  7*P/8  8*P/8
C RAND=  0  5*P/8  2*P/8  7*P/8  4*P/8  1*P/8  6*P/8  3*P/8  0
C Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
C      Several multipliers have been subjected to the spectral test
C (see Knuth, p. 82).  Four suitable multipliers roughly in order of
C goodness according to the spectral test are
C    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
C    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
C    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
C    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
C
C      In the table below LOG10(NU(I)) gives roughly the number of
C random decimal digits in the random numbers considered I at a time.
C C is the primary measure of goodness.  In both cases bigger is better.
C
C                   LOG10 NU(I)              C(I)
C       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
C
C    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
C    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
C    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
C    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
C   Best
C    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
C
C             Input Argument --
C R      If R=0., the next random number of the sequence is generated.
C        If R .LT. 0., the last generated number will be returned for
C          possible use in a restart procedure.
C        If R .GT. 0., the sequence of random numbers will start with
C          the seed R mod 1.  This seed is also returned as the value of
C          RAND provided the arithmetic is done exactly.
C
C             Output Value --
C RAND   a pseudo-random number between 0. and 1.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  RAND
      DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
      DATA IC /1731/
      DATA IX1, IX0 /0, 0/
C***FIRST EXECUTABLE STATEMENT  RAND
      IF (R.LT.0.) GO TO 10
      IF (R.GT.0.) GO TO 20
C
C           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
C                   + IA0*IX0) + IA0*IX0
C
      IY0 = IA0*IX0
      IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0
      IY0 = IY0 + IC
      IX0 = MOD (IY0, 2048)
      IY1 = IY1 + (IY0-IX0)/2048
      IX1 = MOD (IY1, 2048)
C
 10   RAND = IX1*2048 + IX0
      RAND = RAND / 4194304.
      RETURN
C
 20   IX1 = AMOD(R,1.)*4194304. + 0.5
      IX0 = MOD (IX1, 2048)
      IX1 = (IX1-IX0)/2048
      GO TO 10
C
      END

  **/

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

static inline unsigned long int slatec_get (void *vstate);
static double slatec_get_double (void *vstate);
static void slatec_set (void *state, unsigned long int s);

typedef struct
  {
    long int x0, x1;
  }
slatec_state_t;

static const long P = 4194304;
static const long a1 = 1536;
static const long a0 = 1029;
static const long a1ma0 = 507;
static const long c = 1731;

static inline unsigned long int
slatec_get (void *vstate)
{
  long y0, y1;
  slatec_state_t *state = (slatec_state_t *) vstate;

  y0 = a0 * state->x0;
  y1 = a1 * state->x1 + a1ma0 * (state->x0 - state->x1) + y0;
  y0 = y0 + c;
  state->x0 = y0 % 2048;
  y1 = y1 + (y0 - state->x0) / 2048;
  state->x1 = y1 % 2048;

  return state->x1 * 2048 + state->x0;
}

static double
slatec_get_double (void *vstate)
{
  return slatec_get (vstate) / 4194304.0 ;
}

static void
slatec_set (void *vstate, unsigned long int s)
{
  slatec_state_t *state = (slatec_state_t *) vstate;

  /* Only eight seeds are permitted.  This is pretty limiting, but
     at least we are guaranteed that the eight sequences are different */

  s = s % 8;
  s *= P / 8;

  state->x0 = s % 2048;
  state->x1 = (s - state->x0) / 2048;
}

static const gsl_rng_type slatec_type =
{"slatec",			/* name */
 4194303,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (slatec_state_t),
 &slatec_set,
 &slatec_get,
 &slatec_get_double};

const gsl_rng_type *gsl_rng_slatec = &slatec_type;

ext: c


text:
/*********************
 * Time intergral KDK scheme.
 * kick and drifts.
 *
 * This code was initially modified by Jun Koda,
 * from the original serial COLA code
 * by Svetlin Tassev.
 *
 * The kick and drift still supports a COLA compat-mode.
 * Most of the nasty factors are for COLA compat-mode
 * (not needed in PM)
 * We also added a 2LPT mode that does just 2LPT.
 *
 *  Yu Feng <rainwoodman@gmail.com>
 *
 */

#include <math.h>
#include <string.h>
#include <assert.h>
#include <alloca.h>
#include <mpi.h>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_errno.h>

#include <fastpm/libfastpm.h>
#include <fastpm/logging.h>

#include "pmpfft.h"
#include "vpm.h"

static double
Sq(double ai, double af, double aRef, double nLPT, FastPMCosmology * c, int USE_NONSTDDA);

static double
Sphi(double ai, double af, double aRef, double nLPT, FastPMCosmology * c, int USE_NONSTDDA);

static inline void
fastpm_drift_lookup(FastPMDriftFactor * drift, double af, double * dyyy, double * da1, double * da2)
{
    double ind;

    if(af == drift->af) {
        *dyyy = drift->dyyy[drift->nsamples - 1];
        *da1  = drift->da1[drift->nsamples - 1];
        *da2  = drift->da2[drift->nsamples - 1];
        return;
    }
    if(af == drift->ai) {
        *dyyy = drift->dyyy[0];
        *da1  = drift->da1[0];
        *da2  = drift->da2[0];
        return;
    }
    {
        ind = (af - drift->ai) / (drift->af - drift->ai) * (drift->nsamples - 1);
        int l = floor(ind);
        double u = l + 1 - ind;
        double v = ind - l;
        if(l + 1 >= drift->nsamples) {
            fastpm_raise(-1, "drift beyond factor's available range. ");
        }
        *dyyy = drift->dyyy[l] * u + drift->dyyy[l + 1] * v;
        *da1  = drift->da1[l] * u + drift->da1[l + 1] * v;
        *da2  = drift->da2[l] * u + drift->da2[l + 1] * v;
    }
}

inline void
fastpm_drift_one(FastPMDriftFactor * drift, FastPMStore * p, ptrdiff_t i, double xo[3], double af)
{

    double dyyy_f, da1_f, da2_f;
    double dyyy_i, da1_i, da2_i;
    double dyyy, da1, da2;

    fastpm_drift_lookup(drift, af, &dyyy_f, &da1_f, &da2_f);
    fastpm_drift_lookup(drift, p->meta.a_x, &dyyy_i, &da1_i, &da2_i);

    dyyy = dyyy_f - dyyy_i;
    da1 = da1_f - da1_i;
    da2 = da2_f - da2_i;

    int d;
    for(d = 0; d < 3; d ++) {
        double v;
        switch(drift->forcemode) {
            case FASTPM_FORCE_2LPT:
                xo[d] = p->x[i][d] + p->dx1[i][d] * da1 + p->dx2[i][d] * da2;
            break;
            case FASTPM_FORCE_ZA:
                xo[d] = p->x[i][d] + p->dx1[i][d] * da1;
            break;
            case FASTPM_FORCE_FASTPM:
            case FASTPM_FORCE_PM:
                xo[d] = p->x[i][d] + p->v[i][d] * dyyy;
            break;
            case FASTPM_FORCE_COLA:
                /* For cola, remove the lpt velocity to find the residual velocity v*/
                v = p->v[i][d] - (p->dx1[i][d]*drift->Dv1 + p->dx2[i][d]*drift->Dv2);
                xo[d] = p->x[i][d] + v * dyyy;
                xo[d] += p->dx1[i][d] * da1 + p->dx2[i][d] * da2;
            break;
        }
        /* if PGDCorrection is enabled, add it */
        if(p->pgdc) {
            /* no drift; to protect the pgdc line */
            if (drift->ai == drift->af) continue;
            xo[d] += 0.5 * p->pgdc[i][d] * dyyy / drift->dyyy[drift->nsamples-1];
        }
    }
}
static inline void
fastpm_kick_lookup(FastPMKickFactor * kick, double af, double * dda, double * Dv1, double * Dv2)
{
    double ind;

    if(af == kick->af) {
        *dda = kick->dda[kick->nsamples - 1];
        *Dv1 = kick->Dv1[kick->nsamples - 1];
        *Dv2 = kick->Dv2[kick->nsamples - 1];
        return;
    }
    if(af == kick->ai) {
        *dda = kick->dda[0];
        *Dv1 = kick->Dv1[0];
        *Dv2 = kick->Dv2[0];
        return;
    }
    {
        ind = (af - kick->ai) / (kick->af - kick->ai) * (kick->nsamples - 1);
        int l = floor(ind);
        double u = l + 1 - ind;
        double v = ind - l;
        if(l + 1 >= kick->nsamples) {
            fastpm_raise(-1, "kick beyond factor's available range. ");
        }
        *dda = kick->dda[l] * u + kick->dda[l + 1] * v;
        *Dv1 = kick->Dv1[l] * u + kick->Dv1[l + 1] * v;
        *Dv2 = kick->Dv2[l] * u + kick->Dv2[l + 1] * v;
    }
}

inline void
fastpm_kick_one(FastPMKickFactor * kick, FastPMStore * p, ptrdiff_t i, float vo[3], double af)
{
    double dda_i, Dv1_i, Dv2_i;
    double dda_f, Dv1_f, Dv2_f;
    double dda, Dv1, Dv2;

    fastpm_kick_lookup(kick, af, &dda_f, &Dv1_f, &Dv2_f);
    fastpm_kick_lookup(kick, p->meta.a_v, &dda_i, &Dv1_i, &Dv2_i);
    dda = dda_f - dda_i;
    Dv1 = Dv1_f - Dv1_i;
    Dv2 = Dv2_f - Dv2_i;

    int d;
    for(d = 0; d < 3; d++) {
        float ax = p->acc[i][d];       // unlike a_x, which means a at which x is calcd
        if(kick->forcemode == FASTPM_FORCE_COLA) {
            ax += (p->dx1[i][d]*kick->q1 + p->dx2[i][d]*kick->q2);
        }
        vo[d] = p->v[i][d] + ax * dda;
        if(kick->forcemode == FASTPM_FORCE_COLA) {
            vo[d] += (p->dx1[i][d] * Dv1 + p->dx2[i][d] * Dv2);
        }
    }
}

// Leap frog time integration

void
fastpm_kick_store(FastPMKickFactor * kick,
    FastPMStore * pi, FastPMStore * po, double af)
{
    int np = pi->np;

    // Kick using acceleration at a= ac
    // Assume forces at a=ac is in particles->force

    int i;
#pragma omp parallel for
    for(i=0; i<np; i++) {
        int d;
        float vo[3];
        fastpm_kick_one(kick, pi, i, vo, af);
        for(d = 0; d < 3; d++) {
            po->v[i][d] = vo[d];
        }
    }

    //velocity is now at a= avel1
    po->meta.a_v = af;
}

static double G_p(FastPMGrowthInfo * growth_info)
{
    /* integral of G_p */
    return growth_info->D1;
}
static double g_p(FastPMGrowthInfo * growth_info)
{
    return DGrowthFactorDa(growth_info);
}

static double G_f(FastPMGrowthInfo * growth_info)
{
    /* integral of g_f */
    double a = growth_info->a;
    return a * a * a * HubbleEa(a, growth_info->c) * g_p(growth_info);
}

static double g_f(FastPMGrowthInfo * growth_info)
{
    double a = growth_info->a;
    FastPMCosmology * c = growth_info->c;

    double E = HubbleEa(a, c);
    double dEda = DHubbleEaDa(a, c);

    double dDda = g_p(growth_info);
    double d2Dda2 = D2GrowthFactorDa2(growth_info);

    double g_f = 3 * a * a * E * dDda
                   + a * a * a * dEda * dDda
                   + a * a * a * E * d2Dda2;
    return g_f;
}

void fastpm_kick_init(FastPMKickFactor * kick, FastPMSolver * fastpm, double ai, double ac, double af)
{
    FastPMCosmology * c = fastpm->cosmology;
    kick->forcemode = fastpm->config->FORCE_TYPE;

    FastPMGrowthInfo gi_i;
    FastPMGrowthInfo gi_c;
    FastPMGrowthInfo gi_e;

    fastpm_growth_info_init(&gi_i, ai, c);
    fastpm_growth_info_init(&gi_c, ac, c);

    double E_i = HubbleEa(ai, c);
    double E_c = HubbleEa(ac, c);

    double D1_i = gi_i.D1;
    double D2_i = gi_i.D2;
    double f1_i = gi_i.f1;
    double f2_i = gi_i.f2;

    double D1_c = gi_c.D1;
    double D2_c = gi_c.D2;

    double Omega_m0 = Omega_source(1, c);
    double Omega_mc = Omega_source(ac, c);

    // kick->q1,2 are used for the COLA force implementation.
    // growth_mode = ODE and LCDM should match for an LCDM background,
    // but neither is guaranteed accurate for a background with radiaiton.
    // We advise using LCDM mode for forcemode = FASTPM_FORCE_COLA, as in the
    // original implementation of FastPM.
    kick->q1 = D1_c;
    switch (c->growth_mode){
        case FASTPM_GROWTH_MODE_LCDM:
            kick->q2 = D1_c*D1_c * (1.0 + 7.0/3.0 * pow(Omega_mc, 1.0/143.0));
        break;
        case FASTPM_GROWTH_MODE_ODE:
            kick->q2 = D1_c*D1_c * (1 - D1_c*D1_c/D2_c);
        break;
        default:
            fastpm_raise(-1, "Please enter a valid growth mode.\n");
    }

    kick->nsamples = 32;
    int i;

    double Dv1i = D1_i * ai * ai * E_i * f1_i;
    double Dv2i = D2_i * ai * ai * E_i * f2_i;
    for(i = 0; i < kick->nsamples; i ++) {
        double ae = ai * (1.0 * (kick->nsamples - 1 - i) / (kick->nsamples - 1))
                  + af * (1.0 * i / (kick->nsamples - 1));

        fastpm_growth_info_init(&gi_e, ae, c);
        double D1_e = gi_e.D1;
        double f1_e = gi_e.f1;
        double D2_e = gi_e.D2;
        double f2_e = gi_e.f2;
        double E_e = HubbleEa(ae, c);

        if(kick->forcemode == FASTPM_FORCE_FASTPM) {
            kick->dda[i] = -1.5 * Omega_mc * ac
               * E_c
               * (G_f(&gi_e) - G_f(&gi_i)) / g_f(&gi_c);
        } else {
            kick->dda[i] = -1.5 * Omega_m0
                * Sphi(ai, ae, ac, fastpm->config->nLPT, c, kick->forcemode == FASTPM_FORCE_COLA);
        }
        kick->Dv1[i] = D1_e * ae * ae * E_e * f1_e - Dv1i;
        kick->Dv2[i] = D2_e * ae * ae * E_e * f2_e - Dv2i;
    }

    kick->ai = ai;
    kick->ac = ac;
    kick->af = af;

    /* Output growth and FastPM factor at af for reference.
    This is a weird place to put this, but it's convenient because G and g are static */
    fastpm_info("Growth/FastPM factors at a = %6.4f: D1=%g, D2=%g, f1=%g, f2=%g, G_p=%g, G_f=%g, g_p=%g, g_f=%g\n",
               ai,
               gi_i.D1,
               gi_i.D2,
               gi_i.f1,
               gi_i.f2,
               G_p(&gi_i),
               G_f(&gi_i),
               g_p(&gi_i),
               g_f(&gi_i));
}

void
fastpm_drift_init(FastPMDriftFactor * drift, FastPMSolver * fastpm,
                double ai, double ac, double af)
{
    FastPMCosmology * c = fastpm->cosmology;
    drift->forcemode = fastpm->config->FORCE_TYPE;

    FastPMGrowthInfo gi_i;
    FastPMGrowthInfo gi_c;
    FastPMGrowthInfo gi_e;

    fastpm_growth_info_init(&gi_i, ai, c);
    fastpm_growth_info_init(&gi_c, ac, c);

    double E_c = HubbleEa(ac, c);

    double D1_i = gi_i.D1;
    double D2_i = gi_i.D2;

    double D1_c = gi_c.D1;
    double D2_c = gi_c.D2;
    double f1_c = gi_c.f1;
    double f2_c = gi_c.f2;

    drift->nsamples = 32;
    int i;

    for(i = 0; i < drift->nsamples; i ++ ) {
        double ae = ai * (1.0 * (drift->nsamples - 1 - i) / (drift->nsamples - 1))
                  + af * (1.0 * i / (drift->nsamples - 1));

        fastpm_growth_info_init(&gi_e, ae, c);     // overwrite each iteration
        double D1_e = gi_e.D1;
        double D2_e = gi_e.D2;

        if (drift->forcemode == FASTPM_FORCE_FASTPM) {
            drift->dyyy[i] = 1 / (ac * ac * ac * E_c)
                        * (G_p(&gi_e) - G_p(&gi_i)) / g_p(&gi_c);
        } else {
            drift->dyyy[i] = Sq(ai, ae, ac, fastpm->config->nLPT, c, drift->forcemode == FASTPM_FORCE_COLA);
        }
        drift->da1[i] = D1_e - D1_i;    // change in D_1lpt
        drift->da2[i] = D2_e - D2_i;  // change in D_2lpt
    }
    drift->af = af;
    drift->ai = ai;
    drift->ac = ac;
    drift->Dv1 = D1_c * ac * ac * E_c * f1_c;
    drift->Dv2 = D2_c * ac * ac * E_c * f2_c;
}

void
fastpm_drift_store(FastPMDriftFactor * drift,
               FastPMStore * pi, FastPMStore * po,
               double af)
{
    int np = pi->np;

    int i;
    // Drift
#pragma omp parallel for
    for(i=0; i<np; i++) {
        double xo[3] = {0};
        fastpm_drift_one(drift, pi, i, xo, af);
        int d;
        for(d = 0; d < 3; d ++) {
            po->x[i][d] = xo[d];
        }
    }
    po->meta.a_x = af;
}

//
// Functions for our modified time-stepping (used when StdDA=0):
//
struct iparam {
    FastPMCosmology * cosmology;
    double nLPT;
};

double gpQ(double a, double nLPT) {
    return pow(a, nLPT);
}

static double stddriftfunc (double a, struct iparam * iparam) {
    return 1 / (pow(a, 3) * HubbleEa(a, iparam->cosmology));
}

static double nonstddriftfunc (double a, struct iparam * iparam) {
    return gpQ(a, iparam->nLPT)/(pow(a, 3) * HubbleEa(a, iparam->cosmology));
}

static double stdkickfunc (double a, struct iparam * iparam) {
    return 1/ (pow(a, 2) * HubbleEa(a, iparam->cosmology));
}

static double integrand(double a, void * params) {
    void ** p = (void**) params;
    double (*func)(double a, struct iparam * s) = p[0];
    struct iparam * s = p[1];
    return func(a, s);
}

double integrate(double ai, double af,
        struct iparam * iparam,
        double (*func)(double , struct iparam * )) {

    gsl_integration_workspace * w
        = gsl_integration_workspace_alloc (5000);

    gsl_function F;
    double error;
    double result;

    F.params = (void*[]){func, iparam};
    F.function = integrand;

    gsl_integration_qag (&F, ai, af, 0, 1e-8, 5000, 6,
            w, &result, &error);

    gsl_integration_workspace_free (w);
    return result;
}

/*
       When StdDA=0, one needs to set nLPT.
       assumes time dep. for velocity = B a^nLPT
       nLPT is a real number. Sane values lie in the range (-4,3.5). Cannot be 0, but of course can be -> 0 (say 0.001).
       See Section A.3 of TZE.
       */

static double
Sq(double ai, double af, double aRef, double nLPT, FastPMCosmology * c, int USE_NONSTDDA)
{
    double resultstd, result;
    struct iparam iparam[1];
    iparam->cosmology = c;
    iparam->nLPT = nLPT;

    resultstd = integrate(ai, af, iparam, stddriftfunc);

    result = integrate(ai, af, iparam, nonstddriftfunc);
    result /= gpQ(aRef, nLPT);

    /*
    fastpm_info("ref time = %6.4f, std drift =%g, non std drift = %g \n",
        aRef, resultstd, result); */

    if (USE_NONSTDDA)
        return result;
    else
        return resultstd;
}

double DERgpQ(double a, double nLPT) {
    /* This must return d(gpQ)/da */
    return nLPT*pow(a, nLPT-1);
}

static double
Sphi(double ai, double af, double aRef, double nLPT, FastPMCosmology * c, int USE_NONSTDDA)
{
    double result;
    double resultstd;

    struct iparam iparam[1];
    iparam->cosmology = c;
    iparam->nLPT = nLPT;

    result = (gpQ(af, nLPT) - gpQ(ai, nLPT)) * aRef
        / (pow(aRef, 3) * HubbleEa(aRef, c) * DERgpQ(aRef, nLPT));

    resultstd = integrate(ai, af, iparam, stdkickfunc);

    /*
    fastpm_info("ref time = %6.4f, std kick = %g, non std kick = %g\n",
            aRef, resultstd, result); */

    if (USE_NONSTDDA) {
        return result;
    } else {
        return resultstd;
    }
}

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
/* tz_3dgeom_utils.c
 *
 * 17-Jan-2008 Initial write: Ting Zhao
 */
#include <stdio.h>
#include <math.h>
#include "tz_geo3d_utils.h"
#ifdef HAVE_LIBGSL
#  if defined(HAVE_INLINE)
#    undef HAVE_INLINE
#    define INLINE_SUPPRESSED
#  endif
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#  if defined(INLINE_SUPPRESSED)
#    define HAVE_INLINE
#  endif
#endif
#include "tz_error.h"
#include "tz_utilities.h"
#include "tz_constant.h"
#include "tz_geoangle_utils.h"
#include "tz_geometry.h"
#include "tz_3dgeom.h"
#include "tz_darray.h"
#include "tz_math.h"

#define TZ_DIST_3D_EPS 0.0000000000001 /* rounding torelance */

/* Geo3d_Translate_Coordinate(): 3D coordinate translation.
 *
 * Args: x - x coordinate to translate. It stores the result after return;
 *       y - y coordinate to translate. It stores the result after return;
 *       z - z coordinate to translate. It stores the result after return;
 *       dx - x axis translation;
 *       dy - y axis translation;
 *       dz - z axis translation.
 *
 * Return: void.
 */
void Geo3d_Translate_Coordinate(double *x, double *y, double *z,
				double dx, double dy, double dz)
{
  *x += dx;
  *y += dy;
  *z += dz;
}

/* Geo3d_Coordinate_Offset(): 3D coordinates offset.
 *
 * Args: x1, y1, z1 - start coordinate;
 *       x2, y2, z2 - end coordinate;
 *       dx, dy, dz - store the result.
 *
 * Return: void.
 */
void Geo3d_Coordinate_Offset(double x1, double y1, double z1,
			     double x2, double y2, double z2,
			     double *dx, double *dy, double *dz)
{
  *dx = x2 - x1;
  *dy = y2 - y1;
  *dz = z2 - z1;
}

/*
 * Geo3d_Rotate_Coordinate(): 3D  rotation.
 *
 * Args: x, y, z - the point to be rotated;
 *       theta - angle around the X axis;
 *       psi - angle around the Z axis;
 *       reverse - reverse the rotation or not.
 *
 * Return: void.
 */
void Geo3d_Rotate_Coordinate(double *x, double *y, double *z,
			     double theta, double psi, BOOL reverse)
{
  double p[3] = {*x, *y, *z};
  Rotate_XZ(p, p, 1, theta, psi, reverse);
  *x = p[0];
  *y = p[1];
  *z = p[2];
}

/* Geo3d_Orientation_Normal(): calculate the normal vector at a certain
 *                             direction.
 *
 * Args: theta - radian around X axis;
 *       psi - radian around Z axis;
 *       x - resulted x coordinate;
 *       y - resulted y coordinate;
 *       z - resulted z coordinate.
 *
 * Return: void.
 */
void Geo3d_Orientation_Normal(double theta, double psi,
			      double *x, double *y, double *z)
{
  /* rotaion of the vector (0, 0, 1) */
  double sin_theta = sin(theta);
  *x = sin_theta * sin(psi);
  *y = -sin_theta * cos(psi);
  //*z = cos(theta);
  *z = sqrt(1.0 - sin_theta * sin_theta);
}

void Geo3d_Normal_Orientation(double x, double y, double z,
			      double *theta, double *psi)
{
  /*
  *theta = acos(z);
  *psi = Vector_Angle(x, y) + TZ_PI_2;
  */
  //*psi = Vector_Angle(x, y) + TZ_PI_2 * 3.0;

  /* The range of theta will be [-pi, pi] and
   * the range of psi will be [-pi/2, pi/2]. */
  *theta = acos(z);

  if (*theta < GEOANGLE_COMPARE_EPS) {
    *psi = 0.0;
  } else {
    if (y >= 0.0) {
      *theta = -*theta;
      *psi = Vector_Angle(x, y) - TZ_PI_2;
    } else {
      *psi = Vector_Angle(x, y) - TZ_PI_2 * 3.0;
    }
  }
}

void Geo3d_Coord_Orientation(double x, double y, double z,
			      double *theta, double *psi)
{
  double r = Geo3d_Orgdist(x, y, z);
  if (r > TZ_DIST_3D_EPS) {
    x /= r;
    y /= r;
    z /= r;
    Geo3d_Normal_Orientation(x, y, z, theta, psi);
  } else {
    *theta = 0.0;
    *psi = 0.0;
  }
}

void Geo3d_Rotate_Orientation(double rtheta, double rpsi,
			      double *theta, double *psi)
{
  double coord[3];
  Geo3d_Orientation_Normal(*theta, *psi, coord, coord + 1, coord + 2);
  Rotate_XZ(coord, coord, 1, rtheta, rpsi, 0);
  Geo3d_Normal_Orientation(coord[0], coord[1], coord[2], theta, psi);
}

double Geo3d_Dot_Product(double x1, double y1, double z1,
			 double x2, double y2, double z2)
{
  return x1 * x2 + y1 * y2 + z1 * z2;
}

void Geo3d_Cross_Product(double x1, double y1, double z1,
			 double x2, double y2, double z2,
			 double *x, double *y, double *z)
{
  *x = y1 * z2 - y2 * z1;
  *y = x2 * z1 - x1 * z2;
  *z = x1 * y2 - x2 * y1;
}

double Geo3d_Orgdist_Sqr(double x, double y, double z)
{
  return x * x + y * y + z * z;
}

double Geo3d_Orgdist(double x, double y, double z)
{
  return sqrt(x * x + y * y + z * z);
}

double Geo3d_Dist_Sqr(double x1, double y1, double z1,
		      double x2, double y2, double z2)
{
  double dx = x1 - x2;
  double dy = y1 - y2;
  double dz = z1 - z2;

  return dx * dx + dy * dy + dz * dz;
}

double Geo3d_Dist(double x1, double y1, double z1,
		  double x2, double y2, double z2)
{
  double dx = x1 - x2;
  double dy = y1 - y2;
  double dz = z1 - z2;

  return sqrt(dx * dx + dy * dy + dz * dz);
}

double Geo3d_Angle2(double x1, double y1, double z1,
		    double x2, double y2, double z2)
{
  double d1 = Geo3d_Orgdist_Sqr(x1, y1, z1);
  if (d1 == 0.0) {
    return 0.0;
  }

  double d2 = Geo3d_Orgdist_Sqr(x2, y2, z2);
  if (d2 == 0.0) {
    return 0.0;
  }

  double d12 = Geo3d_Dot_Product(x1, y1, z1, x2, y2, z2) / sqrt(d1 * d2);

  /* to avoid invalid value caused by rounding error */
  d12 = (fabs(d12) > 1.0) ? round(d12) : d12;

  return acos(d12);
}

void Geo3d_Lineseg_Break(const double *line_start, const double *line_end,
			 double lambda, double *point)
{
  int i;
  for (i = 0; i < 3; i++) {
    point[i] = (1 - lambda) * line_start[i] + lambda * line_end[i];
  }
}

#define GEO3D_LINESEG_DIST2_EPS TZ_DIST_3D_EPS

/*
 * reference: http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline3d/
 */
static int geo3d_line_line_dist(const double *line1_start,
				const double *line1_end,
				const double *line2_start,
				const double *line2_end,
				double *dc1, double *dc2,
				double *intersect1, double *intersect2,
				double eps)
{
  double ds[3];  /*13*/
  double de[3];  /*24*/

  int i;
  for (i = 0; i < 3; i++) {
    dc1[i] = line1_end[i] - line1_start[i]; /*21*/
    dc2[i] = line2_end[i] - line2_start[i]; /*43*/
    ds[i] = line1_start[i] - line2_start[i];
    de[i] = line1_end[i] - line2_end[i];
  }

  double d2121 = dc1[0] * dc1[0] + dc1[1] * dc1[1] + dc1[2] * dc1[2];
  double d4343 = dc2[0] * dc2[0] + dc2[1] * dc2[1] + dc2[2] * dc2[2];

  if (d2121 < eps) {
    if (d4343 < eps) {
      return -3;
    }

    return -1;
  }

  if (d4343 < eps) {
    return -2;
  }

  double d4321 = dc2[0] * dc1[0] + dc2[1] * dc1[1] + dc2[2] * dc1[2];

  double denom = d2121 * d4343 - d4321 * d4321;

  double d1343 = ds[0] * dc2[0] + ds[1] * dc2[1] + ds[2] * dc2[2];
  double d1321 = ds[0] * dc1[0] + ds[1] * dc1[1] + ds[2] * dc1[2];

  double numer = d1343 * d4321 - d1321 * d4343;

  /* parallel */
  if (denom <= eps) {
    return 0;
  }

  *intersect1 = numer / denom;
  *intersect2 = (d1343 + d4321 * (*intersect1)) / d4343;

  return 1;
}

static int geo3d_point_line_dist(const double *point,
				 const double *line_start,
				 const double *line_end,
				 double *lamda, double eps)
{
  double line_vec[3];
  double point_vec[3];
  int i;
  for (i = 0; i < 3; i++) {
    line_vec[i] = line_end[i] - line_start[i];
    point_vec[i] = point[i] - line_start[i];
  }

  double length_sqr = Geo3d_Orgdist_Sqr(line_vec[0], line_vec[1], line_vec[2]);

  if (length_sqr < eps) {
    return 0;
  }

  /* projection from point to line */
  double dot = point_vec[0] * line_vec[0] + point_vec[1] * line_vec[1] +
    point_vec[2] * line_vec[2];

  *lamda = dot / length_sqr;

  return 1;
}

double Geo3d_Line_Line_Dist(double line1_start[], double line1_end[],
			    double line2_start[], double line2_end[])
{
  double intersect1, intersect2;
  double dc1[3], dc2[3];

  double lamda;

  switch (geo3d_line_line_dist(line1_start, line1_end, line2_start, line2_end,
			       dc1, dc2, &intersect1, &intersect2,
			       TZ_DIST_3D_EPS)) {
  case -3: /* point to point */
    return Geo3d_Dist(line1_start[0], line1_start[1], line1_start[2],
		      line2_start[0], line2_start[1], line2_start[2]);
  case -1: /* point to line2 */
    geo3d_point_line_dist(line1_start, line2_start, line2_end, &lamda,
			  TZ_DIST_3D_EPS);
    return Geo3d_Dist(line1_start[0], line1_start[1], line1_start[2],
		      line2_start[0] + lamda * dc2[0],
		      line2_start[1] + lamda * dc2[1],
		      line2_start[2] + lamda * dc2[2]);
  case -2: /* point to line 1 */
    geo3d_point_line_dist(line2_start, line1_start, line1_end, &lamda,
			  TZ_DIST_3D_EPS);
    return Geo3d_Dist(line2_start[0], line2_start[1], line2_start[2],
		      line1_start[0] + lamda * dc1[0],
		      line1_start[1] + lamda * dc1[1],
		      line1_start[2] + lamda * dc1[2]);
  case 0: /* parallel */
    geo3d_point_line_dist(line1_start, line2_start, line2_end, &lamda,
			  TZ_DIST_3D_EPS);
    return Geo3d_Dist(line1_start[0], line1_start[1], line1_start[2],
		      line2_start[0] + lamda * dc2[0],
		      line2_start[1] + lamda * dc2[1],
		      line2_start[2] + lamda * dc2[2]);
  default: /* line to line */
    return sqrt(Geo3d_Dist_Sqr(line1_start[0] + intersect1 * dc1[0],
			       line1_start[1] + intersect1 * dc1[1],
			       line1_start[2] + intersect1 * dc1[2],
			       line2_start[0] + intersect2 * dc2[0],
			       line2_start[1] + intersect2 * dc2[1],
			       line2_start[2] + intersect2 * dc2[2]));
  }
}

double Geo3d_Lineseg_Lineseg_Dist(double line1_start[], double line1_end[],
				  double line2_start[], double line2_end[],
				  double *intersect1, double *intersect2,
				  int *cond)
{
  double dc1[3]; /*21*/
  double dc2[3]; /*43*/

  double d1, d2;

  switch (geo3d_line_line_dist(line1_start, line1_end, line2_start, line2_end,
			       dc1, dc2, intersect1, intersect2,
			       TZ_DIST_3D_EPS)) {
  case 0: /* parallel */
    *cond = 9;
    d1 = Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				  intersect1);
    d2 = Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				  intersect2);
    if (d1 <= d2) {
      *intersect2 = *intersect1;
      *intersect1 = 1.0;
      return d1;
    } else {
      *intersect1 = 0.0;
      return d2;
    }
  case -1:
    *intersect1 = 0.0;
    *cond = 10;
    return Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				    intersect2);
  case -2:
    *intersect2= 0.0;
    *cond = 10;
    return Geo3d_Point_Lineseg_Dist(line2_start, line1_start, line1_end,
				    intersect1);
  case -3:
    *cond = 10;
    *intersect1 = 0.0;
    *intersect2 = 0.0;
    return Geo3d_Dist(line1_start[0], line1_start[1], line1_start[2],
		      line2_start[0], line2_start[1], line2_start[2]);
  default:
#if defined BREAK_IS_IN_RANGE
#  undef BREAK_IS_IN_RANGE
#endif
#define BREAK_IS_IN_RANGE(mu) ((mu >= 0.0) && (mu <= 1.0))

    if (BREAK_IS_IN_RANGE(*intersect1) && BREAK_IS_IN_RANGE(*intersect2)) {
      *cond = 0;
      return sqrt(Geo3d_Dist_Sqr(line1_start[0] + *intersect1 * dc1[0],
				 line1_start[1] + *intersect1 * dc1[1],
				 line1_start[2] + *intersect1 * dc1[2],
				 line2_start[0] + *intersect2 * dc2[0],
				 line2_start[1] + *intersect2 * dc2[1],
				 line2_start[2] + *intersect2 * dc2[2]));
    } else if ((*intersect1 < 0.0) && BREAK_IS_IN_RANGE(*intersect2)) {
      *cond = 1;
      *intersect1 = 0.0;
      return Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				      intersect2);
    } else if ((*intersect1 > 0.0) && BREAK_IS_IN_RANGE(*intersect2)) {
      *cond = 2;
      *intersect1 = 1.0;
      return Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				      intersect2);
    } else if ((*intersect2 < 0.0) && BREAK_IS_IN_RANGE(*intersect1)) {
      *cond = 3;
      *intersect2 = 0.0;
      return Geo3d_Point_Lineseg_Dist(line2_start, line1_start, line1_end,
				      intersect1);
    } else if ((*intersect2 > 1.0) && BREAK_IS_IN_RANGE(*intersect1)) {
      *cond = 4;
      *intersect2 = 1.0;
      return Geo3d_Point_Lineseg_Dist(line2_end, line1_start, line1_end,
				      intersect1);
    } else if ((*intersect1 < 0.0) && (*intersect2 < 0.0)) {
      *cond = 5;
      d1 = Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				    intersect2);
      d2 = Geo3d_Point_Lineseg_Dist(line2_start, line1_start, line1_end,
				    intersect1);
      if (d1 <= d2) {
	*intersect1 = 0.0;
	return d1;
      } else {
	*intersect2 = 0.0;
	return d2;
      }
    } else if ((*intersect1 > 1.0) && (*intersect2 < 0.0)) {
      *cond = 6;
      d1 = Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				    intersect2);
      d2 = Geo3d_Point_Lineseg_Dist(line2_start, line1_start, line1_end,
				    intersect1);
      if (d1 <= d2) {
	*intersect1 = 1.0;
	return d1;
      } else {
	*intersect2 = 0.0;
	return d2;
      }
    } else if ((*intersect1 < 0.0) && (*intersect2 > 1.0)) {
      *cond = 7;
      d1 = Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				    intersect2);
      d2 = Geo3d_Point_Lineseg_Dist(line2_end, line1_start, line1_end,
				    intersect1);
      if (d1 <= d2) {
	*intersect1 = 0.0;
	return d1;
      } else {
	*intersect2 = 1.0;
	return d2;
      }
    } else if ((*intersect1 > 1.0) && (*intersect2 > 1.0)) {
      *cond = 8;
      d1 = Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				    intersect2);
      d2 = Geo3d_Point_Lineseg_Dist(line2_end, line1_start, line1_end,
				    intersect1);
      if (d1<= d2) {
	*intersect1 = 1.0;
	return d1;
      } else {
	*intersect2 = 1.0;
	return d2;
      }
    }
  }

  return -1.0;
}

double Geo3d_Lineseg_Dist2(double line1_start[], double line1_end[],
			   double line2_start[], double line2_end[],
			   double *intersect1, double *intersect2,
			   int *cond)
{
  double dc1[3]; /*21*/
  double dc2[3]; /*43*/
  double ds[3];  /*13*/
  double de[3];  /*24*/

  int i;
  for (i = 0; i < 3; i++) {
    dc1[i] = line1_end[i] - line1_start[i];
    dc2[i] = line2_end[i] - line2_start[i];
    ds[i] = line1_start[i] - line2_start[i];
    de[i] = line1_end[i] - line2_end[i];
  }

  double d2121 = dc1[0] * dc1[0] + dc1[1] * dc1[1] + dc1[2] * dc1[2];
  double d4343 = dc2[0] * dc2[0] + dc2[1] * dc2[1] + dc2[2] * dc2[2];

  ASSERT(d2121 > GEO3D_LINESEG_DIST2_EPS, "invalid length");
  ASSERT(d4343 > GEO3D_LINESEG_DIST2_EPS, "invalid length");

  double d4321 = dc2[0] * dc1[0] + dc2[1] * dc1[1] + dc2[2] * dc1[2];

  double denom = d2121 * d4343 - d4321 * d4321;

  double d1343 = ds[0] * dc2[0] + ds[1] * dc2[1] + ds[2] * dc2[2];
  double d1321 = ds[0] * dc1[0] + ds[1] * dc1[1] + ds[2] * dc1[2];

  double numer = d1343 * d4321 - d1321 * d4343;

  if (denom <= GEO3D_LINESEG_DIST2_EPS) {
    *cond = 9;
    return
      dmin2(Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				     intersect1),
	    Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				     intersect2)
	    );
  }

  *intersect1 = numer / denom;
  *intersect2 = (d1343 + d4321 * (*intersect1)) / d4343;

#if defined BREAK_IS_IN_RANGE
#  undef BREAK_IS_IN_RANGE
#endif
#define BREAK_IS_IN_RANGE(mu) ((mu >= 0.0) && (mu <= 1.0))

  if (BREAK_IS_IN_RANGE(*intersect1) && BREAK_IS_IN_RANGE(*intersect2)) {
    *cond = 0;
    return sqrt(Geo3d_Dist_Sqr(line1_start[0] + *intersect1 * dc1[0],
			       line1_start[1] + *intersect1 * dc1[1],
			       line1_start[2] + *intersect1 * dc1[2],
			       line2_start[0] + *intersect2 * dc2[0],
			       line2_start[1] + *intersect2 * dc2[1],
			       line2_start[2] + *intersect2 * dc2[2]));
  } else if ((*intersect1 < 0.0) && BREAK_IS_IN_RANGE(*intersect2)) {
    *cond = 1;
    return Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				    intersect2);
  } else if ((*intersect1 > 0.0) && BREAK_IS_IN_RANGE(*intersect2)) {
    *cond = 2;
    return Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				    intersect2);
  } else if ((*intersect2 < 0.0) && BREAK_IS_IN_RANGE(*intersect1)) {
    *cond = 3;
    return Geo3d_Point_Lineseg_Dist(line2_start, line1_start, line1_end,
				    intersect1);
  } else if ((*intersect2 > 1.0) && BREAK_IS_IN_RANGE(*intersect1)) {
    *cond = 4;
    return Geo3d_Point_Lineseg_Dist(line2_end, line1_start, line1_end,
				    intersect1);
  } else if ((*intersect1 < 0.0) && (*intersect2 < 0.0)) {
    *cond = 5;
    return
      dmin2(Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				     intersect2),
	    Geo3d_Point_Lineseg_Dist(line2_start, line1_start, line1_end,
				     intersect1)
	    );
  } else if ((*intersect1 > 1.0) && (*intersect2 < 0.0)) {
    *cond = 6;
    return
      dmin2(Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				     intersect2),
	    Geo3d_Point_Lineseg_Dist(line2_start, line1_start, line1_end,
				     intersect1)
	    );
  } else if ((*intersect1 < 0.0) && (*intersect2 > 1.0)) {
    *cond = 7;
    return
      dmin2(Geo3d_Point_Lineseg_Dist(line1_start, line2_start, line2_end,
				     intersect2),
	    Geo3d_Point_Lineseg_Dist(line2_end, line1_start, line1_end,
				     intersect1)
	    );
  } else if ((*intersect1 > 1.0) && (*intersect2 > 1.0)) {
    *cond = 8;
    return
      dmin2(Geo3d_Point_Lineseg_Dist(line1_end, line2_start, line2_end,
				     intersect2),
	    Geo3d_Point_Lineseg_Dist(line2_end, line1_start, line1_end,
				     intersect1)
	    );
  }

  return -1.0;
}

#define GEO3D_POINT_LINESEG_EPS 0.0000000000001
double Geo3d_Point_Lineseg_Dist(const double *point, const double *line_start,
				const double *line_end, double *lamda)
{
  double tmp_lamda;
  int status = geo3d_point_line_dist(point, line_start, line_end, &tmp_lamda,
				     GEO3D_POINT_LINESEG_EPS);

  double dist;
  if (status == 0) {
    dist = Geo3d_Dist(line_start[0], line_start[1], line_start[2],
		      point[0], point[1], point[2]);
  } else {
    if ((tmp_lamda >= 0.0) && (tmp_lamda <= 1.0)) {
      dist =
	Geo3d_Dist(point[0], point[1], point[2],
		   (1.0 - tmp_lamda) * line_start[0] + tmp_lamda * line_end[0],
		   (1.0 - tmp_lamda) * line_start[1] + tmp_lamda * line_end[1],
		   (1.0 - tmp_lamda) * line_start[2] + tmp_lamda * line_end[2]);
    } else if (tmp_lamda < 0.0) {
      tmp_lamda = 0.0;
      dist = sqrt(Geo3d_Dist_Sqr(line_start[0], line_start[1], line_start[2],
				 point[0], point[1], point[2]));
    } else {
      tmp_lamda = 1.0;
      dist = sqrt(Geo3d_Dist_Sqr(line_end[0], line_end[1], line_end[2],
				 point[0], point[1], point[2]));
    }
  }

  /*
  double line_vec[3];
  double point_vec[3];
  int i;
  for (i = 0; i < 3; i++) {
    line_vec[i] = line_end[i] - line_start[i];
    point_vec[i] = point[i] - line_start[i];
  }
  double length_sqr = Geo3d_Orgdist_Sqr(line_vec[0], line_vec[1], line_vec[2]);

  double tmp_lamda = 0.0;
  double dist = -1.0;

  if (length_sqr <= GEO3D_POINT_LINESEG_EPS) {
    dist = sqrt(Geo3d_Dist_Sqr(line_start[0], line_start[1], line_start[2],
			       point[0], point[1], point[2]));
  } else {
    double dot = point_vec[0] * line_vec[0] + point_vec[1] * line_vec[1] +
      point_vec[2] * line_vec[2];

    double point_vec_lensqr = Geo3d_Orgdist_Sqr(point_vec[0], point_vec[1],
					      point_vec[2]);

    if (point_vec_lensqr <= GEO3D_POINT_LINESEG_EPS) {
      dist = 0.0;
    } else {
      tmp_lamda = dot / length_sqr;
      if ((tmp_lamda >= 0.0) && (tmp_lamda <= 1.0)) {
	dist = sqrt(point_vec_lensqr - dot * dot / length_sqr);
      } else if (tmp_lamda < 0.0) {
	dist = sqrt(Geo3d_Dist_Sqr(line_start[0], line_start[1], line_start[2],
				   point[0], point[1], point[2]));
      } else {
	dist = sqrt(Geo3d_Dist_Sqr(line_end[0], line_end[1], line_end[2],
				   point[0], point[1], point[2]));
      }
    }
  }
  */
  if (lamda != NULL) {
    *lamda = tmp_lamda;
  }

  return dist;
}

void Geo3d_Cov_Orientation(double *cov, double *vec, double *ext)
{
#ifdef HAVE_LIBGSL
  gsl_matrix_view gmv = gsl_matrix_view_array(cov, 3, 3);

  gsl_vector *eval = gsl_vector_alloc(3);
  gsl_matrix *evec = gsl_matrix_alloc(3,3);
  gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(3);

  gsl_eigen_symmv(&(gmv.matrix), eval, evec, w);
  gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);

  if(ext != NULL) {
    darraycpy(ext, gsl_vector_const_ptr(eval, 0), 0, 3);
  }

  gsl_vector_view v = gsl_vector_view_array(vec, 3);
  gsl_matrix_get_col(&(v.vector), evec, 0);

  gsl_vector_free(eval);
  gsl_matrix_free(evec);
  gsl_eigen_symmv_free(w);
#else
  double eigv[3];
  Matrix_Eigen_Value_Cs(cov[0], cov[4], cov[8], cov[1], cov[2], cov[5], eigv);
  Matrix_Eigen_Vector_Cs(cov[0], cov[4], cov[8], cov[1], cov[2], cov[5],
      eigv[0], vec);
  if(ext != NULL) {
    darraycpy(ext, eigv, 0, 3);
  }
#endif
}


text:
// Copyright (C) 2021 Intel Corporation
// SPDX-License-Identifier: Apache-2.0

#ifndef _HEBench_ClearText_MatMult_H_7e5fa8c2415240ea93eff148ed73539b
#define _HEBench_ClearText_MatMult_H_7e5fa8c2415240ea93eff148ed73539b

#include <gsl/gsl>

#include "hebench/api_bridge/cpp/hebench.hpp"

#include "clear_benchmark.h"

template <class T>
class MatMult_Benchmark : public ClearTextBenchmark
{
private:
    HEBERROR_DECLARE_CLASS_NAME(MatMult_Benchmark)

public:
    MatMult_Benchmark(hebench::cpp::BaseEngine &engine,
                      const hebench::APIBridge::BenchmarkDescriptor &bench_desc,
                      const hebench::APIBridge::WorkloadParams &bench_params);
    ~MatMult_Benchmark() override;

    hebench::APIBridge::Handle encode(const hebench::APIBridge::PackedData *p_parameters) override;
    void decode(hebench::APIBridge::Handle encoded_data, hebench::APIBridge::PackedData *p_native) override;

    hebench::APIBridge::Handle load(const hebench::APIBridge::Handle *p_local_data, std::uint64_t count) override;
    void store(hebench::APIBridge::Handle remote_data,
               hebench::APIBridge::Handle *p_local_data, std::uint64_t count) override;

    hebench::APIBridge::Handle operate(hebench::APIBridge::Handle h_remote_packed,
                                       const hebench::APIBridge::ParameterIndexer *p_param_indexers) override;

    //std::int64_t classTag() const override { return BaseBenchmark::classTag() | ClearTextBenchmark::tag; }

protected:
    std::array<std::uint64_t, 2> m_rows;
    std::array<std::uint64_t, 2> m_cols;

private:
    static void matMult(gsl::span<T> &result,
                        const gsl::span<const T> &M0, const gsl::span<const T> &M1,
                        std::size_t columns);
};

#include "inl/bench_matmult.inl"

#endif // defined _HEBench_ClearText_MatMult_H_7e5fa8c2415240ea93eff148ed73539b

max_stars_repo_name: jlakness-intel/backend-cpu-helib

max_issues_repo_name: jlakness-intel/backend-cpu-helib

max_forks_repo_name: jlakness-intel/backend-cpu-helib


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

#include "cosmocalc.h"

double weff(double a)
{
  if(a != 1.0)
    return cosmoData.w0 + cosmoData.wa - cosmoData.wa*(a - 1.0)/log(a);
  else
    return cosmoData.w0;
}

double hubble_noscale(double a)
{
  return sqrt(cosmoData.OmegaM/a/a/a + cosmoData.OmegaK/a/a + cosmoData.OmegaL/pow(a,3.0*(1.0 + weff(a))));
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#pragma once
#include <gsl/gsl_vector.h>
#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include "NodeVisitor.h"
#include "VarStore.h"
#include <map>
#include "FloatSupport.h"
#include <glpk.h>

#include <iostream>
#include <fstream>

using namespace std;

class Line { // m(x - x0) + d < 0 (x0 is fixed for a particular evaluator, hence not stored in the line)
public:
	double d;
	gsl_vector* m;
	Line(double _d, gsl_vector* _m): d(_d), m(_m) {}
	~Line(void) {
		delete m;
	}

	double getOffset(const gsl_vector* x0) { // m.x0 - d
		double r = 0.0;
		for (int i = 0; i < x0->size; i++) {
			r += gsl_vector_get(x0, i) * gsl_vector_get(m, i);
		}
		return r - d;
	}

	double getCoeff(int i) {
		return gsl_vector_get(m, i);
	}

	double getDist() {
		double r = 0;
		for (int i = 0; i < m->size; i++) {
			r += gsl_vector_get(m, i) * gsl_vector_get(m, i);
		}
		return d/sqrt(r);
	}

	string print() {
		stringstream str;
		str << "<";
		for (int i = 0; i < m->size; i++) {
			str << gsl_vector_get(m, i) << " ";
		}
		str << d << ">";
		return str.str();
	}
};

class Region { // Intersection of lines
public:
	vector<Line*> lines;
	int size() {
		return lines.size();
	}
	Line* get(int i) {
		return lines[i];
	}
	void addLine(Line* l) {
		if (find(lines.begin(), lines.end(), l) == lines.end()) {
			lines.push_back(l);
		}
	}

	static Region* intersect(Region* r1, Region* r2) {
		Region* r = new Region();
		for (int i = 0; i < r1->size(); i++) {
			r->addLine(r1->get(i));
		}
		for (int i = 0; i < r2->size(); i++) {
			r->addLine(r2->get(i));
		}
		return r;
	}

	static Region* copy(Region* r1) {
		Region* r = new Region();
		for (int i = 0; i < r1->size(); i++) {
			r->addLine(r1->get(i));
		}
		return r;
	}

	double getMaxDist() {
		double d = - numeric_limits<double>::max();
		for (int i = 0; i < lines.size(); i++) {
			if (lines[i]->getDist() > d) {
				d = lines[i]->getDist();
			}
		}
		return d;
	}

	gsl_vector* getMaxGrad() {
		double d = - numeric_limits<double>::max();
		gsl_vector* g;
		for (int i = 0; i < lines.size(); i++) {
			if (lines[i]->d > d) {
				d = lines[i]->d;
				g = lines[i]->m;
			}
		}
		return g;

	}

	~Region(void) {
		for (int i = 0; i < lines.size(); i++) {
			delete lines[i];
		}
	}

	string print() {
		stringstream str;
		str << "{";
		for (int i = 0; i < lines.size(); i++) {
			str << lines[i]->print() << " ";
		}
		str << "}";
		return str.str();
	}

};

class Value {
public:
	double val;
	gsl_vector* grad;
	Value(double v, gsl_vector* g): val(v), grad(g) {}
	~Value(void) {
		delete grad;
	}

	double getOffset(const gsl_vector* x0) { // d - m.x0
		double r = 0.0;
		for (int i = 0; i < x0->size; i++) {
			r += gsl_vector_get(x0, i) * gsl_vector_get(grad, i);
		}
		return val - r;
	}

	double getCoeff(int i) {
		return gsl_vector_get(grad, i);
	}

	static Value* add(Value* v1, Value* v2);
	static Value* mult(Value* v1, Value* v2);
	static Value* div(Value* v1, Value* v2);
	static Value* neg(Value* v1);
	static Value* copy(Value* v1);
	static Value* arctan(Value* v1);
	static Value* tan(Value* v1);
	static Value* cos(Value* v1);
	static Value* sin(Value* v1);
	static Value* sqrt(Value* v1);
	static Value* square(Value* v1);

	string print() {
		stringstream str;
		str << "<";
		for (int i = 0; i < grad->size; i++) {
			str << gsl_vector_get(grad, i) << " ";
		}
		str << val << ">";
		return str.str();
	}
};

class RegionValuePair {
public:
	Region* r;
	Value* v;
	RegionValuePair(Region* _r, Value* _v): r(_r), v(_v) {}

	string print() {
		stringstream str;
		str << r->print() << " -> " << v->print();
		return str.str();
	}
};

typedef vector<RegionValuePair*> ValuesList;

class LP {
	glp_prob* lp;
	int ia[1+1000], ja[1+1000];
	double ar[1+1000];
	double LOWBND = -32.0; // TODO: don;t hardcode
	double HIGHBND = 32.0;
	double PRECISION = 1e-5;

public:
	double optimize(RegionValuePair* rv, const gsl_vector* x0) {
		lp = glp_create_prob();
		glp_set_obj_dir(lp, GLP_MIN);
		int numrows = rv->r->size();
		if (numrows > 0) {
			glp_add_rows(lp, numrows);
			for (int i = 0; i < numrows; i++) {
				glp_set_row_name(lp, (i+1), "r" + (i+1));
				double c = rv->r->lines[i]->getOffset(x0);
				glp_set_row_bnds(lp, (i+1), GLP_UP, 0.0, c - PRECISION);
			}
		}
		int numcols = x0->size;
		glp_add_cols(lp, numcols);
		for (int i = 0; i < numcols; i++) {
			glp_set_col_name(lp, (i+1), "x" + (i+1));
			glp_set_col_bnds(lp, (i+1), GLP_DB, LOWBND, HIGHBND);
			glp_set_obj_coef(lp, (i+1), rv->v->getCoeff(i));
		}
		int counter = 1;
		for (int i = 0; i < numrows; i++) {
			for (int j = 0; j < numcols; j++) {
				ia[counter] = (i+1);
				ja[counter] = (j+1);
				ar[counter] = rv->r->lines[i]->getCoeff(j);
				counter++;
			}
		}
		Assert(counter-1 <= 1000, "Out of bounds");

		glp_load_matrix(lp, counter-1, ia, ja, ar);
		glp_simplex(lp, NULL);
		int status = glp_get_prim_stat(lp);
		cout << "status: " << status << endl;
		if (status == GLP_FEAS) {
			double z = glp_get_obj_val(lp);
			cout << "z = " << z;
			for (int i = 0; i < numcols; i++) {
				cout << " x" << i << " = " << glp_get_col_prim(lp, (i+1));
			}
			glp_delete_prob(lp);
			return z + rv->v->getOffset(x0);
		} else {
			return 1000;
		}

	}

	bool check(Region* r, const gsl_vector* x0) {
		int numrows = r->size();
		if (numrows == 0) return true;

		lp = glp_create_prob();
		glp_set_obj_dir(lp, GLP_MIN);
		glp_add_rows(lp, numrows);
		for (int i = 0; i < numrows; i++) {
			stringstream str;
			str << "r" << (i+1);
			glp_set_row_name(lp, (i+1), str.str().c_str());
			double c = r->lines[i]->getOffset(x0);
			glp_set_row_bnds(lp, (i+1), GLP_UP, 0.0, c-PRECISION);
		}
		int numcols = x0->size;
		glp_add_cols(lp, numcols);
		for (int i = 0; i < numcols; i++) {
			stringstream str;
			str << "x" << (i+1);
			glp_set_col_name(lp, (i+1), str.str().c_str());
			glp_set_col_bnds(lp, (i+1), GLP_DB, LOWBND, HIGHBND);
			glp_set_obj_coef(lp, (i+1), 0.0);
		}
		int counter = 1;
		for (int i = 0; i < numrows; i++) {
			for (int j = 0; j < numcols; j++) {
				ia[counter] = (i+1);
				ja[counter] = (j+1);
				ar[counter] = r->lines[i]->getCoeff(j);
				counter++;
			}
		}

		glp_load_matrix(lp, counter-1, ia, ja, ar);
		glp_smcp param;
		glp_init_smcp(&param);
		param.msg_lev = GLP_MSG_ERR;
		glp_simplex(lp, &param);
		int status = glp_get_prim_stat(lp);
		glp_delete_prob(lp);
		if (status == GLP_FEAS) {
			return true;
		} else {
			return false;
		}
	}

};

class GlobalEvaluator: NodeVisitor
{
	FloatManager& floats;
	BooleanDAG& bdag;
  map<string, int> floatCtrls; // Maps float ctrl names to indices with grad vectors
  int nctrls; // number of float ctrls
	gsl_vector* ctrls; // Ctrl values
	vector<ValuesList> values; // Keeps track of values for each node for different regions
	map<int, int> inputValues; // Maps node id to values set by the SAT solver
	double error = 0.0;
	gsl_vector* errorGrad;

	int DEFAULT_INP = -1;
	static gsl_vector* tmp;
	LP* lp;

	int MAX_REGIONS = 16;
	int ctr = 0;
public:
  GlobalEvaluator(BooleanDAG& bdag_p, FloatManager& _floats, const map<string, int>& floatCtrls_p);
  ~GlobalEvaluator(void);

  virtual void visit( SRC_node& node );
  virtual void visit( DST_node& node );
  virtual void visit( CTRL_node& node );
  virtual void visit( PLUS_node& node );
  virtual void visit( TIMES_node& node );
  virtual void visit( ARRACC_node& node );
  virtual void visit( DIV_node& node );
  virtual void visit( MOD_node& node );
  virtual void visit( NEG_node& node );
  virtual void visit( CONST_node& node );
  virtual void visit( LT_node& node );
  virtual void visit( EQ_node& node );
	virtual void visit( AND_node& node );
	virtual void visit( OR_node& node );
	virtual void visit( NOT_node& node );
  virtual void visit( ARRASS_node& node );
  virtual void visit( UFUN_node& node );
  virtual void visit( TUPLE_R_node& node );
	virtual void visit( ASSERT_node& node );

  double run(const gsl_vector* ctrls_p, map<int, int>& inputValues_p, gsl_vector* errorGrad_p);
	void truncate(ValuesList& vl, const gsl_vector* x0, int numRegions);
	void truncate(vector<Region*>& vl, const gsl_vector* x0, int numRegions);
	void truncateUnfeasibleRegions(ValuesList& vl, const gsl_vector* x0);
	void truncateUnfeasibleRegions(vector<Region*>& vl, const gsl_vector* x0);
	void truncateFarRegions(ValuesList& vl, const gsl_vector* x0, int numRegions);
	void truncateFarRegions(vector<Region*>& vl, const gsl_vector* x0, int numRegions);

	void setvalues(bool_node& bn, ValuesList& v) {
		values[bn.id] = v;
	}

	ValuesList& getvalues(bool_node& bn) {
		return values[bn.id];
	}

	ValuesList& getvalues(bool_node* bn) {
		return getvalues(*bn);
	}

	void print() {
    for (int i = 0; i < bdag.size(); i++) {
      cout << bdag[i]->lprint() << endl;
		}
  }

	bool isFloat(bool_node& bn) {
		return (bn.getOtype() == OutType::FLOAT);
	}

	bool isFloat(bool_node* bn) {
		return (bn->getOtype() == OutType::FLOAT);
	}

	int getInputValue(bool_node& bn) {
		if (inputValues.find(bn.id) != inputValues.end()) {
			int val = inputValues[bn.id];
			Assert(val == 0 || val == 1, "NYI: Integer values");
			return val;
		} else {
			return DEFAULT_INP;
		}
	}

	int getInputValue(bool_node* bn) {
		return getInputValue(*bn);
	}

	static gsl_vector* default_grad(int nctrls) {
		gsl_vector* g = gsl_vector_alloc(nctrls);
		for (int i = 0; i < nctrls; i++) {
			gsl_vector_set(g, i, 0);
		}
		return g;
	}

	static gsl_vector* add_grad(gsl_vector* g1, gsl_vector* g2);
	static gsl_vector* mult_grad(double v1, double v2, gsl_vector* g1, gsl_vector* g2);
	static gsl_vector* div_grad(double v1, double v2, gsl_vector* g1, gsl_vector* g2);
	static gsl_vector* neg_grad(gsl_vector* g1);
	static gsl_vector* copy_grad(gsl_vector* g1);
	static gsl_vector* sub_grad(gsl_vector* g1, gsl_vector* g2);
	static gsl_vector* arctan_grad(double v1, gsl_vector* g1);
	static gsl_vector* tan_grad(double v1, gsl_vector* g1);
	static gsl_vector* cos_grad(double v1, gsl_vector* g1);
	static gsl_vector* sin_grad(double v1, gsl_vector* g1);
	static gsl_vector* sqrt_grad(double v1, gsl_vector* g1);
	static gsl_vector* square_grad(double v1, gsl_vector* g1);

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <errno.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include "mpi.h"

//data specific parameters
#define KMEANS 3 //number of cluster groups
#define DIMENSIONS 4 //number of dimensions of data
#define DATA_ROWS 150

//run specific parameters
#define THREAD_COUNT 6 //count of thread used by OpenMP
#define FITNESS_FUNCTION fitness5 //specify fitness function used, 5 possibilities (fitness1,fitness2,...,fitness5)
#define SWARM_SIZE (200)
#define MAX_ITERATIONS (10000)
#define GEN_MUTATION 4 //how many parts of genes will be mutated
#define ELITE 4 //how many genes will be preserved without change

#define REAL long double
#define size_t int
#define ROOT 0

long double vectors[DATA_ROWS][DIMENSIONS]; //here are stored data vectors

long double fitness[SWARM_SIZE];
int chrom[SWARM_SIZE][DATA_ROWS];
int new_chromosoms[SWARM_SIZE][DATA_ROWS];

long double fit_data_vectors[KMEANS][DIMENSIONS][DATA_ROWS]; //used to compute fitness
#pragma omp threadprivate(fit_data_vectors)
char dataset[100];

/*
 * FITNEES FUNCTION
 * fitness1 - vseobecne kriterium (VVV)
 * fitness2 - spriemerovane kovariancne matice sa pouziju na fitness (EEE)
 * fitness3 - stopa kovariancnych matic (VII)
 * fitness4 - sucet stvorcov euklidovskych vzdialenosti fitness (KMEANS)
 * fitness5 - spriemerovane stopy kovariancnych matic (EII)
 */

long double my_mean(const long double data[], const size_t stride,
		const size_t size) {
	long double mean = 0;
	size_t i;

	for (i = 0; i < size; i++) {
		mean += (data[i * stride] - mean) / (i + 1);
	}

	return mean;
}

long double _my_covariance(const long double data1[], const size_t stride1,
		const long double data2[], const size_t stride2, const size_t n,
		const long double mean1, const long double mean2) {

	long double covariance = 0;

	size_t i;

	/* find the sum of the squares */
	for (i = 0; i < n; i++) {
		const long double delta1 = (data1[i * stride1] - mean1);
		const long double delta2 = (data2[i * stride2] - mean2);
		covariance += (delta1 * delta2 - covariance) / (i + 1);
	}

	return covariance;
}

long double my_gsl_stats_covariance(const long double data1[],
		const size_t stride1, const long double data2[], const size_t stride2,
		const size_t n) {
	const long double mean1 = my_mean(data1, stride1, n);
	const long double mean2 = my_mean(data2, stride2, n);

	return _my_covariance(data1, stride1, data2, stride2, n, mean1, mean2);
}

long double my_gsl_linalg_LU_det(gsl_matrix_long_double * LU, int signum) {
	size_t i, n = LU->size1;
	long double det = (long double) signum;

	for (i = 0; i < n; i++) {
		det *= gsl_matrix_long_double_get(LU, i, i);
	}

	return det;
}

void my_gsl_linalg_LU_decomp(gsl_matrix_long_double * A, gsl_permutation * p,
		int *signum) {
	if (A->size1 != A->size2) {
		printf("LU decomposition requires square matrix");
	} else if (p->size != A->size1) {
		printf("permutation length must match matrix size");
	} else {
		const size_t N = A->size1;
		size_t i, j, k;

		*signum = 1;
		gsl_permutation_init(p);

		for (j = 0; j < N - 1; j++) {
			/* Find maximum in the j-th column */

			long double ajj, max = fabs(gsl_matrix_long_double_get(A, j, j));
			size_t i_pivot = j;

			for (i = j + 1; i < N; i++) {
				long double aij = fabs(gsl_matrix_long_double_get(A, i, j));

				if (aij > max) {
					max = aij;
					i_pivot = i;
				}
			}

			if (i_pivot != j) {
				gsl_matrix_long_double_swap_rows(A, j, i_pivot);
				gsl_permutation_swap(p, j, i_pivot);
				*signum = -(*signum);
			}

			ajj = gsl_matrix_long_double_get(A, j, j);

			if (ajj != 0.0) {
				for (i = j + 1; i < N; i++) {
					long double aij = gsl_matrix_long_double_get(A, i, j) / ajj;
					gsl_matrix_long_double_set(A, i, j, aij);

					for (k = j + 1; k < N; k++) {
						long double aik = gsl_matrix_long_double_get(A, i, k);
						long double ajk = gsl_matrix_long_double_get(A, j, k);
						gsl_matrix_long_double_set(A, i, k, aik - aij * ajk);
					}
				}
			}
		}
	}
}

int random_int(int min_num, int max_num) {
	int result = 0;
	result = (rand() % (max_num - min_num)) + min_num;
	return result;
}

int my_ceil(float num) {
	int inum = (int) num;
	if (num == (float) inum) {
		return inum;
	}
	return inum + 1;
}

void my_export(char filename[100]) {
	int i;
	FILE *fp;
	fp = fopen(filename, "w");
	for (i = 0; i < DATA_ROWS; i++) {
		fprintf(fp, "%d\n", chrom[0][i] + 1);
	}
	fclose(fp);
}

long double fitness1(int i) {
	int j, k, h, s;
	long double det;
	int fit_cluster_volume[KMEANS];
	gsl_matrix_long_double *fit_matrix[KMEANS];
	gsl_permutation * p = gsl_permutation_alloc(DIMENSIONS);
	long double result = 0;

	for (j = 0; j < KMEANS; j++) {
		fit_cluster_volume[j] = 0;
	}

	for (j = 0; j < DATA_ROWS; j++) {
		for (h = 0; h < DIMENSIONS; h++) {
			fit_data_vectors[chrom[i][j]][h][fit_cluster_volume[chrom[i][j]]] =
					vectors[j][h];
		}

		fit_cluster_volume[chrom[i][j]]++;
	}

	for (j = 0; j < KMEANS; j++) {
		fit_matrix[j] = gsl_matrix_long_double_calloc(DIMENSIONS, DIMENSIONS);
		for (k = 0; k < DIMENSIONS; k++) {
			for (h = 0; h < DIMENSIONS; h++) {
				gsl_matrix_long_double_set(fit_matrix[j], k, h,
						my_gsl_stats_covariance(fit_data_vectors[j][k], 1,
								fit_data_vectors[j][h], 1,
								fit_cluster_volume[j]));
			}
		}
	}

	for (j = 0; j < KMEANS; j++) {
		s = 1;
		my_gsl_linalg_LU_decomp(fit_matrix[j], p, &s);
		det = my_gsl_linalg_LU_det(fit_matrix[j], s);
		if (det != 0.0) {
			result += gsl_sf_log_abs(det) * (double) fit_cluster_volume[j];
		}
		gsl_matrix_long_double_free(fit_matrix[j]);
	}
	return result;
}

long double fitness2(int i) {
	int j, k, h, s;
	long double m, det;
	int fit_cluster_volume[KMEANS];
	gsl_matrix_long_double *fit_matrix[KMEANS];
	gsl_matrix_long_double *fit_reuslt_matrix;
	gsl_permutation * p = gsl_permutation_alloc(DIMENSIONS);
	long double result = 0;

	for (j = 0; j < KMEANS; j++) {
		fit_cluster_volume[j] = 0;
	}

	for (j = 0; j < DATA_ROWS; j++) {
		for (h = 0; h < DIMENSIONS; h++) {
			fit_data_vectors[chrom[i][j]][h][fit_cluster_volume[chrom[i][j]]] =
					vectors[j][h];
		}

		fit_cluster_volume[chrom[i][j]]++;
	}

	for (j = 0; j < KMEANS; j++) {
		fit_matrix[j] = gsl_matrix_long_double_calloc(DIMENSIONS, DIMENSIONS);
		for (k = 0; k < DIMENSIONS; k++) {
			for (h = 0; h < DIMENSIONS; h++) {
				gsl_matrix_long_double_set(fit_matrix[j], k, h,
						my_gsl_stats_covariance(fit_data_vectors[j][k], 1,
								fit_data_vectors[j][h], 1,
								fit_cluster_volume[j]));
			}
		}
	}

	fit_reuslt_matrix = gsl_matrix_long_double_calloc(DIMENSIONS, DIMENSIONS);

	for (j = 0; j < KMEANS; j++) {
		for (k = 0; k < DIMENSIONS; k++) {
			for (h = 0; h < DIMENSIONS; h++) {

				m = gsl_matrix_long_double_get(fit_reuslt_matrix, k, h);
				m += fit_cluster_volume[j]
						* gsl_matrix_long_double_get(fit_matrix[j], k, h);
				gsl_matrix_long_double_set(fit_reuslt_matrix, k, h, m);
			}
		}
		gsl_matrix_long_double_free(fit_matrix[j]);
	}

	for (k = 0; k < DIMENSIONS; k++) {
		for (h = 0; h < DIMENSIONS; h++) {
			m = gsl_matrix_long_double_get(fit_reuslt_matrix, k, h);
			m = m / (double) DATA_ROWS;
			gsl_matrix_long_double_set(fit_reuslt_matrix, k, h, m);
		}
	}

	s = 1;
	my_gsl_linalg_LU_decomp(fit_reuslt_matrix, p, &s);
	det = my_gsl_linalg_LU_det(fit_reuslt_matrix, s);
	if (det != 0.0) {
		result += gsl_sf_log_abs(det);
	}
	gsl_matrix_long_double_free(fit_reuslt_matrix);
	return result;
}

long double fitness3(int i) {
	int j, k, h;
	int fit_cluster_volume[KMEANS];
	gsl_matrix_long_double *fit_matrix[KMEANS];
	long double result = 0;

	for (j = 0; j < KMEANS; j++) {
		fit_cluster_volume[j] = 0;
	}

	for (j = 0; j < DATA_ROWS; j++) {
		for (h = 0; h < DIMENSIONS; h++) {
			fit_data_vectors[chrom[i][j]][h][fit_cluster_volume[chrom[i][j]]] =
					vectors[j][h];
		}

		fit_cluster_volume[chrom[i][j]]++;
	}

	for (j = 0; j < KMEANS; j++) {
		fit_matrix[j] = gsl_matrix_long_double_calloc(DIMENSIONS, DIMENSIONS);
		for (k = 0; k < DIMENSIONS; k++) {
			for (h = 0; h < DIMENSIONS; h++) {
				gsl_matrix_long_double_set(fit_matrix[j], k, h,
						my_gsl_stats_covariance(fit_data_vectors[j][k], 1,
								fit_data_vectors[j][h], 1,
								fit_cluster_volume[j]));
			}
		}
	}

	for (j = 0; j < KMEANS; j++) {
		for (k = 0; k < DIMENSIONS; k++) {
			result += gsl_matrix_long_double_get(fit_matrix[j], k, k);
		}
		gsl_matrix_long_double_free(fit_matrix[j]);
	}
	return result;
}

long double fitness4(int agent) {
	int i, j, v;
	int weight_matrix[DATA_ROWS][KMEANS];
	long double center_matrix[KMEANS][DIMENSIONS];
	long double sum_w;
	long double sum_wx;
	long double result = 0;

	for (i = 0; i < DATA_ROWS; i++) {
		for (j = 0; j < KMEANS; j++) {
			if (chrom[agent][i] == j) {
				weight_matrix[i][j] = 1;
			} else {
				weight_matrix[i][j] = 0;
			}
		}
	}

	for (j = 0; j < KMEANS; j++) {
		for (v = 0; v < DIMENSIONS; v++) {
			sum_w = 0.0;
			for (i = 0; i < DATA_ROWS; i++) {
				sum_w += weight_matrix[i][j];
			}

			sum_wx = 0.0;
			for (i = 0; i < DATA_ROWS; i++) {
				sum_wx += weight_matrix[i][j] * vectors[i][v];
			}

			center_matrix[j][v] = sum_wx / sum_w;
		}
	}

	for (j = 0; j < KMEANS; j++) {
		for (i = 0; i < DATA_ROWS; i++) {
			sum_w = 0.0;
			for (v = 0; v < DIMENSIONS; v++) {
				sum_w += weight_matrix[i][j]
						* ((vectors[i][v] - center_matrix[j][v])
								* (vectors[i][v] - center_matrix[j][v]));
			}
			result += sqrt(sum_w);

		}
	}
	return result;
}

long double fitness5(int i) {
	int j, k, h;
	long double m;
	int fit_cluster_volume[KMEANS];
	gsl_matrix_long_double *fit_matrix[KMEANS];
	gsl_matrix_long_double *fit_reuslt_matrix;
	long double result = 0;

	for (j = 0; j < KMEANS; j++) {
		fit_cluster_volume[j] = 0;
	}

	for (j = 0; j < DATA_ROWS; j++) {
		for (h = 0; h < DIMENSIONS; h++) {
			fit_data_vectors[chrom[i][j]][h][fit_cluster_volume[chrom[i][j]]] =
					vectors[j][h];
		}

		fit_cluster_volume[chrom[i][j]]++;
	}

	for (j = 0; j < KMEANS; j++) {
		fit_matrix[j] = gsl_matrix_long_double_calloc(DIMENSIONS, DIMENSIONS);
		for (k = 0; k < DIMENSIONS; k++) {
			for (h = 0; h < DIMENSIONS; h++) {
				gsl_matrix_long_double_set(fit_matrix[j], k, h,
						my_gsl_stats_covariance(fit_data_vectors[j][k], 1,
								fit_data_vectors[j][h], 1,
								fit_cluster_volume[j]));
			}
		}
	}

	fit_reuslt_matrix = gsl_matrix_long_double_calloc(DIMENSIONS, DIMENSIONS);

	for (j = 0; j < KMEANS; j++) {
		for (k = 0; k < DIMENSIONS; k++) {
			for (h = 0; h < DIMENSIONS; h++) {

				m = gsl_matrix_long_double_get(fit_reuslt_matrix, k, h);
				m += fit_cluster_volume[j]
						* gsl_matrix_long_double_get(fit_matrix[j], k, h);
				gsl_matrix_long_double_set(fit_reuslt_matrix, k, h, m);
			}
		}
		gsl_matrix_long_double_free(fit_matrix[j]);
	}

	for (k = 0; k < DIMENSIONS; k++) {
		for (h = 0; h < DIMENSIONS; h++) {
			m = gsl_matrix_long_double_get(fit_reuslt_matrix, k, h);
			m = m / (double) DATA_ROWS;
			gsl_matrix_long_double_set(fit_reuslt_matrix, k, h, m);
		}
	}

	for (k = 0; k < DIMENSIONS; k++) {
		result += gsl_matrix_long_double_get(fit_reuslt_matrix, k, k);
	}
	gsl_matrix_long_double_free(fit_reuslt_matrix);

	return result;
}

void quicksort(long double list[], int ch[][DATA_ROWS], int n) {
	int i, j;
	long double pivot;
	long double temp;
	int temp_ch[DATA_ROWS];

	if (n < 2)
		return;
	pivot = list[n / 2];
	for (i = 0, j = n - 1;; i++, j--) {
		while (list[i] < pivot)
			i++;
		while (pivot < list[j])
			j--;
		if (i >= j)
			break;
		//swap
		temp = list[i];
		list[i] = list[j];
		list[j] = temp;
		memcpy(&temp_ch, &ch[i][0], sizeof(int) * DATA_ROWS);
		memcpy(&ch[i][0], &ch[j][0], sizeof(int) * DATA_ROWS);
		memcpy(&ch[j][0], &temp_ch, sizeof(int) * DATA_ROWS);

	}
	quicksort(list, ch, i);
	quicksort(list + i, ch + i, n - i);
}

int *int_array_2dto1d(int array[][DATA_ROWS], int first_dimension,
		int second_dimension) {
	int i, j;
	int *subarray = malloc(sizeof(int) * first_dimension * second_dimension);
	for (i = 0; i < first_dimension; i++) {
		for (j = 0; j < second_dimension; j++) {
			subarray[i * second_dimension + j] = array[i][j];
		}
	}
	return subarray;
}

void int_array_1dto2d(int *array, int solver, int elements_per_proc,
		int first_dimension, int second_dimension) {
	int i, j;
	for (i = 0; i < first_dimension; i++) {
		for (j = 0; j < second_dimension; j++) {
			chrom[solver * elements_per_proc + i][j] = array[i
					* second_dimension + j];
		}
	}
}

int main(int argc, char *argv[]) {
	int i, j, k, h, it;

	int temp_ch[DATA_ROWS];
	int *sub_chrom;
	int * chrom1d = NULL;
	long double *sub_fitness;

	unsigned long long int mating_h[SWARM_SIZE];
	unsigned long long int max_mating = 0;
	int mates[SWARM_SIZE];

	int r1, r2;

	time_t t, tt;

	int solver, solvers_count;
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &solver);
	MPI_Comm_size(MPI_COMM_WORLD, &solvers_count);

	MPI_Datatype my_type_scatter;
	MPI_Type_contiguous(DATA_ROWS, MPI_INT, &my_type_scatter);
	MPI_Type_commit(&my_type_scatter);

	if (SWARM_SIZE % 2 == 1) {
		printf("SWARM SIZE need to be dividable/partible by 2\n");
		exit(1);
	}

	if (argc >= 2) {
		strcpy(dataset, argv[1]);
	} else {
		printf("Program argument (dataset to read) not found.\n");
		exit(1);
	}

	if (SWARM_SIZE % solvers_count != 0) {
		printf(
				"SWARM_SIZE need to be divadable/partible by count of mpi solvers\n");
		exit(1);
	}

	srand(time(NULL));
	omp_set_num_threads(THREAD_COUNT);

	if (solver == ROOT) {
		for (i = 0; i < SWARM_SIZE / 2; i++) {
			mating_h[i] = gsl_pow_int((SWARM_SIZE / 2 - i), 5);
			max_mating += mating_h[i];
			mating_h[i + SWARM_SIZE / 2] = 0;
		}
	}

//read data
	FILE *fp;
	char *token;
	fp = fopen(dataset, "r");

	if (fp != NULL) {
		int lineNumber = 0;
		char line[1024];
		while (fgets(line, sizeof line, fp) != NULL) {
			token = strtok(line, ",");
			for (i = 0; i < DIMENSIONS - 1; i++) {
				vectors[lineNumber][i] = atof(token);
				token = strtok(NULL, ",");
			}
			vectors[lineNumber][DIMENSIONS - 1] = atof(token);
			lineNumber++;
		}
		fclose(fp);
	} else {
		printf("%s\n", strerror(errno));
		exit(1);
	}
//end of read data

	time(&t);

//Generate initial population.
	if (solver == ROOT) {
		for (i = 0; i < SWARM_SIZE; i++) {
			for (j = 0; j < DATA_ROWS; j++) {
				chrom[i][j] = random_int(0, KMEANS);
			}
		}
	}

	//Compute fitness of each chromosome.
#pragma omp parallel for shared(i,fitness)
	for (i = 0; i < SWARM_SIZE; i++) {
		fitness[i] = FITNESS_FUNCTION(i);
	}

	quicksort(fitness, chrom, SWARM_SIZE);

	for (it = 0; it < MAX_ITERATIONS; it++) {
		//Natural selection. Select mates.
		if (solver == ROOT) {
			for (i = 0; i < SWARM_SIZE; i++) {
				j = random_int(1, max_mating);
				h = 0;
				k = mating_h[h];
				while (j > k) {
					k += mating_h[++h];
				}
				mates[i] = h;
			}
		}

		//Mating
		if (solver == ROOT) {
			for (i = 0; i < SWARM_SIZE; i++) {
				do {
					r1 = random_int(0, DATA_ROWS);
					r2 = random_int(0, DATA_ROWS);
				} while (r1 == r2);

				if (r1 < r2) {
					j = r1;
					k = r2;
				} else {
					j = r2;
					k = r1;
				}

				//first kid
				memcpy(&temp_ch[0], &chrom[mates[i]][0], j * sizeof(int));
				memcpy(&temp_ch[j], &chrom[mates[i + 1]][j],
						(k - j) * sizeof(int));
				memcpy(&temp_ch[k], &chrom[mates[i]][k],
						(DATA_ROWS - k) * sizeof(int));
				memcpy(new_chromosoms[i], temp_ch, DATA_ROWS * sizeof(int));

				//second kid
				memcpy(&temp_ch[0], &chrom[mates[i + 1]][0], j * sizeof(int));
				memcpy(&temp_ch[j], &chrom[mates[i]][j], (k - j) * sizeof(int));
				memcpy(&temp_ch[k], &chrom[mates[i + 1]][k],
						(DATA_ROWS - k) * sizeof(int));
				memcpy(new_chromosoms[i + 1], temp_ch, DATA_ROWS * sizeof(int));
				i++;
			}
		}

		//copy new population
		if (solver == ROOT) {
			for (i = ELITE; i < SWARM_SIZE; i++) {
				memcpy(&chrom[i], &new_chromosoms[i], sizeof(int) * DATA_ROWS);
			}
		}

		//Mutation.
		if (solver == ROOT) {
			for (i = ELITE; i < SWARM_SIZE; i++) {
				for (j = 0; j < GEN_MUTATION; j++) {
					k = random_int(0, DATA_ROWS);
					chrom[i][k] = random_int(0, KMEANS);
				}
			}
		}

		//mpi section
		if (solver == ROOT) {
			chrom1d = int_array_2dto1d(chrom, SWARM_SIZE, DATA_ROWS);
		}

		sub_chrom = malloc(
				sizeof(int) * (SWARM_SIZE / solvers_count) * DATA_ROWS);

		MPI_Scatter(chrom1d,
		SWARM_SIZE / solvers_count, my_type_scatter, sub_chrom,
		SWARM_SIZE / solvers_count, my_type_scatter,
		ROOT,
		MPI_COMM_WORLD);

		int_array_1dto2d(sub_chrom, solver, SWARM_SIZE / solvers_count,
				SWARM_SIZE / solvers_count, DATA_ROWS);
		sub_fitness = malloc(
				sizeof(long double) * (SWARM_SIZE / solvers_count));

		//Compute fitness of each chromosome.
#pragma omp parallel for shared(i,fitness,sub_fitness)
		for (i = solver * (SWARM_SIZE / solvers_count);
				i
						< solver * (SWARM_SIZE / solvers_count)
								+ (SWARM_SIZE / solvers_count); i++) {
			sub_fitness[i - solver * (SWARM_SIZE / solvers_count)] =
					FITNESS_FUNCTION(i);
		}

		MPI_Gather(sub_fitness,
		SWARM_SIZE / solvers_count,
		MPI_LONG_DOUBLE, fitness,
		SWARM_SIZE / solvers_count,
		MPI_LONG_DOUBLE,
		ROOT,
		MPI_COMM_WORLD);

		free(sub_chrom);
		free(sub_fitness);
		if (solver == ROOT && chrom1d != NULL) {
			free(chrom1d);
		}
		//end of mpi section

		//Order population
		if (solver == ROOT) {
			quicksort(fitness, chrom, SWARM_SIZE);
		}
	}

	if (solver == ROOT) {

		time(&tt);
		printf("Program bezal %.f sekund.\n", difftime(tt, t));

		printf("fitness: %Lf\n", fitness[0]);

		if (argc >= 3) {
			my_export(argv[2]);
		} else {
			my_export("/home/lukas/workspace_parallel/mpi_gen/export.txt");
		}

		printf("SWARM SIZE: %d\n", SWARM_SIZE);
		printf("MAX_ITERATIONS: %d\n", MAX_ITERATIONS);
		printf("GEN_MUTATION: %d\n", GEN_MUTATION);
		printf("ELITE: %d\n", ELITE);
		printf("DIMENSIONS: %d\n", DIMENSIONS);
		printf("KMEANS: %d\n", KMEANS);
		printf("THREAD_COUNT per instance: %d\n", THREAD_COUNT);
		printf("mpi solvers_count: %d\n", solvers_count);
		printf(" FITNESS_FUNCTION ");

	}

	MPI_Finalize();

	return 0;

}

max_stars_repo_path: Impl of clus alg/GenClust/MPI_GenClust.c

max_stars_repo_name: PetoLau/GenClust

max_issues_repo_name: PetoLau/GenClust

max_forks_repo_name: PetoLau/GenClust


text:
#ifndef CONSUMER_H
#define CONSUMER_H

#include <gsl/gsl>
#include <producer.h>

class consumer
{
public:

    explicit consumer(gsl::not_null<producer *> p)
    { p->print_msg(); }

    ~consumer() = default;

public:

    consumer(consumer &&) noexcept = default;
    consumer &operator=(consumer &&) noexcept = default;

    consumer(const consumer &) = delete;
    consumer &operator=(const consumer &) = delete;
};

#endif


text:
#include <gsl/gsl_poly.h>

size_t remove_leading_zeros(double *h, size_t n)
{

    size_t m = 0;
    double small = 1e-5;
    size_t i;
    for(i=n; i>0; i--)
    {

        if (fabs(h[i-1]) > small)
        {
            m = i;
            break;
        }
    }

    return m;

}

size_t find_roots(double *z, double *h, size_t n)
{

    size_t m;
    m = remove_leading_zeros(&h[0], n);

    if (m > 1)
    {

        gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(m);
        gsl_poly_complex_solve(h, m, w, z);
        gsl_poly_complex_workspace_free (w);

    }

    return m;

}


text:
/*
 # Copyright (C) 2009-2013 Thierry Sousbie
 # University of Tokyo / CNRS, 2009
 #
 # This file is part of porject DisPerSE
 #
 #  Author          : Thierry Sousbie
 #  Contact         : tsousbie@gmail.com
 #
 #  Licenses        : This file is 'dual-licensed', you have to choose one
 #                    of the two licenses below to apply.
 #
 #                    CeCILL-C
 #                    The CeCILL-C license is close to the GNU LGPL.
 #                    ( http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html )
 #
 #                or  CeCILL v2.0
 #                    The CeCILL license is compatible with the GNU GPL.
 #                    ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
 #
 #  This software is governed either by the CeCILL or the CeCILL-C license
 #  under French law and abiding by the rules of distribution of free software.
 #  You can  use, modify and or redistribute the software under the terms of
 #  the CeCILL or CeCILL-C licenses as circulated by CEA, CNRS and INRIA
 #  at the following URL : "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and  rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #  same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL and CeCILL-C licenses and that you accept its terms.
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "simplex.h"

#ifdef HAVE_GSL

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

double SimplexVolumed(double **coord,int nvert, int ndims)
{
    int i,j,k;
    double det;
    int s;
    double m_data[(nvert-1)*(nvert-1)];
    double w_data[(nvert-1)*(nvert-1)];
/*
    double w_data[3*3] = {1.,0.,0.,
			  1.,1.,0.,
			  1.,1.,1.};*/
    //double wt_data[(nvert-1)*ndims];
    gsl_permutation * perm;

    gsl_matrix_view m = gsl_matrix_view_array (m_data, nvert-1, nvert-1);
    gsl_matrix_view w = gsl_matrix_view_array (w_data, nvert-1, ndims);
    //gsl_matrix_view wt = gsl_matrix_view_array (wt_data, nvert-1,ndims);

    if (nvert==2)
    {
	for (i=0,det=0;i<ndims;i++)
	    det+=(coord[1][i]-coord[0][i])*(coord[1][i]-coord[0][i]);
	return sqrt(det);
    }

    perm = gsl_permutation_alloc (nvert-1);

    for (i=0,k=0;i<nvert-1;i++)
	for (j=0;j<nvert-1;j++,k++)
	    w_data[k]=coord[i+1][j]-coord[0][j];

    //gsl_matrix_transpose_memcpy (wt,w);

    gsl_blas_dgemm (CblasNoTrans, CblasTrans,
		    1.0, &w.matrix, &w.matrix,
		    0.0, &m.matrix);

    gsl_linalg_LU_decomp(&m.matrix,perm,&s);
    det=gsl_linalg_LU_det(&m.matrix,s);

    gsl_permutation_free(perm);

    for (i=2,j=1;i<=nvert-1;i++) j*=i;

    return sqrt(fabs(det))/j;
}

double SimplexVolumef(float **coord,int nvert, int ndims)
{
    int i,j,k;
    double det;
    int s;
    double m_data[(nvert-1)*(nvert-1)];
    double w_data[(nvert-1)*(nvert-1)];
/*
    double w_data[3*3] = {1.,0.,0.,
			  1.,1.,0.,
			  1.,1.,1.};*/
    //double wt_data[(nvert-1)*ndims];
    gsl_permutation * perm;

    gsl_matrix_view m = gsl_matrix_view_array (m_data, nvert-1, nvert-1);
    gsl_matrix_view w = gsl_matrix_view_array (w_data, nvert-1, ndims);
    //gsl_matrix_view wt = gsl_matrix_view_array (wt_data, nvert-1,ndims);

    if (nvert==2)
    {
	for (i=0,det=0;i<ndims;i++)
	    det+=(double)(coord[1][i]-coord[0][i])*(double)(coord[1][i]-coord[0][i]);
	return sqrt(det);
    }

    perm = gsl_permutation_alloc (nvert-1);

    for (i=0,k=0;i<nvert-1;i++)
	for (j=0;j<nvert-1;j++,k++)
	    w_data[k]=coord[i+1][j]-coord[0][j];

    //gsl_matrix_transpose_memcpy (wt,w);

    gsl_blas_dgemm (CblasNoTrans, CblasTrans,
		    1.0, &w.matrix, &w.matrix,
		    0.0, &m.matrix);

    gsl_linalg_LU_decomp(&m.matrix,perm,&s);
    det=gsl_linalg_LU_det(&m.matrix,s);

    gsl_permutation_free(perm);

    for (i=2,j=1;i<=nvert-1;i++) j*=i;

    return sqrt(fabs(det))/j;
}

// center must be allocated to ndims vals at least.
// coord are the vertice coords
// center is the center of the resulting ndims-sphere
// function returns its radius squared
double SimplexSphere(double **coord, int ndims,double *center)
{
    double a_data[ndims*ndims];
    double b_data[ndims];

    gsl_matrix_view m=gsl_matrix_view_array (a_data, ndims, ndims);
    gsl_vector_view b=gsl_vector_view_array (b_data, ndims);
    gsl_vector_view x = gsl_vector_view_array (center,ndims);
    gsl_permutation * p = gsl_permutation_alloc (ndims);

    int s;
    int i,j,k;
    double d;

    for (i=0,k=0;i<ndims;i++)
    {
	b_data[i]=0;
	for (j=0;j<ndims;j++,k++)
	{
	    a_data[k]=((double)coord[0][j]-(double)coord[i+1][j]);
	    b_data[i]+=((double)coord[0][j]*(double)coord[0][j]-(double)coord[i+1][j]*(double)coord[i+1][j]);
	}
	b_data[i]*=0.5;
    }

    gsl_linalg_LU_decomp (&m.matrix, p, &s);
    gsl_linalg_LU_solve (&m.matrix, p, &b.vector, &x.vector);
    gsl_permutation_free (p);

    for (i=0,d=0;i<ndims;i++)
	d+=((double)center[i]-(double)coord[0][i])*((double)center[i]-(double)coord[0][i]);

    return d;
}

#else

double SimplexVolumef(float **coord,int nvert, int ndims)
{
    int i;
    double det;

    if (nvert==2)
    {
	for (i=0,det=0;i<ndims;i++)
	    det+=(double)(coord[1][i]-coord[0][i])*(double)(coord[1][i]-coord[0][i]);
	return sqrt(det);
    }

    fprintf(stderr,"ERROR: Function 'SimplexVolume' needs GSL.\n");
    fprintf(stderr,"Please, recompile with GSL.\n");
    return -1;
}

double SimplexSphere(double **coords, int ndims,double *center)
{
    fprintf(stderr,"ERROR: Function 'SimplexSphere' needs GSL.\n");
    fprintf(stderr,"Please, recompile with GSL.\n");

    return -1;
}

#endif