import os
import numpy as np
from scipy.optimize import minimize
outputfile='temporary.csv'


# initial parameters
sm=0.95
fsn=0.5
fss=0.5
fsv=0.5

x0=[sm, fsn, fss, fsv]

# parameters obey the following bounds
#bnds = ((0, 1), (0, 1), (0, 1))

#bnds = ((0.5, 1.0), (0.1, 0.9), (0.1, 1.5), (0.1, 0.5))
bnds = ((0.25, 1.0), (0.1, 1.5), (0.1, 1.5), (0.1, 0.75))

# possible target quantities
# --------------------------
theta0=1 #38
# region deformation
am0=1.672631   # Mesophyl
an0=1.608788   # Nectary
as0=1.790434   # Side
av0=1.294615   # Vasculature
# strain anisotropy
anim=0.6247   # Mesophyl
anin=0.6172   # Nectary
anis=0.5489   # Side
aniv=0.6958   # Vasculature

# experimental standard deviation
stdtheta=0.21 #38
stdam=0.23   # Mesophyl
stdan=0.11   # Nectary
stdas=0.21   # Side
stdav=0.06   # Vasculature
#
stdanim=0.0532   # Mesophyl
stdanin=0.0378   # Nectary
stdanis=0.0741  # Side
stdaniv=0.0333   # Vasculature

a0=[theta0, am0, an0, as0, av0, anim, anin, anis, aniv]
stda=[stdtheta, stdam, stdan, stdas, stdav, stdanim, stdanin, stdanis, stdaniv]

# target quantities included in the objective function
var=np.array([1, 2, 3, 4]) # here theta is not included
#var=np.array([0, 1, 2, 3, 4]) # here theta (side angle) is included
#var=np.array([0, 1, 2, 3, 4, 5, 6, 7, 8]) # side angle, deformation and anisotropy are included
#var=np.array([5, 6, 7, 8]) # side angle, deformation and anisotropy are included

'''
# objectiv function 1
def f(x):
	res=0
	for i in range(len(x)):
		res+=(x[i]-a0[var[i]])**2/(x[i]+a0[var[i]])**2
	return res
'''

# objectiv function taking into account the experimental standard deviation for the weight
def f(x):
	res=0
	for i in range(len(x)):
		res+=(x[i]-a0[var[i]])**2/(stda[var[i]])**2
	return res



def objective(x):
	[sm, fsn, fss, fsv] = x
	os.system("FreeFem++ pulvinus4optimize.cpp -fsn "+str(fsn)+" -fss "+str(fss)+" -fsv "+str(fsv)+" -sm "+str(sm)+" -outfile "+outputfile)
	data=np.loadtxt(outputfile, delimiter=';', dtype=float)
	xvar=data[var+4]
	print(xvar)
	return f(xvar)

#
# GLOBAL OPTIMIZATION
#

from scipy.optimize import differential_evolution #New in version 0.15.0.

res = differential_evolution(objective, bounds=bnds)
res


'''

###############################################################
###############################################################
# fitting the observed expansions        K=lambda+2*mu/3 (d=3 formula)
###############################################################
###############################################################

Ok: Normal End
[1.67263 1.60879 1.79043 1.29461]
Out[2]: 
     fun: 7.656738339961133e-09
 message: 'Optimization terminated successfully.'
    nfev: 4745
     nit: 78
 success: True
       x: array([0.54006051, 0.95448193, 1.21647743, 0.50956303])


Results in the results file:

/*
sm;fsn;fss;fsv;sideangle;am;an;as;av;sam;san;sas;sav;
*/
-0.540061;0.954482;1.21648;0.509563;0.477428;1.67263;1.60879;1.79043;1.29461;0.963544;0.973407;0.942094;0.935191


###############################################################
###############################################################
# fitting the observed expansions        K=lambda+mu (d=2 formula)
###############################################################
###############################################################
Ok: Normal End
[1.67263 1.60879 1.79043 1.29461]
Out[3]: 
     fun: 7.656738339961133e-09
 message: 'Optimization terminated successfully.'
    nfev: 5825
     nit: 96
 success: True
       x: array([0.46445277, 0.95448186, 1.21647529, 0.5095638 ])


Results in the results file:

/*
sm;fsn;fss;fsv;sideangle;am;an;as;av;sam;san;sas;sav;
*/
-0.464453;0.954482;1.21648;0.509564;0.477427;1.67263;1.60879;1.79043;1.29461;0.963544;0.973407;0.942095;0.935191



'''