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 <stdlib.h>
#include <time.h>
#include <cblas.h>
#include <execinfo.h>

//<<PAD>>//
//#include <mpi.h>
//#include "UTIL/FIFO/pad_list.h"
#include "discrete_event_generator.h"
#include "threadpool.h"
#include <pthread.h>
#include "oaisim_functions.h"
//<<PAD>>//

#include "SIMULATION/RF/defs.h"
#include "PHY/types.h"
#include "PHY/defs.h"
#include "PHY/vars.h"
#include "MAC_INTERFACE/vars.h"

//#ifdef OPENAIR2
#include "LAYER2/MAC/defs.h"
#include "LAYER2/MAC/vars.h"
#include "RRC/LITE/vars.h"
#include "PHY_INTERFACE/vars.h"
//#endif

#include "ARCH/CBMIMO1/DEVICE_DRIVER/vars.h"

#ifdef IFFT_FPGA
//#include "PHY/LTE_REFSIG/mod_table.h"
#endif //IFFT_FPGA

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

#include "oaisim.h"
#include "oaisim_config.h"
#include "UTIL/OCG/OCG_extern.h"
#include "cor_SF_sim.h"
#include "UTIL/OMG/omg_constants.h"

//#include "UTIL/LOG/vcd_signal_dumper.h"

#define RF

//#define DEBUG_SIM

#define MCS_COUNT 24//added for PHY abstraction
#define N_TRIALS 1

/*
  DCI0_5MHz_TDD0_t          UL_alloc_pdu;
  DCI1A_5MHz_TDD_1_6_t      CCCH_alloc_pdu;
  DCI2_5MHz_2A_L10PRB_TDD_t DLSCH_alloc_pdu1;
  DCI2_5MHz_2A_M10PRB_TDD_t DLSCH_alloc_pdu2;
 */

#define UL_RB_ALLOC computeRIV(lte_frame_parms->N_RB_UL,0,24)
#define CCCH_RB_ALLOC computeRIV(lte_frame_parms->N_RB_UL,0,3)
#define RA_RB_ALLOC computeRIV(lte_frame_parms->N_RB_UL,0,3)
#define DLSCH_RB_ALLOC 0x1fff

#define DECOR_DIST 100
#define SF_VAR 10

//constant for OAISIM soft realtime calibration
#define SF_DEVIATION_OFFSET_NS 100000 //= 0.1ms : should be as a number of UE
#define SLEEP_STEP_US       100 //  = 0.01ms could be adaptive, should be as a number of UE
#define K 2                  // averaging coefficient
#define TARGET_SF_TIME_NS 1000000       // 1ms = 1000000 ns

//#ifdef OPENAIR2
//uint16_t NODE_ID[1];
//uint8_t NB_INST = 2;
//#endif //OPENAIR2
extern int otg_times;
extern int for_times;
extern int if_times;
int for_main_times = 0;

frame_t frame=0;
char stats_buffer[16384];
channel_desc_t *eNB2UE[NUMBER_OF_eNB_MAX][NUMBER_OF_UE_MAX];
channel_desc_t *UE2eNB[NUMBER_OF_UE_MAX][NUMBER_OF_eNB_MAX];
Signal_buffers *signal_buffers_g;
//Added for PHY abstraction
node_desc_t *enb_data[NUMBER_OF_eNB_MAX];
node_desc_t *ue_data[NUMBER_OF_UE_MAX];
//double sinr_bler_map[MCS_COUNT][2][16];
//double sinr_bler_map_up[MCS_COUNT][2][16];
//extern double SINRpost_eff[301];
extern int mcsPost;
extern int  nrbPost;
extern int frbPost;
extern void kpi_gen();

extern uint16_t Nid_cell;
extern uint8_t target_dl_mcs;
extern uint8_t rate_adaptation_flag;
extern double snr_dB, sinr_dB;
extern uint8_t set_seed;
extern uint8_t cooperation_flag;          // for cooperative communication
extern uint8_t abstraction_flag, ethernet_flag;
extern uint8_t ue_connection_test;
extern int map1,map2;
extern double **ShaF;
// pointers signal buffers (s = transmit, r,r0 = receive)
extern double **s_re, **s_im, **r_re, **r_im, **r_re0, **r_im0;
extern Node_list ue_node_list;
extern Node_list enb_node_list;
extern int pdcp_period, omg_period;
extern LTE_DL_FRAME_PARMS *frame_parms;
// time calibration for soft realtime mode
extern struct timespec time_spec;
extern unsigned long time_last, time_now;
extern int td, td_avg, sleep_time_us;

int eMBMS_active = 0;

threadpool_t * pool;

#ifdef OPENAIR2
extern int pfd[2];
#endif

// this should reflect the channel models in openair1/SIMULATION/TOOLS/defs.h
mapping small_scale_names[] = {
  {"custom", custom},
  {"SCM_A", SCM_A},
  {"SCM_B", SCM_B},
  {"SCM_C", SCM_C},
  {"SCM_D", SCM_D},
  {"EPA", EPA},
  {"EVA", EVA},
  {"ETU", ETU},
  {"Rayleigh8", Rayleigh8},
  {"Rayleigh1", Rayleigh1},
  {"Rayleigh1_800", Rayleigh1_800},
  {"Rayleigh1_corr", Rayleigh1_corr},
  {"Rayleigh1_anticorr", Rayleigh1_anticorr},
  {"Rice8", Rice8},
  {"Rice1", Rice1},
  {"Rice1_corr", Rice1_corr},
  {"Rice1_anticorr", Rice1_anticorr},
  {"AWGN", AWGN},
  {NULL, -1}
};

//static void *sigh(void *arg);
void terminate(void);

void
help (void)
{
  printf
  ("Usage: oaisim -h -a -F -C tdd_config -V -R N_RB_DL -e -x transmission_mode -m target_dl_mcs -r(ate_adaptation) -n n_frames -s snr_dB -k ricean_factor -t max_delay -f forgetting factor -A channel_model -z cooperation_flag -u nb_local_ue -U UE mobility -b nb_local_enb -B eNB_mobility -M ethernet_flag -p nb_master -g multicast_group -l log_level -c ocg_enable -T traffic model -D multicast network device\n");

  printf ("-h provides this help message!\n");
  printf ("-a Activates PHY abstraction mode\n");
  printf ("-F Activates FDD transmission (TDD is default)\n");
  printf ("-C [0-6] Sets TDD configuration\n");
  printf ("-R [6,15,25,50,75,100] Sets N_RB_DL\n");
  printf ("-e Activates extended prefix mode\n");
  printf ("-m Gives a fixed DL mcs\n");
  printf ("-r Activates rate adaptation (DL for now)\n");
  printf ("-n Set the number of frames for the simulation\n");
  printf ("-s snr_dB set a fixed (average) SNR, this deactivates the openair channel model generator (OCM)\n");
  printf ("-S snir_dB set a fixed (average) SNIR, this deactivates the openair channel model generator (OCM)\n");
  printf ("-k Set the Ricean factor (linear)\n");
  printf ("-t Set the delay spread (microseconds)\n");
  printf ("-f Set the forgetting factor for time-variation\n");
  printf ("-A set the multipath channel simulation,  options are: SCM_A, SCM_B, SCM_C, SCM_D, EPA, EVA, ETU, Rayleigh8, Rayleigh1, Rayleigh1_corr,Rayleigh1_anticorr, Rice8,, Rice1, AWGN \n");
  printf ("-b Set the number of local eNB\n");
  printf ("-u Set the number of local UE\n");
  printf ("-M Set the machine ID for Ethernet-based emulation\n");
  printf ("-p Set the total number of machine in emulation - valid if M is set\n");
  printf ("-g Set multicast group ID (0,1,2,3) - valid if M is set\n");
  printf ("-l Set the global log level (8:trace, 7:debug, 6:info, 4:warn, 3:error) \n");
  printf
  ("-c [1,2,3,4] Activate the config generator (OCG) to process the scenario descriptor, or give the scenario manually: -c template_1.xml \n");
  printf ("-x Set the transmission mode (1,2,5,6 supported for now)\n");
  printf ("-z Set the cooperation flag (0 for no cooperation, 1 for delay diversity and 2 for distributed alamouti\n");
  printf ("-T activate the traffic generator: 0 for NONE, 1 for CBR, 2 for M2M, 3 for FPS Gaming, 4 for mix\n");
  printf ("-B Set the mobility model for eNB, options are: STATIC, RWP, RWALK, \n");
  printf ("-U Set the mobility model for UE, options are: STATIC, RWP, RWALK \n");
  printf ("-E Random number generator seed\n");
  printf ("-P enable protocol analyzer : 0 for wireshark interface, 1: for pcap , 2 : for tshark \n");
  printf ("-I Enable CLI interface (to connect use telnet localhost 1352)\n");
  printf ("-V Enable VCD dump, file = openair_vcd_dump.vcd\n");
  printf ("-G Enable background traffic \n");
  printf ("-O [mme ipv4 address] Enable MME mode\n");
  printf ("-Z Reserved\n");
}

#ifdef OPENAIR2
void omv_end (int pfd, Data_Flow_Unit omv_data);
int omv_write (int pfd,  Node_list enb_node_list, Node_list ue_node_list, Data_Flow_Unit omv_data);
#endif

//<<<< PAD >>>>//
#define PAD 1
//#define PAD_FINE 1
//#define PAD_SYNC 1
#define JOB_REQUEST_TAG 246
#define JOB_REPLY_TAG 369
#define FRAME_END 888
#define NO_JOBS_TAG 404
#define JOB_DIST_DEBUG 33

//Global Variables
int worker_number;
int frame_number = 1;
//<<<< PAD >>>>//

//<<<< DEG >>>>//
extern End_Of_Sim_Event end_event; //Could later be a list of condition_events
extern Event_List event_list;
//<<<< DEG >>>>//

extern Packet_OTG_List *otg_pdcp_buffer;

void run(int argc, char *argv[]);

#ifdef PAD
void pad_init()
{

  int UE_id, i;

  pool = threadpool_create(PAD);

  if (pool == NULL) {
    printf("ERROR threadpool allocation\n");
    return;
  }

  signal_buffers_g = malloc(NB_UE_INST * sizeof(Signal_buffers));

  if (abstraction_flag == 0) {
    for (UE_id = 0; UE_id < NB_UE_INST; UE_id++) {
      signal_buffers_g[UE_id].s_re = malloc(2*sizeof(double*));
      signal_buffers_g[UE_id].s_im = malloc(2*sizeof(double*));
      signal_buffers_g[UE_id].r_re = malloc(2*sizeof(double*));
      signal_buffers_g[UE_id].r_im = malloc(2*sizeof(double*));
      signal_buffers_g[UE_id].r_re0 = malloc(2*sizeof(double*));
      signal_buffers_g[UE_id].r_im0 = malloc(2*sizeof(double*));

      for (i=0; i<2; i++) {

        signal_buffers_g[UE_id].s_re[i] = malloc(FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        bzero(signal_buffers_g[UE_id].s_re[i],FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        signal_buffers_g[UE_id].s_im[i] = malloc(FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        bzero(signal_buffers_g[UE_id].s_im[i],FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        signal_buffers_g[UE_id].r_re[i] = malloc(FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        bzero(signal_buffers_g[UE_id].r_re[i],FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        signal_buffers_g[UE_id].r_im[i] = malloc(FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        bzero(signal_buffers_g[UE_id].r_im[i],FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        signal_buffers_g[UE_id].r_re0[i] = malloc(FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        bzero(signal_buffers_g[UE_id].r_re0[i],FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        signal_buffers_g[UE_id].r_im0[i] = malloc(FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
        bzero(signal_buffers_g[UE_id].r_im0[i],FRAME_LENGTH_COMPLEX_SAMPLES*sizeof(double));
      }
    }
  }
}

void pad_finalize()
{

  int ret, i;
  module_id_t UE_id;

  ret = threadpool_destroy(pool);

  if (ret)
    printf("ERROR threadpool destroy = %d\n", ret);

  if (abstraction_flag == 0) {

    for (UE_id = 0; UE_id < NB_UE_INST; UE_id++) {
      for (i = 0; i < 2; i++) {
        free(signal_buffers_g[UE_id].s_re[i]);
        free(signal_buffers_g[UE_id].s_im[i]);
        free(signal_buffers_g[UE_id].r_re[i]);
        free(signal_buffers_g[UE_id].r_im[i]);
      }

      free(signal_buffers_g[UE_id].s_re);
      free(signal_buffers_g[UE_id].s_im);
      free(signal_buffers_g[UE_id].r_re);
      free(signal_buffers_g[UE_id].r_im);
    }

    //free node by node here same pattern as below
  }

  free(signal_buffers_g);
}

void pad_inject_job(int eNB_flag, int nid, int frame, int next_slot, int last_slot, enum Job_type type, int ctime)
{

  int ret;
  Job_elt *job_elt;

  job_elt = malloc(sizeof(Job_elt));
  job_elt->next = NULL;
  (job_elt->job).eNB_flag = eNB_flag;
  (job_elt->job).nid = nid;
  (job_elt->job).frame = frame;
  (job_elt->job).next_slot = next_slot;
  (job_elt->job).last_slot = last_slot;
  (job_elt->job).type = type;
  (job_elt->job).ctime = ctime;

  ret = threadpool_add(pool, job_elt);

  if (ret) {
    printf("ERROR threadpool_add %d\n", ret);
    return;
  }
}

void pad_synchronize()
{
  pthread_mutex_lock(&(pool->sync_lock));

  while(pool->active > 0) {
    pthread_cond_wait(&(pool->sync_notify), &(pool->sync_lock));
  }

  pthread_mutex_unlock(&(pool->sync_lock));
}

#endif
//<<PAD(DEG_MAIN)>>//
int main (int argc, char *argv[])
{

  //Mobility *mobility_frame_10;
  //Application_Config *application_frame_20;

  //Here make modifications on the mobility and traffic new models
  //mob_frame_10 -> ...
  //application_frame_30 -> ...

  //schedule(ET_OMG, 10, NULL, mobility_frame_10);
  //schedule(ET_OTG, 15, NULL, application_frame_20);

  //event_list_display(&event_list);

  schedule_end_of_simulation(FRAME, 100);

  run(argc, argv);

  return 0;
}
//<<PAD>>//

//<<PAD(RUN)>>//
void run(int argc, char *argv[])
{

  int32_t i;
  module_id_t UE_id, eNB_id;
  Job_elt *job_elt;
  int ret;

  clock_t t;

  Event_elt *user_defined_event;
  Event event;

  // Framing variables
  int32_t slot, last_slot, next_slot;

  FILE *SINRpost;
  char SINRpost_fname[512];
  sprintf(SINRpost_fname,"postprocSINR.m");
  SINRpost = fopen(SINRpost_fname,"w");
  // variables/flags which are set by user on command-line
  double snr_direction,snr_step=1.0;//,sinr_direction;

  lte_subframe_t direction;
  char fname[64],vname[64];

#ifdef XFORMS
  // current status is that every UE has a DL scope for a SINGLE eNB (eNB_id=0)
  // at eNB 0, an UL scope for every UE
  FD_lte_phy_scope_ue  *form_ue[NUMBER_OF_UE_MAX];
  FD_lte_phy_scope_enb *form_enb[NUMBER_OF_UE_MAX];
  char title[255];
#endif

#ifdef PROC
  int node_id;
  int port,Process_Flag=0,wgt,Channel_Flag=0,temp;
#endif

  // uint8_t awgn_flag = 0;

#ifdef PRINT_STATS
  int len;
  FILE *UE_stats[NUMBER_OF_UE_MAX], *UE_stats_th[NUMBER_OF_UE_MAX], *eNB_stats, *eNB_avg_thr, *eNB_l2_stats;
  char UE_stats_filename[255];
  char UE_stats_th_filename[255];
  char eNB_stats_th_filename[255];
#endif

#ifdef SMBV
  uint8_t config_smbv = 0;
  char smbv_ip[16];
  strcpy(smbv_ip,DEFAULT_SMBV_IP);
#endif

#ifdef OPENAIR2
  Data_Flow_Unit omv_data;
#endif

  //time_t t0,t1;
  //clock_t start, stop;

  //double **s_re2[MAX_eNB+MAX_UE], **s_im2[MAX_eNB+MAX_UE], **r_re2[MAX_eNB+MAX_UE], **r_im2[MAX_eNB+MAX_UE], **r_re02, **r_im02;
  //double **r_re0_d[MAX_UE][MAX_eNB], **r_im0_d[MAX_UE][MAX_eNB], **r_re0_u[MAX_eNB][MAX_UE],**r_im0_u[MAX_eNB][MAX_UE];
  //default parameters

  //{
  /* INITIALIZATIONS */

  target_dl_mcs = 0;
  rate_adaptation_flag = 0;
  oai_emulation.info.n_frames = 0xffff;//1024;          //10;
  oai_emulation.info.n_frames_flag = 0;//fixme
  snr_dB = 30;
  cooperation_flag = 0;         // default value 0 for no cooperation, 1 for Delay diversity, 2 for Distributed Alamouti

  //Default values if not changed by the user in get_simulation_options();
  pdcp_period = 1;
  omg_period = 10;

  mRAL_init_default_values(); //Default values
  eRAL_init_default_values(); //Default values

  init_oai_emulation(); //Default values

  get_simulation_options(argc, argv); //Command-line options

  oaisim_config(); // config OMG and OCG, OPT, OTG, OLG

  //To fix eventual conflict on the value of n_frames
  if (oai_emulation.info.n_frames_flag) {
    schedule_end_of_simulation(FRAME, oai_emulation.info.n_frames);
  }

  VCD_SIGNAL_DUMPER_INIT(); // Initialize VCD LOG module

#ifdef OPENAIR2
  init_omv();
#endif

  check_and_adjust_params(); //Before this call, NB_UE_INST and NB_eNB_INST are not set correctly

  init_otg_pdcp_buffer();

#ifdef PRINT_STATS

  for (UE_id=0; UE_id<NB_UE_INST; UE_id++) {
    sprintf(UE_stats_filename,"UE_stats%d.txt",UE_id);
    UE_stats[UE_id] = fopen (UE_stats_filename, "w");
  }

  eNB_stats = fopen ("eNB_stats.txt", "w");
  printf ("UE_stats=%p, eNB_stats=%p\n", UE_stats, eNB_stats);

  eNB_avg_thr = fopen ("eNB_stats_th.txt", "w");

#endif

  LOG_I(EMU,"total number of UE %d (local %d, remote %d) mobility %s \n", NB_UE_INST,oai_emulation.info.nb_ue_local,oai_emulation.info.nb_ue_remote,
        oai_emulation.topology_config.mobility.eNB_mobility.eNB_mobility_type.selected_option);
  LOG_I(EMU,"Total number of eNB %d (local %d, remote %d) mobility %s \n", NB_eNB_INST,oai_emulation.info.nb_enb_local,oai_emulation.info.nb_enb_remote,
        oai_emulation.topology_config.mobility.UE_mobility.UE_mobility_type.selected_option);
  LOG_I(OCM,"Running with frame_type %d, Nid_cell %d, N_RB_DL %d, EP %d, mode %d, target dl_mcs %d, rate adaptation %d, nframes %d, abstraction %d, channel %s\n",
        oai_emulation.info.frame_type, Nid_cell, oai_emulation.info.N_RB_DL, oai_emulation.info.extended_prefix_flag, oai_emulation.info.transmission_mode,target_dl_mcs,rate_adaptation_flag,
        oai_emulation.info.n_frames,abstraction_flag,oai_emulation.environment_system_config.fading.small_scale.selected_option);

  init_seed(set_seed);

  init_openair1();

  init_openair2();

  init_ocm();

#ifdef XFORMS
  init_xforms();
#endif

  printf ("before L2 init: Nid_cell %d\n", PHY_vars_eNB_g[0]->lte_frame_parms.Nid_cell);
  printf ("before L2 init: frame_type %d,tdd_config %d\n",
          PHY_vars_eNB_g[0]->lte_frame_parms.frame_type,
          PHY_vars_eNB_g[0]->lte_frame_parms.tdd_config);

  init_time();

#ifdef PAD
  pad_init();
#endif

  if (ue_connection_test == 1) {
    snr_direction = -snr_step;
    snr_dB=20;
    sinr_dB=-20;
  }

  frame = 0;
  slot = 0;

  LOG_I(EMU,">>>>>>>>>>>>>>>>>>>>>>>>>>> OAIEMU initialization done <<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");

  printf ("after init: Nid_cell %d\n", PHY_vars_eNB_g[0]->lte_frame_parms.Nid_cell);
  printf ("after init: frame_type %d,tdd_config %d\n",
          PHY_vars_eNB_g[0]->lte_frame_parms.frame_type,
          PHY_vars_eNB_g[0]->lte_frame_parms.tdd_config);

  t = clock();

  while (!end_of_simulation()) {

    last_slot = (slot - 1)%20;

    if (last_slot <0)
      last_slot+=20;

    next_slot = (slot + 1)%20;

    oai_emulation.info.time_ms = frame * 10 + (next_slot>>1);
    oai_emulation.info.frame = frame;

    if (slot == 0) { //Frame's Prologue

      //Run the aperiodic user-defined events
      while ((user_defined_event = event_list_get_head(&event_list)) != NULL) {

        event = user_defined_event->event;

        if (event.frame == frame) {
          switch (event.type) {
          case ET_OMG:
            update_omg_model(event.key, event.value); //implement it with assigning the new values to that of oai_emulation & second thing is to ensure mob model is always read from oai_emulation
            user_defined_event = event_list_remove_head(&event_list);
            break;

          case ET_OTG:
            update_otg_model(event.key, event.value);
            user_defined_event = event_list_remove_head(&event_list);
            break;
          }
        } else {
          break;
        }
      }

      //Comment (handle cooperation flag) deleted here. Look at oaisim.c to see it
      if (ue_connection_test==1) {
        if ((frame%20) == 0) {
          snr_dB += snr_direction;
          sinr_dB -= snr_direction;
        }

        if (snr_dB == -20) {
          snr_direction=snr_step;
        } else if (snr_dB==20) {
          snr_direction=-snr_step;
        }
      }

      update_omg(); // frequency is defined in the omg_global params configurable by the user

      update_omg_ocm();

#ifdef OPENAIR2

      // check if pipe is still open
      if ((oai_emulation.info.omv_enabled == 1) ) {
        omv_write(pfd[1], enb_node_list, ue_node_list, omv_data);
      }

#endif

#ifdef DEBUG_OMG

      if ((((int) oai_emulation.info.time_s) % 100) == 0) {
        for (UE_id = oai_emulation.info.first_ue_local; UE_id < (oai_emulation.info.first_ue_local + oai_emulation.info.nb_ue_local); UE_id++) {
          get_node_position (UE, UE_id);
        }
      }

#endif
      update_ocm();
    }

    direction = subframe_select(frame_parms,next_slot>>1);

    if((next_slot %2) ==0)
      clear_eNB_transport_info(oai_emulation.info.nb_enb_local);

    for (eNB_id=oai_emulation.info.first_enb_local;
         (eNB_id<(oai_emulation.info.first_enb_local+oai_emulation.info.nb_enb_local)) && (oai_emulation.info.cli_start_enb[eNB_id]==1);
         eNB_id++) {
      for_main_times += 1;
      //printf ("debug: Nid_cell %d\n", PHY_vars_eNB_g[eNB_id]->lte_frame_parms.Nid_cell);
      //printf ("debug: frame_type %d,tdd_config %d\n", PHY_vars_eNB_g[eNB_id]->lte_frame_parms.frame_type,PHY_vars_eNB_g[eNB_id]->lte_frame_parms.tdd_config);
      LOG_D(EMU,"PHY procedures eNB %d for frame %d, slot %d (subframe TX %d, RX %d) TDD %d/%d Nid_cell %d\n",
            eNB_id, frame, slot, next_slot >> 1,last_slot>>1,
            PHY_vars_eNB_g[eNB_id]->lte_frame_parms.frame_type,
            PHY_vars_eNB_g[eNB_id]->lte_frame_parms.tdd_config,PHY_vars_eNB_g[eNB_id]->lte_frame_parms.Nid_cell);

      //Appliation
#ifdef PAD_FINE
      pad_inject_job(1, eNB_id, frame, next_slot, last_slot, JT_OTG, oai_emulation.info.time_ms);
#else
      update_otg_eNB(eNB_id, oai_emulation.info.time_ms);
#endif

      //Access layer
      if (frame % pdcp_period == 0) {
#ifdef PAD_FINE
        pad_inject_job(1, eNB_id, frame, next_slot, last_slot, JT_PDCP, oai_emulation.info.time_ms);
#else
        pdcp_run(frame, 1, 0, eNB_id);//PHY_vars_eNB_g[eNB_id]->Mod_id
#endif
      }

      //Phy/Mac layer
#ifdef PAD_FINE
      pad_inject_job(1, eNB_id, frame, next_slot, last_slot, JT_PHY_MAC, oai_emulation.info.time_ms);
#else
      phy_procedures_eNB_lte (last_slot, next_slot, PHY_vars_eNB_g[eNB_id], abstraction_flag, no_relay, NULL);
#endif

#ifdef PRINT_STATS

      if(last_slot==9 && frame%10==0)
        if(eNB_avg_thr)
          fprintf(eNB_avg_thr,"%d %d\n",PHY_vars_eNB_g[eNB_id]->frame,(PHY_vars_eNB_g[eNB_id]->total_system_throughput)/((PHY_vars_eNB_g[eNB_id]->frame+1)*10));

      if (eNB_stats) {
        len = dump_eNB_stats(PHY_vars_eNB_g[eNB_id], stats_buffer, 0);
        rewind (eNB_stats);
        fwrite (stats_buffer, 1, len, eNB_stats);
        fflush(eNB_stats);
      }

#ifdef OPENAIR2

      if (eNB_l2_stats) {
        len = dump_eNB_l2_stats (stats_buffer, 0);
        rewind (eNB_l2_stats);
        fwrite (stats_buffer, 1, len, eNB_l2_stats);
        fflush(eNB_l2_stats);
      }

#endif
#endif
    }

#ifdef PAD_SYNC

    if ((direction == SF_DL) || ((direction == SF_S) && (next_slot%2==0)) )
      pad_synchronize();

#endif

    // Call ETHERNET emulation here
    //emu_transport (frame, last_slot, next_slot, direction, oai_emulation.info.frame_type, ethernet_flag);

    if ((next_slot % 2) == 0)
      clear_UE_transport_info (oai_emulation.info.nb_ue_local);

    for (UE_id = oai_emulation.info.first_ue_local;
         (UE_id < (oai_emulation.info.first_ue_local+oai_emulation.info.nb_ue_local)) && (oai_emulation.info.cli_start_ue[UE_id]==1);
         UE_id++)
      if (frame >= (UE_id * 20)) {    // activate UE only after 20*UE_id frames so that different UEs turn on separately

        LOG_D(EMU,"PHY procedures UE %d for frame %d, slot %d (subframe TX %d, RX %d)\n",
              UE_id, frame, slot, next_slot >> 1,last_slot>>1);

        if (PHY_vars_UE_g[UE_id]->UE_mode[0] != NOT_SYNCHED) {
          if (frame>0) {
            PHY_vars_UE_g[UE_id]->frame = frame;

            //Application UE
#ifdef PAD_FINE
            pad_inject_job(0, UE_id, frame, next_slot, last_slot, JT_OTG, oai_emulation.info.time_ms);
#else
            update_otg_UE(UE_id + NB_eNB_INST, oai_emulation.info.time_ms);
#endif

            //Access layer UE
            if (frame % pdcp_period == 0) {
#ifdef PAD_FINE
              pad_inject_job(0, UE_id, frame, next_slot, last_slot, JT_PDCP, oai_emulation.info.time_ms);
#else
              pdcp_run(frame, 0, UE_id, 0);
#endif
            }

            //Phy/Mac layer UE
#ifdef PAD_FINE
            pad_inject_job(0, UE_id, frame, next_slot, last_slot, JT_PHY_MAC, oai_emulation.info.time_ms);
#else
            phy_procedures_UE_lte (last_slot, next_slot, PHY_vars_UE_g[UE_id], 0, abstraction_flag, normal_txrx, no_relay, NULL);
            ue_data[UE_id]->tx_power_dBm = PHY_vars_UE_g[UE_id]->tx_power_dBm;
#endif
          }
        } else {
          if (abstraction_flag==1) {
            LOG_E(EMU, "sync not supported in abstraction mode (UE%d,mode%d)\n", UE_id, PHY_vars_UE_g[UE_id]->UE_mode[0]);
            exit(-1);
          }

          if ((frame>0) && (last_slot == (LTE_SLOTS_PER_FRAME-2))) {
#ifdef PAD_FINE
            pad_inject_job(0, UE_id, frame, next_slot, last_slot, JT_INIT_SYNC, oai_emulation.info.time_ms);
#else
            initial_sync(PHY_vars_UE_g[UE_id],normal_txrx);
#endif
            /* LONG write output comment DELETED here */
          }
        }

#ifdef PRINT_STATS

        if(last_slot==2 && frame%10==0)
          if (UE_stats_th[UE_id])
            fprintf(UE_stats_th[UE_id],"%d %d\n",frame, PHY_vars_UE_g[UE_id]->bitrate[0]/1000);

        if (UE_stats[UE_id]) {
          len = dump_ue_stats (PHY_vars_UE_g[UE_id], stats_buffer, 0, normal_txrx, 0);
          rewind (UE_stats[UE_id]);
          fwrite (stats_buffer, 1, len, UE_stats[UE_id]);
          fflush(UE_stats[UE_id]);
        }

#endif
      }

#ifdef PAD_SYNC

    if ((direction == SF_UL) || ((direction == SF_S) && (next_slot%2==1)) )
      pad_synchronize();

#endif

    emu_transport (frame, last_slot, next_slot,direction, oai_emulation.info.frame_type, ethernet_flag);

    if ((direction  == SF_DL)|| (frame_parms->frame_type==0)) {
      for (UE_id=0; UE_id<NB_UE_INST; UE_id++) {

#ifdef PAD
        pad_inject_job(0, UE_id, frame, next_slot, last_slot, JT_DL, oai_emulation.info.time_ms);
#else
        do_DL_sig(r_re0,r_im0,r_re,r_im,s_re,s_im,eNB2UE,enb_data,ue_data,next_slot,abstraction_flag,frame_parms,UE_id);
#endif
      }
    }

    if ((direction  == SF_UL)|| (frame_parms->frame_type==0)) { //if ((subframe<2) || (subframe>4))
      do_UL_sig(r_re0,r_im0,r_re,r_im,s_re,s_im,UE2eNB,enb_data,ue_data,next_slot,abstraction_flag,frame_parms,frame);

      /*
            int ccc;
            fprintf(SINRpost,"SINRdb For eNB New Subframe : \n ");
            for(ccc = 0 ; ccc<301; ccc++)
            {
              fprintf(SINRpost,"_ %f ", SINRpost_eff[ccc]);
            }
            fprintf(SINRpost,"SINRdb For eNB : %f \n ", SINRpost_eff[ccc]);
            */
    }

    if ((direction == SF_S)) {//it must be a special subframe
      if (next_slot%2==0) {//DL part
        for (UE_id=0; UE_id<NB_UE_INST; UE_id++) {
#ifdef PAD
          pad_inject_job(0, UE_id, frame, next_slot, last_slot, JT_DL, oai_emulation.info.time_ms);
#else
          do_DL_sig(r_re0,r_im0,r_re,r_im,s_re,s_im,eNB2UE,enb_data,ue_data,next_slot,abstraction_flag,frame_parms,UE_id);
#endif
        }

        /*
              for (aarx=0;aarx<UE2eNB[1][0]->nb_rx;aarx++)
              for (aatx=0;aatx<UE2eNB[1][0]->nb_tx;aatx++)
              for (k=0;k<UE2eNB[1][0]->channel_length;k++)
              printf("SB(%d,%d,%d)->(%f,%f)\n",k,aarx,aatx,UE2eNB[1][0]->ch[aarx+(aatx*UE2eNB[1][0]->nb_rx)][k].r,UE2eNB[1][0]->ch[aarx+(aatx*UE2eNB[1][0]->nb_rx)][k].i);
            */
      } else { // UL part
        /*#ifdef PAD
        pthread_mutex_lock(&(pool->sync_lock));
        while(pool->active != 0) {
        pthread_cond_wait(&(pool->sync_notify), &(pool->sync_lock));
        }
        pthread_mutex_unlock(&(pool->sync_lock));
        #endif*/
        do_UL_sig(r_re0,r_im0,r_re,r_im,s_re,s_im,UE2eNB,enb_data,ue_data,next_slot,abstraction_flag,frame_parms,frame);
        /*
                int ccc;
                fprintf(SINRpost,"SINRdb For eNB New Subframe : \n ");
                for(ccc = 0 ; ccc<301; ccc++)
                {
                  fprintf(SINRpost,"_ %f ", SINRpost_eff[ccc]);
                }
                fprintf(SINRpost,"SINRdb For eNB : %f \n ", SINRpost_eff[ccc]);
                */
      }
    }

    if ((last_slot == 1) && (frame == 0)
        && (abstraction_flag == 0) && (oai_emulation.info.n_frames == 1)) {

      write_output ("dlchan0.m", "dlch0",
                    &(PHY_vars_UE_g[0]->lte_ue_common_vars.common_vars_rx_data_per_thread[subframe&0x1].dl_ch_estimates[0][0][0]),
                    (6 * (PHY_vars_UE_g[0]->lte_frame_parms.ofdm_symbol_size)), 1, 1);
      write_output ("dlchan1.m", "dlch1",
                    &(PHY_vars_UE_g[0]->lte_ue_common_vars.common_vars_rx_data_per_thread[subframe&0x1].dl_ch_estimates[1][0][0]),
                    (6 * (PHY_vars_UE_g[0]->lte_frame_parms.ofdm_symbol_size)), 1, 1);
      write_output ("dlchan2.m", "dlch2",
                    &(PHY_vars_UE_g[0]->lte_ue_common_vars.common_vars_rx_data_per_thread[subframe&0x1].dl_ch_estimates[2][0][0]),
                    (6 * (PHY_vars_UE_g[0]->lte_frame_parms.ofdm_symbol_size)), 1, 1);
      write_output ("pbch_rxF_comp0.m", "pbch_comp0",
                    PHY_vars_UE_g[0]->lte_ue_pbch_vars[0]->rxdataF_comp[0], 6 * 12 * 4, 1, 1);
      write_output ("pbch_rxF_llr.m", "pbch_llr",
                    PHY_vars_UE_g[0]->lte_ue_pbch_vars[0]->llr, (frame_parms->Ncp == 0) ? 1920 : 1728, 1, 4);
    }

    if (next_slot %2 == 0) {
      clock_gettime (CLOCK_REALTIME, &time_spec);
      time_last = time_now;
      time_now = (unsigned long) time_spec.tv_nsec;
      td = (int) (time_now - time_last);

      if (td>0) {
        td_avg = (int)(((K*(long)td) + (((1<<3)-K)*((long)td_avg)))>>3); // in us
        LOG_T(EMU,"sleep frame %d, time_now %ldus,time_last %ldus,average time difference %ldns, CURRENT TIME DIFF %dus, avgerage difference from the target %dus\n",
              frame, time_now,time_last,td_avg, td/1000,(td_avg-TARGET_SF_TIME_NS)/1000);
      }

      if (td_avg<(TARGET_SF_TIME_NS - SF_DEVIATION_OFFSET_NS)) {
        sleep_time_us += SLEEP_STEP_US;
        LOG_D(EMU,"Faster than realtime increase the avg sleep time for %d us, frame %d\n",
              sleep_time_us,frame);
        // LOG_D(EMU,"Faster than realtime increase the avg sleep time for %d us, frame %d, time_now %ldus,time_last %ldus,average time difference %ldns, CURRENT TIME DIFF %dus, avgerage difference from the target %dus\n",    sleep_time_us,frame, time_now,time_last,td_avg, td/1000,(td_avg-TARGET_SF_TIME_NS)/1000);
      } else if (td_avg > (TARGET_SF_TIME_NS + SF_DEVIATION_OFFSET_NS)) {
        sleep_time_us-= SLEEP_STEP_US;
        LOG_D(EMU,"Slower than realtime reduce the avg sleep time for %d us, frame %d, time_now\n",
              sleep_time_us,frame);
        //LOG_T(EMU,"Slower than realtime reduce the avg sleep time for %d us, frame %d, time_now %ldus,time_last %ldus,average time difference %ldns, CURRENT TIME DIFF %dus, avgerage difference from the target %dus\n",     sleep_time_us,frame, time_now,time_last,td_avg, td/1000,(td_avg-TARGET_SF_TIME_NS)/1000);
      }
    } // end if next_slot%2

    slot++;

    if (slot == 20) { //Frame's Epilogue
      frame++;
      slot = 0;

      // if n_frames not set by the user or is greater than max num frame then set adjust the frame counter
      if ( (oai_emulation.info.n_frames_flag == 0) || (oai_emulation.info.n_frames >= 0xffff) ) {
        frame %=(oai_emulation.info.n_frames-1);
      }

      oai_emulation.info.time_s += 0.01;

      if ((frame>=1)&&(frame<=9)&&(abstraction_flag==0)) {
        write_output("UEtxsig0.m","txs0", PHY_vars_UE_g[0]->lte_ue_common_vars.txdata[0],FRAME_LENGTH_COMPLEX_SAMPLES,1,1);
        sprintf(fname,"eNBtxsig%d.m",frame);
        sprintf(vname,"txs%d",frame);
        write_output(fname,vname, PHY_vars_eNB_g[0]->lte_eNB_common_vars.txdata[0][0],FRAME_LENGTH_COMPLEX_SAMPLES,1,1);
        write_output("eNBtxsigF0.m","txsF0",PHY_vars_eNB_g[0]->lte_eNB_common_vars.txdataF[0][0],PHY_vars_eNB_g[0]->lte_frame_parms.symbols_per_tti*PHY_vars_eNB_g[0]->lte_frame_parms.ofdm_symbol_size,1,1);

        write_output("UErxsig0.m","rxs0", PHY_vars_UE_g[0]->lte_ue_common_vars.rxdata[0],FRAME_LENGTH_COMPLEX_SAMPLES,1,1);
        write_output("eNBrxsig0.m","rxs0", PHY_vars_eNB_g[0]->lte_eNB_common_vars.rxdata[0][0],FRAME_LENGTH_COMPLEX_SAMPLES,1,1);
      }

#ifdef XFORMS
      do_xforms();
#endif

      // calibrate at the end of each frame if there is some time  left
      if((sleep_time_us > 0)&& (ethernet_flag ==0)) {
        LOG_I(EMU,"[TIMING] Adjust average frame duration, sleep for %d us\n",sleep_time_us);
        usleep(sleep_time_us);
        sleep_time_us=0; // reset the timer, could be done per n SF
      }

#ifdef SMBV

      if ((frame == config_frames[0]) || (frame == config_frames[1]) || (frame == config_frames[2]) || (frame == config_frames[3])) {
        smbv_frame_cnt++;
      }

#endif
    }
  }

  t = clock() - t;
  printf("rrc Duration of the simulation: %f seconds\n",((float)t)/CLOCKS_PER_SEC);

  fclose(SINRpost);
  LOG_I(EMU,">>>>>>>>>>>>>>>>>>>>>>>>>>> OAIEMU Ending <<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");

  free(otg_pdcp_buffer);

#ifdef SMBV

  if (config_smbv) {
    smbv_send_config (smbv_fname,smbv_ip);
  }

#endif

  //Perform KPI measurements
  if (oai_emulation.info.otg_enabled==1)
    kpi_gen();

#ifdef PAD
  pad_finalize();
#endif

  // relase all rx state
  if (ethernet_flag == 1) {
    emu_transport_release ();
  }

  if (abstraction_flag == 0) {
    /*
           #ifdef IFFT_FPGA
           free(txdataF2[0]);
           free(txdataF2[1]);
           free(txdataF2);
           free(txdata[0]);
           free(txdata[1]);
           free(txdata);
           #endif
         */

    for (i = 0; i < 2; i++) {
      free (s_re[i]);
      free (s_im[i]);
      free (r_re[i]);
      free (r_im[i]);
    }

    free (s_re);
    free (s_im);
    free (r_re);
    free (r_im);
    lte_sync_time_free ();
  }

  //  pthread_join(sigth, NULL);

  // added for PHY abstraction
  if (oai_emulation.info.ocm_enabled == 1) {
    for (eNB_id = 0; eNB_id < NUMBER_OF_eNB_MAX; eNB_id++)
      free(enb_data[eNB_id]);

    for (UE_id = 0; UE_id < NUMBER_OF_UE_MAX; UE_id++)
      free(ue_data[UE_id]);
  } //End of PHY abstraction changes

#ifdef OPENAIR2
  mac_top_cleanup();
#endif

#ifdef PRINT_STATS

  for(UE_id=0; UE_id<NB_UE_INST; UE_id++) {
    if (UE_stats[UE_id])
      fclose (UE_stats[UE_id]);

    if(UE_stats_th[UE_id])
      fclose (UE_stats_th[UE_id]);
  }

  if (eNB_stats)
    fclose (eNB_stats);

  if (eNB_avg_thr)
    fclose (eNB_avg_thr);

  if (eNB_l2_stats)
    fclose (eNB_l2_stats);

#endif

  // stop OMG
  stop_mobility_generator(oai_emulation.info.omg_model_ue);//omg_param_list.mobility_type
#ifdef OPENAIR2

  if (oai_emulation.info.omv_enabled == 1)
    omv_end(pfd[1],omv_data);

#endif

  if ((oai_emulation.info.ocm_enabled == 1) && (ethernet_flag == 0) && (ShaF != NULL))
    destroyMat(ShaF,map1, map2);

  if (opt_enabled == 1 )
    terminate_opt();

  if (oai_emulation.info.cli_enabled)
    cli_server_cleanup();

  //bring oai if down
  terminate();
  logClean();
  VCD_SIGNAL_DUMPER_CLOSE();
  //printf("FOR MAIN TIMES = %d &&&& OTG TIMES = %d <-> FOR TIMES = %d <-> IF TIMES = %d\n", for_main_times, otg_times, for_times, if_times);

}
//<<PAD>>//

void terminate(void)
{
  int i;
  char interfaceName[8];

  for (i=0; i < NUMBER_OF_eNB_MAX+NUMBER_OF_UE_MAX; i++)
    if (oai_emulation.info.oai_ifup[i]==1) {
      sprintf(interfaceName, "oai%d", i);
      bringInterfaceUp(interfaceName,0);
    }
}

#ifdef OPENAIR2
int omv_write (int pfd,  Node_list enb_node_list, Node_list ue_node_list, Data_Flow_Unit omv_data)
{
  int i,j;
  omv_data.end=0;

  //omv_data.total_num_nodes = NB_UE_INST + NB_eNB_INST;
  for (i=0; i<NB_eNB_INST; i++) {
    if (enb_node_list != NULL) {
      omv_data.geo[i].x = (enb_node_list->node->X_pos < 0.0)? 0.0 : enb_node_list->node->X_pos;
      omv_data.geo[i].y = (enb_node_list->node->Y_pos < 0.0)? 0.0 : enb_node_list->node->Y_pos;
      omv_data.geo[i].z = 1.0;
      omv_data.geo[i].mobility_type = oai_emulation.info.omg_model_enb;
      omv_data.geo[i].node_type = 0; //eNB
      enb_node_list = enb_node_list->next;
      omv_data.geo[i].Neighbors=0;

      for (j=NB_eNB_INST; j< NB_UE_INST + NB_eNB_INST ; j++) {
        if (is_UE_active(i,j - NB_eNB_INST ) == 1) {
          omv_data.geo[i].Neighbor[omv_data.geo[i].Neighbors]=  j;
          omv_data.geo[i].Neighbors++;
          LOG_D(OMG,"[eNB %d][UE %d] is_UE_active(i,j) %d geo (x%d, y%d) num neighbors %d\n", i,j-NB_eNB_INST, is_UE_active(i,j-NB_eNB_INST),
                omv_data.geo[i].x, omv_data.geo[i].y, omv_data.geo[i].Neighbors);
        }
      }
    }
  }

  for (i=NB_eNB_INST; i<NB_UE_INST+NB_eNB_INST; i++) {
    if (ue_node_list != NULL) {
      omv_data.geo[i].x = (ue_node_list->node->X_pos < 0.0) ? 0.0 : ue_node_list->node->X_pos;
      omv_data.geo[i].y = (ue_node_list->node->Y_pos < 0.0) ? 0.0 : ue_node_list->node->Y_pos;
      omv_data.geo[i].z = 1.0;
      omv_data.geo[i].mobility_type = oai_emulation.info.omg_model_ue;
      omv_data.geo[i].node_type = 1; //UE
      //trial
      omv_data.geo[i].state = 1;
      omv_data.geo[i].rnti = 88;
      omv_data.geo[i].connected_eNB = 0;
      omv_data.geo[i].RSRP = 66;
      omv_data.geo[i].RSRQ = 55;
      omv_data.geo[i].Pathloss = 44;
      omv_data.geo[i].RSSI[0] = 33;
      omv_data.geo[i].RSSI[1] = 22;
      omv_data.geo[i].RSSI[2] = 11;

      ue_node_list = ue_node_list->next;
      omv_data.geo[i].Neighbors=0;

      for (j=0; j< NB_eNB_INST ; j++) {
        if (is_UE_active(j,i-NB_eNB_INST) == 1) {
          omv_data.geo[i].Neighbor[ omv_data.geo[i].Neighbors]=j;
          omv_data.geo[i].Neighbors++;
          LOG_D(OMG,"[UE %d][eNB %d] is_UE_active  %d geo (x%d, y%d) num neighbors %d\n", i-NB_eNB_INST,j, is_UE_active(j,i-NB_eNB_INST),
                omv_data.geo[i].x, omv_data.geo[i].y, omv_data.geo[i].Neighbors);
        }
      }
    }
  }

  if( write( pfd, &omv_data, sizeof(struct Data_Flow_Unit) ) == -1 )
    perror( "write omv failed" );

  return 1;
}

void omv_end (int pfd, Data_Flow_Unit omv_data)
{
  omv_data.end=1;

  if( write( pfd, &omv_data, sizeof(struct Data_Flow_Unit) ) == -1 )
    perror( "write omv failed" );
}
#endif


text:
#pragma once

#include <set>
#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#endif
#include "Util.h"
using namespace std;

class Interface;

#ifdef _NOGSL
class gsl_vector;
#endif

class LocalState {
public:
	vector<gsl_vector*> startStates;
	vector<gsl_vector*> localSols; // local sol for different betas
	vector<double> errors;

	LocalState(int ncontrols, int nbetas) {
		for (int i = 0; i < nbetas; i++) {
			startStates.push_back(gsl_vector_alloc(ncontrols));
			localSols.push_back(gsl_vector_alloc(ncontrols));
		}
		errors.resize(nbetas);
	}
	~LocalState() {
		for (size_t i = 0; i < localSols.size(); i++) {
			gsl_vector_free(localSols[i]);
		}
		for (size_t i = 0; i < startStates.size(); i++) {
			gsl_vector_free(startStates[i]);
		}
	}

	void print() {
		cout << "beta = -1:  " << Util::print(localSols[0]) << endl;
		cout << "beta = -10: " << Util::print(localSols[1]) << endl;
		cout << "beta = -50: " << Util::print(localSols[2]) << endl;
	}
};

class OptimizationWrapper {
public:
	virtual bool optimize(Interface* inputs, gsl_vector* initState, const set<int>& constraints, int minimizeNode, bool suppressPrint, int MAX_TRIES, LocalState* localState, int level = -1) = 0;
	virtual gsl_vector* getMinState() = 0;
	virtual void randomizeCtrls(gsl_vector* x, Interface* inputs, const set<int>& constraints, int minimizeNode) = 0;
	virtual double getObjectiveVal() = 0;
	//virtual bool maximize(Interface* inputs, const gsl_vector* initState, const set<int>& assertConstraints, int minimizeNode, float beta, Predicate* pred, int predVal, int level) = 0;
};

ext: h

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
/*
Copyright 2010-2011, D. E. Shaw Research.
All rights reserved.

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

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

* 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.

* Neither the name of D. E. Shaw Research nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER 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.
*/
#include <gsl/gsl_randist.h>
#include <stdio.h>
#include "random123/philox.h"
#include "random123/threefry.h"
#include "random123/gsl_microrng.h"

/* Compute pi, using the gsl_ran_flat distribution with
   an underlying threefry4x64 counter-based rng (cbrng).
   We can call cbrng 8 times between calls to cbrng_reset */

GSL_MICRORNG(cbrng, threefry4x64); /* creates gsl_rng_cbrng */

#include "pi_check.h"

int main(int argc, char **argv){
    unsigned long hits = 0, tries = 0;
    gsl_rng *r;
    (void)argc; (void)argv; /* unused */

    threefry4x64_ctr_t c = {{0}};
    threefry4x64_key_t k = {{0}};
    r = gsl_rng_alloc(gsl_rng_cbrng);
    printf("%lu uniforms from %s\n", NTRIES, gsl_rng_name(r));
    while (tries < NTRIES) {
            double x, y;
            c.v[0]++; /* increment the counter */
            cbrng_reset(r, c, k); /* reset the rng to the new counter */
            x = gsl_ran_flat (r, -1.0, 1.0);
            y = gsl_ran_flat (r, -1.0, 1.0);
            if( x*x + y*y < 1.0 )
                hits++;
	    tries++;
    }
    gsl_rng_free (r);
    return pi_check(hits, tries);
}

ext: c

max_stars_repo_name: kabicm/arbor

max_issues_repo_name: kabicm/arbor

max_forks_repo_name: kabicm/arbor


text:
/**
 *
 * @file qwrapper_zlanhe.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlanhe(Quark *quark, Quark_Task_Flags *task_flags,
                       int norm, PLASMA_enum uplo, int N,
                       const PLASMA_Complex64_t *A, int LDA, int szeA,
                       int szeW, double *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANHE;
    QUARK_Insert_Task(quark, CORE_zlanhe_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(double)*szeW,                 NULL,          SCRATCH,
        sizeof(double),                     result,         OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlanhe_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, PLASMA_enum uplo, int N,
                          const PLASMA_Complex64_t *A, int LDA, int szeA,
                          int szeW, double *result,
                          double *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANHE;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_zlanhe_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double)*szeF,                 result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_zlanhe_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm, VALUE,
            sizeof(PLASMA_enum),                &uplo, VALUE,
            sizeof(int),                        &N,    VALUE,
            sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,  VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT,
            sizeof(double)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlanhe_quark = PCORE_zlanhe_quark
#define CORE_zlanhe_quark PCORE_zlanhe_quark
#endif
void CORE_zlanhe_quark(Quark *quark)
{
    double *normA;
    int norm;
    PLASMA_enum uplo;
    int N;
    PLASMA_Complex64_t *A;
    int LDA;
    double *work;

    quark_unpack_args_7(quark, norm, uplo, N, A, LDA, work, normA);
    *normA = LAPACKE_zlanhe_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm), lapack_const(uplo),
        N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlanhe_f1_quark = PCORE_zlanhe_f1_quark
#define CORE_zlanhe_f1_quark PCORE_zlanhe_f1_quark
#endif
void CORE_zlanhe_f1_quark(Quark *quark)
{
    double *normA;
    int norm;
    PLASMA_enum uplo;
    int N;
    PLASMA_Complex64_t *A;
    int LDA;
    double *work;
    double *fake;

    quark_unpack_args_8(quark, norm, uplo, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_zlanhe_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm), lapack_const(uplo),
        N, A, LDA, work);
}


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

#include <stdexcept>
#include <gsl/gsl>

namespace winlogo::utils {

/**
 * This exception is raised when failing to run VirtualProtect.
 */
class VirtualProtectFailedError : public std::runtime_error {
public:
    VirtualProtectFailedError() noexcept;
};

/**
 * This is a RAII class for temporarily changing the protection of a memory region.
 */
class MemoryProtectionGuard {
public:
    MemoryProtectionGuard(void* startAddress, size_t size, DWORD protection);

    NO_COPY(MemoryProtectionGuard);
    NO_MOVE(MemoryProtectionGuard);

    ~MemoryProtectionGuard() noexcept;

private:
    /// The start address of the memory region for which the protection should be changed.
    void* m_startAddress;
    /// The size of the memory region for which the protection should be changed.
    size_t m_size;
    /// The old protection to be restored in the end.
    DWORD m_oldProtection;

    /**
     * Run VirtualProtect for the stored address and store the previous protection in
     * m_oldProtection.
     */
    bool doVirtualProtect(DWORD protection) noexcept;
};

}  // namespace winlogo::utils


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

#include <config.h>
#include <math.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_dft_complex.h>
#include <gsl/gsl_errno.h>

#include "complex_internal.h"

int
gsl_dft_complex_forward (const double data[],
                         const size_t stride, const size_t n,
                         double result[])
{
  gsl_fft_direction sign = forward;
  int status = gsl_dft_complex (data, stride, n, result, sign);
  return status;
}

int
gsl_dft_complex_backward (const double data[],
                         const size_t stride, const size_t n,
                         double result[])
{
  gsl_fft_direction sign = backward;
  int status = gsl_dft_complex (data, stride, n, result, sign);
  return status;
}

int
gsl_dft_complex_inverse (const double data[],
                         const size_t stride, const size_t n,
                         double result[])
{
  gsl_fft_direction sign = backward;
  int status = gsl_dft_complex (data, stride, n, result, sign);

  /* normalize inverse fft with 1/n */

  {
    const ATOMIC norm = 1.0 / n;
    size_t i;
    for (i = 0; i < n; i++)
      {
        REAL(result,stride,i) *= norm;
        IMAG(result,stride,i) *= norm;
      }
  }
  return status;
}

int
gsl_dft_complex (const double data[],
                 const size_t stride, const size_t n,
                 double result[],
                 const gsl_fft_direction sign)
{

  size_t i, j, exponent;

  const double d_theta = 2.0 * ((int) sign) * M_PI / (double) n;

  /* FIXME: check that input length == output length and give error */

  for (i = 0; i < n; i++)
    {
      ATOMIC sum_real = 0;
      ATOMIC sum_imag = 0;

      exponent = 0;

      for (j = 0; j < n; j++)
        {
          double theta = d_theta * (double) exponent;
          /* sum = exp(i theta) * data[j] */

          ATOMIC w_real = cos (theta);
          ATOMIC w_imag = sin (theta);

          ATOMIC data_real = REAL(data,stride,j);
          ATOMIC data_imag = IMAG(data,stride,j);

          sum_real += w_real * data_real - w_imag * data_imag;
          sum_imag += w_real * data_imag + w_imag * data_real;

          exponent = (exponent + i) % n;
        }
      REAL(result,stride,i) = sum_real;
      IMAG(result,stride,i) = sum_imag;
    }

  return 0;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#pragma once

#include <gsl/gsl>

struct ScanWordResult;
struct GlyphVertex2d;
struct TigFont;
struct TigTextStyle;
struct TigRect;
struct TigFontMetrics;
using namespace gsl;

namespace gfx {
	class RenderingDevice;
	class ShapeRenderer2d;
	class TextEngine;
}

class FontRenderer {
public:
	explicit FontRenderer(gfx::RenderingDevice& g);
	~FontRenderer();

	void RenderRun(cstring_span<> text,
		int x,
		int y,
		const TigRect& bounds,
		TigTextStyle& style,
		const TigFont& font);

private:

	void RenderGlyphs(const GlyphVertex2d* vertices2d, int textureId, int glyphCount);
	static void Rotate2d(float& x, float& y,
		float rotCos, float rotSin,
		float centerX, float centerY);

	struct Impl;
	std::unique_ptr<Impl> mImpl;

};

/*
Separates a block of text given flags into words split up
on lines and renders them.
*/
class TextLayouter {
	friend class FontRenderer;
public:
	TextLayouter(gfx::RenderingDevice& device, gfx::ShapeRenderer2d& shapeRenderer);
	~TextLayouter();

	void LayoutAndDraw(gsl::cstring_span<> text, const TigFont &font, TigRect& extents, TigTextStyle& style);

	void Measure(const TigFont &font, const TigTextStyle &style, TigFontMetrics &metrics);

private:
	void DrawBackgroundOrOutline(const TigRect& rect, const TigTextStyle& style);
	static int GetGlyphIdx(char ch, const char *text);
	ScanWordResult ScanWord(const char* text,
		int firstIdx,
		int textLength,
		int tabWidth,
		bool lastLine,
		const TigFont& font,
		const TigTextStyle& style,
		int remainingSpace);
	std::pair<int, int> TextLayouter::MeasureCharRun(cstring_span<> text,
		const TigTextStyle& style,
		const TigRect& extents,
		int extentsWidth,
		const TigFont& font,
		int linePadding,
		bool lastLine);
	bool HasMoreText(cstring_span<> text, int tabWidth);

	void LayoutAndDrawVanilla(gsl::cstring_span<> text,
		const TigFont &font,
		TigRect& extents,
		TigTextStyle& style);

	void MeasureVanilla(const TigFont &font,
		const TigTextStyle &style,
		TigFontMetrics &metrics) const;

	uint32_t MeasureVanillaLine(const TigFont &font, const TigTextStyle &style, const char *text) const;
	uint32_t MeasureVanillaParagraph(const TigFont &font, const TigTextStyle &style, const char *text) const;
	uint32_t CountLinesVanilla(uint32_t maxWidth, uint32_t maxLines, const char *text, const TigFont &font, const TigTextStyle &style) const;

	static const char *sEllipsis;
	gfx::TextEngine &mTextEngine;
	FontRenderer mRenderer;
	std::unique_ptr<class FontsMapping> mMapping;
	gfx::ShapeRenderer2d& mShapeRenderer;
};

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
#include "linearwaves.h"
#include <gsl/gsl_integration.h>
typedef struct GSL_params {
    double x;
    int m;
} GSL_params;

void init_params(void) {

    params.mp = 1e-5;
    params.softening = .6*.05;
    params.softening2 = params.softening*params.softening;
    params.a = 1;
    params.indirect = TRUE;
    return;

}

double potential(double phi, double x) {
    double res = -pow(x*x + params.a*params.a + params.softening2 - 2*params.a*x*cos(phi),-.5);
    if (params.indirect) {
        res += cos(phi)*params.a/(x*x);
    }
    return res;
}

double dr_potential(double phi, double x) {
    double res = ( x - params.a*cos(phi)) * pow(x*x + params.a*params.a + params.softening2 - 2*params.a*x*cos(phi),-1.5);

    if (params.indirect) {
        res -= 2* cos(phi)*params.a/(x*x*x);
    }
    return res;
}

double dp_potential(double phi, double x) {
    double res = params.a*x*sin(phi)* pow(x*x + params.a*params.a + params.softening2 - 2*params.a*x*cos(phi),-1.5);

    if (params.indirect) {
        res +=  sin(phi)*params.a/(x*x);
    }
    return res;
}

double pfunc(double phi, void *params) {
    GSL_params p = *(GSL_params *)params;
    double x  = p.x;
    int m = p.m;
    return cos(m*phi) * potential(phi,x);
}
double dpfunc(double phi, void *params) {
    GSL_params p = *(GSL_params *)params;
    double x  = p.x;
    int m = p.m;
    return cos(m*phi) * dr_potential(phi,x);
}

void force(double x, int m, double complex *res) {
    int nspace = 1000;
    int gsl_order = 5;
    gsl_integration_workspace *w = gsl_integration_workspace_alloc(nspace);

    gsl_function F1, F2;
    F1.function = &pfunc;
    F2.function = &dpfunc;
    GSL_params p;
    p.x = x;
    p.m = m;
    F1.params = &p;
    F2.params = &p;

    double tol = params.tol;
    double error;

    double res_r, dr_res_r;
    gsl_integration_qag(&F1,0,M_PI, tol, tol, nspace, gsl_order , w, &res_r, &error);
    gsl_integration_qag(&F2,0,M_PI, tol, tol, nspace, gsl_order , w, &dr_res_r, &error);
    res_r *= params.mp/M_PI;
    dr_res_r *= params.mp/M_PI;

    res[0] = dr_res_r;// *-1;
    res[1] = res_r ;//* -I * m/x;
    res[2] = 0;
    gsl_integration_workspace_free(w);
    return;
}


text:
static char help[] =
"Solves time-dependent heat equation in 2D using TS.  Option prefix -ht_.\n"
"Equation is  u_t = D_0 laplacian u + f.  Domain is (0,1) x (0,1).\n"
"Boundary conditions are non-homogeneous Neumann in x and periodic in y.\n"
"Energy is conserved (for these particular conditions/source) and an extra\n"
"'monitor' is demonstrated.  Discretization is by centered finite differences.\n"
"Converts the PDE into a system  X_t = G(t,X) (PETSc type 'nonlinear') by\n"
"method of lines.  Uses backward Euler time-stepping by default.\n";

#include <petsc.h>

typedef struct {
  double D0;    // conductivity
} HeatCtx;

static double f_source(double x, double y) {
    return 3.0 * exp(-25.0 * (x-0.6) * (x-0.6)) * sin(2.0*PETSC_PI*y);
}

static double gamma_neumann(double y) {
    return sin(6.0 * PETSC_PI * y);
}

extern PetscErrorCode Spacings(DMDALocalInfo*, double*, double*);
extern PetscErrorCode EnergyMonitor(TS, int, double, Vec, void*);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, double, double**,
                                           double**, HeatCtx*);
extern PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo*, double, double**,
                                           Mat, Mat, HeatCtx*);

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  HeatCtx        user;
  TS             ts;
  Vec            u;
  DM             da;
  DMDALocalInfo  info;
  double         t0, tf;
  PetscBool      monitorenergy = PETSC_FALSE;

  PetscInitialize(&argc,&argv,(char*)0,help);

  user.D0  = 1.0;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "ht_", "options for heat", ""); CHKERRQ(ierr);
  ierr = PetscOptionsReal("-D0","constant thermal diffusivity",
           "heat.c",user.D0,&user.D0,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-monitor","also display total heat energy at each step",
           "heat.c",monitorenergy,&monitorenergy,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

//STARTDMDASETUP
  ierr = DMDACreate2d(PETSC_COMM_WORLD,
      DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR,
      5,4,PETSC_DECIDE,PETSC_DECIDE,  // default to hx=hx=0.25 grid
      1,1,                            // degrees of freedom, stencil width
      NULL,NULL,&da); CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(da,&u); CHKERRQ(ierr);
//ENDDMDASETUP

//STARTTSSETUP
  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetDM(ts,da); CHKERRQ(ierr);
  ierr = TSSetApplicationContext(ts,&user); CHKERRQ(ierr);
  ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
           (DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDATSSetRHSJacobianLocal(da,
           (DMDATSRHSJacobianLocal)FormRHSJacobianLocal,&user); CHKERRQ(ierr);
  if (monitorenergy) {
      ierr = TSMonitorSet(ts,EnergyMonitor,&user,NULL); CHKERRQ(ierr);
  }
  ierr = TSSetType(ts,TSBDF); CHKERRQ(ierr);
  ierr = TSSetTime(ts,0.0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,0.1); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,0.001); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
//ENDTSSETUP

  // report on set up
  ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
  ierr = TSGetMaxTime(ts,&tf); CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
           "solving on %d x %d grid for t0=%g to tf=%g ...\n",
           info.mx,info.my,t0,tf); CHKERRQ(ierr);

  // solve
  ierr = VecSet(u,0.0); CHKERRQ(ierr);   // initial condition
  ierr = TSSolve(ts,u); CHKERRQ(ierr);

  VecDestroy(&u);  TSDestroy(&ts);  DMDestroy(&da);
  return PetscFinalize();
}

PetscErrorCode Spacings(DMDALocalInfo *info, double *hx, double *hy) {
    if (hx)  *hx = 1.0 / (double)(info->mx-1);
    if (hy)  *hy = 1.0 / (double)(info->my);   // periodic direction
    return 0;
}

//STARTMONITOR
PetscErrorCode EnergyMonitor(TS ts, int step, double time, Vec u,
                             void *ctx) {
    PetscErrorCode ierr;
    HeatCtx        *user = (HeatCtx*)ctx;
    double         lenergy = 0.0, energy, dt, hx, hy, **au;
    int            i,j;
    MPI_Comm       com;
    DM             da;
    DMDALocalInfo  info;

    ierr = TSGetDM(ts,&da); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(da,u,&au); CHKERRQ(ierr);
    for (j = info.ys; j < info.ys + info.ym; j++) {
        for (i = info.xs; i < info.xs + info.xm; i++) {
            if ((i == 0) || (i == info.mx-1))
                lenergy += 0.5 * au[j][i];
            else
                lenergy += au[j][i];
        }
    }
    ierr = DMDAVecRestoreArrayRead(da,u,&au); CHKERRQ(ierr);
    ierr = Spacings(&info,&hx,&hy); CHKERRQ(ierr);
    lenergy *= hx * hy;
    ierr = PetscObjectGetComm((PetscObject)(da),&com); CHKERRQ(ierr);
    ierr = MPI_Allreduce(&lenergy,&energy,1,MPI_DOUBLE,MPI_SUM,com); CHKERRQ(ierr);
    ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"  energy = %9.2e     nu = %8.4f\n",
                energy,user->D0*dt/(hx*hy)); CHKERRQ(ierr);
    return 0;
}
//ENDMONITOR

//STARTRHSFUNCTION
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info,
                                    double t, double **au,
                                    double **aG, HeatCtx *user) {
  PetscErrorCode ierr;
  int      i, j, mx = info->mx;
  double   hx, hy, x, y, ul, ur, uxx, uyy;

  ierr = Spacings(info,&hx,&hy); CHKERRQ(ierr);
  for (j = info->ys; j < info->ys + info->ym; j++) {
      y = hy * j;
      for (i = info->xs; i < info->xs + info->xm; i++) {
          x = hx * i;
          // apply Neumann b.c.s
          ul = (i == 0) ? au[j][i+1] + 2.0 * hx * gamma_neumann(y)
                        : au[j][i-1];
          ur = (i == mx-1) ? au[j][i-1] : au[j][i+1];
          uxx = (ul - 2.0 * au[j][i]+ ur) / (hx*hx);
          // DMDA is periodic in y
          uyy = (au[j-1][i] - 2.0 * au[j][i]+ au[j+1][i]) / (hy*hy);
          aG[j][i] = user->D0 * (uxx + uyy) + f_source(x,y);
      }
  }
  return 0;
}
//ENDRHSFUNCTION

//STARTRHSJACOBIAN
PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo *info,
                                    double t, double **au,
                                    Mat J, Mat P, HeatCtx *user) {
    PetscErrorCode ierr;
    int            i, j, ncols;
    const double   D = user->D0;
    double         hx, hy, hx2, hy2, v[5];
    MatStencil     col[5],row;

    ierr = Spacings(info,&hx,&hy); CHKERRQ(ierr);
    hx2 = hx * hx;  hy2 = hy * hy;
    for (j = info->ys; j < info->ys+info->ym; j++) {
        row.j = j;  col[0].j = j;
        for (i = info->xs; i < info->xs+info->xm; i++) {
            // set up a standard 5-point stencil for the row
            row.i = i;
            col[0].i = i;
            v[0] = - 2.0 * D * (1.0 / hx2 + 1.0 / hy2);
            col[1].j = j-1;  col[1].i = i;    v[1] = D / hy2;
            col[2].j = j+1;  col[2].i = i;    v[2] = D / hy2;
            col[3].j = j;    col[3].i = i-1;  v[3] = D / hx2;
            col[4].j = j;    col[4].i = i+1;  v[4] = D / hx2;
            ncols = 5;
            // if at the boundary, edit the row back to 4 nonzeros
            if (i == 0) {
                ncols = 4;
                col[3].j = j;  col[3].i = i+1;  v[3] = 2.0 * D / hx2;
            } else if (i == info->mx-1) {
                ncols = 4;
                col[3].j = j;  col[3].i = i-1;  v[3] = 2.0 * D / hx2;
            }
            ierr = MatSetValuesStencil(P,1,&row,ncols,col,v,INSERT_VALUES); CHKERRQ(ierr);
        }
    }

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}
//ENDRHSJACOBIAN


text:
#ifndef _EM_UTILS_H_
#define _EM_UTILS_H_ 1

#include <petsc.h>

#include <string>

std::string string_format(const char*, ...);
std::string parse_string(const std::string &s);

class LogEventHelper {
public:
  LogEventHelper(PetscLogEvent ple) : ple_(ple) { PetscLogEventBegin(ple_, 0, 0, 0, 0); }
  ~LogEventHelper() { PetscLogEventEnd(ple_, 0, 0, 0, 0); }

private:
  PetscLogEvent ple_;
};

class LogStageHelper {
public:
  LogStageHelper(const std::string &name) {
    PetscLogStage pls;
    PetscLogStageRegister(name.c_str(), &pls);
    PetscLogStagePush(pls);
  }
  ~LogStageHelper() { PetscLogStagePop(); }
};

#endif

max_stars_repo_name: emfem/emfem

max_issues_repo_name: emfem/emfem

max_forks_repo_name: emfem/emfem


text:
/**
 *
 * @file qwrapper_zlag2c.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions mixed zc -> ds
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlag2c(Quark *quark, Quark_Task_Flags *task_flags,
                       int m, int n, int nb,
                       const PLASMA_Complex64_t *A, int lda,
                       PLASMA_Complex32_t *B, int ldb,
                       PLASMA_sequence *sequence, PLASMA_request *request)
{
    DAG_CORE_LAG2C;
    QUARK_Insert_Task(quark, CORE_zlag2c_quark, task_flags,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    B,                 OUTPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_sequence*),           &sequence,  VALUE,
        sizeof(PLASMA_request*),            &request,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlag2c_quark = PCORE_zlag2c_quark
#define CORE_zlag2c_quark PCORE_zlag2c_quark
#endif
void CORE_zlag2c_quark(Quark *quark)
{
    int m;
    int n;
    PLASMA_Complex64_t *A;
    int lda;
    PLASMA_Complex32_t *B;
    int ldb;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    int info;

    quark_unpack_args_8(quark, m, n, A, lda, B, ldb, sequence, request);
    info = LAPACKE_zlag2c_work(LAPACK_COL_MAJOR, m, n, A, lda, B, ldb);
    if (sequence->status == PLASMA_SUCCESS && info != 0)
        plasma_sequence_flush(quark, sequence, request, info);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_clag2z(Quark *quark, Quark_Task_Flags *task_flags,
                       int m, int n, int nb,
                       const PLASMA_Complex32_t *A, int lda,
                       PLASMA_Complex64_t *B, int ldb)
{
    QUARK_Insert_Task(quark, CORE_clag2z_quark, task_flags,
        sizeof(int),                        &m,     VALUE,
        sizeof(int),                        &n,     VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    A,             INPUT,
        sizeof(int),                        &lda,   VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    B,             INOUT,
        sizeof(int),                        &ldb,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_clag2z_quark = PCORE_clag2z_quark
#define CORE_clag2z_quark PCORE_clag2z_quark
#endif
void CORE_clag2z_quark(Quark *quark)
{
    int m;
    int n;
    PLASMA_Complex32_t *A;
    int lda;
    PLASMA_Complex64_t *B;
    int ldb;

    quark_unpack_args_6(quark, m, n, A, lda, B, ldb);
    LAPACKE_clag2z_work(LAPACK_COL_MAJOR, m, n, A, lda, B, ldb);
}


text:
/* Provides the Lua Linear module. See LICENSE for license terms. */

#ifndef LUALINEAR_INCLUDED
#define LUALINEAR_INCLUDED

#include <assert.h>
#include <stdlib.h>
#include <lua.h>
#include <lauxlib.h>
#include <cblas.h>

/**
 * Meta tables.
 */
#define LUALINEAR_VECTOR_METATABLE "linear.vector"
#define LUALINEAR_MATRIX_METATABLE "linear.matrix"

/**
 * OpenMP.
 */
#define LUALINEAR_OMP_MINSIZE 2048

/**
 * Vector.
 */
struct vector {
        int size;
        int inc;
        float *values;
        int ref;
};

/**
 * Pushes a new vector onto the stack.
 */
inline struct vector *lualinear_newvector (lua_State *L, int size) {
	struct vector *vector;

	assert(size >= 1);
	vector = (struct vector *)lua_newuserdata(L, sizeof(struct vector));
	vector->size = size;
	vector->inc = 1;
	vector->values = NULL;
	vector->ref = LUA_NOREF;
	luaL_getmetatable(L, LUALINEAR_VECTOR_METATABLE);
	lua_setmetatable(L, -2);
	vector->values = (float *)calloc(size, sizeof(float));
	if (vector->values == NULL) {
		luaL_error(L, "cannot allocate values");
	}
        return vector;
}

/**
 * Pushes an existing vector onto the stack.
 */
inline struct vector *lualinear_wrapvector (lua_State *L, int size,
		float *values) {
	struct vector *vector;

	assert(size >= 1);
	vector = (struct vector *)lua_newuserdata(L, sizeof(struct vector));
	vector->size = size;
	vector->inc = 1;
	vector->values = values;
	vector->ref = LUA_REFNIL;
	luaL_getmetatable(L, LUALINEAR_VECTOR_METATABLE);
	lua_setmetatable(L, -2);
	return vector;
}

/**
 * Matrix.
 */
struct matrix {
        int rows;
        int cols;
        int ld;
        CBLAS_ORDER order;
        float *values;
        int ref;
};

/**
 * Pushes a new matrix onto the stack.
 */
inline struct matrix *lualinear_newmatrix (lua_State *L, int rows, int cols,
                CBLAS_ORDER order) {
	struct matrix *matrix;

	assert(rows >= 1 && cols >= 1);
	matrix = (struct matrix *)lua_newuserdata(L, sizeof(struct matrix));
	matrix->rows = rows;
	matrix->cols = cols;
	matrix->ld = order == CblasRowMajor ? cols : rows;
	matrix->order = order;
	matrix->values = NULL;
	matrix->ref = LUA_NOREF;
	luaL_getmetatable(L, LUALINEAR_MATRIX_METATABLE);
	lua_setmetatable(L, -2);
	matrix->values = (float *)calloc((size_t)rows * cols, sizeof(float));
	if (matrix->values == NULL) {
		luaL_error(L, "cannot allocate values");
	}
	return matrix;
}

/**
 * Pushes an existing matrix onto the stack-
 */
inline struct matrix *lualinear_wrapmatrix (lua_State *L, int rows, int cols,
		CBLAS_ORDER order, float *values) {
	struct matrix *matrix;

	assert(rows >= 1 && cols >= 1);
	matrix = (struct matrix *)lua_newuserdata(L, sizeof(struct matrix));
	matrix->rows = rows;
	matrix->cols = cols;
	matrix->ld = order == CblasRowMajor ? cols : rows;
	matrix->order = order;
	matrix->values = values;
	matrix->ref = LUA_REFNIL;
	luaL_getmetatable(L, LUALINEAR_MATRIX_METATABLE);
	lua_setmetatable(L, -2);
	return matrix;
}

/**
 * Opens the Linear module in a Lua state.
 */
int luaopen_linear (lua_State *L);

#endif /* LUALINEAR_INCLUDED */


text:
/* $Id$  */
/*--------------------------------------------------------------------*/
/*;  Copyright (C) 2013                                               */
/*;  Associated Universities, Inc. Washington DC, USA.                */
/*;                                                                   */
/*;  This program is free software; you can redistribute it and/or    */
/*;  modify it under the terms of the GNU General Public License as   */
/*;  published by the Free Software Foundation; either version 2 of   */
/*;  the License, or (at your option) any later version.              */
/*;                                                                   */
/*;  This program is distributed in the hope that it will be useful,  */
/*;  but WITHOUT ANY WARRANTY; without even the implied warranty of   */
/*;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    */
/*;  GNU General Public License for more details.                     */
/*;                                                                   */
/*;  You should have received a copy of the GNU General Public        */
/*;  License along with this program; if not, write to the Free       */
/*;  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,     */
/*;  MA 02139, USA.                                                   */
/*;                                                                   */
/*;Correspondence about this software should be addressed as follows: */
/*;         Internet email: bcotton@nrao.edu.                         */
/*;         Postal address: William Cotton                            */
/*;                         National Radio Astronomy Observatory      */
/*;                         520 Edgemont Road                         */
/*;                         Charlottesville, VA 22903-2475 USA        */
/*--------------------------------------------------------------------*/

#include <time.h>
#include <gsl/gsl_randist.h>
#include "ObitThread.h"
#include "ObitOTFUtil.h"

/*----------------Obit: Merx mollis mortibus nuper ------------------*/
/**
 * \file ObitVEGASUtil.c
 * ObitOTF class utility function definitions.
 */

/*---------------Private structures----------------*/

/*---------------Private function prototypes----------------*/
/*----------------------Public functions---------------------------*/
/**
 * Average the frequencies in a GBT/VEGAS OTF
 * \param inVEGAS  Input VEGAS, any calibration and flagging should be applied.
 * \param outVEGAS Output VEGAS, must already be defined
 * \param chAvg  Number of channels to average, -1 => all
 * \param err    Error stack
 */
void ObitVEGASUtilAverage(ObitOTF *inOTF, ObitOTF *outOTF, olong chAvg,
			  ObitErr *err)
{
  const ObitClassInfo *ParentClass;
  ObitIOCode retCode;
  gboolean doCalSelect, done;
  olong firstRec, navg, nchanIn, nchanOut;
  olong ichan, ochan, ochn, ifeed, nfeed, istok, nstok, indx, ondx;
  ObitInfoType type;
  ObitIOAccess access;
  gint32 dim[MAXINFOELEMDIM] = {1,1,1,1,1};
  ofloat sum, sumWt, *iData, *oData;
  olong NPIO, iran, irec;
  ObitOTFDesc *iDesc, *oDesc;
  ObitHistory *inHist=NULL, *outHist=NULL;
  /* Don't copy Cal and Soln or data or flag tables */
  gchar *exclude[]={"OTFSoln", "OTFCal", "OTFScanData", "OTFFlag", NULL};
  gchar *routine = "ObitVEGASUtilAverage";

    /* error checks */
  g_assert (ObitErrIsA(err));
  if (err->error) return;
  g_assert (ObitOTFIsA(inOTF));
  g_assert (ObitOTFIsA(outOTF));
  /* Input and output must be different */
  Obit_return_if_fail ((!ObitOTFSame (inOTF, outOTF, err)), err,
		       "%s: Output cannot be the same as the input", routine);

  /* deep copy any base class members */
  ParentClass = ((ObitClassInfo*)inOTF->ClassInfo)->ParentClass;
  g_assert ((ParentClass!=NULL) && (ParentClass->ObitCopy!=NULL));
  ParentClass->ObitCopy (inOTF, outOTF, err);
  outOTF->mySel     = newObitOTFSel (outOTF->name);
  outOTF->tableList = newObitTableList(outOTF->name);
  outOTF->geom      = newObitOTFArrayGeom(outOTF->name);

  /* Calibration wanted? */
  doCalSelect = FALSE;
  ObitInfoListGetTest(inOTF->info, "doCalSelect", &type, (gint32*)dim,
		      &doCalSelect);
  if (doCalSelect) access = OBIT_IO_ReadCal;
  else access = OBIT_IO_ReadWrite;

  /* Make sure NPIO the same */
  NPIO = 1000;
  dim[0] = 1;
  ObitInfoListGetTest (inOTF->info, "nRecPIO", &type, dim,  &NPIO);
  ObitInfoListAlwaysPut (inOTF->info, "nRecPIO", OBIT_long, dim,  &NPIO);
  ObitInfoListAlwaysPut (outOTF->info, "nRecPIO", OBIT_long, dim,  &NPIO);

  /* Open Input Data */
  retCode = ObitOTFOpen (inOTF, access, err);
  if ((retCode != OBIT_IO_OK) || (err->error>0)) goto cleanup;

  /* How many channels to average? Defaults to all. */
  iDesc = inOTF->myDesc;
  nchanIn = iDesc->inaxes[iDesc->jlocf];
  if (chAvg>0) navg = chAvg;
  else         navg = nchanIn;
  navg     = MAX (1, MIN (navg, nchanIn));
  nchanOut = MAX (1, (olong)(((ofloat)(nchanIn)/navg)+0.9999));

  /* Number of feeds */
  nfeed = iDesc->inaxes[iDesc->jlocfeed];
  /* NUmber of Stokes */
  nstok = iDesc->inaxes[iDesc->jlocs];

  /* Copy descriptor */
  outOTF->myDesc = (gpointer)ObitOTFDescCopy(iDesc, outOTF->myDesc, err);
  outOTF->myDesc->nrecord = 0;     /* may not copy all */
  oDesc = outOTF->myDesc;

  /* Averaging in frequency */
  oDesc->inaxes[oDesc->jlocf]      = nchanOut;
  oDesc->colRepeat[oDesc->ncol-1] /= navg;
  oDesc->cdelt[oDesc->jlocf]      *= navg;
  oDesc->crpix[oDesc->jlocf]      /= navg;
  ObitOTFDescIndex (oDesc);

  /* copy/average Array Geometry - this time with full information
     NO, VEGAS doesn't have one "detector" per frequency
  outOTF->geom = ObitOTFArrayGeomAver(inOTF->geom, iDesc, outOTF->geom, oDesc, err);
  if (err->error) goto cleanup; */

  /* Open Output Data */
  retCode = ObitOTFOpen (outOTF, OBIT_IO_WriteOnly, err) ;
  if ((retCode != OBIT_IO_OK) || (err->error>0)) goto cleanup;

  /* Copy any history */
  inHist  = newObitHistoryValue("in history", inOTF->info, err);
  outHist = newObitHistoryValue("out history", outOTF->info, err);
  outHist = ObitHistoryCopy (inHist, outHist, err);
  if (err->error)  goto cleanup;
  inHist  = ObitHistoryUnref(inHist);
  outHist = ObitHistoryUnref(outHist);

   /* Copy tables before data  */
  retCode = ObitOTFCopyTables (inOTF, outOTF, exclude, NULL, err);
  if (err->error) goto cleanup;

  /* Close and reopen input to init calibration which will have
     been disturbed by the table copy */
  retCode = ObitOTFClose (inOTF, err);
  if (err->error) goto cleanup;

  retCode = ObitOTFOpen (inOTF, access, err);
  if ((retCode != OBIT_IO_OK) || (err->error>0)) goto cleanup;

  /* Loop over data */
  done = (retCode != OBIT_IO_OK);
  while (!done) {

    /* read buffer */#
    retCode = ObitOTFRead (inOTF, NULL, err);
    if (err->error) goto cleanup;
    done = (retCode == OBIT_IO_EOF); /* done? */
    if (done) break;
    firstRec = inOTF->myDesc->firstRec;

    iData = inOTF->buffer;   /* Input data pointer */
    oData = outOTF->buffer;  /* Output data pointer */

    /* How many? */
    outOTF->myDesc->numRecBuff = inOTF->myDesc->numRecBuff;

    /* Loop over buffer */
    for (irec=0; irec<iDesc->numRecBuff; irec++) {

      /* Copy random (descriptive parameters */
      for (iran=0; iran<iDesc->numDesc; iran++) oData[iran] = iData[iran];

      /* Average input to output */
      /* Loop over Stokes */
      for (istok=0; istok<nstok; istok++) {
	/* Loop over Feed */
	for (ifeed=0; ifeed<nfeed; ifeed++) {
	  /* Outer loop over channel */
	  ochn = 0;  /* output channel 0-rel index */
	  for (ochan=0; ochan<nchanIn; ochan+=navg) {
	    /* Inner frequency loop summing */
	    sum = sumWt = 0.0;
	    indx = iDesc->ilocdata + istok*iDesc->incs + ifeed*iDesc->incfeed;
	    ondx = oDesc->ilocdata + istok*oDesc->incs + ifeed*oDesc->incfeed;
	    for (ichan=ochan; ichan<ochan+navg; ichan++) {
	      sum   += iData[indx+ichan*iDesc->incf] * iData[indx+ichan*iDesc->incf+1];
	      sumWt += iData[indx+ichan*iDesc->incf+1];
	    } /* end inner loop */
	    /* Save average to outout */
	    if (sumWt>0.0) {
	      oData[ondx+ochn*oDesc->incf]   = sum/sumWt;
	      oData[ondx+ochn*oDesc->incf+1] = sumWt;
	    } else { /* bad */
	      oData[ondx+ochn*oDesc->incf]   = 0.0;
	      oData[ondx+ochn*oDesc->incf+1] = 0.0;
	    }
	    ochn++;
	  } /* end outer channel loop */
	} /* end feed loop */
      } /* end Stokes loop */
      iData += iDesc->lrec;  /* Update buffer pointers */
      oData += oDesc->lrec;
    } /* end loop over buffer */

    /* Write buffer */
    if (outOTF->myDesc->numRecBuff>0) retCode = ObitOTFWrite (outOTF, NULL, err);
    if (err->error) goto cleanup;
  } /* end loop over file */

 cleanup:
  /* unset output buffer  */
  outOTF->buffer = NULL;
  outOTF->bufferSize = 0;
  retCode = ObitOTFClose (outOTF, err); /* Close output */

  /* Close input */
  retCode = ObitOTFClose (inOTF, err);

  /* cleanup */
  if (err->error) Obit_traceback_msg (err, routine, inOTF->name);

} /* end ObitOTFUtilAverage */

max_stars_repo_name: sarrvesh/Obit

max_issues_repo_name: sarrvesh/Obit

max_forks_repo_name: sarrvesh/Obit


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_3.h"

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

/**
 * Function to obtain the coefficients of a 6 steps 3rd order Runge-Kutta
 * method.
 */
int
rk_tb_6_3 (Optimize * optimize) ///< Optimize struct.
{
  long double A[3], B[3], C[3], D[3];
  long double *tb, *r;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3: 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];
  b51 (tb) = r[11];
  b52 (tb) = r[12];
  b53 (tb) = r[13];
  b54 (tb) = r[14];
  b65 (tb) = r[15];
  b64 (tb) = r[16];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = 0.5L - b64 (tb) * t4 (tb) - b65 (tb) * t5 (tb);
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = 1.L / 3.L - b64 (tb) * sqr (t4 (tb)) - b65 (tb) * sqr (t5 (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 / 6.L - b64 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                                 + b43 (tb) * t3 (tb))
    - b65 (tb) * (b51 (tb) * t1 (tb) + b52 (tb) * t2 (tb)
                  + b53 (tb) * t3 (tb) + b54 (tb) * t4 (tb));
  solve_3 (A, B, C, D);
  if (isnan (D[0]) || isnan (D[1]) || isnan (D[2]))
    return 0;
  b63 (tb) = D[2];
  b62 (tb) = D[1];
  b61 (tb) = D[0];
  rk_b_6 (tb);
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 3rd order, 4th order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_6_3t (Optimize * optimize)        ///< Optimize struct.
{
  long double A[4], B[4], C[4], D[4], E[4];
  long double *tb, *r;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3t: 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];
  b51 (tb) = r[11];
  b52 (tb) = r[12];
  b53 (tb) = r[13];
  b54 (tb) = r[14];
  b65 (tb) = r[15];
  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);
  C[3] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[3] = b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb) + b43 (tb) * t3 (tb);
  E[3] = 1.L / 6.L - b65 (tb) * (b51 (tb) * t1 (tb) + b52 (tb) * t2 (tb)
                                 + b53 (tb) * t3 (tb) + b54 (tb) * t4 (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];
  rk_b_6 (tb);
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 2nd-3rd order Runge-Kutta
 * pair.
 */
int
rk_tb_6_3p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3p: start\n");
#endif
  if (!rk_tb_6_3 (optimize))
    return 0;
  tb = optimize->coefficient;
  e51 (tb) = 0.5L / t1 (tb);
  e52 (tb) = e53 (tb) = 0.L;
  rk_e_6 (tb);
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 6 steps 2nd-3rd order, 3rd-4th order
 * in equations depending only in time, Runge-Kutta pair.
 */
int
rk_tb_6_3tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3tp: start\n");
#endif
  if (!rk_tb_6_3t (optimize))
    return 0;
  tb = optimize->coefficient;
  e63 (tb) = e64 (tb) = 0.L;
  e62 (tb) = (1.L / 3.L - 0.5L * t1 (tb)) / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (e62 (tb)))
    return 0;
  e61 (tb) = (0.5L - e52 (tb) * t2 (tb)) / t1 (tb);
  if (isnan (e61 (tb)))
    return 0;
  rk_e_6 (tb);
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_tb_6_3tp: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 6 steps 3rd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_3 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_objective_tb_6_3: 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 (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_3
  fprintf (stderr, "rk_objective_tb_6_3: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_3: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 6 steps 3rd order, 4th
 * order in equations depending only in time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_3t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_objective_tb_6_3t: 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 (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_3
  fprintf (stderr, "rk_objective_tb_6_3t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_3t: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 6 steps 2nd-3rd order
 * Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_3p (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_objective_tb_6_3p: 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 (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 (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_3
  fprintf (stderr, "rk_objective_tb_6_3p: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_3p: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 6 steps 2nd-3rd order,
 * 3rd-4th order in equations depending only in time, Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_6_3tp (RK * rk) ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_6_3
  fprintf (stderr, "rk_objective_tb_6_3tp: 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 (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 (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_3
  fprintf (stderr, "rk_objective_tb_6_3tp: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_6_3tp: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
/*
 *  parameters.h
 *  Spike
 *
 *  Created by Ben Evans on 18/08/2008.
 *  Copyright 2008 University of Oxford. All rights reserved.
 *
 */

#ifndef _PARAMETERS_H
#define _PARAMETERS_H

#include <stdbool.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_version.h>

typedef int tstep; // Signed to allow spikeTimes[0] = -BIG //long long

typedef unsigned char uchar;

struct SIMULATION {
    bool Xgrid;
    tstep tally;
    tstep ptTS;
    tstep trainTS;
    tstep testTS;
    tstep totTS;
    float minTau;
    double start;
    double elapsed;
    double realSecPerSimSec;
} SIM;

typedef struct DIMENSIONS {
    int nRows;
    int nCols;
    int nFilt;
} DIM;

//extern SIMULATION SIM;

/*typedef enum {
    Off,
    //Uniform,
    Gaussian
} NOISE;*/

typedef enum {
    MinD,
    ConstD,
    UniformD,
    GaussD,
    SOMD
} DELAY;

typedef enum {
    //Zero,
    Constant, // Zero set using appropriate Dg modifier
    Uniform,
    Gaussian,
    SOM
} INITIALISATION;

typedef enum {
    None,
    MaintainLength,
    MaintainSum
} NORMALISATION;

/********************** Main Network Parmeters *************************/
// Update read_parameters and defaults.m when a new parameter is added
typedef struct {
    // Simulation
    float DT;
    bool initialised; // Internal
    int loops;
    bool pretrain;// = true;
    bool train;// = true;
    bool priorPhases;
    bool isolateEfE;
    bool trainPause;
    bool noise; // int
    float noiseScale;
    float SigmaE;
    float SigmaI;
    NORMALISATION normalise;
    int nRecordsPL;// = 0;
    int nRecords;
    int * vRecords;
    //float TotalTime;
    float MaxTime;
    float EpochTime;
    float TestTime;
    int EpochMS;
    int TestMS;
    int TotalMS;
    int RecordMS;
    //tstep TotalTS;
    int TSperMS;
    int spkBuffer;
    //int inpSpkBuff;
    bool printConnections;
    bool saveInputSpikes;
    bool probConnect;
    bool loadWeights;

    // Stimuli
    char * imgDir;
    char * imgList;
    bool randStimOrder;// = true;
    bool randTransOrder;
    bool randTransDirection;
    bool interleaveTrans;
    bool localRep;// = true;
    float current;// = 1.25e-9;
    float currentSpread;
    bool loadStimuli;
    bool stimGroups;
    int nBG; // # shared between groups
    int nWG; // # specific within each group
    int nGroups;
    int nStimuli;
    int nTransPS;
    bool newTestSet;
    int M; // Degree of training i.e. train with M of the NCK
    int K; // Number of stimuli presented simultaneously
    //int nTestGroups;
    int nTestStimuli;
    int nTestTransPS;
    float transP_Train;
    float transP_Test;
    int shift;
    int nFiringNeurons;
    float a; // Sparseness of stimuli
    bool useFilteredImages;
    bool gabor;
    int nScales;
    int * vScales;
    int nOrients;
    int * vOrients;
    int nPhases;
    int * vPhases;
    int sInputs;
    int nRows;
    int nCols;

    // Network
    int nLayers;// = 2;
    int nWLayers;
    bool inputInhib;
    int nExcit;// = 120;
    int * vExcit;
    int LvExcit;
    int nSynEfE;
    float * pCnxEfE;
    int LpEfE;
    int nSynElE;
    float * pCnxElE;
    int LpElE;
    int nSynIE;
    float * pCnxIE;
    int LpIE;
    int nInhib;// = 40;
    float rInhib;
    int * vInhib;
    int LvInhib;
    int nSynEI;
    float * pCnxEI;
    int LpEI;
    int nSynII;
    float * pCnxII;
    int LpII;
    INITIALISATION initEfE;
    float iEfE;
    bool axonDelay; //DELAY axonDelay;
    DELAY delayEfE;
    DELAY delayElE;
    DELAY delayEI;
    float d_const;
    float d_min;
    float d_max;
    float d_mean;
    float d_sd;
    float spatialScale;
    float condSpeed;
    float maxDelay;
    bool SOM;
    bool SOMinput;
    float SOMsigE;
    float SOMsigI;
    float SOMclip;
    //float SOMstrE;
    //float SOMstrI;
    bool trainEfE; // Internal var modified through isolateEfE
    bool trainElE;
    INITIALISATION initElE;
    float iElE;
    DIM * layDim;
    int * vSquare;

    // Cell bodies
    float capE;// = 2.0e-10;
    float capI;// = 10.0e-12;
    float gLeakE;// = 10.0e-9;
    float gLeakI;// = 5.0e-9;
    float VrestE;
    float VrestI;
    float VhyperE;
    float VhyperI;
    float ThreshE;
    float ThreshI;
    float VrevE;
    float VrevI;
    float refract;// = 0.002
    bool adaptation;
    float alphaCa;
    float tauCa;
    float gAHP;
    float VK;

    // Synapses (afferent axons)
    bool noSTDPdelay;
    float alphaC;// = 0.5;
    float tauC;// = 0.005;
    float alphaD;// = 0.5;
    float tauD;// = 0.009;
    float learnR;// = 0.1;
    float DgEfE; //modEf; // alter the strength of Excit f-f synapses
    //float tauEfE;// tauEE = 0.01;
    float DgElE;//Dg_ElE
    //float tauElE;
    float tauEE;
    float DgIE;//Dg_IE = 0.5;
    float tauIE;// = 0.001
    float DgEI;//Dg_EI = 0.5;
    float tauEI;// = 0.001;
    float DgII;//Dg_II = 0.5;
    float tauII;// = 0.001;
    float gMax;// = 48.0e-10;
} PARAMS;

//PARAMS * mp;
extern PARAMS * mp;
extern gsl_rng * mSeed;
extern gsl_rng ** states;

#endif

max_stars_repo_name: bdevans/spike

max_issues_repo_name: bdevans/spike

max_forks_repo_name: bdevans/spike


text:
#pragma once

#include <imageview/ContinuousImageView.h>
#include <imageview/ImageView.h>

#include <gsl/assert>

namespace imageview {

template <class PixelFormat, bool Mutable>
constexpr ImageView<PixelFormat, Mutable> crop(ImageView<PixelFormat, Mutable> image, unsigned int first_row,
                                               unsigned int first_column, unsigned int num_rows,
                                               unsigned int num_columns) {
  Expects(first_row <= image.height());
  Expects(first_column <= image.width());
  Expects(first_row + num_rows <= image.height());
  Expects(first_column + num_columns <= image.width());
  const std::size_t data_offset = (first_row * image.stride() + first_column) * PixelFormat::kBytesPerPixel;
  const std::size_t new_data_size =
      (num_rows == 0) ? 0 : ((num_rows - 1) * image.stride() + num_columns) * PixelFormat::kBytesPerPixel;
  const auto data_new = image.data().subspan(data_offset, new_data_size);
  return ImageView<PixelFormat>(num_rows, num_columns, image.stride(), data_new, image.pixelFormat());
}

template <class PixelFormat, bool Mutable>
constexpr ImageView<PixelFormat, Mutable> crop(ContinuousImageView<PixelFormat, Mutable> image, unsigned int first_row,
                                               unsigned int first_column, unsigned int num_rows,
                                               unsigned int num_columns) {
  return crop(ImageView<PixelFormat, Mutable>(image), first_column, first_column, num_rows, num_columns);
}

}  // namespace imageview

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview