text:
/**
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF licenses this file
 * to you under the Apache License, Version 2.0 (the
 * "License"); you may not use this file except in compliance
 * with the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#ifndef SINGA_CORE_TENSOR_TENSOR_MATH_CPP_H_
#define SINGA_CORE_TENSOR_TENSOR_MATH_CPP_H_

#include "./tensor_math.h"
#include <cfloat>
#include "singa/core/common.h"
#include <math.h>

#ifdef USE_CBLAS
#include <cblas.h>
#endif

namespace singa {

template <>
void Abs<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                           Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = fabs(inPtr[i]);
  }
}

template <>
void Add<float, lang::Cpp>(const size_t num, const Block *in, const float x,
                           Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = inPtr[i] + x;
  }
}

template <>
void Add<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                           Block *out, Context *ctx) {
  // CHECK_EQ(ctx->stream, nullptr);
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = in1Ptr[i] + in2Ptr[i];
  }
}

template <>
void Clamp<float, lang::Cpp>(const size_t num, const float low,
                             const float high, const Block *in, Block *out,
                             Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    if (inPtr[i] > high) {
      outPtr[i] = high;
    } else if (inPtr[i] < low) {
      outPtr[i] = low;
    } else {
      outPtr[i] = inPtr[i];
    }
  }
}

template <>
void Div<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                           Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    CHECK_NE(in2Ptr[i], 0.f);
    outPtr[i] = in1Ptr[i] / in2Ptr[i];
  }
}

template <>
void Div<float, lang::Cpp>(const size_t num, const float x, const Block *in,
                           Block *out, Context *ctx) {
  const float *inPtr = static_cast<const float *>(in->data());
  float *outPtr = static_cast<float *>(out->mutable_data());
  for (size_t i = 0; i < num; i++) {
    CHECK_NE(inPtr[i], 0.f);
    outPtr[i] = x / inPtr[i];
  }
}

template <>
void EltwiseMult<float, lang::Cpp>(const size_t num, const Block *in,
                                   const float x, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = inPtr[i] * x;
  }
}

template <>
void EltwiseMult<float, lang::Cpp>(const size_t num, const Block *in1,
                                   const Block *in2, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = in1Ptr[i] * in2Ptr[i];
  }
}
template <>
void Exp<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                           Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = exp(inPtr[i]);
  }
}

template <>
void GE<float, lang::Cpp>(const size_t num, const Block *in, const float x,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr[i] >= x) ? 1.f : 0.f;
  }
}

template <>
void GE<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr1 = static_cast<const float *>(in1->data());
  const float *inPtr2 = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr1[i] >= inPtr2[i]) ? 1.f : 0.f;
  }
}
template <>
void GT<float, lang::Cpp>(const size_t num, const Block *in, const float x,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr[i] > x) ? 1.f : 0.f;
  }
}
template <>
void GT<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr1 = static_cast<const float *>(in1->data());
  const float *inPtr2 = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr1[i] > inPtr2[i]) ? 1.f : 0.f;
  }
}

template <>
void LE<float, lang::Cpp>(const size_t num, const Block *in, const float x,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr[i] <= x) ? 1.f : 0.f;
  }
}
template <>
void LE<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr1 = static_cast<const float *>(in1->data());
  const float *inPtr2 = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr1[i] <= inPtr2[i]) ? 1.f : 0.f;
  }
}
template <>
void Log<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                           Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    CHECK_GT(inPtr[i], 0.f);
    outPtr[i] = log(inPtr[i]);
  }
}
template <>
void LT<float, lang::Cpp>(const size_t num, const Block *in, const float x,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr[i] < x) ? 1.f : 0.f;
  }
}
template <>
void LT<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                          Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr1 = static_cast<const float *>(in1->data());
  const float *inPtr2 = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr1[i] < inPtr2[i]) ? 1.f : 0.f;
  }
}

template <>
void Pow<float, lang::Cpp>(const size_t num, const Block *in, const float x,
                           Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = pow(inPtr[i], x);
  }
}

template <>
void Pow<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                           Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = pow(in1Ptr[i], in2Ptr[i]);
  }
}
template <>
void ReLU<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                            Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr[i] >= 0.f) ? inPtr[i] : 0.f;
  }
}
template <>
void Set<float, lang::Cpp>(const size_t num, const float x, Block *out,
                           Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  for (size_t i = 0; i < num; i++) outPtr[i] = x;
}
template <>
void Set<int, lang::Cpp>(const size_t num, const int x, Block *out,
                           Context *ctx) {
  int *outPtr = static_cast<int *>(out->mutable_data());
  for (size_t i = 0; i < num; i++) outPtr[i] = x;
}

template <>
void Sigmoid<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                               Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = 1.f / (1.f + exp(-inPtr[i]));
  }
}

template <>
void Sign<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                            Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = (inPtr[i] > 0) - (inPtr[i] < 0);
  }
}

template <>
void Sqrt<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                            Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    CHECK_GE(inPtr[i], 0.f);
    outPtr[i] = sqrt(inPtr[i]);
  }
}
/*
template <>
void Square<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                              Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = inPtr[i] * inPtr[i];
  }
}
*/

template <>
void Sub<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                           Block *out, Context *ctx) {
  // CHECK_EQ(ctx->stream, nullptr);
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = in1Ptr[i] - in2Ptr[i];
  }
}

// sum all elements of input into out
// TODO(wangwei) optimize using omp
template <>
void Sum<float, lang::Cpp>(const size_t num, const Block *in, float *out,
                           Context *ctx) {
  float s = 0.f;
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    s += inPtr[i];
  }
  *out = s;
}

template <>
void Tanh<float, lang::Cpp>(const size_t num, const Block *in, Block *out,
                            Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = tanh(inPtr[i]);
  }
}

// ===============Random operations==========================================
template <>
void Bernoulli<float, lang::Cpp>(const size_t num, const float p, Block *out,
                                 Context *ctx) {
  std::bernoulli_distribution distribution(p);
  float *outPtr = static_cast<float *>(out->mutable_data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = distribution(ctx->random_generator) ? 1.0f : 0.0f;
  }
}

template <>
void Gaussian<float, lang::Cpp>(const size_t num, const float mean,
                                const float std, Block *out, Context *ctx) {
  std::normal_distribution<float> distribution(mean, std);
  float *outPtr = static_cast<float *>(out->mutable_data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = static_cast<float>(distribution(ctx->random_generator));
  }
}
template <>
void Uniform<float, lang::Cpp>(const size_t num, const float low,
                               const float high, Block *out, Context *ctx) {
  std::uniform_real_distribution<float> distribution(low, high);
  float *outPtr = static_cast<float *>(out->mutable_data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] = static_cast<float>(distribution(ctx->random_generator));
  }
}

// ====================Blas operations======================================

template <>
void DGMM<float, lang::Cpp>(const bool side_right, const size_t nrow,
                            const size_t ncol, const Block *M, const Block *v,
                            Block *out, Context *ctx) {
  const float *MPtr = static_cast<const float *>(M->data());
  const float *vPtr = static_cast<const float *>(v->data());
  float *outPtr = static_cast<float *>(out->mutable_data());
  if (side_right) {
    for (size_t r = 0; r < nrow; r++) {
      size_t offset = r * ncol;
      for (size_t c = 0; c < ncol; c++) {
        outPtr[offset + c] = MPtr[offset + c] * vPtr[c];
      }
    }
  } else {
    for (size_t r = 0; r < nrow; r++) {
      size_t offset = r * ncol;
      for (size_t c = 0; c < ncol; c++) {
        outPtr[offset + c] = MPtr[offset + c] * vPtr[r];
      }
    }
  }
}

#ifdef USE_CBLAS
template <>
void Amax<float, lang::Cpp>(const size_t num, const Block *in, size_t *out,
                            Context *ctx) {
  const float *inPtr = static_cast<const float *>(in->data());
  *out = cblas_isamax(num, inPtr, 1);
}

template <>
void Asum<float, lang::Cpp>(const size_t num, const Block *in, float *out,
                            Context *ctx) {
  const float *inPtr = static_cast<const float *>(in->data());
  *out = cblas_sasum(num, inPtr, 1);
}

template <>
void Axpy<float, lang::Cpp>(const size_t num, const float alpha,
                            const Block *in, Block *out, Context *ctx) {
  const float *inPtr = static_cast<const float *>(in->data());
  float *outPtr = static_cast<float *>(out->mutable_data());
  cblas_saxpy(num, alpha, inPtr, 1, outPtr, 1);
}

template <>
void Dot<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                           float *out, Context *ctx) {
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  *out = cblas_sdot(num, in1Ptr, 1, in2Ptr, 1);
}
template <>
void Scale<float, lang::Cpp>(const size_t num, const float x, Block *out,
                             Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  cblas_sscal(num, x, outPtr, 1);
}
template <>
void Nrm2<float, lang::Cpp>(const size_t num, const Block *in, float *out,
                            Context *ctx) {
  const float *inPtr = static_cast<const float *>(in->data());
  *out = cblas_snrm2(num, inPtr, 1);
}

template <>
void GEMV<float, lang::Cpp>(bool trans, const size_t m, const size_t n,
                            const float alpha, const Block *A, const Block *v,
                            const float beta, Block *out, Context *ctx) {
  const float *APtr = static_cast<const float *>(A->data());
  const float *vPtr = static_cast<const float *>(v->data());
  float *outPtr = static_cast<float *>(out->mutable_data());
  if (!trans) {
    cblas_sgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, APtr, n, vPtr, 1,
                beta, outPtr, 1);
  } else {
    cblas_sgemv(CblasRowMajor, CblasTrans, n, m, alpha, APtr, m, vPtr, 1, beta,
                outPtr, 1);
  }
}

template <>
void GEMM<float, lang::Cpp>(const bool transA, const bool transB,
                            const size_t nrowA, const size_t ncolB,
                            const size_t ncolA, const float alpha,
                            const Block *A, const Block *B, const float beta,
                            Block *C, Context *ctx) {
  auto transa = transA ? CblasTrans : CblasNoTrans;
  auto transb = transB ? CblasTrans : CblasNoTrans;
  auto lda = transA ? nrowA : ncolA;
  auto ldb = transB ? ncolA : ncolB;
  auto ldc = ncolB;
  const float *APtr = static_cast<const float *>(A->data());
  const float *BPtr = static_cast<const float *>(B->data());
  float *CPtr = static_cast<float *>(C->mutable_data());
  cblas_sgemm(CblasRowMajor, transa, transb, nrowA, ncolB, ncolA, alpha, APtr,
	  lda, BPtr, ldb, beta, CPtr, ldc);
}

#else

template <>
void Amax<float, lang::Cpp>(const size_t num, const Block *in, size_t *out,
                            Context *ctx) {
  size_t maxPos = 0;
  float maxVal = 0;
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    if (i == 0) {
      maxVal = inPtr[i];
    } else if (inPtr[i] > maxVal) {
      maxVal = inPtr[i];
      maxPos = i;
    }
  }
  *out = maxPos;
}
template <>
void Amin<float, lang::Cpp>(const size_t num, const Block *in, size_t *out,
                            Context *ctx) {
  size_t minPos = 0;
  float minVal = 0;
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    if (i == 0) {
      minVal = inPtr[i];
    } else if (inPtr[i] > minVal) {
      minVal = inPtr[i];
      minPos = i;
    }
  }
  *out = minPos;
}

template <>
void Asum<float, lang::Cpp>(const size_t num, const Block *in, float *out,
                            Context *ctx) {
  float sum = 0;
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    sum += fabs(inPtr[i]);
  }
}

template <>
void Axpy<float, lang::Cpp>(const size_t num, const float alpha,
                            const Block *in, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] += alpha * inPtr[i];
  }
}

template <>
void Scale<float, lang::Cpp>(const size_t num, const float x, Block *out,
                             Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  for (size_t i = 0; i < num; i++) {
    outPtr[i] *= x;
  }
}

template <>
void Dot<float, lang::Cpp>(const size_t num, const Block *in1, const Block *in2,
                           float *out, Context *ctx) {
  float sum = 0;
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  for (size_t i = 0; i < num; i++) {
    sum += in1Ptr[i] * in2Ptr[i];
  }
}

template <>
void GEMV<float, lang::Cpp>(bool trans, const size_t m, const size_t n,
                            const float alpha, const Block *A, const Block *v,
                            const float beta, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *APtr = static_cast<const float *>(A->data());
  const float *vPtr = static_cast<const float *>(v->data());
  for (size_t r = 0; r < m; r++) {
    float sum = 0;
    for (size_t c = 0; c < n; c++) {
      size_t idx = trans ? c * m + r : r * n + c;
      sum += APtr[idx] * vPtr[c];
    }
    outPtr[r] = alpha * sum + beta * outPtr[r];
  }
}

#endif  // USE_CBLAS
template <>
void ComputeCrossEntropy<float, lang::Cpp>(bool int_target,
                                           const size_t batchsize,
                                           const size_t dim, const Block *p,
                                           const Block *t, Block *loss,
                                           Context *ctx) {
  const float *pPtr = static_cast<const float *>(p->data());
  const int *tPtr = static_cast<const int *>(t->data());
  float *lossPtr = static_cast<float *>(loss->mutable_data());
  if (int_target) {
    for (size_t i = 0; i < batchsize; i++) {
      int truth_idx = tPtr[i];
      CHECK_GE(truth_idx, 0);
      float prob_of_truth = pPtr[i * dim + truth_idx];
      lossPtr[i] = -std::log((std::max)(prob_of_truth, FLT_MIN));
    }
  } else {
    for (size_t i = 0;i < batchsize; i++) {
      float sum = 0.f;
      for (size_t j = 0; j < dim; j++) {
        sum += tPtr[i * dim + j];
      }
      float loss = 0.f;
      for (size_t j = 0, offset = i * dim; j < dim; j++, offset++) {
        loss -= tPtr[offset] / sum * std::log((std::max)(pPtr[offset], FLT_MIN));
      }
      lossPtr[i] = loss;
    }
  }
}

template <>
void SoftmaxCrossEntropyBwd<float, lang::Cpp>(bool int_target,
                                              const size_t batchsize,
                                              const size_t dim, const Block *p,
                                              const Block *t, Block *grad,
                                              Context *ctx) {
  CHECK_EQ(p, grad) << "Use the same pointer to optimize performance";
  // const float* pPtr = static_cast<const float*>(p->data());
  const int *tPtr = static_cast<const int *>(t->data());
  float *gradPtr = static_cast<float *>(grad->mutable_data());

  if (int_target) {
    for (size_t i = 0; i < batchsize; i++) {
      int truth_idx = static_cast<int>(tPtr[i]);
      CHECK_GE(truth_idx, 0);
      gradPtr[i * dim + truth_idx] -= 1.0;
    }
  } else {
    for (size_t i = 0; i < batchsize; i++) {
      float sum = 0.f;
      for (size_t j = 0; j < dim; j++) {
        sum += tPtr[i * dim + j];
      }
      for (size_t j = 0, offset = i * dim; j < dim; j++, offset++) {
        gradPtr[offset] -= tPtr[offset] / sum;
      }
    }
  }
}

template <>
void RowMax<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                              const Block *in, Block *out, Context *ctx) {
  const float *inPtr = static_cast<const float *>(in->data());
  float *outPtr = static_cast<float *>(out->mutable_data());
  for (size_t r = 0; r < nrow; r++) {
    int offset = (int)(r * ncol);
    float maxval = inPtr[offset];
    for (size_t c = 1; c < ncol; c++)
      maxval = (std::max)(maxval, inPtr[offset + c]);
    outPtr[r] = maxval;
  }
}

// =========Matrix operations ================================================
/*
template <>
void AddCol<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                              const Block *A, const Block *v, Block *out,
                              Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *APtr = static_cast<const float *>(A->data());
  const float *vPtr = static_cast<const float *>(v->data());
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[offset + c] = APtr[offset + c] + vPtr[r];
    }
  }
}

template <>
void AddRow<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                              const Block *A, const Block *v, Block *out,
                              Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *APtr = static_cast<const float *>(A->data());
  const float *vPtr = static_cast<const float *>(v->data());
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[offset + c] = APtr[offset + c] + vPtr[c];
    }
  }
}
template <>
void Outer<float, lang::Cpp>(const size_t m, const size_t n, const Block *in1,
                             const Block *in2, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *in1Ptr = static_cast<const float *>(in1->data());
  const float *in2Ptr = static_cast<const float *>(in2->data());
  for (size_t r = 0; r < m; r++) {
    size_t offset = r * n;
    for (size_t c = 0; c < n; c++) {
      outPtr[offset + c] = in1Ptr[r] * in2Ptr[c];
    }
  }
}
template <>
void Softmax<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                               const Block *in, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  float *bPtr = new float[ncol];
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    float denom = 0.f;
    for (size_t c = 0; c < ncol; c++) {
      bPtr[c] = exp(inPtr[offset + c]);
      denom += bPtr[c];
    }
    for (size_t c = 0; c < ncol; c++) {
      size_t idx = offset + c;
      outPtr[idx] = bPtr[c] / denom;
    }
  }
  delete bPtr;
}

template <>
void SumColumns<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                                  const Block *in, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t c = 0; c < ncol; c++) {
    outPtr[c] = 0.f;
  }
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[c] += inPtr[offset + c];
    }
  }
}

template <>
void SumRows<float, lang::Cpp>(const size_t nrow, const size_t ncol,
                               const Block *in, Block *out, Context *ctx) {
  float *outPtr = static_cast<float *>(out->mutable_data());
  const float *inPtr = static_cast<const float *>(in->data());
  for (size_t r = 0; r < nrow; r++) {
    size_t offset = r * ncol;
    outPtr[r] = 0.f;
    for (size_t c = 0; c < ncol; c++) {
      outPtr[r] += inPtr[offset + c];
    }
  }
}
*/
}  // namespace singa

#endif  // SINGA_CORE_TENSOR_TENSOR_MATH_CPP_H_

max_issues_repo_name: TsinghuaDatabaseGroup/incubator-singa

max_forks_repo_name: TsinghuaDatabaseGroup/incubator-singa


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <sys/stat.h>
#include <sys/types.h>

#include "allvars.h"
#include "proto.h"

#ifdef OUTPUTLINEOFSIGHT

#ifdef OUTPUTLINEOFSIGHT_SPECTRUM
#define  PIXELS 512		/* number of bins along line of sight */
#else
#define  PIXELS 1
#endif

#define  N_LOS  10		/* number of lines of sight selected  */

static double H_a, Wmax;

struct line_of_sight
{
  int xaxis, yaxis, zaxis;
  double Xpos, Ypos;		/* relative position of line-of-sight on face of box */
  double BoxSize, Wmax, Time;

  /* total gas density */
  double Rho[PIXELS];
  double Vpec[PIXELS];
  double Temp[PIXELS];
  double Metallicity[PIXELS];

  /* neutral hydrogen */
  double RhoHI[PIXELS];
  double NHI[PIXELS];
  double VpecHI[PIXELS];
  double TempHI[PIXELS];
  double TauHI[PIXELS];

  /* HeII quantities */
  double RhoHeII[PIXELS];
  double NHeII[PIXELS];
  double VpecHeII[PIXELS];
  double TempHeII[PIXELS];
  double TauHeII[PIXELS];
}
 *Los, *LosGlobal;

struct line_of_sight_particles
{
  MyFloat Pos[3];
  MyFloat Hsml;
  MyFloat Vz;
  MyFloat Utherm;
  MyFloat Mass;
  MyFloat Metallicity;
}
 *particles;

void lineofsight_output(void)
{
  char buf[500];
  int n, s, next;
  double ti;

  next = find_next_lineofsighttime(All.Ti_nextlineofsight);

  ti = All.TimeBegin * exp(next * All.Timebase_interval);

  if(ThisTask == 0)
    {
      printf("Line of sight output! ThisTask=%d Time=%g  NextTime=%g\n", ThisTask, All.Time, ti);
      fflush(stdout);
    }

  H_a = hubble_function(All.Time);
  Wmax = All.Time * H_a * All.BoxSize;

  if(ThisTask == 0)
    {
      sprintf(buf, "%s/los", All.OutputDir);
      mkdir(buf, 02755);
    }

  Los = mymalloc(sizeof(struct line_of_sight));
  LosGlobal = mymalloc(sizeof(struct line_of_sight));

  for(n = 0, s = 0; n < N_LOS; n++)
    {
      if(s + 3 >= RNDTABLE)
	{
	  set_random_numbers();
	  s = 0;
	}

      Los->zaxis = (int) (3.0 * get_random_number(s++));
      switch (Los->zaxis)
	{
	case 2:
	  Los->xaxis = 0;
	  Los->yaxis = 1;
	  break;
	case 0:
	  Los->xaxis = 1;
	  Los->yaxis = 2;
	  break;
	case 1:
	  Los->xaxis = 2;
	  Los->yaxis = 0;
	  break;
	}

      Los->Xpos = All.BoxSize * get_random_number(s++);
      Los->Ypos = All.BoxSize * get_random_number(s++);

#ifdef OUTPUTLINEOFSIGHT_SPECTRUM
      add_along_lines_of_sight();
      sum_over_processors_and_normalize();
      absorb_along_lines_of_sight();
      output_lines_of_sight(n);
#endif

#ifdef OUTPUTLINEOFSIGHT_PARTICLES
      find_particles_and_save_them(n);
#endif
    }

  myfree(LosGlobal);
  myfree(Los);
}

