#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jun  5 14:41:33 2018

@author: mirjam
"""
import numpy as np

def linmin(f,a,c,epsilon=1e-5):     
    
    r = (-1+np.sqrt(5))/2 #golden-cut position (~61.8%)
    # Compute points at 38.2 and 61.8% of [a,c] interval
    b = a + (1-r)*(c-a) 
    d = a + r*(c-a)
    
    while abs(d-b) > epsilon:
        if f(b) < f(d):
            c = d
        else:
            a = b
        
        #Recompute b, d
        b = a + (1-r)*(c-a)
        d = a + r*(c-a)

    return (c+a)/2


def f(x): 
    return(x**2)

print("linmin(f,-1,1): "+str(linmin(f,-1,1)))
print("linmin(f,-2,100): "+str(linmin(f,-2,100)))
print("linmin(f,10,20): "+str(linmin(f,10,20)))
print("linmin(f,-20,-10): "+str(linmin(f,-20,-10)))

raw_input("Press Enter to continue...")

#######################################################################
# part 2
#######################################################################

print("--- part 2 ---")

def my_paratrafo(x = np.array([1])):
    f = 1./(1+np.exp(-1*x))
    p = f/np.sum(f)
    return(p)

def entropy(x):
    p = my_paratrafo(x)
    entropy = np.sum(p*np.log(p));
    return(entropy)
    

# Copy this loop into your solution! From here until the end

# Number of dimensions
N = 10
# Draw intial probabilities and normalize them
x0 = np.random.rand(N)
p0 = my_paratrafo(x0)

print("start:")
print(np.around(p0,decimals=3))

p_old = p0
p_new = 1e5
# Matrix with initial unit vectors as columns
U = np.diag(np.ones(N))
counter = 0

# Here, the actual loop of Powell's method starts
while np.linalg.norm(p_old - p_new) > 1e-5:
    counter = counter + 1
    x0_old = x0
    # Normalize the current probabilities
    p_old = my_paratrafo(x0)
    i_del = 0
    norm_del = 0
    # loop over the current vectors ui and determine the minimum of lambda for each x0+lambda*ui
    for i in range(0,N-1):
        ui = U[:,i]
        def entropy_fun(mylambda):
            x_lambda = x0 + mylambda*ui
            return(entropy(x_lambda))
        mylambda = linmin(entropy_fun, -10, 10)
        x0 = x0 + mylambda*ui
        # Save i if mylambda * norm(ui) the largest (this one will be deleted afterwards)
        if mylambda*np.linalg.norm(ui)>norm_del:
            i_del = i
        
    p_new = my_paratrafo(x0)
    
    # Powells method: update search vectors by deleting one and adding the x0 direction
    x_dir = x0-x0_old
    #U = np.append(np.delete(U,0,1),np.reshape(x_dir, (N,1)), axis=1)
    U = np.append(np.delete(U,i_del,1),np.reshape(x_dir, (N,1)), axis=1)
    
    # go once more in direction of the new search vector
    def entropy_fun(mylambda):
        x_lambda = x0 + mylambda*x_dir
        return(entropy(x_lambda))
    mylambda = linmin(entropy_fun, -10, 10)
    x0 = x0 + mylambda*x_dir

# Normalize and print the final result
p_final = my_paratrafo(x0)

print("result after " + str(counter) + " iterations:")
print(p_final)
print("result rounded:")
print(np.around(p_final,decimals=3))