#!/usr/bin/env python

# Above statement tells your terminal that the script should be run in 
# a python enviroment.

# --------------------------------------------------------- #
#  Script illustrating some basic python functionality
#   
#  Authors: Lars Egholm Pedersen (Niels Bohr Institute)
#           Troels C. Petersen (Niels Bohr Institute)
#  Date:    10-11-2017 (latest update)
#
#  Run by: ./PythonIntro.py
#
#  Since output is rather long it can by recommeded to run by
#  ./PythonIntro.py | less 
# --------------------------------------------------------- #


# Tell python which modules it should load. You can now access the 
# functionality of e.g. numpy by np.sqrt( number ) to get the square root etc.

from __future__ import division, print_function
import numpy as np



# -------------------------------------------------------------------------------
# Variable declaration
# -------------------------------------------------------------------------------

# Declare an integer
myint = 1

# Declare a float
myfloat = 1.0

# Notice the difference in the declaration.
# Be carefull of the difference. One common error is if you want to compare 
# your variables to a number:
#   myint == 1       # Is true no matter what
#   myfloat == 1.0   # Could be. myfloat is a number with ??? digits, it could be either
#                      0.9999999999999999999998 or 1.000000000000000001

# Some simple rules
another_float = myint * 2.0        # int comined with float gives float

another_int   = int(myfloat + 0.5) # Rounding floats to int
                                   # When rounding, you will always get the closest int lower
                                   # than the float. Thus if you add 0.5 before rounding, you
                                   # will always get the closest int to the original float

# One important operater you will use time and again is the modulus operator %.
# a % b returns the remainder when deviding the integer a with the integer b.
# So if you wish to find out if 9312312 is divisible by 17, you can check this
# by asking if 93112312 % 17 is zero. More on this below.

# Strings
mystring = "This is a string"      # Use strings to store strings of letters



# -------------------------------------------------------------------------------
# Printing output
# -------------------------------------------------------------------------------

# You can use the print statement to print the output of your script to the terminal:
print("This is an int", myint)
print("This is a float", myfloat)
print(mystring)

print("This is an int", another_int, ". This is a float", another_float)
print("It is also possible (and a good way!) to write",)           # Note that the "," makes the line continue!

print("ints {0}, floats {1} and strings: '{2}' like this".format(another_int, another_float, mystring))
# In the above statement, the numbers in {} each point to the n'th argument in format(). This way you decide
# on the order yourself. You can also re-use an argument multiple times ({0} {0}). If no formatting is chosen,
# Python will print each type clearly (though max 6 decimals for Python 2.7).
# NOTE: You don't even need the numbers. The n'th {} just refers to the n'th argument. (Does not work for Python 2.6)

# or like this:
print("ints {0:d}, floats {1:8.3f} and strings: '{2:s}' like this".format(another_int, 333.333, mystring))
# In this statement, they are now formatted by adding a colon (:) and the formatting string.
# Here you tell the print statement that it is going to fill an integer {:d} (digit), a float with 5 characters,
# where three are after the decimal point {:5.3f} and a string {:s} into the line.
# NOTE: Even if you print e.g. 456.78 with a {:4.2f} statement, it will give 456.78 (not 6.78).

# If you are used to the old notation with a percent sign (%) instead of .format(), you can still do it this way,
# print("A string: %s, an int: %d, a float: %.12f"%("Hello",123,0.1232311231231231312323))
# This also uses % instead of {} for formatting like in other programming languages (e.g. C++).
# NOTE: .format() is much more powerful than % which is why we'll be using it.

print("") # Print an empty line
print("New line \n line starts here \t and is tabulated") 
print("") # Print an empty line


# -------------------------------------------------------------------------------
# Lists 
# -------------------------------------------------------------------------------

# Often you want to store a large number of variables together. Intead of having
# to declare them one by one. The list offers a simple way of organising these:

