Laboratoire RDP » Biophysics Team » Dandelion pappus morphing

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

'''