# -*- coding: utf-8 -*-
"""NNratio_1D.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/13zC6nTSvjEFa_-FFy6-Rpy5zH8bUUlMJ
"""

from __future__ import print_function
from __future__ import division
import keras
import os
from keras.models import Sequential
from keras import layers
from keras import backend as K
from keras.layers import Activation, Dense
from keras.constraints import Constraint
from keras.constraints import max_norm

from math import *
import numpy as np
import scipy.integrate as integrate
import scipy.stats as st
import time
import concurrent.futures as cf

##!cat /proc/cpuinfo

##!cat /proc/meminfo

#!apt-get -qq install -y graphviz && pip install -q pydot
#import pydot


"""# Main"""

#sigmoid function
def sigmoid(x):
  return 1 / (1 + exp(-x))

#@title Glabal Parameters
global epochs = 100000 #@param {type:"integer"}
global batch_size = 1000 #@param {type:"integer"}
#Geteps = 0.005 #@param {type:"number"}
global Geteps = 0.0015811 #@param {type:"number"}
global Getmu = 0.8 #@param {type:"number"}
global Getsigma = 0.02 #@param {type:"number"}
global cut = 0. #@param {type:"number"}
global gcut = 0. #@param {type:"number"}

def SM(x):
  return exp(-8*x)

def BSM(x):
  return exp(-(x-Getmu)**2/(2*Getsigma**2))
#Normalize distribution
SM_norm = integrate.quad(lambda y :SM(y),cut,1)
BSM_norm = integrate.quad(lambda y :BSM(y),cut,1)
#SM_norm_c = integrate.quad(lambda y :SMn(y),gcut,1)

#normalized distribution
def SMn(x):
  return exp(-8*x)/SM_norm[0]

def BSMn(x):
  return (SMn(x)+Geteps*BSM(x)/BSM_norm[0])/(1+Geteps)

SM_norm_c = integrate.quad(lambda y :SM(y),gcut,1)

def SMnc(x):
  return exp(-8*x)/SM_norm_c[0]


#define probability distribution function
class P_SM(st.rv_continuous):
    def _pdf(self,x):
        return SMn(x)
      
class P_BSM(st.rv_continuous):
    def _pdf(self,x):
        return BSMn(x)        
      
class P_SMc(st.rv_continuous):
    def _pdf(self,x):
        return SMnc(x)
      
      
      
SM_gen = P_SM(a=cut,b=1,name='sm_sample')
BSM_gen = P_BSM(a=cut,b=1,name='bsm_sample')
SMc_gen = P_SMc(a=gcut,b=1,name='smc_sample')

global NRef = 200000
global NR = 20000
global Nbsm = NR*(1+Geteps)

def bsm2smRatio(x):
  return log(BSMn(x)/SMn(x)*Nbsm/NR)

def sample_generator(i):
    sample2k["sm_sample{0}".format(i)] = SM_gen.rvs(size=np.random.poisson(NR, 1)[0])
    sample2k["bsm_sample{0}".format(i)] = BSM_gen.rvs(size=np.random.poisson(NR * (1 + Geteps), 1)[0])


#start = time.time()
#with cf.ProcessPoolExecutor(max_workers=6) as executor:
#    executor.map(sample_generator, range(1, k))

#end = time.time()
#print("took {0} secs".format(end - start))
k=500
sample20k = {}
for i in range(1,k):
  sample20k["sm_sample{0}".format(i)]= SM_gen.rvs(size=np.random.poisson(NR,1)[0])
  sample20k["bsm_sample{0}".format(i)]= BSM_gen.rvs(size=np.random.poisson(NR*(1+Geteps),1)[0])

sfilename = 'SM_and_BSM_samples_2k_poisson_1_500'
np.save(sfilename+'.npy',sample20k)