simplelist = [0.0, 0.0, 0.0, 0.0, 0.0]   # List of 5 floating point zeros
simplelist[2] = 13.0   # Assign the '2' entry to 13.0.

# NOTE: As in most modern programming languages, you count from zero:
# Entry number : 0    1    2    3    4             
# simplelist = [0.0, 0.0, 0.0, 0.0, 0.0]

# If you now want to append some number to your list you can use the append function:
simplelist.append( 1.0 ) 
simplelist.append( 5.0 )
simplelist.append( 3.0 )

# One good thing to remember with python lists is that they are cyclical, which means
# you can also access the numbers in the following way:

# Entry number  -8   -7   -6   -5   -4   -3   -2   -1
# simplelist = [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 5.0, 3.0]

# The reason this is smart is that if you use the append function and want to access
# The last appended value, you can use simplelist[-1] instead of simplelist[len(simplelist)-1],
# where len( list ) is the function that gives you the length of a list.


# A lot of times it is cumbersome to declare each entry in the list by hand.
# If it is for instance necessary to store 100 numbers in your analysis, you
# can declare a list to contain them in the following way (if you are a python
# expert you probably know 10 other ways of doing this).

alotofzeros = [0.0]*100

# In case you want to have have a list of only numbers, you can take advantage of a 
# numpy array. These are more similar to Matlab lists in a lot of ways.

numpyarray = np.array([0., 21., 10., -4., 1/5.])

# numpy arrays allow you to do the following:
# numpyarray_2 = numpyarray_1 + 1
# numpyarray_4 = numpyarray_3*12.4

# this is an easy way to create for example a range:
numpyrange = np.arange(5)

# or other evenly spaced ranges from 0 to 5 in 10 intervals:
numpylinspace = np.linspace(0, 5, 10+1) # can you see why we need a +1 here?

#we can also initialize an empty list:
numpyzeros = np.zeros(10000)

print("My simple list     : ", simplelist) 
print("Which has          : ", len(simplelist), "Entries")
print("The first entry is : ", simplelist[1])
print("The last entry is  : ", simplelist[-1], simplelist[len(simplelist)-1])
print("")
print("A large list   : ", alotofzeros)    # Not so nice output, will be adressed below.
print("")
print("A numpy array   : ", numpyarray)   
print("A numpy range   : ", numpyrange)   
print("A numpy linearly spaced range   : ", numpylinspace)   
print("An array of zeroes   : ", numpyzeros)   
print("")


# A great strenght of the python list is that you can store just about everything you want to
# in them. Say if you want to store a name and a score for some persons : 

mixedlist = [ ["Max"  ,  7 ] , # This is a list of lists, each containing a string and an integer 
              ["John" ,  12] , 
              ["Nick" , -3 ] ,
              ["Niels",  10] ]

# If you want to sort by score you can use the sorting function : 
sortedlist = sorted( mixedlist, key=lambda name : name[1] )
print(sortedlist)

# This is a bit more advanced since, you have to tell python that it should
# sort it by the 'second entry' in each list (i.e. column 1)
# If you only have a list with numbers thelist = [0, -214, 35, 12, 343, 1224],
# you can simply use sorted( thelist ) or thelist.sort()

# If instead you wanted to sort by name and not score, you would do:
sortedlist_v2 = sorted( mixedlist, key=lambda name : name[0] )
print(sortedlist_v2)
print("")

# The printed output is not very nice/readable, but more on that just below!




# -------------------------------------------------------------------------------
# Looping and logical statemants
# -------------------------------------------------------------------------------

# A lot of times you want to repeat a calculation with some different initial parameters
# which is where the loop comes in handy:

# This is a simple loop
for i in range( 10 ) : 
    # This is the block within the loop
    print("Counting to ten : ", i)

# This will tell the intepreter that it should make a list with the numbers 0, ..., 9
# For each of these numbers, it will assign i to the given value and perform whatever
# is inside the 'loop block' with that specific value of i.

# If you want to "control" the loop range more, options are:
for i in range( 4, 10 ) : print("Counting from four to ten : ", i)
for i in range( 4, 11, 2 ) : print("Counting every second number from four to ten : ", i)

print("")
# We can now return to how you can print list more nicely : 
for ientry in mixedlist : # ientry will now be the ith entry in the mixed list,
                          # i.e ientry = "["max", 7] the first time etc.
    print("Name {0:s} \t -> \t Score {1:3d}".format( ientry[0], ientry[1] ))


# If you both want the counter 'i' and the entry value 'ientry' of the list, the enumerate
# function provides a simple way to do this:
print("")
for i, ientry in enumerate( sortedlist ) : 
    print("Counter {0:d} \t Name {1:s} \t -> \t Score {2:3d}".format( i, ientry[0], ientry[1] ))
print("")



# -------------------------------------------------------------------------------
# A note on programming blocks
# -------------------------------------------------------------------------------

# As is the case with programming, you will structure your code as blocks, just like above

# statement1 : 
#     
#    This is coding block1
#
#    statement2 : 
#        
#        This is the coding block2
#
#    ...
#
# Note that all the variables you declare in statement2 is only known in coding block2
# But everything you define in statement1 is also known in statement2
# Most programming languages you have probably seen incloses the statements
# with brackets or curly braces. In python this is done solely by indention 
# (The recommeded pratice is to use four spaces). This means that the interpreter
# knows where e.g. coding block2 is (8 spaces) and when it stops (when you are back
# to 4 spaces).
# To learn more about the python development teams stance on braces, try to import:
#   from __future__ import braces 
# in your terminal interpreter or the very beginning of this program!



# -------------------------------------------------------------------------------
# Back to looping and logical statemants
# -------------------------------------------------------------------------------

# To compare values, you can use the if statement
anint   = 3
afloat  = 314.0
astring = "string"

if anint == 3 : 
    print("The integer value is three ")

if afloat > 1000.0 : 
    print("The float value is greater than 1000")

elif 100.0 < afloat < 500.0 : 
    print("The float value is somewhere between 100. and 500.")

else : 
    print("The float is neither > 1000 nor in the range [100,500]")

if astring != "blah" : 
    print("'astring' does not say 'blah'")

if astring == "string" : 
    print("'astring' says", astring)
print("")


# -------------------------------------------------------------------------------

# One very nice thing with python is you can combine lists, for loops and if statements
# in an intuitive way. If you e.g. want a list containing the first 20 square numbers, you
# can use what is called a list comprehension: 
sqrnumbers = [isqr*isqr for isqr in range( 20 ) ]
for isqr in sqrnumbers :
    print("This is a square number ", isqr)
print("")

# If you only want the even square numbers, you can use the modulus operator : 
evensqrnumbers = [ isqr*isqr for isqr in range( 20 ) if isqr % 2 == 0 ]
for isqr in evensqrnumbers :
    print("This is an even square number", isqr)
print("")


# -------------------------------------------------------------------------------

# Another form of loop is the while statement
whilevar = 0.0
while ( np.tan( whilevar ) < 1000.0 ) : 
    whilevar += 3.1415926535897932385 / 10000
#    whilevar += np.pi / 10000
    
print("{0:6.3f} is the first value whose tangent is greater than 1000 ({1:8.3f}).".format(whilevar, np.tan( whilevar )),)
print("(in steps of pi/10000)")
print("")

# Above while loop test weather tan( whilevar ) is less than 1000.
# If it is greater than 1000 it will exit the loop. In the block
# inside the loop it will increment whilevar with pi / 10000, and 
# run the test again.



# -------------------------------------------------------------------------------
# Random numbers
# -------------------------------------------------------------------------------

# To generate some randoms numbers, we will use the numpys class random, which is both
# very efficient, and produces high quality random numbers (think about what that means!).

# Get the object from the numpy random module
r = np.random
r.seed(42)       # If you don't set any value you get a new (series of) random number(s) every time you run.
                   # With all other seeds, you get the same (series of) random number(s).


# Print some random numbers uniformely distributed between 0 and 1:
for i in range( 10 ) : 
    print("This is a uniformly generated number", r.uniform())   # Calls the uniform function
print("")


r.seed(42)
# Print some random numbers uniformely distributed between 0 and 1:
for i in range( 10 ) : 
    print("This is a uniformly generated number", r.uniform())   # Calls the uniform function
print("")



# Create some lists containing random numbers:

# List with 1000 uniformly distributed random numbers:
uniflist = [r.uniform() for i in range( 1000 ) ] 
uniflist = r.uniform(size=1000)

# List with 1000 unit normally (Gaussianly) distributed random numbers:
gauslist = r.normal(0.0, 1.0, 1000)

# Print a couple of numbers from the list
for i in range( 10 ) : 
    print("Uniform number {0:5.3f} \t Gaussian number {1:5.3f}".format( uniflist[i], gauslist[i] ))
print("")


# -------------------------------------------------------------------------------
# Functions
# -------------------------------------------------------------------------------

# To define a function, use the following structure
#
# def ( input ) : 
#     function block
#     return value
#
# Normally, it is good pratice to put this in the beginning of your script,
# but can be done anywhere, and for the sake of readability we will just do
# it here in this case:

# Define a function that returns the square of a number
def sqr( a ) : 
    return a**2 # or a*a

# Print the square of some of your random numbers : 
for i in range( 10 ) : 
    print("Uniform : {0:5.3f}^2 = {1:5.3f} \t Gauss : {2:6.3f}^2 = {3:5.3f}".format(
            uniflist[i], sqr(uniflist[i]), gauslist[i], sqr(gauslist[i]) ))
print("")

# Print the square of some of your random numbers, now using numpy: 
for i in range( 10 ) : 
    print("Uniform : {0:5.3f}^2 = {1:5.3f} \t Gauss : {2:6.3f}^2 = {3:5.3f}".format(
            uniflist[i], np.square(uniflist[i]), gauslist[i], np.square(gauslist[i]) ))
print("")


# Define a more complicated function that calculates the mean and RMS of the values in a list.

# Can you define it yourself?

def MeanAndRMS( inlist = [] ) : 

    if len( inlist ) == 0  : 
        print("Ups, called function with an empty list")
        print("If I don't exit now, I will divide by zero when calculating the mean")
        return [-9999999999.9, -99999999999.9]

    elif len( inlist ) == 1  : 
        print("Ups, called function with a list of length one")
        print("If I don't exit now, I will divide by zero when calculating the RMS")
        return [inlist[0], -99999999999.9]

    # Find the mean:
    mu = 0.0

    # Find the standard deviation (RMS):
    rms = 0.0

    return [mu, rms]                           # Return a list with two values: [mean, rms]


# Calculate the mean and rms of the values in your two lists with random numbers :
mu_sigma_unif = MeanAndRMS( uniflist )
mu_sigma_gaus = MeanAndRMS( gauslist )

# Print the results. Do these values make sense?
print("For 1000 uniformely distributed numbers:   mu = {0:5.3f}   rms = {1:5.3f}".format(mu_sigma_unif[0], mu_sigma_unif[1]))
print("For 1000 Gaussianly distributed numbers:   mu = {0:5.3f}   rms = {1:5.3f}".format(mu_sigma_gaus[0], mu_sigma_gaus[1]))
print("Is this what you would expect? \n")


# -------------------------------------------------------------------------------
# Writing and reading files
# -------------------------------------------------------------------------------

# Often you will have to read and write ascii files (i.e. human readable text files)
# Start with writing a file randomnumbers.dat containing some random numbers:

# Open 'randomnumbers.dat' in write mode. Note that the file is only open in the
# block opened by the 'with ... as file' statement.
with open( 'randomnumbers.dat', 'w' ) as file : 

    for i in range( 25 ) : 

        # Numbers distributed according to the following pdfs:
        #   uniform, gaussian, poissonian, breitwigner, landau
        outline = "{0:9.4f} \t {1:9.4f} \t {2:9.0f} \t {3:9.4f} \t {4:9.4f} \n".format(r.uniform(), r.normal(), r.poisson( 10.0 ), r.exponential(), r.power(1) )
        print(outline,)

        file.write( outline )


# Now we will read the same file in again.

# Define lists that can contain the table of random numbers:
unif = []
gaus = []
pois = []
exp = []
power = []

with open( 'randomnumbers.dat', 'r' ) as file : 
    
    # Loop through each line in the file
    for iline in file : 

        # Now we have the whole line in string format, but want it back in floating point format. 
        # The command below does the following:
        # Take iline and strip it of '\n and \t' (i.e. newlines and tabulations)
        # '1.0 \t 2.0 \t 3.0 \n' -> '1.0    2.0    3.0'
        # (strip('somestring') will in general remove 'somestring' but default is '\n \t'
        # Take the result of this and split it into a list :
        # '1.0    2.0    3.0' -> ['1.0', '2.0', '3.0']
        # Take this an convert it to floating point numbers
        # ['1.0', '2.0', '3.0'] -> [1.0, 2.0, 3.0]

        # format_line = [float( ientry ) for ientry in iline.strip().split() ]

        strip_line = iline.strip()
        format_line = strip_line.split()

        unif.append( float(format_line[0]) )        # Now "unif" gets the first number in the line.
        gaus.append( float(format_line[1]) )        # Now "gaus" gets the second numb...
        pois.append( float(format_line[2]) )
        exp.append( float(format_line[3]) )
        power.append( float(format_line[4]) )

#could also be done in Numpy by:
#unif, gaus, pois, exp, power = np.loadtxt('randomnumbers.dat', unpack=True)

# To check that the above worked, print the lines 10 to 19: 
print("\nRandom numbers (range 10-19 read from file): ")
for i in range( 10, 20 ) : 
    print("{0:8.4f} \t {1:8.4f} {2:8.4f} \t {3:8.4f} \t {4:8.4f}".format(unif[i], gaus[i], pois[i], exp[i], power[i] ))
print("\n")


# Lastly create a latex table you can copy paste into your paper/thesis
# with the mean and rms of the numbers you have genereated:
mu_sigma_unif = MeanAndRMS( unif )
mu_sigma_gaus = MeanAndRMS( gaus )
mu_sigma_pois = MeanAndRMS( pois )
mu_sigma_exp = MeanAndRMS( exp )
mu_sigma_power = MeanAndRMS( power )

print("\\begin{tabular}{|l|cc|}")
print("\t \\hline")
print("\t \\multicolumn{3}{|c|}{PDF Means and RMS}")
print("\t \\hline")
print("\t {0:15s} & {1:7s} & {2:7s}\\\\".format( "pdf", "$\\mu$", "$\\sigma$" ))
print("\t \\hline")
print("\t {0:15s} & {1:7.3f} & {2:7.3f} \\\\".format( "Uniform"    , mu_sigma_unif[0],  mu_sigma_unif[1]  ))
print("\t {0:15s} & {1:7.3f} & {2:7.3f} \\\\".format( "Gaussian"   , mu_sigma_gaus[0],  mu_sigma_gaus[1]  ))
print("\t {0:15s} & {1:7.3f} & {2:7.3f} \\\\".format( "Poissonian" , mu_sigma_pois[0],  mu_sigma_pois[1]  ))
print("\t {0:15s} & {1:7.3f} & {2:7.3f} \\\\".format( "Exponential", mu_sigma_exp[0],   mu_sigma_exp[1]   ))
print("\t {0:15s} & {1:7.3f} & {2:7.3f} \\\\".format( "Power"      , mu_sigma_power[0], mu_sigma_power[1] ))
print("\t \\hline")
print("\\end{tabular}")

print("")

# Done