void find_particles_and_save_them(int num)
{
  int n, k, count_local, *countlist, counttot, rep;
  double a3inv;
  MyFloat dx, dy, r2;
  char fname[1000];
  MPI_Status status;
  FILE *fd = 0;

  countlist = mymalloc(sizeof(int) * NTask);
  particles = mymalloc(sizeof(struct line_of_sight_particles) * N_gas);

  a3inv = 1.0 / (All.Time * All.Time * All.Time);

  for(n = 0, count_local = 0; n < N_gas; n++)
    {
      if(P[n].Type == 0)
	{
	  dx = los_periodic(P[n].Pos[Los->xaxis] - Los->Xpos);
	  dy = los_periodic(P[n].Pos[Los->yaxis] - Los->Ypos);

	  r2 = dx * dx + dy * dy;

	  if(r2 < PPP[n].Hsml * PPP[n].Hsml)
	    {
	      for(k = 0; k < 3; k++)
		particles[count_local].Pos[k] = P[n].Pos[k];

	      particles[count_local].Hsml = PPP[n].Hsml;
	      particles[count_local].Vz = P[n].Vel[Los->zaxis];
	      particles[count_local].Utherm = SphP[n].Entropy / GAMMA_MINUS1 * pow(SphP[n].d.Density *
										   a3inv, GAMMA_MINUS1);
	      particles[count_local].Mass = P[n].Mass;
	      particles[count_local].Metallicity = P[n].Metallicity;

	      count_local++;
	    }
	}
    }

  MPI_Gather(&count_local, 1, MPI_INT, countlist, 1, MPI_INT, 0, MPI_COMM_WORLD);

  if(ThisTask == 0)
    {
      sprintf(fname, "%s/los/part_los_z%05.3f_%03d.dat", All.OutputDir, 1 / All.Time - 1, num);

      if(!(fd = fopen(fname, "w")))
	{
	  printf("can't open file `%s`\n", fname);
	  endrun(112);
	}

      fwrite(&count_local, sizeof(int), 1, fd);	/* will be overwritten later */
      fwrite(&LosGlobal->xaxis, sizeof(int), 1, fd);
      fwrite(&LosGlobal->yaxis, sizeof(int), 1, fd);
      fwrite(&LosGlobal->zaxis, sizeof(int), 1, fd);
      fwrite(&LosGlobal->Xpos, sizeof(double), 1, fd);
      fwrite(&LosGlobal->Ypos, sizeof(double), 1, fd);
      fwrite(&LosGlobal->BoxSize, sizeof(double), 1, fd);
      fwrite(&LosGlobal->Wmax, sizeof(double), 1, fd);
      fwrite(&LosGlobal->Time, sizeof(double), 1, fd);
    }

  for(rep = 0, counttot = 0; rep < NTask; rep++)
    {
      if(ThisTask != 0 && rep == ThisTask && count_local > 0)
	MPI_Ssend(particles, sizeof(struct line_of_sight_particles) * count_local, MPI_BYTE, 0,
		  TAG_PDATA, MPI_COMM_WORLD);

      if(ThisTask == 0)
	{
	  if(rep > 0 && countlist[rep] > 0)
	    MPI_Recv(particles, sizeof(struct line_of_sight_particles) * countlist[rep],
		     MPI_BYTE, rep, TAG_PDATA, MPI_COMM_WORLD, &status);

	  fwrite(particles, sizeof(struct line_of_sight_particles), countlist[rep], fd);

	  counttot += countlist[rep];
	}
    }

  if(ThisTask == 0)
    {
      fclose(fd);
      if(!(fd = fopen(fname, "r+")))
	{
	  printf("can't open file `%s'\n", fname);
	  endrun(113);
	}

      fseek(fd, 0, SEEK_CUR);
      fwrite(&counttot, sizeof(int), 1, fd);
      fclose(fd);
    }

  myfree(particles);
  myfree(countlist);
}

void add_along_lines_of_sight(void)
{
  int n, bin, i, iz0, iz1, iz;
  double dx, dy, dz, r, r2, ne, nh0, nHeII, utherm, temp, meanWeight;
  double u, wk, weight, a3inv, h3inv;
  double z0, z1, dmax1, dmax2;

  for(i = 0; i < PIXELS; i++)
    {
      Los->Rho[i] = 0;
      Los->Vpec[i] = 0;
      Los->Temp[i] = 0;
      Los->Metallicity[i] = 0;

      Los->RhoHI[i] = 0;
      Los->NHI[i] = 0;
      Los->VpecHI[i] = 0;
      Los->TempHI[i] = 0;
      Los->TauHI[i] = 0;

      Los->RhoHeII[i] = 0;
      Los->NHeII[i] = 0;
      Los->VpecHeII[i] = 0;
      Los->TempHeII[i] = 0;
      Los->TauHeII[i] = 0;
    }

  a3inv = 1.0 / (All.Time * All.Time * All.Time);

  for(n = 0; n < N_gas; n++)
    {
      if(P[n].Type == 0)
	{
	  dx = los_periodic(P[n].Pos[Los->xaxis] - Los->Xpos);
	  dy = los_periodic(P[n].Pos[Los->yaxis] - Los->Ypos);

	  r2 = dx * dx + dy * dy;

	  if(r2 < PPP[n].Hsml * PPP[n].Hsml)
	    {
	      z0 = (P[n].Pos[Los->zaxis] - PPP[n].Hsml) / All.BoxSize * PIXELS;
	      z1 = (P[n].Pos[Los->zaxis] + PPP[n].Hsml) / All.BoxSize * PIXELS;
	      iz0 = (int) z0;
	      iz1 = (int) z1;
	      if(z0 < 0)
		iz0 -= 1;

	      for(iz = iz0; iz <= iz1; iz++)
		{
		  dz = los_periodic((iz + 0.5) / PIXELS * All.BoxSize - P[n].Pos[Los->zaxis]);
		  r = sqrt(r2 + dz * dz);

		  if(PPP[n].Hsml > All.BoxSize)
		    {
		      printf("Here:%d  n=%d %g\n", ThisTask, n, PPP[n].Hsml);
		      endrun(89);
		    }

		  if(r < PPP[n].Hsml)
		    {
		      u = r / PPP[n].Hsml;
		      h3inv = 1.0 / (PPP[n].Hsml * PPP[n].Hsml * PPP[n].Hsml);

		      if(u < 0.5)
			wk = h3inv * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		      else
			wk = h3inv * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);

		      bin = iz;
		      while(bin >= PIXELS)
			bin -= PIXELS;
		      while(bin < 0)
			bin += PIXELS;

		      ne = SphP[n].Ne;
		      utherm = DMAX(All.MinEgySpec,
				    SphP[n].Entropy / GAMMA_MINUS1 * pow(SphP[n].d.Density *
									 a3inv, GAMMA_MINUS1));

		      AbundanceRatios(utherm, SphP[n].d.Density * a3inv, &ne, &nh0, &nHeII);

		      meanWeight = 4.0 / (3 * HYDROGEN_MASSFRAC + 1 + 4 * HYDROGEN_MASSFRAC * ne);

		      temp = meanWeight * PROTONMASS / BOLTZMANN * GAMMA_MINUS1 * utherm
			* All.UnitEnergy_in_cgs / All.UnitMass_in_g;

		      /* do total gas */
		      weight = P[n].Mass * wk;
		      Los->Rho[bin] += weight;
		      Los->Metallicity[bin] += P[n].Metallicity * weight;
		      Los->Temp[bin] += temp * weight;
		      Los->Vpec[bin] += P[n].Vel[Los->zaxis] * weight;

		      /* do neutral hydrogen */
		      weight = nh0 * HYDROGEN_MASSFRAC * P[n].Mass * wk;
		      Los->RhoHI[bin] += weight;
		      Los->TempHI[bin] += temp * weight;
		      Los->VpecHI[bin] += P[n].Vel[Los->zaxis] * weight;

		      /* do HeII */
		      weight = 4 * nHeII * HYDROGEN_MASSFRAC * P[n].Mass * wk;
		      Los->RhoHeII[bin] += weight;
		      Los->TempHeII[bin] += temp * weight;
		      Los->VpecHeII[bin] += P[n].Vel[Los->zaxis] * weight;
		    }
		}
	    }
	}
    }
}

void sum_over_processors_and_normalize(void)
{
  int bin;

  MPI_Reduce(Los->Rho, LosGlobal->Rho, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(Los->Metallicity, LosGlobal->Metallicity, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(Los->Temp, LosGlobal->Temp, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(Los->Vpec, LosGlobal->Vpec, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

  MPI_Reduce(Los->RhoHI, LosGlobal->RhoHI, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(Los->TempHI, LosGlobal->TempHI, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(Los->VpecHI, LosGlobal->VpecHI, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

  MPI_Reduce(Los->RhoHeII, LosGlobal->RhoHeII, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(Los->TempHeII, LosGlobal->TempHeII, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(Los->VpecHeII, LosGlobal->VpecHeII, PIXELS, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

  if(ThisTask == 0)
    {
      /* normalize results by the weights */
      for(bin = 0; bin < PIXELS; bin++)
	{
	  /* total gas density */
	  LosGlobal->Metallicity[bin] /= LosGlobal->Rho[bin];
	  LosGlobal->Temp[bin] /= LosGlobal->Rho[bin];
	  LosGlobal->Vpec[bin] /= (All.Time * LosGlobal->Rho[bin]);

	  /* neutral hydrogen quantities */
	  LosGlobal->VpecHI[bin] /= LosGlobal->RhoHI[bin];
	  LosGlobal->TempHI[bin] /= LosGlobal->RhoHI[bin];
	  LosGlobal->NHI[bin] = LosGlobal->RhoHI[bin] * (All.UnitMass_in_g / PROTONMASS);

	  /* HeII quantities */
	  LosGlobal->VpecHeII[bin] /= (All.Time * LosGlobal->RhoHeII[bin]);
	  LosGlobal->TempHeII[bin] /= LosGlobal->RhoHeII[bin];
	  LosGlobal->NHeII[bin] = LosGlobal->RhoHeII[bin] * (All.UnitMass_in_g / (4 * PROTONMASS));
	}
    }
}

void absorb_along_lines_of_sight(void)
{
  double dz, dv, b, fac, fac_HeII;
  int bin, k;

  if(ThisTask == 0)
    {
      dz = All.BoxSize / PIXELS;

      for(bin = 0; bin < PIXELS; bin++)
	{
	  LosGlobal->TauHI[bin] = 0;
	  LosGlobal->TauHeII[bin] = 0;

	  for(k = 0; k < PIXELS; k++)
	    {
	      dv = (k - bin);

	      while(dv < -PIXELS / 2)
		dv += PIXELS;
	      while(dv > PIXELS / 2)
		dv -= PIXELS;

	      dv = (dv * Wmax / PIXELS + LosGlobal->VpecHI[k]) * All.UnitVelocity_in_cm_per_s;

	      b = sqrt(2 * BOLTZMANN * LosGlobal->TempHI[k] / PROTONMASS);

	      LosGlobal->TauHI[bin] += LosGlobal->NHI[k] * exp(-dv * dv / (b * b)) / b * dz;

	      /* now HeII */
	      dv = (k - bin);

	      while(dv < -PIXELS / 2)
		dv += PIXELS;
	      while(dv > PIXELS / 2)
		dv -= PIXELS;

	      dv = (dv * Wmax / PIXELS + LosGlobal->VpecHeII[k]) * All.UnitVelocity_in_cm_per_s;

	      b = sqrt(2 * BOLTZMANN * LosGlobal->TempHeII[k] / (4 * PROTONMASS));

	      LosGlobal->TauHeII[bin] += LosGlobal->NHeII[k] * exp(-dv * dv / (b * b)) / b * dz;
	    }
	}

      /* multiply with correct prefactors */

      /*  to get things into cgs units */
      fac = 1 / pow(All.UnitLength_in_cm, 2);

      fac *= All.HubbleParam * All.HubbleParam;

      fac *= OSCILLATOR_STRENGTH * M_PI * LYMAN_ALPHA * sqrt(3 * THOMPSON / (8 * M_PI));	/* Ly-alpha cross section */

      fac *= C / (All.Time * All.Time) / sqrt(M_PI);

      /* Note: For HeII, the oscillator strength is equal to that of HI,
         and the Lyman-alpha wavelength is 4 times shorter */

      fac_HeII = fac * (OSCILLATOR_STRENGTH_HeII / OSCILLATOR_STRENGTH) * (LYMAN_ALPHA_HeII / LYMAN_ALPHA);

      for(bin = 0; bin < PIXELS; bin++)
	{
	  LosGlobal->TauHI[bin] *= fac;
	  LosGlobal->TauHeII[bin] *= fac_HeII;
	}

      LosGlobal->BoxSize = All.BoxSize;
      LosGlobal->Wmax = Wmax;
      LosGlobal->Time = All.Time;
    }

}

void output_lines_of_sight(int num)
{
  FILE *fd;
  int dummy;
  char fname[400];

  sprintf(fname, "%s/los/spec_los_z%05.3f_%03d.dat", All.OutputDir, 1 / All.Time - 1, num);

  if(!(fd = fopen(fname, "w")))
    {
      printf("can't open file `%s`\n", fname);
      exit(1);
    }

  dummy = PIXELS;
  fwrite(&dummy, sizeof(int), 1, fd);
  fwrite(&LosGlobal->BoxSize, sizeof(double), 1, fd);
  fwrite(&LosGlobal->Wmax, sizeof(double), 1, fd);
  fwrite(&LosGlobal->Time, sizeof(double), 1, fd);
  fwrite(&LosGlobal->Xpos, sizeof(double), 1, fd);
  fwrite(&LosGlobal->Ypos, sizeof(double), 1, fd);
  fwrite(&LosGlobal->xaxis, sizeof(int), 1, fd);
  fwrite(&LosGlobal->yaxis, sizeof(int), 1, fd);
  fwrite(&LosGlobal->zaxis, sizeof(int), 1, fd);

  fwrite(LosGlobal->TauHI, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->TempHI, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->VpecHI, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->NHI, sizeof(double), PIXELS, fd);

  fwrite(LosGlobal->TauHeII, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->TempHeII, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->VpecHeII, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->NHeII, sizeof(double), PIXELS, fd);

  fwrite(LosGlobal->Rho, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->Vpec, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->Temp, sizeof(double), PIXELS, fd);
  fwrite(LosGlobal->Metallicity, sizeof(double), PIXELS, fd);

  fclose(fd);
}

int find_next_lineofsighttime(int time0)
{
  double a1, a2, df1, df2, u1, u2;
  double logTimeBegin, logTimeMax;
  int i1, i2, im, time1;

  logTimeBegin = log(All.TimeBegin);
  logTimeMax = log(All.TimeMax);

  a1 = logTimeBegin + time0 * All.Timebase_interval;

  u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  i1 = (int) u1;
  if(i1 >= DRIFT_TABLE_LENGTH)
    i1 = DRIFT_TABLE_LENGTH - 1;

  if(i1 <= 1)
    df1 = u1 * GravKickTable[0];
  else
    df1 = GravKickTable[i1 - 1] + (GravKickTable[i1] - GravKickTable[i1 - 1]) * (u1 - i1);

  df2 = df1 + All.BoxSize / (C / All.UnitVelocity_in_cm_per_s);

  i2 = DRIFT_TABLE_LENGTH - 1;

  while(i2 - i1 > 0)
    {
      im = (i1 - 1 + i2) / 2;

      if(GravKickTable[im] > df2)
	i2 = im;
      else
	i1 = im + 1;
    }

  u2 = (df2 - GravKickTable[i2 - 1]) / (GravKickTable[i2] - GravKickTable[i2 - 1]) + i2;

  a2 = u2 * (logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH + logTimeBegin;

  time1 = (int) ((a2 - logTimeBegin) / All.Timebase_interval + 0.5);

  return time1;
}

double los_periodic(double x)
{
  if(x >= 0.5 * All.BoxSize)
    x -= All.BoxSize;
  if(x < -0.5 * All.BoxSize)
    x += All.BoxSize;

  return x;
}

#endif

ext: c

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/* fracbm-fpt-mc (2019)
 * Authors: Benjamin Walter, Kay Wiese
 * Header File */

// LIBRARIES
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <unistd.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf.h>
#include <fftw3.h>
#include <lapacke.h>
#include <cblas.h>

// MACROS
#define IJ2K(a,b) (a+b*(b+1)/2) // Converts matrix indices
#define ARRAY_REALLOC_FACTOR 2.0 // Factor for realloc
#define MIN(a,b) ( (a < b) ? (a) : (b))
#define MAX(a,b) ( (a > b) ? (a) : (b))
#define ABS(a) ((a > 0) ? (a): (-a) )
#define ALLOC(p,n)  (p)=malloc( (n) * sizeof(*(p))); if( (p) == NULL){printf("Allocation of '%s' failed. Terminate. \n", #p); exit(2); }
#define FFT_ALLOC(p,n)  (p)=fftw_malloc( (n) * sizeof(*(p))); if( (p) == NULL){printf("Allocation of '%s' failed. Terminate. \n", #p); exit(2); }
#define REALLOC(p,n)  (p)=realloc( (p) , (n) * sizeof(*(p))); if( (p) == NULL){printf("Allocation of '%s' failed. Terminate. \n", #p); exit(2); }

// STRUCT
// Complex numbers
typedef struct complex_z
{
	 double r;
	 double i;
} complex_z;

typedef struct bridge_process // This is the bridge process spanned by two midpoints. In its center, a midpoint will be generated
{
	double right_time;	// Time of right endpoint
	double right_value;	// Height of right endpoint
	double left_time;	// Time of left endpoint
	double left_value;	// Height of left endpoint
	int generation;		// Generation
	double threshold;	// They know where the upper threshold is
	double critical_strip;	// And the critical strip;
	int centre_critical;	// 1 - Center is in critical strip -> Sub-divide; 0 - Off the critical strip -> Stop
	int crossing_threshold; // 1 - The connection of both endpoints crosses the threshold
	struct bridge_process * left_sub_bridge;
	struct bridge_process * right_sub_bridge;
	struct bridge_process * parental_bridge;
	struct bridge_process * root_bridge;
	struct bridge_process * previous_root;
	struct bridge_process ** pointer_stack; // Only root node has this. List of all bridges attached to the root. Free every pointer in this list.
	int stack_top;
} bridge_process;

typedef struct triag_matrix
{
	/* This struct stores a sequence of points X_1, X_2, ... and their inverse correlation matrix. It can be dynamically managed as points are added. Convention for triagonal matrix: column major form and 'upper' triagonal form.*/
	long size; // This is the size of the matrix itself, that is 0, 1, ..., size - 1 are array indices for both vectors 'trajectory_x' and 'trajectory_t', and 0, 1, ..., size*(size+1)/2 - 1 are indices for upper triagonal matrix. Note that size=N, meanwhile fracbm has N+1 entries. So size is the number of 'free' points, the first one being fixed.
	long array_length; // This is the length of the array. If size gets to large, realloc. If size == array_length, this array is full.
	double * inv_corr_matrix; // ( size * ( size + 1) / 2) entries in symmetric inverse correlation matrix of all points already known
	double * trajectory_x; // All points of the trajectory already known. Length = size + 1 (X_0 = 0 doesn't count, and is neglected in inverse correlation matrix (null mode)).
	double * trajectory_t; // And the corresponding time points. Length = size + 1
	double hurst_parameter;
} triag_matrix;

// FUNCTIONS
void initialise( fftw_complex** ,  fftw_complex** , fftw_complex** ,  fftw_complex** ,  double** , complex_z** , long N, gsl_rng**, const gsl_rng_type**, int);
void initialise_trajectory ( double**, long);
void initialise_inverse_correlation_matrix(double***, long );
double erfcinv(double);
void write_correlation_exponents(double*, long, double, double);
void write_correlation(fftw_complex*, double*, long);
void write_inverse_correlation_matrix(double **, long, double);
void generate_random_vector(complex_z*, fftw_complex*, fftw_complex*,long, double);
void set_to_zero(double*, long);
void integrate_noise(double*, fftw_complex*, double, double, long, double, int*, double);
void find_fpt(double*, double*, double, long, double, double**, double, int);
void fpt_to_zvar(double, double, double);
void initialise_critical_bridge(bridge_process**, double, double, double, double, double, double, bridge_process*);
void split_and_search_bridge(bridge_process*, int*, double*,  double);

bridge_process* check_this_bridge(bridge_process*, double*, double);
bridge_process* split_bridge(bridge_process*);
double generate_random_conditional_midpoint(double, double, double, double);
double crossing_time_of_bridge(bridge_process*);
double time_time_correlation(double, double, double);

void enlarge_QI(void);
void print_QI(void);
void print_bridge(bridge_process*);
void copy_QI(double**, int, double, double*, double );
void free_tree(bridge_process**);
void free_bridge(bridge_process**);

// GLOBAL VARIABLES
extern int max_generation;
extern double *gamma_N_vec;
extern double *g_vec;
extern double lin_drift, frac_drift; // Additional linear and fractional drift constants: Z_t = X_t + lin_drift * t + frac_drift * t^(2*hurst). FPT is searched for Z_t.
extern double *xfracbm; // The fBM trajectory with *no* drift is 'xfracbm'. The process with drift is labelled 'fracbm'.

// GSL RNG
extern const gsl_rng_type *T;
extern gsl_rng *r;
extern int seed;
extern triag_matrix *QI; //malloc somewhere

ext: h

max_stars_repo_name: benjamin-w/fracbm-fpt-mc

max_issues_repo_name: benjamin-w/fracbm-fpt-mc

max_forks_repo_name: benjamin-w/fracbm-fpt-mc


text:
/*
 * SpanDSP - a series of DSP components for telephony
 *
 * makecss.c - Create the composite source signal (CSS) for G.168 testing.
 *
 * Written by Steve Underwood <steveu@coppice.org>
 *
 * Copyright (C) 2003 Steve Underwood
 *
 * All rights reserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License version 2, as
 * published by the Free Software Foundation.
 *
 * 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.
 */

/*! \page makecss_page CSS construction for G.168 testing
\section makecss_page_sec_1 What does it do?
???.

\section makecss_page_sec_2 How does it work?
???.
*/

#if defined(HAVE_CONFIG_H)
#include "config.h"
#endif

#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <time.h>
#include <stdio.h>
#include <fcntl.h>
#include <sndfile.h>
#if defined(HAVE_FFTW3_H)
#include <fftw3.h>
#else
#include <fftw.h>
#endif

#include "spandsp.h"
#include "spandsp/g168models.h"

#if !defined(NULL)
#define NULL (void *) 0
#endif

#define FAST_SAMPLE_RATE    44100.0

#define C1_VOICED_SAMPLES   2144    /* 48.62ms at 44100 samples/second => 2144.142 */
#define C1_NOISE_SAMPLES    8820    /* 200ms at 44100 samples/second => 8820.0 */
#define C1_SILENCE_SAMPLES  4471    /* 101.38ms at 44100 samples/second => 4470.858 */

#define C3_VOICED_SAMPLES   3206    /* 72.69ms at 44100 samples/second => 3205.629 */
#define C3_NOISE_SAMPLES    8820    /* 200ms at 44100 samples/second => 8820.0 */
#define C3_SILENCE_SAMPLES  5614    /* 127.31ms at 44100 samples/second => 5614.371 */

static double scaling(double f, double start, double end, double start_gain, double end_gain)
{
    double scale;

    scale = start_gain + (f - start)*(end_gain - start_gain)/(end - start);
    return scale;
}
/*- End of function --------------------------------------------------------*/

static double peak(const int16_t amp[], int len)
{
    int16_t peak;
    int i;

    peak = 0;
    for (i = 0;  i < len;  i++)
    {
        if (abs(amp[i]) > peak)
            peak = abs(amp[i]);
    }
    return peak/32767.0;
}
/*- End of function --------------------------------------------------------*/

static double rms(const int16_t amp[], int len)
{
    double ms;
    int i;

    ms = 0.0;
    for (i = 0;  i < len;  i++)
        ms += amp[i]*amp[i];
    return sqrt(ms/len)/32767.0;
}
/*- End of function --------------------------------------------------------*/

static double rms_to_dbm0(double rms)
{
    return 20.0*log10(rms) + DBM0_MAX_POWER;
}
/*- End of function --------------------------------------------------------*/

static double rms_to_db(double rms)
{
    return 20.0*log10(rms);
}
/*- End of function --------------------------------------------------------*/

int main(int argc, char *argv[])
{
#if defined(HAVE_FFTW3_H)
    double in[8192][2];
    double out[8192][2];
#else
    fftw_complex in[8192];
    fftw_complex out[8192];
#endif
    fftw_plan p;
    int16_t voiced_sound[8192];
    int16_t noise_sound[8830];
    int16_t silence_sound[8192];
    int i;
    int voiced_length;
    int randy;
    double f;
    double pk;
    double ms;
    double scale;
    SNDFILE *filehandle;
    SF_INFO info;
    awgn_state_t *noise_source;

    memset(&info, 0, sizeof(info));
    info.frames = 0;
    info.samplerate = FAST_SAMPLE_RATE;
    info.channels = 1;
    info.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
    info.sections = 1;
    info.seekable = 1;
    if ((filehandle = sf_open("sound_c1.wav", SFM_WRITE, &info)) == NULL)
    {
        fprintf(stderr, "    Failed to open result file\n");
        exit(2);
    }

    printf("Generate C1\n");
    /* The sequence is
            48.62ms  of voiced sound from table C.1 of G.168
            200.0ms  of pseudo-noise
            101.38ms of silence
       The above is then repeated phase inverted.

       The voice comes straight from C.1, repeated enough times to
       fill out the 48.62ms period - i.e. 16 copies of the sequence.

       The pseudo noise section is random numbers filtered by the spectral
       pattern in Figure C.2 */

    /* The set of C1 voice samples is ready for use in the output file. */
    voiced_length = sizeof(css_c1)/sizeof(css_c1[0]);
    for (i = 0;  i < voiced_length;  i++)
        voiced_sound[i] = css_c1[i];
    pk = peak(voiced_sound, voiced_length);
    ms = rms(voiced_sound, voiced_length);
    printf("Voiced level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));

#if defined(HAVE_FFTW3_H)
    p = fftw_plan_dft_1d(8192, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
#else
    p = fftw_create_plan(8192, FFTW_BACKWARD, FFTW_ESTIMATE);
#endif
    for (i = 0;  i < 8192;  i++)
    {
#if defined(HAVE_FFTW3_H)
        in[i][0] = 0.0;
        in[i][1] = 0.0;
#else
        in[i].re = 0.0;
        in[i].im = 0.0;
#endif
    }
    for (i = 1;  i <= 3715;  i++)
    {
        f = FAST_SAMPLE_RATE*i/8192.0;

#if 1
        if (f < 50.0)
            scale = -60.0;
        else if (f < 100.0)
            scale = scaling(f, 50.0, 100.0, -25.8, -12.8);
        else if (f < 200.0)
            scale = scaling(f, 100.0, 200.0, -12.8, 17.4);
        else if (f < 215.0)
            scale = scaling(f, 200.0, 215.0, 17.4, 17.8);
        else if (f < 500.0)
            scale = scaling(f, 215.0, 500.0, 17.8, 12.2);
        else if (f < 1000.0)
            scale = scaling(f, 500.0, 1000.0, 12.2, 7.2);
        else if (f < 2850.0)
            scale = scaling(f, 1000.0, 2850.0, 7.2, 0.0);
        else if (f < 3600.0)
            scale = scaling(f, 2850.0, 3600.0, 0.0, -2.0);
        else if (f < 3660.0)
            scale = scaling(f, 3600.0, 3660.0, -2.0, -20.0);
        else if (f < 3680.0)
            scale = scaling(f, 3600.0, 3680.0, -20.0, -30.0);
        else
            scale = -60.0;
#else
        scale = 0.0;
#endif
        randy = ((rand() >> 10) & 0x1)  ?  1.0  :  -1.0;
#if defined(HAVE_FFTW3_H)
        in[i][0] = randy*pow(10.0, scale/20.0)*35.0;
        in[8192 - i][0] = -in[i][0];
#else
        in[i].re = randy*pow(10.0, scale/20.0)*35.0;
        in[8192 - i].re = -in[i].re;
#endif
    }
#if defined(HAVE_FFTW3_H)
    fftw_execute(p);
#else
    fftw_one(p, in, out);
#endif
    for (i = 0;  i < 8192;  i++)
    {
#if defined(HAVE_FFTW3_H)
        noise_sound[i] = out[i][1];
#else
        noise_sound[i] = out[i].im;
#endif
    }
    pk = peak(noise_sound, 8192);
    ms = rms(noise_sound, 8192);
    printf("Filtered noise level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));

    for (i = 0;  i < 8192;  i++)
        silence_sound[i] = 0.0;

    for (i = 0;  i < 16;  i++)
        sf_writef_short(filehandle, voiced_sound, voiced_length);
    printf("%d samples of voice\n", 16*voiced_length);
    sf_writef_short(filehandle, noise_sound, 8192);
    sf_writef_short(filehandle, noise_sound, C1_NOISE_SAMPLES - 8192);
    printf("%d samples of noise\n", C1_NOISE_SAMPLES);
    sf_writef_short(filehandle, silence_sound, C1_SILENCE_SAMPLES);
    printf("%d samples of silence\n", C1_SILENCE_SAMPLES);

    /* Now phase invert the C1 set of samples. */
    voiced_length = sizeof(css_c1)/sizeof(css_c1[0]);
    for (i = 0;  i < voiced_length;  i++)
        voiced_sound[i] = -css_c1[i];
    for (i = 0;  i < 8192;  i++)
        noise_sound[i] = -noise_sound[i];

    for (i = 0;  i < 16;  i++)
        sf_writef_short(filehandle, voiced_sound, voiced_length);
    printf("%d samples of voice\n", 16*voiced_length);
    sf_writef_short(filehandle, noise_sound, 8192);
    sf_writef_short(filehandle, noise_sound, C1_NOISE_SAMPLES - 8192);
    printf("%d samples of noise\n", C1_NOISE_SAMPLES);
    sf_writef_short(filehandle, silence_sound, C1_SILENCE_SAMPLES);
    printf("%d samples of silence\n", C1_SILENCE_SAMPLES);

    if (sf_close(filehandle))
    {
        fprintf(stderr, "    Cannot close speech file '%s'\n", "sound_c1.wav");
        exit(2);
    }

    memset(&info, 0, sizeof(info));
    info.frames = 0;
    info.samplerate = FAST_SAMPLE_RATE;
    info.channels = 1;
    info.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
    info.sections = 1;
    info.seekable = 1;
    if ((filehandle = sf_open("sound_c3.wav", SFM_WRITE, &info)) == NULL)
    {
        fprintf(stderr, "    Failed to open result file\n");
        exit(2);
    }

    printf("Generate C3\n");
    /* The sequence is
            72.69ms  of voiced sound from table C.3 of G.168
            200.0ms  of pseudo-noise
            127.31ms of silence
       The above is then repeated phase inverted.

       The voice comes straight from C.3, repeated enough times to
       fill out the 72.69ms period - i.e. 14 copies of the sequence.

       The pseudo noise section is AWGN filtered by the spectral
       pattern in Figure C.2. Since AWGN has the quality of being its
       own Fourier transform, we can use an approach like the one above
       for the C1 signal, using AWGN samples instead of randomly alternating
       ones and zeros. */

    /* Take the supplied set of C3 voice samples. */
    voiced_length = (sizeof(css_c3)/sizeof(css_c3[0]));
    for (i = 0;  i < voiced_length;  i++)
        voiced_sound[i] = css_c3[i];
    pk = peak(voiced_sound, voiced_length);
    ms = rms(voiced_sound, voiced_length);
    printf("Voiced level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));

    noise_source = awgn_init_dbm0(NULL, 7162534, rms_to_dbm0(ms));
    for (i = 0;  i < 8192;  i++)
        noise_sound[i] = awgn(noise_source);
    pk = peak(noise_sound, 8192);
    ms = rms(noise_sound, 8192);
    printf("Unfiltered noise level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));

    /* Now filter them */
#if defined(HAVE_FFTW3_H)
    p = fftw_plan_dft_1d(8192, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
#else
    p = fftw_create_plan(8192, FFTW_BACKWARD, FFTW_ESTIMATE);
#endif
    for (i = 0;  i < 8192;  i++)
    {
#if defined(HAVE_FFTW3_H)
        in[i][0] = 0.0;
        in[i][1] = 0.0;
#else
        in[i].re = 0.0;
        in[i].im = 0.0;
#endif
    }
    for (i = 1;  i <= 3715;  i++)
    {
        f = FAST_SAMPLE_RATE*i/8192.0;

#if 1
        if (f < 50.0)
            scale = -60.0;
        else if (f < 100.0)
            scale = scaling(f, 50.0, 100.0, -25.8, -12.8);
        else if (f < 200.0)
            scale = scaling(f, 100.0, 200.0, -12.8, 17.4);
        else if (f < 215.0)
            scale = scaling(f, 200.0, 215.0, 17.4, 17.8);
        else if (f < 500.0)
            scale = scaling(f, 215.0, 500.0, 17.8, 12.2);
        else if (f < 1000.0)
            scale = scaling(f, 500.0, 1000.0, 12.2, 7.2);
        else if (f < 2850.0)
            scale = scaling(f, 1000.0, 2850.0, 7.2, 0.0);
        else if (f < 3600.0)
            scale = scaling(f, 2850.0, 3600.0, 0.0, -2.0);
        else if (f < 3660.0)
            scale = scaling(f, 3600.0, 3660.0, -2.0, -20.0);
        else if (f < 3680.0)
            scale = scaling(f, 3600.0, 3680.0, -20.0, -30.0);
        else
            scale = -60.0;
#else
        scale = 0.0;
#endif
#if defined(HAVE_FFTW3_H)
        in[i][0] = noise_sound[i]*pow(10.0, scale/20.0)*0.0106;
        in[8192 - i][0] = -in[i][0];
#else
        in[i].re = noise_sound[i]*pow(10.0, scale/20.0)*0.0106;
        in[8192 - i].re = -in[i].re;
#endif
    }
#if defined(HAVE_FFTW3_H)
    fftw_execute(p);
#else
    fftw_one(p, in, out);
#endif
    for (i = 0;  i < 8192;  i++)
    {
#if defined(HAVE_FFTW3_H)
        noise_sound[i] = out[i][1];
#else
        noise_sound[i] = out[i].im;
#endif
    }
    pk = peak(noise_sound, 8192);
    ms = rms(noise_sound, 8192);
    printf("Filtered noise level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));

    for (i = 0;  i < 14;  i++)
        sf_writef_short(filehandle, voiced_sound, voiced_length);
    printf("%d samples of voice\n", 14*voiced_length);
    sf_writef_short(filehandle, noise_sound, 8192);
    sf_writef_short(filehandle, noise_sound, C3_NOISE_SAMPLES - 8192);
    printf("%d samples of noise\n", C3_NOISE_SAMPLES);
    sf_writef_short(filehandle, silence_sound, C3_SILENCE_SAMPLES);
    printf("%d samples of silence\n", C3_SILENCE_SAMPLES);

    /* Now phase invert the C3 set of samples. */
    voiced_length = (sizeof(css_c3)/sizeof(css_c3[0]));
    for (i = 0;  i < voiced_length;  i++)
        voiced_sound[i] = -css_c3[i];
    for (i = 0;  i < 8192;  i++)
        noise_sound[i] = -noise_sound[i];

    for (i = 0;  i < 14;  i++)
        sf_writef_short(filehandle, voiced_sound, voiced_length);
    printf("%d samples of voice\n", 14*voiced_length);
    sf_writef_short(filehandle, noise_sound, 8192);
    sf_writef_short(filehandle, noise_sound, C3_NOISE_SAMPLES - 8192);
    printf("%d samples of noise\n", C3_NOISE_SAMPLES);
    sf_writef_short(filehandle, silence_sound, C3_SILENCE_SAMPLES);
    printf("%d samples of silence\n", C3_SILENCE_SAMPLES);

    if (sf_close(filehandle))
    {
        fprintf(stderr, "    Cannot close speech file '%s'\n", "sound_c3.wav");
        exit(2);
    }

    fftw_destroy_plan(p);
    return 0;
}
/*- End of function --------------------------------------------------------*/
/*- End of file ------------------------------------------------------------*/

max_stars_repo_name: zengfanmao/mpds

max_issues_repo_name: zengfanmao/mpds

max_forks_repo_name: zengfanmao/mpds


text:
#pragma once
#include "IntervalGrad.h"

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "CustomSolver.h"
#endif

#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include "NodeVisitor.h"
#include "VarStore.h"
#include <map>
#include "FloatSupport.h"

#include <iostream>

using namespace std;

#define FINE_GRAIN_RANGES 1

class RangeDiff: 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<IntervalGrad*> ranges; // Keeps track of ranges for each node
	vector<DistanceGrad*> distances; // Keeps track of distance metric for boolean nodes
	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;
	int assertCtr;
	bool foundFailure;

	static constexpr float ASSERT_PENALTY = 1.0;

public:
	int failedAssert;

  RangeDiff(BooleanDAG& bdag_p, FloatManager& _floats, const map<string, int>& floatCtrls_p);
  ~RangeDiff(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 setrange(bool_node& bn, IntervalGrad* r) {
		ranges[bn.id] = r;
	}

	IntervalGrad* r(bool_node& bn) {
		IntervalGrad* interval = ranges[bn.id];
		if (interval == NULL) {
			gsl_vector* l = gsl_vector_alloc(nctrls);
			gsl_vector* h = gsl_vector_alloc(nctrls);
			interval = new IntervalGrad(0, 0, l, h);
			setrange(bn, interval);
		}
		return interval;
	}

	IntervalGrad* r(bool_node* bn) {
		return r(*bn);
	}

	void setdistance(bool_node& bn, DistanceGrad* d) {
		distances[bn.id] = d;
	}

	DistanceGrad* d(bool_node& bn) {
		DistanceGrad* dist = distances[bn.id];
		if (dist == NULL) {
			gsl_vector* g = gsl_vector_alloc(nctrls);
			dist = new DistanceGrad(0, g);
			setdistance(bn, dist);
		}
		return dist;
	}

	DistanceGrad* d(bool_node* bn) {
		return d(*bn);
	}

	void print() {
		for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << " ";
			if (bdag[i]->getOtype() == OutType::FLOAT) {
				cout << r(bdag[i])->print() << endl;
			} else {
				cout << d(bdag[i])->print() << endl;
			}
		}
	}
	void printFull() {
    for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << endl;
			if (bdag[i]->getOtype() == OutType::FLOAT) {
				cout << r(bdag[i])->printFull() << endl;
			} else {
				cout << d(bdag[i])->printFull() << 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);
	}
	void computeError(double dist, int expected, gsl_vector* dg, bool_node& node, bool relax = false);

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend