MAJ terminée. Nous sommes passés en version 14.6.2 . Pour consulter les "releases notes" associées c'est ici :

https://about.gitlab.com/releases/2022/01/11/security-release-gitlab-14-6-2-released/
https://about.gitlab.com/releases/2022/01/04/gitlab-14-6-1-released/

Commit 55080b74 authored by HERASIMENKA Alesia's avatar HERASIMENKA Alesia
Browse files

cleaning #2

parent 5a00bd8e
%% Cell type:markdown id: tags:
## Initializations
%% Cell type:code id: tags:
``` python
import numpy as np
import scipy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
from nutopy import nle
from nutopy import tools
from nutopy import ocp
from nutopy import path
from scipy.linalg import null_space
from gve import *
from os import system
# Definition of optical parameters
rho = 0.6 # Specular reflection coefficient
a = 1. - rho # Absorption coefficient
s = 0.94 # Diffuse reflection coefficient
eps_f = 0.05 # Emissivity front coeff 0.05
eps_b = 0.55 # Emissivity back coeff 0.55
Bf = 0.79 # Front Lambertian coeff 0.79
Bb = 0.55 # Back Lambertian coeff 0.55
eps = (Bf * eps_f - Bb * eps_b) / (eps_b + eps_f)
# Initial guess of pI taken from Matlab Convex Programming
pI = np.array([ - 0.1577, -0.3003, -0.1011, -0.2680, 1.0000])
M1_guess = 1.2252
M2_guess = 5.0039
I1_guess = np.array([-0.9623, 0.1885, 0.6944, -0.0214, 0.0558])
I2_guess = np.array([-0.9633, 0.1899, 0.6946, -0.0220, 0.0559])
z = np.hstack((pI, M1_guess, M2_guess, I1_guess, I2_guess))
def deg2rad(x):
y = x * np.pi / 180.
return y
# Initial state
cont = 0. # 0 for triangular cone, 1 for real drop
mu = 1.
I = np.array([deg2rad(10), deg2rad(20), deg2rad(30), 1., 0.5]) # Initial state
d = np.array([ 0., 0., 0., 0., 1.]) # Direction of the displacement
b = np.array([1. - rho * s, 2. * rho * s, Bf * rho * (1. - s) + (1. - rho) * eps ])
I0 = np.array([0., 0., 0., 0., 0]) # delta I (0)
M0 = 0
Mf = 2*np.pi
# Real cone force definition
cBeta = (- b[0] * b[2] - 2. * b[1] * b[2] + \
np.sqrt(b[0] ** 2 * b[2] ** 2 - 4. * b[0] * b[2] ** 2 * b[1] + 8. * b[0] ** 2 * b[1] ** 2 + 4. * b[1] ** 3 * b[0])) / \
(4. * b[0] * b[1] + 2. * b[1] ** 2);
sBeta = np.sqrt(1 - cBeta ** 2);
fs = b[0] * cBeta + (b[1] * cBeta ** 2 + b[2] * cBeta) * cBeta;
fperp = (b[1] * cBeta ** 2 + b[2] * cBeta) * sBeta;
fCone = np.array([fs, fperp])
pars = np.hstack((mu, I, fCone, b))
d_orth = scipy.linalg.null_space(np.array([[ 0., 0., 0., 0., 1.]])) # Orthogonal space of b
d_orth = d_orth.transpose()
```
%% Cell type:markdown id: tags:
## Hamiltonian \& flow
%% Cell type:code id: tags:
``` python
!python -m numpy.f2py -c hfun.f90 hfun.pyf -m hfun > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d.f90 hfun_d.pyf -m hfun_d > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d_d.f90 hfun_d_d.pyf -m hfun_d_d > /dev/null 2>&1
```
%% Cell type:code id: tags:
``` python
from hfun import hfun
from hfun_d import hfun_d
from hfun_d_d import hfun_d_d
hfun_u = lambda M, q, p, pars, cont : hfun.hfun_u(M, q, p, pars, cont)
dhfun_u = lambda M, dM, q, dq, p, dp, pars, cont, dcont : hfun_d.hfun_u_d(M, dM, q, dq, p, dp, pars, cont, dcont)
d2hfun_u = lambda M, dM, d2M, q, dq, d2q, p, dp, d2p, pars, cont, dcont, d2cont : hfun_d_d.hfun_u_d_d(M, dM, d2M, q, dq, d2q, p, dp, d2p, pars, cont, dcont, d2cont)
```
%% Cell type:code id: tags:
``` python
def dhfun_0(M, dM, q, dq, p, dp, pars, cont, dcont):
h = 0.0
dh = 0.0
return h, dh
def d2hfun_0(M, dM, d2M, q, dq, d2q, p, dp, d2p, pars, cont, dcont, d2cont):
h = 0.0
dh = 0.0
d2h = 0.0
return h, dh, d2h
def hfun_0(M, q, p, pars, cont):
h = 0.0
return h
```
%% Cell type:code id: tags:
``` python
hfun_u = tools.tensorize(dhfun_u, d2hfun_u, tvars=(1, 2, 3, 5), full=True)(hfun_u)
hfun_0 = tools.tensorize(dhfun_0, d2hfun_0, tvars=(1, 2, 3, 5), full=True)(hfun_0)
Hu = ocp.Hamiltonian(hfun_u)
H0 = ocp.Hamiltonian(hfun_0)
fu = ocp.Flow(Hu)
f0 = ocp.Flow(H0)
```
%% Cell type:markdown id: tags:
## Shooting function
%% Cell type:code id: tags:
``` python
def dshoot(z, dz, pars, d, d_orth, cont, dcont):
I0 = np.zeros(5)
M0 = 0.
Mf = 2. * np.pi
pI = z[0 : 5]
M1 = z[5]
M2 = z[6]
I1 = z[7 : 12]
I2 = z[12 : 17]
dpI = dz[0 : 5]
dM1 = dz[5]
dM2 = dz[6]
dI1 = dz[7 : 12]
dI2 = dz[12 : 17]
dM0 = 0.
dI0 = np.zeros(5)
dMf = 0.
#----------------------------------
(I1sol, dI1sol), (pI1, dpI1) = fu((M0, dM0), (I0, dI0), (pI, dpI), (M1, dM1), pars, cont)
(I2sol, dI2sol), (pI2, dpI2) = f0((M1, dM1), (I1, dI1), (pI, dpI), (M2, dM2), pars, cont)
(If, dIf), (pIf, dpIf) = fu((M2, dM2), (I2, dI2), (pI, dpI), (Mf, dMf), pars, cont)
#----------------------------------
s = np.zeros(17)
s[0] = pIf[0] * d[0] + pIf[1] * d[1] + pIf[2] * d[2] + pIf[3] * d[3] + pIf[4] * d[4] - np.linalg.norm(d)**2
s[1 : 5] = np.dot(d_orth, If)
s[7 : 12] = I1 - I1sol
s[12 : 17] = I2 - I2sol
ds = np.zeros(17)
ds[0] = dpIf[0] * d[0] + dpIf[1] * d[1] + dpIf[2] * d[2] + dpIf[3] * d[3] + dpIf[4] * d[4]
ds[1 : 5] = np.dot(d_orth, dIf)
ds[7 : 12] = dI1 - dI1sol
ds[12 : 17] = dI2 - dI2sol
s[5], ds[5] = Hu((M1, dM1), (I1, dI1), (pI, dpI), pars, cont)
s[6], ds[6] = Hu((M2, dM2), (I2, dI2), (pI, dpI), pars, cont)
# on impose pas la continuité de pI car il est constant
# il faut imposer la continuité de l'etat deltaI entre les switchs entre 7 et 17
return s, ds
@tools.tensorize(dshoot, tvars=(1,5), full=True)
def shoot(z, pars, d, d_orth, cont):
I0 = np.zeros(5)
M0 = 0.
Mf = 2. * np.pi
pI = z[0 : 5]
M1 = z[5]
M2 = z[6]
I1 = z[7 : 12]
I2 = z[12 : 17]
#----------------------------------
I1sol, pI1 = fu(M0, I0, pI, M1, pars, cont)
I2sol, pI2 = f0(M1, I1, pI, M2, pars, cont)
If, pIf = fu(M2, I2, pI, Mf, pars, cont)
#----------------------------------
s = np.zeros(17)
s[0] = pIf[0] * d[0] + pIf[1] * d[1] + pIf[2] * d[2] + pIf[3] * d[3] + pIf[4] * d[4] - np.linalg.norm(d)**2
s[1 : 5] = np.dot(d_orth, If)
# theta1 = atan2(pI1*sPerp/(pI1*sDir))
# theta2 = atan2(pI2*sPerp/(pI2*sDir))
# s[5] = theta1 - cone_alpha - deg2rad(90.)
# s[6] = theta2 - cone_alpha - deg2rad(90.)
s[5] = Hu(M1, I1, pI, pars, cont)
s[6] = Hu(M2, I2, pI, pars, cont)
s[7 : 12] = I1 - I1sol
s[12 : 17] = I2 - I2sol
# on impose pas la continuité de pI car il est constant
# il faut imposer la continuité de l'etat deltaI entre les switchs
return s
```
%% Cell type:markdown id: tags:
## Solve initial shooting
%% Cell type:code id: tags:
``` python
dfoo = lambda z, dz, cont: shoot((z, dz), pars, d, d_orth, (cont, 0.))
foo = lambda z, cont: shoot(z, pars, d, d_orth, cont)
foo = tools.tensorize(dfoo, tvars=(1,), full=True)(foo)
```
%% Cell type:code id: tags:
``` python
nleopt = nle.Options(SolverMethod='hybrj', Display='on', TolX=1e-8)
et = time.time(); sol = nle.solve(foo, z, df=foo, options=nleopt, args=cont);
#et = time.time(); sol = nle.solve(foo, z, options=nleopt, args=cont);
```
%% Cell type:code id: tags:
``` python
z_sol = sol.x
print('z_sol =', z_sol)
print('foo =', foo(z_sol, cont))
```
%% Cell type:markdown id: tags:
## Function of continuous integration
%% Cell type:code id: tags:
``` python
# Definition of the homotopic function and its first order derivative
# This function is used to solve S=0 for different values of e=epsilon and tf.
dhomfun = lambda z, dz, cont, dcont: shoot((z, dz), pars, d, d_orth, (cont, dcont))
homfun = lambda z, cont: shoot(z, pars, d, d_orth, cont)
homfun = tools.tensorize(dhomfun, tvars=(1,2), full=True)(homfun)
#def dhomfun(z, dz, lam, dlam):
# s, ds_z = foo((z, dz), lam)
#
# dx = 1.0e-8
# ds_lam = (foo(z, lam - dx) - foo(z, lam + dx)) / dx / 2 * dlam
#
# ds = ds_z + ds_lam
# return s, ds
#@tools.tensorize(dhomfun, tvars=(1, 2), full=True)
#def homfun(z, lam):
# s = foo(z, lam)
# return s
homfun((z_sol, np.zeros(17)), (cont, 1.))
```
%%%% Output: execute_result
(array([ 0.00000000e+00, 1.90181874e-01, 6.96537407e-01, -2.24707161e-02,
9.66402564e-01, -5.57283109e-11, -6.48865406e-11, -9.66402564e-01,
1.90181874e-01, 6.96537407e-01, -2.24707161e-02, 5.64359354e-02,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00]),
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]))
%% Cell type:markdown id: tags:
## Homotopic function
%% Cell type:code id: tags:
``` python
# Making the penalization smaller: homotopy on e
z0 = z_sol
sol_path = path.solve(homfun, z0, 0., 1., df=homfun)
#sol_path = path.solve(homfun, z0, 0., 1.)
```
%% Cell type:code id: tags:
``` python
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment