text:
#pragma once

#include <cassert>
#include <functional>
#include <map>
#include <type_traits>
#include <vector>

#include <gsl/assert>

namespace xmol::utils {

class DeadObserverAccessError : public std::runtime_error {
public:
  using std::runtime_error::runtime_error;
};

template <typename Observer>
class DeadObserverAccessErrorT : public DeadObserverAccessError { /// todo: parametrize on error message instead
public:
  using DeadObserverAccessError::DeadObserverAccessError;
};

enum class ObserverState { ANY, ACTIVE, INVALID };

/// @brief Implements base primitives for observable entity
template <typename Observer> class Observable {
  static_assert(!std::is_reference<Observer>::value);
  static_assert(!std::is_pointer<Observer>::value);

public:
  Observable() = default;
  Observable(Observable&& rhs) noexcept = default;
  Observable(const Observable& rhs) = default;
  Observable& operator=(Observable&& rhs) noexcept = default;
  Observable& operator=(const Observable& rhs) = default;

protected:
  template <ObserverState apply_to = ObserverState::ANY, typename... Args, typename Func = void (Observer::*)(Args...)>
  void notify(Func func, Args&&... args) const {
    static_assert(apply_to != ObserverState::INVALID);
    for (auto& [observer, state] : observers) {
      if (state == ObserverState::ACTIVE) {
        std::invoke(func, observer, std::forward<Args>(args)...);
      } else {
        if (GSL_UNLIKELY(apply_to == ObserverState::ANY)) {
          throw DeadObserverAccessErrorT<Observer>("");
        }
      }
    }
  }

  void add_observer(Observer& ptr) const { observers.emplace(&ptr, ObserverState::ACTIVE); }

  void remove_observer(Observer& ptr) const {
    auto count = observers.erase(&ptr);
    static_cast<void>(count);
    assert(count == 1);
  }

  void clear_observers() const { observers.clear(); }

  /// @brief marks observer as invalid
  /// next broadcast notify would fire an exception
  void invalidate_observer(Observer& ptr) const {
    auto it = observers.find(&ptr);
    assert(it != observers.end());
    it->second = ObserverState::INVALID;
  }

  void move_observer(Observer& from, Observer& to) const {
    remove_observer(from);
    add_observer(to);
  }

public:

  void on_move(Observer& from, Observer& to) {
    remove_observer(from);
    add_observer(to);
  }

  void on_delete(Observer& o) { remove_observer(o); }
  void on_copy(Observer& o) { add_observer(o); }

protected:
  mutable std::map<Observer*, ObserverState> observers;
};

} // namespace xmol


text:
/*
   Copyright [2017-2020] [IBM Corporation]
   Licensed 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 __FABRIC_CONNECTION_BASE_H__
#define __FABRIC_CONNECTION_BASE_H__

#include "mcas_config.h"

#include "buffer_manager.h" /* Buffer_manager */

#include <api/fabric_itf.h> /* IFabric_memory_region, IFabric_server */
#include <gsl/pointers>
#include <algorithm>
#include <list>
#include <queue>

namespace mcas
{

struct open_connection_construct_key;

struct open_connection
{
private:
  gsl::not_null<component::IFabric_server_factory *> _factory;
  gsl::not_null<component::IFabric_server *>         _transport;
public:
  open_connection(
    gsl::not_null<component::IFabric_server_factory *> factory_
    , gsl::not_null<component::IFabric_server *> transport_
    , const open_connection_construct_key &
    )
    : _factory(factory_)
    , _transport(transport_)
  {}

  open_connection(const open_connection &) = delete;
  open_connection &operator=(const open_connection &) = delete;

  auto transport() const { return _transport; }

  ~open_connection()
  {
    /* close connection */
    _factory->close_connection(_transport);
  }
};

class Fabric_connection_base : protected common::log_source {

 public:
  friend class Connection;
  friend class Shard;
  using memory_region_t = component::IFabric_memory_region *;

 protected:
  using buffer_t = Buffer_manager<component::IFabric_memory_control>::buffer_internal;
  using pool_t   = component::IKVStore::pool_t;

  enum class action_type {
        ACTION_NONE = 0, /* unused */
        ACTION_RELEASE_VALUE_LOCK_SHARED,
#if 0
        ACTION_RELEASE_VALUE_LOCK_EXCLUSIVE, /* unused */
        ACTION_POOL_DELETE, /* unused */
#endif
  };

  /* deferred actions */
  struct action_t {
    action_type op;
    void *parm;
  };

  enum class Completion_state {
    COMPLETIONS,
    CLIENT_DISCONNECT,
    NONE,
  };

 private:
  std::unique_ptr<component::IFabric_endpoint_unconnected_server> _preconnection;
  Buffer_manager<component::IFabric_memory_control> _bm;

  /* xx_buffer_outstanding is the signal for completion,
     xx_buffer is the buffer pointer that needs to be freed (and set to null)
  */
protected:
  unsigned              _recv_buffer_posted_count;
private:
  unsigned              _send_buffer_posted_count;

protected:
  /* values for two-phase get
   */
  unsigned _send_value_posted_count;
private:

  std::list<buffer_t *> _completed_recv_buffers;

  open_connection _oc;
 protected:
  size_t                             _max_message_size;

  /* Filled by base class     check_for_posted_value_complete
   * Drained by derived class check_network_completions
   */
  std::queue<action_t> _deferred_unlock;

  /**
   * Ctor
   *
   * @param factory
   * @param fabric_connection
   */
  explicit Fabric_connection_base( //
    unsigned                           debug_level_,
    gsl::not_null<component::IFabric_server_factory *>factory,
    std::unique_ptr<component::IFabric_endpoint_unconnected_server> && preconnection);

  Fabric_connection_base(const Fabric_connection_base &) = delete;
  Fabric_connection_base &operator=(const Fabric_connection_base &) = delete;

  virtual ~Fabric_connection_base();

  static void static_send_callback(void *cnxn, buffer_t *iob) noexcept
  {
    static_cast<Fabric_connection_base *>(cnxn)->send_callback(iob);
  }

  void send_callback(buffer_t *iob) noexcept
  {
    --_send_buffer_posted_count;
    posted_count_log();
    CPLOG(2, "Completed send (%p) freeing buffer", common::p_fmt(iob));
    free_buffer(iob);
  }

  static void static_recv_callback(void *cnxn, buffer_t *iob) noexcept
  {
    static_cast<Fabric_connection_base *>(cnxn)->recv_callback(iob);
  }

  void recv_callback(buffer_t *iob) noexcept
  {
    --_recv_buffer_posted_count;
    posted_count_log();
    _completed_recv_buffers.push_front(iob);
    CPLOG(2, "Completed recv (%p) (complete %zu)", common::p_fmt(iob), _completed_recv_buffers.size());
  }

  static void completion_callback(void *   context,
                                  status_t st,
                                  std::uint64_t,  // completion_flags,
                                  std::size_t len,
                                  void *      error_data,
                                  void *      cnxn) noexcept
  {
    if (LIKELY(st == S_OK)) {
      auto iob = static_cast<buffer_t *>(context);
      iob->completion_cb(cnxn, iob);
    }
    else {
      PERR("Fabric_connection_base: fabric operation failed st != S_OK (st=%d, "
           "context=%p, len=%lu)",
           st, context, len);
      PERR("Error: %s", static_cast<char *>(error_data));
    }
  }

  bool check_for_posted_recv_complete()
  {
    /* don't free buffer (such as above); it will be used for response */
    //    return _recv_buffer_posted_outstanding == false;
    return _completed_recv_buffers.size() > 0;
  }

  void free_recv_buffer()
  {
    for (auto b : _completed_recv_buffers) {
      free_buffer(b);
    }
  }

  size_t recv_buffer_posted_count() const {
    return _recv_buffer_posted_count;
  }

  void post_recv_buffer(buffer_t *buffer)
  {
    _preconnection->post_recv(buffer->iov, buffer->iov + 1, buffer->desc, buffer);
    ++_recv_buffer_posted_count;
    posted_count_log();
    CPLOG(2, "Posted recv (%p) (complete %zu)", common::p_fmt(buffer), _completed_recv_buffers.size());
  }

  void post_send_buffer(gsl::not_null<buffer_t *> buffer)
  {
    const auto iov = buffer->iov;

    /* if packet is small enough use inject */
    if (iov->iov_len <= transport()->max_inject_size()) {
      CPLOG(2, "Fabric_connection_base: posting send with inject (iob %p %p,len=%lu)",
             common::p_fmt(buffer), iov->iov_base, iov->iov_len);

      transport()->inject_send(iov->iov_base, iov->iov_len);
      CPLOG(2, "%s: buffer %p", __func__, common::p_fmt(buffer));
      free_buffer(buffer); /* buffer is immediately complete fi_inject */
    }
    else {
      CPLOG(2, "Fabric_connection_base: posting send (%p, %p)", common::p_fmt(buffer), iov->iov_base);

      transport()->post_send(iov, iov + 1, buffer->desc, buffer);
      ++_send_buffer_posted_count;
      posted_count_log();
      CPLOG(2, "%s buffer (%p)", __func__, common::p_fmt(buffer));
    }
  }

  void post_send_buffer2(gsl::not_null<buffer_t *> buffer, const ::iovec &val_iov, void *val_desc)
  {
    buffer->iov[1] = val_iov;
    buffer->desc[1] = val_desc;

    ++_send_buffer_posted_count;
    posted_count_log();
    CPLOG(2, "Posted send (%p) ... value (%.*s) (len=%lu,ptr=%p)", common::p_fmt(buffer),
         int(val_iov.iov_len), static_cast<char *>(val_iov.iov_base), val_iov.iov_len, val_iov.iov_base);

    transport()->post_send(buffer->iov, buffer->iov + 2, buffer->desc, buffer);
  }

  buffer_t *posted_recv()
  {
    if (_completed_recv_buffers.size() == 0) return nullptr;

    auto iob = _completed_recv_buffers.back();
    _completed_recv_buffers.pop_back();
    CPLOG(2, "Presented recv (%p) (complete %zu)", common::p_fmt(iob), _completed_recv_buffers.size());
    return iob;
  }

  void posted_count_log() const
  {
    CPLOG(2, "POSTs recv %u send %u value %u", _recv_buffer_posted_count, _send_buffer_posted_count,
           _send_value_posted_count);
  }

  Completion_state poll_completions()
  {
    if (_recv_buffer_posted_count != 0 || _send_buffer_posted_count != 0 || _send_value_posted_count != 0) {
#if 0
      PLOG("%s posted_recv %u posted_send %u posted_value %u", __func__, _recv_buffer_posted_count, _send_buffer_posted_count, _send_value_posted_count);
#endif
      try {
        transport()->poll_completions(&Fabric_connection_base::completion_callback, this);
        return Completion_state::COMPLETIONS;
      }
      catch (const std::logic_error &e) {
        return Completion_state::CLIENT_DISCONNECT;
      }
    }
    return Completion_state::NONE;
  }

  /**
   * Forwarders that allow us to avoid exposing transport() and _bm
   *
   */
 public:
  inline auto register_memory(const void *base, size_t len, std::uint64_t key, std::uint64_t flags)
  {
    return transport()->register_memory(base, len, key, flags); /* flags not supported for verbs */
  }

  inline void deregister_memory(memory_region_t region) { return transport()->deregister_memory(region); }

  inline void *get_memory_descriptor(memory_region_t region) { return transport()->get_memory_descriptor(region); }

  inline uint64_t get_memory_remote_key(memory_region_t region) { return transport()->get_memory_remote_key(region); }

 protected:
  inline auto allocate(buffer_t::completion_t c) { return _bm.allocate(c); }

  inline void free_buffer(buffer_t *buffer) { _bm.free(buffer); }

  inline size_t IO_buffer_size() const { return Buffer_manager<component::IFabric_memory_control>::BUFFER_LEN; }

  gsl::not_null<component::IFabric_server *> transport() const { return _oc.transport(); }

  inline std::string get_local_addr() { return transport()->get_local_addr(); }
};

struct open_connection_construct_key
{
private:
  open_connection_construct_key() {};
public:
  friend class Fabric_connection_base;
};

}  // namespace mcas

#endif  // __FABRIC_CONNECTION_BASE_H__

ext: h

max_stars_repo_name: fQuinzan/mcas

max_issues_repo_name: fQuinzan/mcas

max_forks_repo_name: fQuinzan/mcas


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk_6_4.c
 * \brief Source file to optimize Runge-Kutta 6 steps 4th order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_6_4.h"

#define DEBUG_RK_6_4 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 6 steps 4th order Runge-Kutta
 * method.
 */
int
rk_tb_6_4 (Optimize * optimize) ///< Optimize struct.
{
  long double A[4], B[4], C[4], D[4], E[4];
  long double *tb, *r;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_tb_6_4: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t6 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  b41 (tb) = r[7];
  b42 (tb) = r[8];
  b43 (tb) = r[9];
  t5 (tb) = r[10];
  b54 (tb) = r[11];
  b65 (tb) = r[12];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = t4 (tb);
  E[0] = 0.5L - b65 (tb) * t5 (tb);
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = D[0] * t4 (tb);
  E[1] = 1.L / 3.L - b65 (tb) * sqr (t5 (tb));
  A[2] = A[1] * t1 (tb);
  B[2] = B[1] * t2 (tb);
  C[2] = C[1] * t3 (tb);
  D[2] = D[1] * t4 (tb);
  E[2] = 0.25L - b65 (tb) * sqr (t5 (tb)) * t5 (tb);
  A[3] = 0.L;
  B[3] = b21 (tb) * t1 (tb) * (t2 (tb) - t5 (tb));
  C[3] = (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb)) * (t3 (tb) - t5 (tb));
  D[3] = (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb) + b43 (tb) * t3 (tb))
    * (t4 (tb) - t5 (tb));
  E[3] = 0.125L - 1.L / 6.L * t5 (tb);
  solve_4 (A, B, C, D, E);
  if (isnan (E[0]) || isnan (E[1]) || isnan (E[2]) || isnan (E[3]))
    return 0;
  b64 (tb) = E[3];
  b63 (tb) = E[2];
  b62 (tb) = E[1];
  b61 (tb) = E[0];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = (1.L / 6.L - b62 (tb) * b21 (tb) * t1 (tb)
          - b63 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))
          - b64 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                        + b43 (tb) * t3 (tb))) / b65 (tb) - b54 (tb) * t4 (tb);
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = (1.L / 12.L - b62 (tb) * b21 (tb) * sqr (t1 (tb))
          - b63 (tb) * (b31 (tb) * sqr (t1 (tb)) + b32 (tb) * sqr (t2 (tb)))
          - b64 (tb) * (b41 (tb) * sqr (t1 (tb)) + b42 (tb) * sqr (t2 (tb))
                        + b43 (tb) * sqr (t3 (tb)))) / b65 (tb)
    - b54 (tb) * sqr (t4 (tb));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = (1.L / 24.L - b63 (tb) * b32 (tb) * b21 (tb) * t1 (tb)
          - b64 (tb) * (b42 (tb) * b21 (tb) * t1 (tb)
                        + b43 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))))
    / b65 (tb) - b54 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                             + b43 (tb) * t3 (tb));
  solve_3 (A, B, C, D);
  if (isnan (D[0]) || isnan (D[1]) || isnan (D[2]))
    return 0;
  b53 (tb) = D[2];
  b52 (tb) = D[1];
  b51 (tb) = D[0];
  rk_b_6 (tb);
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_tb_6_4: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 4th order, 5th order in
 * equations depending only on time, Runge-Kutta method.
 */
int
rk_tb_6_4t (Optimize * optimize)        ///< Optimize struct.
{
  long double A[5], B[5], C[5], D[5], E[5], F[5];
  long double *tb, *r;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_tb_6_4t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t6 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  b41 (tb) = r[7];
  b42 (tb) = r[8];
  b43 (tb) = r[9];
  t5 (tb) = r[10];
  b54 (tb) = r[11];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = t4 (tb);
  E[0] = t5 (tb);
  F[0] = 0.5L;
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = D[0] * t4 (tb);
  E[1] = E[0] * t5 (tb);
  F[1] = 1.L / 3.L;
  A[2] = A[1] * t1 (tb);
  B[2] = B[1] * t2 (tb);
  C[2] = C[1] * t3 (tb);
  D[2] = D[1] * t4 (tb);
  E[2] = E[1] * t5 (tb);
  F[2] = 0.25L;
  A[3] = A[2] * t1 (tb);
  B[3] = B[2] * t2 (tb);
  C[3] = C[2] * t3 (tb);
  D[3] = D[2] * t4 (tb);
  E[3] = E[2] * t5 (tb);
  F[3] = 0.2L;
  A[4] = A[3] * t1 (tb);
  B[4] = B[3] * t2 (tb);
  C[4] = C[3] * t3 (tb);
  D[4] = D[3] * t4 (tb);
  E[4] = E[3] * t5 (tb);
  F[4] = 1.L / 6.L;
  solve_5 (A, B, C, D, E, F);
  if (isnan (F[0]) || isnan (F[1]) || isnan (F[2]) || isnan (F[3])
      || isnan (F[4]))
    return 0;
  b65 (tb) = F[4];
  b64 (tb) = F[3];
  b63 (tb) = F[2];
  b62 (tb) = F[1];
  b61 (tb) = F[0];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = (1.L / 6.L - b62 (tb) * b21 (tb) * t1 (tb)
          - b63 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))
          - b64 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                        + b43 (tb) * t3 (tb))) / b65 (tb) - b54 (tb) * t4 (tb);
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = (1.L / 12.L - b62 (tb) * b21 (tb) * sqr (t1 (tb))
          - b63 (tb) * (b31 (tb) * sqr (t1 (tb)) + b32 (tb) * sqr (t2 (tb)))
          - b64 (tb) * (b41 (tb) * sqr (t1 (tb)) + b42 (tb) * sqr (t2 (tb))
                        + b43 (tb) * sqr (t3 (tb)))) / b65 (tb)
    - b54 (tb) * sqr (t4 (tb));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = (1.L / 24.L - b63 (tb) * b32 (tb) * b21 (tb) * t1 (tb)
          - b64 (tb) * (b42 (tb) * b21 (tb) * t1 (tb)
                        + b43 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))))
    / b65 (tb) - b54 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                             + b43 (tb) * t3 (tb));
  solve_3 (A, B, C, D);
  if (isnan (D[0]) || isnan (D[1]) || isnan (D[2]))
    return 0;
  b53 (tb) = D[2];
  b52 (tb) = D[1];
  b51 (tb) = D[0];
  rk_b_6 (tb);
#if DEBUG_RK_6_4
  rk_print_tb (optimize, "rk_tb_6_4t", stderr);
  fprintf (stderr, "rk_tb_6_4t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 3th-4th order Runge-Kutta
 * pair.
 */
int
rk_tb_6_4p (Optimize * optimize)        ///< Optimize struct.
{
  long double A[4], B[4], C[4], D[4], E[4], AA[4], BB[4], CC[4], DD[4], EE[4];
  long double *tb, *r;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_tb_6_4p: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t6 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  b41 (tb) = r[7];
  b42 (tb) = r[8];
  b43 (tb) = r[9];
  t5 (tb) = r[10];
  b54 (tb) = r[11];
  b65 (tb) = r[12];
  A[0] = AA[0] = t1 (tb);
  B[0] = BB[0] = t2 (tb);
  C[0] = CC[0] = t3 (tb);
  D[0] = DD[0] = t4 (tb);
  EE[0] = 0.5L;
  E[0] = 0.5L - b65 (tb) * t5 (tb);
  A[1] = AA[1] = A[0] * t1 (tb);
  B[1] = BB[1] = B[0] * t2 (tb);
  C[1] = CC[1] = C[0] * t3 (tb);
  D[1] = DD[1] = D[0] * t4 (tb);
  EE[1] = 1.L / 3.L;
  E[1] = 1.L / 3.L - b65 (tb) * sqr (t5 (tb));
  A[2] = AA[2] = A[1] * t1 (tb);
  B[2] = BB[2] = B[1] * t2 (tb);
  C[2] = CC[2] = C[1] * t3 (tb);
  D[2] = DD[2] = D[1] * t4 (tb);
  EE[2] = 0.25L;
  E[2] = 0.25L - b65 (tb) * sqr (t5 (tb)) * t5 (tb);
  A[3] = AA[3] = 0.L;
  BB[3] = b21 (tb) * t1 (tb);
  B[3] = BB[3] * (t2 (tb) - t5 (tb));
  CC[3] = (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb));
  C[3] = CC[3] * (t3 (tb) - t5 (tb));
  DD[3] = (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb) + b43 (tb) * t3 (tb));
  D[3] = DD[3] * (t4 (tb) - t5 (tb));
  EE[3] = 1.L / 6.L;
  E[3] = 0.125L - 1.L / 6.L * t5 (tb);
  solve_4 (A, B, C, D, E);
  if (isnan (E[0]) || isnan (E[1]) || isnan (E[2]) || isnan (E[3]))
    return 0;
  b64 (tb) = E[3];
  b63 (tb) = E[2];
  b62 (tb) = E[1];
  b61 (tb) = E[0];
  solve_4 (AA, BB, CC, DD, EE);
  if (isnan (EE[0]) || isnan (EE[1]) || isnan (EE[2]) || isnan (EE[3]))
    return 0;
  e64 (tb) = EE[3];
  e63 (tb) = EE[2];
  e62 (tb) = EE[1];
  e61 (tb) = EE[0];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = (1.L / 6.L - b62 (tb) * b21 (tb) * t1 (tb)
          - b63 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))
          - b64 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                        + b43 (tb) * t3 (tb))) / b65 (tb) - b54 (tb) * t4 (tb);
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = (1.L / 12.L - b62 (tb) * b21 (tb) * sqr (t1 (tb))
          - b63 (tb) * (b31 (tb) * sqr (t1 (tb)) + b32 (tb) * sqr (t2 (tb)))
          - b64 (tb) * (b41 (tb) * sqr (t1 (tb)) + b42 (tb) * sqr (t2 (tb))
                        + b43 (tb) * sqr (t3 (tb)))) / b65 (tb)
    - b54 (tb) * sqr (t4 (tb));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = (1.L / 24.L - b63 (tb) * b32 (tb) * b21 (tb) * t1 (tb)
          - b64 (tb) * (b42 (tb) * b21 (tb) * t1 (tb)
                        + b43 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))))
    / b65 (tb) - b54 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                             + b43 (tb) * t3 (tb));
  solve_3 (A, B, C, D);
  if (isnan (D[0]) || isnan (D[1]) || isnan (D[2]))
    return 0;
  b53 (tb) = D[2];
  b52 (tb) = D[1];
  b51 (tb) = D[0];
  rk_b_6 (tb);
  rk_e_6 (tb);
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_tb_6_4p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 3rd-4th order, 4th-5th order
 * in equations depending only on time, Runge-Kutta method.
 */
int
rk_tb_6_4tp (Optimize * optimize)       ///< Optimize struct.
{
  long double A[5], B[5], C[5], D[5], E[5], F[5], AA[4], BB[4], CC[4], DD[4],
    EE[4];
  long double *tb, *r;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_tb_6_4tp: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t6 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  t3 (tb) = r[2];
  b31 (tb) = r[3];
  b32 (tb) = r[4];
  t4 (tb) = r[5];
  b41 (tb) = r[6];
  b42 (tb) = r[7];
  b43 (tb) = r[8];
  t5 (tb) = r[9];
  b54 (tb) = r[10];
  b65 (tb) = r[11];
  A[0] = AA[0] = t1 (tb);
  B[0] = BB[0] = t2 (tb);
  C[0] = CC[0] = t3 (tb);
  D[0] = DD[0] = t4 (tb);
  E[0] = t5 (tb);
  F[0] = EE[0] = 0.5L;
  A[1] = AA[1] = A[0] * t1 (tb);
  B[1] = BB[1] = B[0] * t2 (tb);
  C[1] = CC[1] = C[0] * t3 (tb);
  D[1] = DD[1] = D[0] * t4 (tb);
  E[1] = E[0] * t5 (tb);
  F[1] = EE[1] = 1.L / 3.L;
  A[2] = AA[2] = A[1] * t1 (tb);
  B[2] = BB[2] = B[1] * t2 (tb);
  C[2] = CC[2] = C[1] * t3 (tb);
  D[2] = DD[2] = D[1] * t4 (tb);
  E[2] = E[1] * t5 (tb);
  F[2] = EE[2] = 0.25L;
  A[3] = AA[3] = A[2] * t1 (tb);
  B[3] = BB[3] = B[2] * t2 (tb);
  C[3] = CC[3] = C[2] * t3 (tb);
  D[3] = DD[3] = D[2] * t4 (tb);
  E[3] = E[2] * t5 (tb);
  F[3] = EE[3] = 0.2L;
  A[4] = A[3] * t1 (tb);
  B[4] = B[3] * t2 (tb);
  C[4] = C[3] * t3 (tb);
  D[4] = D[3] * t4 (tb);
  E[4] = E[3] * t5 (tb);
  F[4] = 1.L / 6.L;
  solve_4 (AA, BB, CC, DD, EE);
  if (isnan (EE[0]) || isnan (EE[1]) || isnan (EE[2]) || isnan (EE[3]))
    return 0;
  e64 (tb) = EE[3];
  e63 (tb) = EE[2];
  e62 (tb) = EE[1];
  e61 (tb) = EE[0];
  solve_5 (A, B, C, D, E, F);
  if (isnan (F[0]) || isnan (F[1]) || isnan (F[2]) || isnan (F[3])
      || isnan (F[4]))
    return 0;
  b65 (tb) = F[4];
  b64 (tb) = F[3];
  b63 (tb) = F[2];
  b62 (tb) = F[1];
  b61 (tb) = F[0];
  b21 (tb) = (1.L / 6.L - e63 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))
              - e64 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                            + b43 (tb) * t3 (tb))) / (e62 (tb) * t1 (tb));
  if (isnan (b21 (tb)))
    return 0;
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = (1.L / 6.L - b62 (tb) * b21 (tb) * t1 (tb)
          - b63 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))
          - b64 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                        + b43 (tb) * t3 (tb))) / b65 (tb) - b54 (tb) * t4 (tb);
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = (1.L / 12.L - b62 (tb) * b21 (tb) * sqr (t1 (tb))
          - b63 (tb) * (b31 (tb) * sqr (t1 (tb)) + b32 (tb) * sqr (t2 (tb)))
          - b64 (tb) * (b41 (tb) * sqr (t1 (tb)) + b42 (tb) * sqr (t2 (tb))
                        + b43 (tb) * sqr (t3 (tb)))) / b65 (tb)
    - b54 (tb) * sqr (t4 (tb));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = (1.L / 24.L - b63 (tb) * b32 (tb) * b21 (tb) * t1 (tb)
          - b64 (tb) * (b42 (tb) * b21 (tb) * t1 (tb)
                        + b43 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb))))
    / b65 (tb) - b54 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                             + b43 (tb) * t3 (tb));
  solve_3 (A, B, C, D);
  if (isnan (D[0]) || isnan (D[1]) || isnan (D[2]))
    return 0;
  b53 (tb) = D[2];
  b52 (tb) = D[1];
  b51 (tb) = D[0];
  rk_b_6 (tb);
#if DEBUG_RK_6_4
  rk_print_tb (optimize, "rk_tb_6_4t", stderr);
  fprintf (stderr, "rk_tb_6_4t: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 6 steps 4th order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_4 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b60 (tb) < 0.L)
    o += b60 (tb);
  if (b61 (tb) < 0.L)
    o += b61 (tb);
  if (b62 (tb) < 0.L)
    o += b62 (tb);
  if (b63 (tb) < 0.L)
    o += b63 (tb);
  if (b64 (tb) < 0.L)
    o += b64 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L,
             fmaxl (t1 (tb),
                    fmaxl (t2 (tb),
                           fmaxl (t3 (tb), fmaxl (t4 (tb), t5 (tb))))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_4: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 6 steps 4th order, 5th
 * order in equations depending only on time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_4t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4t: start\n");
#endif
  tb = rk->tb->coefficient;
#if DEBUG_RK_6_4
  rk_print_tb (optimize, "rk_objective_tb_6_4t", stderr);
#endif
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b60 (tb) < 0.L)
    o += b60 (tb);
  if (b61 (tb) < 0.L)
    o += b61 (tb);
  if (b62 (tb) < 0.L)
    o += b62 (tb);
  if (b63 (tb) < 0.L)
    o += b63 (tb);
  if (b64 (tb) < 0.L)
    o += b64 (tb);
  if (b65 (tb) < 0.L)
    o += b65 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L,
             fmaxl (t1 (tb),
                    fmaxl (t2 (tb),
                           fmaxl (t3 (tb), fmaxl (t4 (tb), t5 (tb))))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_4t: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 6 steps 3rd-4th order
 * Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_4p (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4p: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b60 (tb) < 0.L)
    o += b60 (tb);
  if (b61 (tb) < 0.L)
    o += b61 (tb);
  if (b62 (tb) < 0.L)
    o += b62 (tb);
  if (b63 (tb) < 0.L)
    o += b63 (tb);
  if (b64 (tb) < 0.L)
    o += b64 (tb);
  if (e60 (tb) < 0.L)
    o += e60 (tb);
  if (e61 (tb) < 0.L)
    o += e61 (tb);
  if (e62 (tb) < 0.L)
    o += e62 (tb);
  if (e63 (tb) < 0.L)
    o += e63 (tb);
  if (e64 (tb) < 0.L)
    o += e64 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L,
             fmaxl (t1 (tb),
                    fmaxl (t2 (tb),
                           fmaxl (t3 (tb), fmaxl (t4 (tb), t5 (tb))))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4p: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_4p: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 6 steps 3rd-4th order,
 * 4th-5th order in equations depending only in time, Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_4tp (RK * rk) ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4tp: start\n");
#endif
  tb = rk->tb->coefficient;
#if DEBUG_RK_6_4
  rk_print_tb (optimize, "rk_objective_tb_6_4tp", stderr);
#endif
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b60 (tb) < 0.L)
    o += b60 (tb);
  if (b61 (tb) < 0.L)
    o += b61 (tb);
  if (b62 (tb) < 0.L)
    o += b62 (tb);
  if (b63 (tb) < 0.L)
    o += b63 (tb);
  if (b64 (tb) < 0.L)
    o += b64 (tb);
  if (b65 (tb) < 0.L)
    o += b65 (tb);
  if (e60 (tb) < 0.L)
    o += e60 (tb);
  if (e61 (tb) < 0.L)
    o += e61 (tb);
  if (e62 (tb) < 0.L)
    o += e62 (tb);
  if (e63 (tb) < 0.L)
    o += e63 (tb);
  if (e64 (tb) < 0.L)
    o += e64 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L,
             fmaxl (t1 (tb),
                    fmaxl (t2 (tb),
                           fmaxl (t3 (tb), fmaxl (t4 (tb), t5 (tb))))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_6_4
  fprintf (stderr, "rk_objective_tb_6_4tp: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_4tp: end\n");
#endif
  return o;
}

ext: c

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
#pragma once

#include <string>
#include <vector>

#include <boost/random.hpp>
#include <cblas.h>

namespace con {

  using std::vector;
  using std::cout;
  using std::cerr;
  using std::endl;
  using std::string;

  static boost::mt19937 rng(0);

  #define BUG(x) std::cout<<#x<<" = "<<(x)<<std::endl;

  typedef double Real;

  typedef vector<Real> Vec;

  Real sigmoid(const Real &z) {
    return 1.0 / (1.0 + std::exp(-z));
  }

  Real derivativeSigmoid(const Real &v) {
    return v * (1.0 - v);
  }

  Real sqr(const Real &x) {
    return x * x;
  }

  Real randomize(const Real &min, const Real &max) {
    boost::uniform_real<Real> dst(min, max);
    return dst(rng);
  }

  void randomizeVec(const Real &min, const Real &max, Vec *a) {
    for (auto it = a->begin(); it != a->end(); it++) {
      (*it) = randomize(min, max);
    }
  }

  void gaussianRng(const Real &mean, const Real &std, Vec *a) {
    boost::normal_distribution<Real> random_distribution(mean, std);
    boost::variate_generator<boost::mt19937, boost::normal_distribution<Real>>
      variate_generator(rng, random_distribution);

    for (auto it = a->begin(); it != a->end(); it++) {
      *it = variate_generator();
    }
  }

  void clear(Vec *a) {
    for (auto it = a->begin(); it != a->end(); it++) {
      (*it) = 0;
    }
  }

  void clear(vector<Vec> *a) {
    for (int i = 0; i < a->size(); i++) {
      clear(&a->at(i));
    }
  }

  void print(const Vec &a) {
    for (auto x : a) {
      cout << x << " ";
    }
    cout << endl;
  }

  void reshape(const int num, const int width, const int height, const int depth, vector<Vec> *a) {
    a->resize(num);
    for (int i = 0; i < num; i++) {
      a->at(i).resize(width * height * depth);
    }
  }

  int ceilDiv(const int &x, const int &y) {
    return (x + y - 1) / y;
  }

  void gemm(
      const CBLAS_TRANSPOSE &TransA, const CBLAS_TRANSPOSE &TransB,
      const int &M, const int &N, const int &K,
      const Real &alpha, const Vec &vecA, const Vec &vecB,
      const Real &beta, Vec *vecC) {

    const Real *A = &vecA[0];
    const Real *B = &vecB[0];
    Real *C = &vecC->at(0);

    int lda = (TransA == CblasNoTrans) ? K : M;
    int ldb = (TransB == CblasNoTrans) ? N : K;

    cblas_dgemm(CblasRowMajor, TransA, TransB, M, N, K, alpha, A, lda, B,
        ldb, beta, C, N);
  }

  void gemv(
      const CBLAS_TRANSPOSE &TransA,
      const int &M, const int &N,
      const Real &alpha, const Vec &matA, const Vec &vecX,
      const double beta, Vec *vecY) {

    const Real *A = &matA[0];
    const Real *x = &vecX[0];
    Real *y = &vecY->at(0);

    cblas_dgemv(CblasRowMajor, TransA, M, N, alpha, A, N, x, 1, beta, y, 1);
  }

  void vexp(const int &n, const Vec &input, Vec *output) {
    // vsExp(n, &input[0], &output->at(0));
    for (int i = 0; i < n; i++) {
      output->at(i) = exp(input[i]);
    }
  }

  void vdiv(const int &n, const Vec &a, const Vec &b, Vec *c) {
    // vsDiv(n, &a[0], &b[0], &c->at(0));
    for (int i = 0; i < n; i++) {
      c->at(i) = a[i] / b[i];
    }
  }

  void copy(const int &n, const Vec &input, Vec *output) {
    for (int i = 0; i < n; i++) {
      output->at(i) = input[i];
    }
  }

  void ones(const int &n, Vec *v) {
    v->resize(n);
    std::fill(v->begin(), v->end(), 1.0);
  }
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_core.h>
#include <gbpCosmo_NFW_etc.h>
#include <gsl/gsl_sf_expint.h>

// From Alam et al '02
double R_half_V_max_NFW(double M_vir, double z, int mode, cosmo_info **cosmo) {
    double c_vir;
    double R_vir;
    double r_o;
    double r_val = 0.;

    set_NFW_params(M_vir, z, mode, cosmo, &c_vir, &R_vir);

    r_o   = R_vir / c_vir;
    r_val = 0.13 * r_o;

    return (r_val);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/*
** LISA, 2nd level using design matrix
**
** G.Lohmann, July 2018
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>

#define ABS(x) ((x) > 0 ? (x) : -(x))
#define SQR(x) ((x)*(x))

extern gsl_matrix *PseudoInv(gsl_matrix *A,gsl_matrix *B);
extern gsl_matrix *XRead2ndLevel(VString);
extern double t2z(double t,double df);

double EstimateVariance(gsl_matrix *XInv,gsl_vector *con)
{
  int nbeta = XInv->size1;
  double var=0;

  gsl_matrix *bcov = gsl_matrix_calloc(nbeta,nbeta);
  gsl_blas_dgemm(CblasNoTrans,CblasTrans,1.0,XInv,XInv,0.0,bcov);

  gsl_vector *tmp  = gsl_vector_alloc(nbeta);
  gsl_blas_dgemv(CblasTrans,1.0,bcov,con,0.0,tmp);
  gsl_blas_ddot (tmp,con,&var);
  double sigma = sqrt(var);
  gsl_matrix_free(bcov);
  gsl_vector_free(tmp);

  return sigma;
}

/*
** get effective degrees of freedom
*/
double DoF(gsl_matrix *X,gsl_matrix *XInv,double *xtrace)
{
  int i,nimages = X->size1;

  gsl_matrix *R = gsl_matrix_calloc (nimages,nimages);
  gsl_matrix *P = gsl_matrix_calloc (nimages,nimages);
  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,X,XInv,0.0,P);
  gsl_matrix_set_identity(R);
  gsl_matrix_sub(R,P);

  double trace = 0;
  for (i=0; i<nimages; i++)
    trace += gsl_matrix_get(R,i,i);
  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,R,R,0.0,P);
  gsl_matrix_free(R);

  double trace2 = 0;
  for (i=0; i<nimages; i++) {
    trace2 += gsl_matrix_get(P,i,i);
  }
  gsl_matrix_free(P);

  (*xtrace) = trace;
  double df = (trace*trace) / trace2;
  return df;
}

/*
** general linear regression, 2nd level
*/
void GLM2(VImage *src,gsl_matrix *X,gsl_vector *contrast,int *permflag,int *permtable,int *signtable,int signswitch,VImage dest)
{
  int i,j,ip,b,r,c;
  double t=0,z=0,sum=0,tsigma=0,d=0,err=0,var=0,u=0,tiny=TINY;

  int ncols   = VImageNColumns(src[0]);
  int nrows   = VImageNRows(src[0]);
  int nslices = VImageNBands(src[0]);
  int nimages = X->size1;
  VFillImage(dest,VAllBands,0);

  /* permutation or sign switching applied to design matrix */
  gsl_matrix *XP = gsl_matrix_calloc(X->size1,X->size2);
  gsl_matrix_memcpy(XP,X);

  if (signswitch >= 0) {   /* sign switching  (one-sample test) */
    for (i=0; i<X->size1; i++) {
      if (signtable[i] < 0) {
	u = gsl_matrix_get(X,i,signswitch);
	gsl_matrix_set(XP,i,signswitch,-u);
      }
    }
  }

  else {       /*  permutations  */
    for (i=0; i<X->size1; i++) {
      ip = permtable[i];
      for (j=0; j<X->size2; j++) {
	if (permflag[j] > 0) {  /* do not permute columns containing nuisance covariates */
	  u = gsl_matrix_get(X,i,j);
	  gsl_matrix_set(XP,ip,j,u);
	}
      }
    }
  }

  /* pseudo inverse */
  gsl_matrix *XInv = PseudoInv(XP,NULL);

  /* get DoF and variance */
  double trace = 0;
  double df = DoF(XP,XInv,&trace);
  double sigma = EstimateVariance(XInv,contrast);
  if (df < tiny) VError(" Zero degrees of freedom");
  if (sigma < tiny) VError(" No variance in design/contrast");

  /* main loop */
  gsl_set_error_handler_off();
  gsl_vector *y    = gsl_vector_alloc (nimages);
  gsl_vector *yz   = gsl_vector_alloc (nimages);
  gsl_vector *beta = gsl_vector_alloc (X->size2);

  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {

	int nonzero=0;
	for (j=0; j<nimages; j++) {
	  y->data[j] = VGetPixel(src[j],b,r,c);
	  if (fabs(y->data[j]) > 0) nonzero++;
	}
	if (nonzero < nimages-2) continue;

	/* compute beta's */
	gsl_blas_dgemv(CblasNoTrans,1.0,XInv,y,0.0,beta);

	/* residuals */
	gsl_blas_dgemv(CblasNoTrans,1.0,XP,beta,0.0,yz);
	err = 0;
	for (j=0; j<nimages; j++) {
	  d = y->data[j] - yz->data[j];
	  err += d*d;
	}

	/* get z-value */
	t = z = sum = 0;
	gsl_blas_ddot (beta,contrast,&sum);
	if (fabs(sum) < tiny) continue;
	var = err / trace;
	tsigma = sqrt(var) * sigma;
	if (tsigma > tiny) {
	  t = sum / tsigma;
	  z = t2z(t,df);
	  if (sum < 0) z = -z;
	}
	VPixel(dest,b,r,c,VFloat) = z;
      }
    }
  }

  gsl_vector_free(y);
  gsl_vector_free(yz);
  gsl_vector_free(beta);
  gsl_matrix_free(XP);
  gsl_matrix_free(XInv);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <glib.h>
#include <omp.h>

// libgsl0-dev
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <sys/time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_double.h>

double pi(int n, int batch) {
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    int end = (n / batch);
    double sum = 0;

    int i, j;
    // GRand * grand = g_rand_new ();
    //g_rand_double(grand);
    #pragma omp parallel for default(shared)  firstprivate(r) private(j)  reduction (+:sum)
    for (i = 0; i < end; i++) {
        for (j=0; j<batch; j++) {
            double a = gsl_rng_uniform_pos (r);
            double b = gsl_rng_uniform_pos (r);
            double c = (a * a) + (b * b);

            if (c <= 1) {
                sum += c;
            }
        }
    }
    return 8 * sum / n;
}

int main () {
    printf("%0.15f\n", pi(100000000, 1000));
}


text:
// Authors: David Blei (blei@cs.princeton.edu)
//          Sean Gerrish (sgerrish@cs.princeton.edu)
//
// Copyright 2011 Sean Gerrish and David Blei
// All Rights Reserved.
//
// See the README for this package for details about modifying or
// distributing this software.

#ifndef LDA_H
#define LDA_H

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_lambert.h>
#include <gsl/gsl_rng.h>

#include "param.h"
#include "data.h"
#include "gsl-wrappers.h"

/*
 * functions for posterior inference in the latent dirichlet
 * allocation model.
 *
 */

#define LDA_INFERENCE_CONVERGED 1e-8
#define LDA_SEED_INIT 1
#define LDA_INIT_SMOOTH 1.0
#define LDA_EM_CONVERGED 5e-5
#define LDA_USE_VAR_BAYES 0
#define LDA_TOPIC_DIR_PARAM 0.001

// lda model

typedef struct lda {
    int ntopics;         // number of topics
    int nterms;          // vocabulary size
    gsl_matrix* topics;  // each column is a topic (V X K)
    gsl_vector* alpha;   // dirichlet parameters
} lda;

// lda posterior

typedef struct lda_post {
    doc_t* doc;          // document associated to this posterior
    lda* model;          // lda model
    gsl_matrix* phi;     // variational mult parameters (nterms x K)
    gsl_matrix* log_phi; // convenient for computation (nterms x K)
    gsl_vector* gamma;   // variational dirichlet parameters (K)
    gsl_vector* lhood;   // a K+1 vector, sums to the lhood bound
    gsl_vector* doc_weight;  // Not owned by this structure.
    gsl_vector* renormalized_doc_weight;  // Not owned by this structure.
} lda_post;

// lda sufficient statistics

typedef struct lda_suff_stats {
  gsl_matrix* topics_ss;
} lda_suff_stats;

// new lda model and suff stats

lda* new_lda_model(int ntopics, int nterms);
void free_lda_model(lda* m);
lda_suff_stats* new_lda_suff_stats(lda* model);
void reset_lda_suff_stats(lda_suff_stats* ss);
lda_post* new_lda_post(int ntopics, int max_length);
void free_lda_post(lda_post* p);
void initialize_lda_ss_from_data(corpus_t* data, lda_suff_stats* ss);

// posterior inference

double fit_lda_post(int doc_number, int time,
		    lda_post* p, lda_seq* var,
		    gsl_matrix* g,
		    gsl_matrix* g3,
		    gsl_matrix* g4,
		    gsl_matrix* g5);
void init_lda_post(lda_post* p);
void update_gamma(lda_post* p);
void update_phi(int doc_number, int time,
		lda_post* p, lda_seq* var,
		gsl_matrix* g);
void update_phi_dim(int doc_number, int time,
		    lda_post* p, lda_seq* var,
		    gsl_matrix* g);
void update_phi_fixed(int doc_number, int time,
		      lda_post* p, lda_seq* var,
		      gsl_matrix* g3_matrix,
		      gsl_matrix* g4_matrix,
		      gsl_matrix* g5_matrix);
void update_phi_multiple(int doc_number, int time,
			 lda_post* p, lda_seq* var,
			 gsl_matrix* g);

// compute the likelihood bound

double compute_lda_lhood(lda_post* p);

// EM algorithm

double lda_e_step(lda* model, corpus_t* data, lda_suff_stats* ss);
double lda_m_step(lda* model, lda_suff_stats* ss);
void lda_em(lda* model,
	    lda_suff_stats* ss,
	    corpus_t* data,
	    int max_iter,
	    char* outname);

// reading and writing

lda_suff_stats* read_lda_suff_stats(char* filename, int ntopics, int nterms);
void write_lda(lda* model, char* name);
void write_lda_suff_stats(lda_suff_stats* ss, char* name);
lda* read_lda(int ntopics, int nterms, char* name);

void initialize_lda_ss_from_random(corpus_t* data, lda_suff_stats* ss);

#endif

max_stars_repo_name: iwangjian/dtm-lab

max_issues_repo_name: iwangjian/topic-extractor

max_forks_repo_name: iwangjian/topic-extractor


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

/* Driver routine for odeiv2. This is a wrapper for low level GSL
   functions that allows a simple interface to step, control and
   evolve layers.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_machine.h>

static gsl_odeiv2_driver *
driver_alloc (const gsl_odeiv2_system * sys, const double hstart,
              const gsl_odeiv2_step_type * T)
{
  /* Allocates and initializes an ODE driver system. Step and evolve
     objects are allocated here, but control object is allocated in
     another function.
   */

  gsl_odeiv2_driver *state =
    (gsl_odeiv2_driver *) malloc (sizeof (gsl_odeiv2_driver));

  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for driver state",
                      GSL_ENOMEM);
    }

  if (sys == NULL)
    {
      GSL_ERROR_NULL ("gsl_odeiv2_system must be defined", GSL_EINVAL);
    }

  {
    const size_t dim = sys->dimension;

    if (dim == 0)
      {
        GSL_ERROR_NULL
          ("gsl_odeiv2_system dimension must be a positive integer",
           GSL_EINVAL);
      }

    state->sys = sys;

    state->s = gsl_odeiv2_step_alloc (T, dim);

    if (state->s == NULL)
      {
        free (state);
        GSL_ERROR_NULL ("failed to allocate step object", GSL_ENOMEM);
      }

    state->e = gsl_odeiv2_evolve_alloc (dim);
  }

  if (state->e == NULL)
    {
      gsl_odeiv2_step_free (state->s);
      free (state);
      GSL_ERROR_NULL ("failed to allocate evolve object", GSL_ENOMEM);
    }

  if (hstart > 0.0 || hstart < 0.0)
    {
      state->h = hstart;
    }
  else
    {
      GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL);
    }

  state->h = hstart;
  state->hmin = 0.0;
  state->hmax = GSL_DBL_MAX;
  state->nmax = 0;
  state->n = 0;
  state->c = NULL;

  return state;
}

int
gsl_odeiv2_driver_set_hmin (gsl_odeiv2_driver * d, const double hmin)
{
  /* Sets minimum allowed step size fabs(hmin) for driver. It is
     required that hmin <= fabs(h) <= hmax. */

  if ((fabs (hmin) > fabs (d->h)) || (fabs (hmin) > d->hmax))
    {
      GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
    }

  d->hmin = fabs (hmin);

  return GSL_SUCCESS;
}

int
gsl_odeiv2_driver_set_hmax (gsl_odeiv2_driver * d, const double hmax)
{
  /* Sets maximum allowed step size fabs(hmax) for driver. It is
     required that hmin <= fabs(h) <= hmax. */

  if ((fabs (hmax) < fabs (d->h)) || (fabs (hmax) < d->hmin))
    {
      GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
    }

  if (hmax > 0.0 || hmax < 0.0)
    {
      d->hmax = fabs (hmax);
    }
  else
    {
      GSL_ERROR_NULL ("invalid hmax", GSL_EINVAL);
    }

  return GSL_SUCCESS;
}

int
gsl_odeiv2_driver_set_nmax (gsl_odeiv2_driver * d,
                            const unsigned long int nmax)
{
  /* Sets maximum number of allowed steps (nmax) for driver */

  d->nmax = nmax;

  return GSL_SUCCESS;
}

gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_y_new (const gsl_odeiv2_system * sys,
                               const gsl_odeiv2_step_type * T,
                               const double hstart,
                               const double epsabs, const double epsrel)
{
  /* Initializes an ODE driver system with control object of type y_new. */

  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);

  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
    }

  if (epsabs >= 0.0 && epsrel >= 0.0)
    {
      state->c = gsl_odeiv2_control_y_new (epsabs, epsrel);

      if (state->c == NULL)
        {
          gsl_odeiv2_driver_free (state);
          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
        }
    }
  else
    {
      gsl_odeiv2_driver_free (state);
      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
    }

  /* Distribute pointer to driver object */

  gsl_odeiv2_step_set_driver (state->s, state);
  gsl_odeiv2_evolve_set_driver (state->e, state);
  gsl_odeiv2_control_set_driver (state->c, state);

  return state;
}

gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_yp_new (const gsl_odeiv2_system * sys,
                                const gsl_odeiv2_step_type * T,
                                const double hstart,
                                const double epsabs, const double epsrel)
{
  /* Initializes an ODE driver system with control object of type yp_new. */

  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);

  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
    }

  if (epsabs >= 0.0 && epsrel >= 0.0)
    {
      state->c = gsl_odeiv2_control_yp_new (epsabs, epsrel);

      if (state->c == NULL)
        {
          gsl_odeiv2_driver_free (state);
          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
        }
    }
  else
    {
      gsl_odeiv2_driver_free (state);
      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
    }

  /* Distribute pointer to driver object */

  gsl_odeiv2_step_set_driver (state->s, state);
  gsl_odeiv2_evolve_set_driver (state->e, state);
  gsl_odeiv2_control_set_driver (state->c, state);

  return state;
}

gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_standard_new (const gsl_odeiv2_system * sys,
                                      const gsl_odeiv2_step_type * T,
                                      const double hstart,
                                      const double epsabs,
                                      const double epsrel, const double a_y,
                                      const double a_dydt)
{
  /* Initializes an ODE driver system with control object of type
     standard_new.
   */

  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);

  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
    }

  if (epsabs >= 0.0 && epsrel >= 0.0)
    {
      state->c =
        gsl_odeiv2_control_standard_new (epsabs, epsrel, a_y, a_dydt);

      if (state->c == NULL)
        {
          gsl_odeiv2_driver_free (state);
          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
        }
    }
  else
    {
      gsl_odeiv2_driver_free (state);
      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
    }

  /* Distribute pointer to driver object */

  gsl_odeiv2_step_set_driver (state->s, state);
  gsl_odeiv2_evolve_set_driver (state->e, state);
  gsl_odeiv2_control_set_driver (state->c, state);

  return state;
}

gsl_odeiv2_driver *
gsl_odeiv2_driver_alloc_scaled_new (const gsl_odeiv2_system * sys,
                                    const gsl_odeiv2_step_type * T,
                                    const double hstart,
                                    const double epsabs, const double epsrel,
                                    const double a_y, const double a_dydt,
                                    const double scale_abs[])
{
  /* Initializes an ODE driver system with control object of type
     scaled_new.
   */

  gsl_odeiv2_driver *state = driver_alloc (sys, hstart, T);

  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate driver object", GSL_ENOMEM);
    }

  if (epsabs >= 0.0 && epsrel >= 0.0)
    {
      state->c = gsl_odeiv2_control_scaled_new (epsabs, epsrel, a_y, a_dydt,
                                                scale_abs, sys->dimension);

      if (state->c == NULL)
        {
          gsl_odeiv2_driver_free (state);
          GSL_ERROR_NULL ("failed to allocate control object", GSL_ENOMEM);
        }
    }
  else
    {
      gsl_odeiv2_driver_free (state);
      GSL_ERROR_NULL ("epsabs and epsrel must be positive", GSL_EINVAL);
    }

  /* Distribute pointer to driver object */

  gsl_odeiv2_step_set_driver (state->s, state);
  gsl_odeiv2_evolve_set_driver (state->e, state);
  gsl_odeiv2_control_set_driver (state->c, state);

  return state;
}

int
gsl_odeiv2_driver_apply (gsl_odeiv2_driver * d, double *t,
                         const double t1, double y[])
{
  /* Main driver function that evolves the system from t to t1. In
     beginning vector y contains the values of dependent variables at
     t. This function returns values at t=t1 in y. In case of
     unrecoverable error, y and t contains the values after the last
     successful step.
   */

  int sign = 0;
  d->n = 0;

  /* Determine integration direction sign */

  if (d->h > 0.0)
    {
      sign = 1;
    }
  else
    {
      sign = -1;
    }

  /* Check that t, t1 and step direction are sensible */

  if (sign * (t1 - *t) < 0.0)
    {
      GSL_ERROR_NULL
        ("integration limits and/or step direction not consistent",
         GSL_EINVAL);
    }

  /* Evolution loop */

  while (sign * (t1 - *t) > 0.0)
    {
      int s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, d->sys,
                                       t, t1, &(d->h), y);

      if (s != GSL_SUCCESS)
        {
          return s;
        }

      /* Check for maximum allowed steps */

      if ((d->nmax > 0) && (d->n > d->nmax))
        {
          return GSL_EMAXITER;
        }

      /* Set step size if maximum size is exceeded */

      if (fabs (d->h) > d->hmax)
        {
          d->h = sign * d->hmax;
        }

      /* Check for too small step size */

      if (fabs (d->h) < d->hmin)
        {
          return GSL_ENOPROG;
        }

      d->n++;
    }

  return GSL_SUCCESS;
}

int
gsl_odeiv2_driver_apply_fixed_step (gsl_odeiv2_driver * d, double *t,
                                    const double h, const unsigned long int n,
                                    double y[])
{
  /* Alternative driver function that evolves the system from t using
   * n steps of size h. In the beginning vector y contains the values
   * of dependent variables at t. This function returns values at t =
   * t + n * h in y. In case of an unrecoverable error, y and t
   * contains the values after the last successful step.
   */

  unsigned long int i;
  d->n = 0;

  /* Evolution loop */

  for (i = 0; i < n; i++)
    {
      int s = gsl_odeiv2_evolve_apply_fixed_step (d->e, d->c, d->s, d->sys,
                                                  t, h, y);

      if (s != GSL_SUCCESS)
        {
          return s;
        }

      d->n++;
    }

  return GSL_SUCCESS;
}

int
gsl_odeiv2_driver_reset (gsl_odeiv2_driver * d)
{
  /* Reset the driver. Resets evolve and step objects. */

  {
    int s = gsl_odeiv2_evolve_reset (d->e);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = gsl_odeiv2_step_reset (d->s);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  return GSL_SUCCESS;
}

int
gsl_odeiv2_driver_reset_hstart (gsl_odeiv2_driver * d, const double hstart)
{
  /* Resets current driver and sets initial step size to hstart */

  gsl_odeiv2_driver_reset (d);

  if ((d->hmin > fabs (hstart)) || (fabs (hstart) > d->hmax))
    {
      GSL_ERROR_NULL ("hmin <= fabs(h) <= hmax required", GSL_EINVAL);
    }

  if (hstart > 0.0 || hstart < 0.0)
    {
      d->h = hstart;
    }
  else
    {
      GSL_ERROR_NULL ("invalid hstart", GSL_EINVAL);
    }

  return GSL_SUCCESS;
}

void
gsl_odeiv2_driver_free (gsl_odeiv2_driver * state)
{
  if (state->c != NULL)
    {
      gsl_odeiv2_control_free (state->c);
    }

  gsl_odeiv2_evolve_free (state->e);
  gsl_odeiv2_step_free (state->s);
  free (state);
}

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
#ifdef HAVE_CONFIG_H
#include "../config.h"
#endif /* HAVE_CONFIG_H */

#ifndef DO_WITH_GSL

/* This test is not performed */
int main()
{
    return 0;
}

#else

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>

int main(void)
{
  FILE* out_file; /* output stream*/
  int i;
  double cut_off;
  double sigma;
  const int sample_no=100000;   /* sample numbers */
  const int bin_no=100;      /* number of bins in histogram */
  gsl_histogram* my_hist;
  gsl_rng* my_rng;

  out_file=stdout;
  my_hist=gsl_histogram_calloc_uniform(bin_no,-10,10);
  my_rng=gsl_rng_alloc(gsl_rng_default);

  sigma=1;
  cut_off=-4;

  fprintf(out_file,"plot '-' using 1:3 title '1st',");
  fprintf(out_file,"'-' using 1:3 title '2nd',");
  fprintf(out_file,"'-' using 1:3 title '3rd'\n");
  while (cut_off<=4)
    {
      fprintf(out_file,"# cut_off %f, sigma %f\n",cut_off,sigma);
      gsl_histogram_reset(my_hist);
      for (i=0;i<sample_no;i++)
	gsl_histogram_increment(my_hist,gsl_ran_gaussian_tail(my_rng,cut_off,sigma));
      gsl_histogram_fprintf(out_file,my_hist,"%f","%f");
      fprintf(out_file,"e\n");
      cut_off+=4;
    }
  return 0;
}

#endif


text:
#ifndef __GSL_VECTOR_H__
#define __GSL_VECTOR_H__

#include <gsl/vector/gsl_vector_double.h>

#endif /* __GSL_VECTOR_H__ */

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/**
 *
 * @file pclarft_blgtrd.c
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
  @date 2011-05-15
 * @generated c Tue Jan  7 11:45:13 2014
 *
 **/
#include "common.h"
#include <lapacke.h>

#define V(m)     &(V[(m)])
#define TAU(m)   &(TAU[(m)])
#define T(m)   &(T[(m)])
/***************************************************************************/
/**
 *  Parallel compute T2 from bulgechasing of Symetric matrix - static scheduling
 *  Lower case is supported
 **/
/***************************************************************************/
void plasma_pclarft_blgtrd(plasma_context_t *plasma)
{
    int my_core_id = PLASMA_RANK;
    int cores_num  = plasma->world_size;
    /*===========================*/
    int N, NB,Vblksiz;
    PLASMA_Complex32_t *V;
    PLASMA_Complex32_t *T;
    PLASMA_Complex32_t *TAU;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    /*===========================
     *   local variables
     *===========================*/
    int LDT, LDV;
    int Vm, Vn, mt, nt;
    int myrow, mycol, blkj, blki;
    int firstrow;
    int blkid,vpos,taupos,tpos;
    int blkpercore,blkcnt, myid;

    plasma_unpack_args_8(N, NB, Vblksiz, V, T, TAU, sequence, request);
    if (sequence->status != PLASMA_SUCCESS)
        return;

    /* Quick return */
    if (N == 0){
        return;
    }
    if (NB == 0){
        return;
    }
    if (NB == 1){
        return;
    }

    findVTsiz(N, NB, Vblksiz, &blkcnt, &LDV);
    blkpercore = blkcnt/cores_num;
    blkpercore = blkpercore==0 ? 1:blkpercore;
    LDT        = Vblksiz;
    LDV        = NB+Vblksiz-1;

    /*========================================
     * compute the T's in parallel.
     * The Ts are independent so each core pick
     * a T and compute it. The loop is based on
     * the version 113 of the pcunmqr_blgtrd.c
     * which go over the losange block_column
     * by block column. but it is not important
     * here the order because Ts are independent.
     * ========================================
    */
    nt  = plasma_ceildiv((N-1),Vblksiz);
    for (blkj=nt-1; blkj>=0; blkj--) {
        /* the index of the first row on the top of block (blkj) */
        firstrow = blkj * Vblksiz + 1;
        /*find the number of tile for this block */
        if( blkj == nt-1 )
            mt = plasma_ceildiv( N -  firstrow,    NB);
        else
            mt = plasma_ceildiv( N - (firstrow+1), NB);
        /*loop over the tiles find the size of the Vs and apply it */
        for (blki=mt; blki>0; blki--) {
            /*calculate the size of each losange of Vs= (Vm,Vn)*/
            myrow     = firstrow + (mt-blki)*NB;
            mycol     = blkj*Vblksiz;
            Vm = min( NB+Vblksiz-1, N-myrow);
            if( ( blkj == nt-1 ) && ( blki == mt ) ){
                Vn = min (Vblksiz, Vm);
            } else {
                Vn = min (Vblksiz, Vm-1);
            }
            /*calculate the pointer to the Vs and the Ts.
             * Note that Vs and Ts have special storage done
             * by the bulgechasing function*/
            findVTpos(N,NB,Vblksiz,mycol,myrow, &vpos, &taupos, &tpos, &blkid);
            myid = blkid/blkpercore;
            if( my_core_id==(myid%cores_num) ){
                if( ( Vm > 0 ) && ( Vn > 0 ) ){
                    LAPACKE_clarft_work(LAPACK_COL_MAJOR,
                                  lapack_const(PlasmaForward),
                                  lapack_const(PlasmaColumnwise),
                                  Vm, Vn, V(vpos), LDV, TAU(taupos), T(tpos), LDT);
                }
            }
        }
    }
}
#undef V
#undef TAU
#undef T
/***************************************************************************/


text:
/*
 * Licensed to the OpenAirInterface (OAI) Software Alliance under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The OpenAirInterface Software Alliance licenses this file to You under
 * the OAI Public License, Version 1.1  (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.openairinterface.org/?page_id=698
 *
 * 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.
 *-------------------------------------------------------------------------------
 * For more information about the OpenAirInterface (OAI) Software Alliance:
 *      contact@openairinterface.org
 */

#include <string.h>
#include <math.h>
#include <unistd.h>
#include <stdint.h>
#include <stdio.h>
#include <time.h>
#include <cblas.h>

#include "SIMULATION/TOOLS/defs.h"
#include "SIMULATION/RF/defs.h"
#include "PHY/types.h"
#include "PHY/defs.h"
#include "PHY/extern.h"
#include "oaisim_config.h"

#ifdef OPENAIR2
#include "LAYER2/MAC/defs.h"
#include "LAYER2/MAC/extern.h"
#include "UTIL/LOG/log_if.h"
#include "UTIL/LOG/log_extern.h"
#include "RRC/LITE/extern.h"
#include "PHY_INTERFACE/extern.h"
#include "UTIL/OCG/OCG.h"
#include "UTIL/OMG/omg.h"
#include "UTIL/OPT/opt.h" // to test OPT
#endif

#include "SCHED/defs.h"
#include "SCHED/extern.h"

#include "oaisim.h"

#define PI 3.1416
#define Am 20

#define MCL (-70) /*minimum coupling loss (MCL) in dB*/
//double sinr[NUMBER_OF_eNB_MAX][2*25];
/*
extern double sinr_bler_map[MCS_COUNT][2][16];
extern double sinr_bler_map_up[MCS_COUNT][2][16];
double SINRpost_eff[301];
extern double MI_map_4qam[3][162];
extern double MI_map_16qam[3][197];
extern double MI_map_64qam[3][227];
*/
// Extract the positions of UE and ENB from the mobility model

void extract_position (node_list* input_node_list, node_desc_t **node_data, int nb_nodes)
{

  int i;

  for (i=0; i<nb_nodes; i++) {
    if ((input_node_list != NULL) &&  (node_data[i] != NULL)) {

      node_data[i]->x = input_node_list->node->x_pos;

      if (node_data[i]->x <0.0)
        node_data[i]->x = 0.0;

      node_data[i]->y = input_node_list->node->y_pos;

      if (node_data[i]->y <0.0)
        node_data[i]->y = 0.0;

      LOG_D(OCM, "extract_position: added node_data %d with position X: %f and Y: %f \n", i,input_node_list->node->x_pos, input_node_list->node->y_pos );
      input_node_list = input_node_list->next;
    } else {
      LOG_E(OCM, "extract_position: Null pointer!!!\n");
      //exit(-1);
    }
  }
}
void extract_position_fixed_enb  (node_desc_t **node_data, int nb_nodes, frame_t frame)
{
  int i;

  for (i=0; i<nb_nodes; i++) {
    if (i==0) {
      node_data[i]->x = 0;
      node_data[i]->y = 500;
    } else if (i == 1 ) {
      node_data[i]->x = 866;//
      node_data[i]->y = 1000;
    } else if (i == 2 ) {
      node_data[i]->x = 866;
      node_data[i]->y = 0;
    }
  }
}

void extract_position_fixed_ue  (node_desc_t **node_data, int nb_nodes, frame_t frame)
{
  int i;

  if(frame<50)
    for (i=0; i<nb_nodes; i++) {
      if (i==0) {
        node_data[i]->x = 2050;
        node_data[i]->y = 1500;
      } else {
        node_data[i]->x = 2150;
        node_data[i]->y = 1500;
      }
    }
  else {
    for (i=0; i<nb_nodes; i++) {
      if (i==0) {
        node_data[i]->x = 1856 - (frame - 49);
        // if(node_data[i]->x > 2106)
        //  node_data[i]->x = 2106;
        node_data[i]->y = 1813 + (frame - 49);
        // if(node_data[i]->y < 1563)
        //  node_data[i]->y = 1563;
        // if( node_data[i]->x == 2106)
        //   node_data[i]->x = 2106 - (frame - 49);
      } else {
        node_data[i]->x = 2106 - (frame - 49);
        // if(node_data[i]->x < 1856)
        //  node_data[i]->x = 1856;
        node_data[i]->y = 1563 + (frame - 49);
        // if(node_data[i]->y < 1813)
        //  node_data[i]->y = 1813;
      }
    }
  }

}

void init_ue(node_desc_t  *ue_data, UE_Antenna ue_ant)  //changed from node_struct
{

  ue_data->n_sectors = 1;
  ue_data->phi_rad = 2 * PI;
  ue_data->ant_gain_dBi = ue_ant.antenna_gain_dBi;
  ue_data->tx_power_dBm = ue_ant.tx_power_dBm;
  ue_data->rx_noise_level = ue_ant.rx_noise_level_dB; //value in db

}

void init_enb(node_desc_t  *enb_data, eNB_Antenna enb_ant)  //changed from node_struct
{

  int i;
  double sect_angle[3]= {0,2*PI/3,4*PI/3};

  enb_data->n_sectors = enb_ant.number_of_sectors;

  for (i=0; i<enb_data->n_sectors; i++)
    enb_data->alpha_rad[i] = sect_angle[i]; //enb_ant.alpha_rad[i];

  enb_data->phi_rad = enb_ant.beam_width_dB;
  enb_data->ant_gain_dBi = enb_ant.antenna_gain_dBi;
  enb_data->tx_power_dBm = enb_ant.tx_power_dBm;
  enb_data->rx_noise_level = enb_ant.rx_noise_level_dB;

}

void calc_path_loss(node_desc_t* enb_data, node_desc_t* ue_data, channel_desc_t *ch_desc, Environment_System_Config env_desc, double **Shad_Fad)
{

  double dist;
  double path_loss;
  double gain_max;
  double gain_sec[3];
  double alpha, theta;

  int count;

  dist = sqrt(pow((enb_data->x - ue_data->x), 2) + pow((enb_data->y - ue_data->y), 2));

  path_loss = env_desc.fading.free_space_model_parameters.pathloss_0_dB -
              10*env_desc.fading.free_space_model_parameters.pathloss_exponent * log10(dist/1000);
  LOG_D(OCM,"dist %f, Path loss %f\n",dist,ch_desc->path_loss_dB);

  /* Calculating the angle in the range -pi to pi from the slope */
  alpha = atan2((ue_data->y - enb_data->y),(ue_data->x - enb_data->x));

  if (alpha < 0)
    alpha += 2*PI;

  //printf("angle in radians is %lf\n", ue_data[UE_id]->alpha_rad[eNB_id]);
  ch_desc->aoa = alpha;
  ch_desc->random_aoa = 0;

  if (enb_data->n_sectors==1) //assume omnidirectional antenna
    gain_max = 0;
  else {
    gain_max = -1000;

    for(count = 0; count < enb_data->n_sectors; count++) {
      theta = enb_data->alpha_rad[count] - alpha;
      /* gain = -min(Am , 12 * (theta/theta_3dB)^2) */
      gain_sec[count] = -(Am < (12 * pow((theta/enb_data->phi_rad),2)) ? Am : (12 * pow((theta/enb_data->phi_rad),2)));

      if (gain_sec[count]>gain_max)  //take the sector that we are closest too (or where the gain is maximum)
        gain_max = gain_sec[count];
    }
  }

  path_loss += enb_data->ant_gain_dBi + gain_max + ue_data->ant_gain_dBi;

  if (Shad_Fad!=NULL)
    path_loss += Shad_Fad[(int)ue_data->x][(int)ue_data->y];

  ch_desc->path_loss_dB = MCL < path_loss ?  MCL : path_loss;
  //LOG_D(OCM,"x_coordinate\t%f\t,y_coordinate\t%f\t, path_loss %f\n",ue_data->x,ue_data->y,ch_desc->path_loss_dB);
}

void init_snr(channel_desc_t* eNB2UE, node_desc_t *enb_data, node_desc_t *ue_data, double* sinr_dB, double* N0, uint8_t transmission_mode, uint16_t q, uint8_t dl_power_off, uint16_t nb_rb)
{

  double thermal_noise,abs_channel,channelx, channely,channelx_i, channely_i ;
  int count;
  int aarx,aatx;
  uint8_t qq;

  /* Thermal noise is calculated using 10log10(K*T*B) K = Boltzmann's constant T = room temperature B = bandwidth*/
  thermal_noise = -174 + 10*log10(15000); //per RE; value in dBm

  //for (aarx=0; aarx<eNB2UE->nb_rx; aarx++)
  *N0 = thermal_noise + ue_data->rx_noise_level;

  LOG_D(OCM,"Path loss %lf, noise (N0) %lf, signal %lf, snr %lf\n",
        eNB2UE->path_loss_dB,
        thermal_noise + ue_data->rx_noise_level,
        enb_data->tx_power_dBm + eNB2UE->path_loss_dB,
        enb_data->tx_power_dBm + eNB2UE->path_loss_dB - (thermal_noise + ue_data->rx_noise_level));

  if(transmission_mode == 5 && dl_power_off==1)
    transmission_mode = 6;

  switch(transmission_mode) {
  case 1:

    //printf ("coupling factor is %lf\n", coupling);
    for (count = 0; count < (12 * nb_rb); count++) {
      sinr_dB[count] = enb_data->tx_power_dBm
                       + eNB2UE->path_loss_dB
                       - (thermal_noise + ue_data->rx_noise_level)
                       + 10 * log10 (pow(eNB2UE->chF[0][count].x, 2)
                                     + pow(eNB2UE->chF[0][count].y, 2));
      //printf("sinr_dB[%d]: %1f\n",count,sinr_dB[count]);
      //printf("Dl_link SNR for res. block %d is %lf\n", count, sinr[eNB_id][count]);
    }

    break;

  case 2:

    for (count = 0; count < (12 * nb_rb); count++) {
      abs_channel=0;

      for (aarx=0; aarx<eNB2UE->nb_rx; aarx++) {
        for (aatx=0; aatx<eNB2UE->nb_tx; aatx++) {
          abs_channel += (pow(eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x, 2) + pow(eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y, 2));
        }
      }

      sinr_dB[count] = enb_data->tx_power_dBm
                       + eNB2UE->path_loss_dB
                       - (thermal_noise + ue_data->rx_noise_level)
                       + 10 * log10 (abs_channel/2);
      // printf("sinr_dB[%d]: %1f\n",count,sinr_dB[count]);
    }

    break;

  case 5:
    for (count = 0; count < (12 * nb_rb); count++) {
      channelx=0;
      channely=0;
      channelx_i=0;
      channely_i=0;
      qq = (q>>(((count/12)>>2)<<1))&3;
      //printf("pmi_alloc %d: rb %d, pmi %d\n",q,count/12,qq);

      //      qq = q;
      for (aarx=0; aarx<eNB2UE->nb_rx; aarx++) {
        for (aatx=0; aatx<eNB2UE->nb_tx; aatx++) {
          switch(qq) {
          case 0:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channelx_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channelx_i -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely_i -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            }

            break;

          case 1:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channelx_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channelx_i += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely_i += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            }

            break;

          case 2:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channelx_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channely += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channelx_i += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channely_i -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;

            }

            break;

          case 3:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channelx_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely_i = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channely -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channelx_i -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channely_i += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
            }

            break;

          default:
            LOG_E(EMU,"Problem in SINR Calculation for TM5 \n");
            break;

          }//switch(q)

        }//aatx
      }//aarx

      /*  sinr_dB[count] = enb_data->tx_power_dBm
         + eNB2UE->path_loss_dB
        - (thermal_noise + ue_data->rx_noise_level)
        + 10 * log10 ((pow(channelx,2) + pow(channely,2))/2) - 10 * log10 ((pow(channelx_i,2) + pow(channely_i,2))/2);
      */
      sinr_dB[count] = enb_data->tx_power_dBm
                       + eNB2UE->path_loss_dB
                       - (thermal_noise + ue_data->rx_noise_level)
                       + 10 * log10 ((pow(channelx,2) + pow(channely,2))) - 10 * log10 ((pow(channelx_i,2) + pow(channely_i,2))) - 3; // 3dB is subtracted as the tx_power_dBm is to be adjusted on per user basis
      // printf("sinr_dB[%d]: %1f\n",count,sinr_dB[count]);
    }

    break;

  case 6:
    for (count = 0; count < (12 * nb_rb); count++) {
      channelx=0;
      channely=0;
      qq = (q>>(((count/12)>>2)<<1))&3;
      //printf("pmi_alloc %d: rb %d, pmi %d\n",q,count/12,qq);

      //      qq = q;
      for (aarx=0; aarx<eNB2UE->nb_rx; aarx++) {
        for (aatx=0; aatx<eNB2UE->nb_tx; aatx++) {
          switch(qq) {
          case 0:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            }

            break;

          case 1:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            }

            break;

          case 2:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channely += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
            }

            break;

          case 3:
            if (channelx==0 || channely==0) {
              channelx = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
              channely = eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
            } else {
              channelx += eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].y;
              channely -= eNB2UE->chF[aarx+(aatx*eNB2UE->nb_rx)][count].x;
            }

            break;

          default:
            LOG_E(EMU,"Problem in SINR Calculation for TM6 \n");
            break;

          }//switch(q)

        }//aatx
      }//aarx

      sinr_dB[count] = enb_data->tx_power_dBm
                       + eNB2UE->path_loss_dB
                       - (thermal_noise + ue_data->rx_noise_level)
                       + 10 * log10 ((pow(channelx,2) + pow(channely,2))/2);

      // printf("sinr_dB[%d]: %1f\n",count,sinr_dB[count]);
    }

    break;

  default:
    LOG_E(EMU,"Problem in SINR Initialization in sinr_sim.c\n");
    break;
  }//switch
}//function ends

void calculate_sinr(channel_desc_t* eNB2UE, node_desc_t *enb_data, node_desc_t *ue_data, double *sinr_dB, uint16_t nb_rb)
{

  double sir, thermal_noise;
  short count;

  /* Thermal noise is calculated using 10log10(K*T*B) K = Boltzmann's constant T = room temperature B = bandwidth */
  thermal_noise = -174 + 10*log10(15000); //per RE, value in dBm

  for (count = 0; count < 12 * nb_rb; count++) {
    sir = enb_data->tx_power_dBm
          + eNB2UE->path_loss_dB
          - (thermal_noise + ue_data->rx_noise_level)
          + 10 * log10 (pow(eNB2UE->chF[0][count].x, 2)
                        + pow(eNB2UE->chF[0][count].y, 2));

    if (sir > 0)
      sinr_dB[count] -= sir;

    //printf("*****sinr% lf \n",sinr_dB[count]);
  }
}
void get_beta_map()
{
  char *file_path = NULL;
  //int table_len = 0;
  int t,u;
  int mcs = 0;
  char *sinr_bler;
  char buffer[1000];
  FILE *fp;
  double perf_array[13];

  file_path = (char*) malloc(512);

  for (mcs = 0; mcs < MCS_COUNT; mcs++) {
    snprintf( file_path, 512, "%s/SIMULATION/LTE_PHY/BLER_SIMULATIONS/AWGN/AWGN_results/bler_tx1_chan18_nrx1_mcs%d.csv", getenv("OPENAIR1_DIR"), mcs );
    fp = fopen(file_path,"r");

    if (fp == NULL) {
      LOG_E(OCM,"ERROR: Unable to open the file %s! Exitng.\n", file_path);
      exit(-1);
    }

    // else {
    if (fgets (buffer, 1000, fp) != NULL) {
      if (fgets (buffer, 1000, fp) != NULL) {
        table_length[mcs] = 0;

        while (!feof (fp)) {
          u = 0;
          sinr_bler = strtok (buffer, ";");

          while (sinr_bler != NULL) {
            perf_array[u] = atof (sinr_bler);
            u++;
            sinr_bler = strtok (NULL, ";");
          }

          if ((perf_array[4] / perf_array[5]) < 1) {
            sinr_bler_map[mcs][0][table_length[mcs]] = perf_array[0];
            sinr_bler_map[mcs][1][table_length[mcs]] = (perf_array[4] / perf_array[5]);
            table_length[mcs]++;

            if (table_length[mcs]>MCS_TABLE_LENGTH_MAX) {
              LOG_E(OCM,"Error reading MCS table. Increase MCS_TABLE_LENGTH_MAX (mcs %d)!\n",mcs);
              exit(-1);
            }
          }

          if (fgets (buffer, 1000, fp) != NULL) {
          }
        }
      }
    }

    fclose(fp);
    //   }
    LOG_D(OCM,"Print the table for mcs %d\n",mcs);

    for (t = 0; t<table_length[mcs]; t++)
      LOG_D(OCM,"%lf  %lf \n ",sinr_bler_map[mcs][0][t],sinr_bler_map[mcs][1][t]);
  }

  free(file_path);
}

//this function reads and stores the Mutual information tables for the MIESM abstraction.
void get_MIESM_param()
{
  char *file_path = NULL;
  char buffer[10000];
  FILE *fp;
  int qam[3] = {4,16,64};
  int q,cnt;
  char *result = NULL;
  int table_len=0;
  int t;
  file_path = (char*) malloc(512);

  for (q=0; q<3; q++) {
    sprintf(file_path,"%s/SIMU/USER/files/MI_%dqam.csv",getenv("OPENAIR_TARGETS"),qam[q]);
    fp = fopen(file_path,"r");

    if (fp == NULL) {
      printf("ERROR: Unable to open the file %s\n", file_path);
      exit(-1);
    } else {
      cnt=-1;

      switch(qam[q]) {
      case 4:
        while (!feof(fp)) {
          table_len =0;
          cnt++;

          if (fgets(buffer, 10000, fp) != NULL) {
            result = strtok (buffer, ",");

            while (result != NULL) {
              MI_map_4qam[cnt][table_len] = atof (result);
              result = strtok (NULL, ",");
              table_len++;
            }
          }
        }

        for (t = 0; t < 162; t++) {
          // MI_map_4Qam[0][t] = pow(10,0.1*(MI_map_4Qam[0][t]));
          LOG_D(OCM, "MIESM 4QAM Table: %lf  %lf  %1f\n ",MI_map_4qam[0][t],MI_map_4qam[1][t], MI_map_4qam[2][t]);
          // printf("MIESM 4QAM Table: %lf  %lf  %1f\n ",MI_map_4qam[0][t],MI_map_4qam[1][t], MI_map_4qam[2][t]);
        }

        break;

      case 16:
        while (!feof(fp)) {
          table_len =0;
          cnt++;

          if (fgets (buffer, 10000, fp) != NULL) {
            result = strtok (buffer, ",");

            while (result != NULL) {
              MI_map_16qam[cnt][table_len] = atof (result);
              result = strtok (NULL, ",");
              table_len++;
            }
          }
        }

        for (t = 0; t < 197; t++) {
          // MI_map_16Qam[0][t] = pow(10,0.1*(MI_map_16Qam[0][t]));
          LOG_D(OCM, "MIESM 16 QAM Table: %lf  %lf  %1f\n ",MI_map_16qam[0][t],MI_map_16qam[1][t], MI_map_16qam[2][t]);
          // printf("MIESM 16 QAM Table: %lf  %lf  %1f\n ",MI_map_16qam[0][t],MI_map_16qam[1][t], MI_map_16qam[2][t]);
        }

        break;

      case 64:
        while (!feof(fp)) {
          table_len=0;
          cnt++;

          if(cnt==3)
            break;

          if (fgets (buffer, 10000, fp) != NULL) {
            result = strtok(buffer, ",");

            while (result != NULL) {
              MI_map_64qam[cnt][table_len]= atof(result);
              result = strtok(NULL, ",");
              table_len++;
            }
          }
        }

        for (t = 0; t < 227; t++) {
          //MI_map_64Qam[0][t] = pow(10,0.1*(MI_map_64Qam[0][t]));
          LOG_D(OCM, "MIESM 64QAM Table: %lf  %lf  %1f\n ",MI_map_64qam[0][t],MI_map_64qam[1][t], MI_map_64qam[2][t]);
          // printf("MIESM 64QAM Table: %lf  %lf  %1f\n ",MI_map_64qam[0][t],MI_map_64qam[1][t], MI_map_64qam[2][t]);
        }

        break;

      default:
        LOG_E(EMU,"Error, bad input, quitting\n");
        break;
      }

    }

    fclose(fp);
  }

  free(file_path);
}

max_stars_repo_name: davidraditya/OAI-Powder

max_issues_repo_name: davidraditya/OAI-Powder

max_forks_repo_name: davidraditya/OAI-Powder