MultiOpt Namespace Reference

Functions
def	ComputeInertiaProfiles (nodal_values, duration, N_t)
def	RealFun (x)
def	ToyFun (x)
def	main ()

Variables
string	SourceDir = "/home/iostanin/Desktop/MercuryGit/mercurydpm"
string	BuildDir = "/home/iostanin/Desktop/MercuryGit/MercuryBuild"
float	Timestep = 5.56833e-05 * 50
int	N = 5
bool	TEST_MODE = False
dictionary	TEST_OPTIONS = {'xtol':1e-3, 'ftol':1e-4, 'maxiter':3}
dictionary	REAL_OPTIONS = {'xtol':1e-3, 'ftol':1e-3, 'maxiter':20}
float	QMin = 0.5
float	QMax = 1.5
int	I0 = 10
int	ProgDuration = 100
int	SymDuration = 105
int	BaseAngvel = 5
int	NTimesteps = ProgDuration / Timestep
	QInitValues = np.array([1, 1])
	QFinalValues = np.array([1, 1])
	x0 = np.ones(2*N)
	xc = np.random.rand(3*N)
int	delta = np.pi/1000.
	InitialTheta = np.arccos(1/3**0.5)
int	InitialPhi = np.pi/4
int	FinalTheta = delta
int	FinalPhi = delta
int	FunCount = 0
bool	FunBatch = False

Function Documentation

ComputeInertiaProfiles()

def MultiOpt.ComputeInertiaProfiles	(		nodal_values,
			duration,
			N_t
	)

def ComputeInertiaProfiles(nodal_values, duration, N_t):
    # pre-computes time evolution of tensor of inertia for MercuryDPM driver file
    # nodal_values is an array N X 3  equispaced values of three principal moments of inertia
    N = len(nodal_values)//3
    n_I1 = nodal_values[:N]
    n_I2 = nodal_values[N:2*N]
    n_I3 = nodal_values[2*N:]
    x = np.arange(0, duration * (N/(N-1)), duration/(N-1))
    cs1 = CubicSpline(x, n_I1, bc_type = 'clamped')
    cs2 = CubicSpline(x, n_I2, bc_type = 'clamped')
    cs3 = CubicSpline(x, n_I3, bc_type = 'clamped')
    xc = np.arange(0, duration * (N_t/(N_t-1)), duration/(N_t-1))
    I1 = cs1(xc); dI1 = cs1(xc, 1)
    I2 = cs2(xc); dI2 = cs2(xc, 1)
    I3 = cs3(xc); dI3 = cs3(xc, 1)
    return x, n_I1, n_I2, n_I3, xc, I1, dI1, I2, dI2, I3, dI3
 

Referenced by RealFun().

main()

def MultiOpt.main

(

)

def main():
    global SourceDir, BuildDir, initial_phi, InitialTheta, FinalPhi, FinalTheta
 
    # Set up terminal output
    clr = colorClass()
    clr.set_use_colors(True)
    clr.def_colors()
 
    # Hello tag
    print(clr.END + clr.BOLD + clr.VIOLET + "-------------------------------------------")
    print(clr.END + clr.BOLD + clr.VIOLET + "     MercuryDPM MultiOpt experiment        ")
    print(clr.END + clr.BOLD + clr.VIOLET + "-------------------------------------------")
 
 
    # If the values for sourceDir and Build dir are provided in the command, then overwrite the default values
    if (len(sys.argv) == 3):  # Full format: "run source_dir build_dir"
        SourceDir = sys.argv[1]
        BuildDir = sys.argv[2]
    print("Source dir: ", SourceDir)
    print("Build dir: ", BuildDir)
 
    # Path to optimization directory
    OptPath = BuildDir + "/Drivers/Clump/ChangingTOI/opt"
    if ( not os.path.exists( OptPath )):
        os.mkdir( OptPath )
    
    # Let Mercury code know what is the goal orientation
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/init_orientation.txt', (InitialTheta, InitialPhi), delimiter =' ')
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/final_orientation.txt', (FinalTheta, FinalPhi), delimiter =' ')
 
    print(clr.END)
 
    if (TEST_MODE):
        res = minimize(ToyFun, x0, method='Powell', options = TEST_OPTIONS)
    else:
        res = minimize(RealFun, x0, method='Powell', options = REAL_OPTIONS)
 
        # Alternative optimization tools
        #res = minimizeSPSA(real_fun, x0=x0, niter=1000, paired=False)
        #res = minimizeCompass(real_fun, x0=x0, deltatol=0.01, paired=False)
 
    print(res)
 
    if (TEST_MODE): print(clr.END + "Converged to the known solution: ", np.allclose(res["x"], xc))
 
# Main script that runs the maneuver optimization procedure

References Eigen::internal.print().

RealFun()

def MultiOpt.RealFun

(

x

)

def RealFun(x):
    import os
    import subprocess
    global SourceDir, BuildDir, N
    global ProgDuration, SymDuration, NTimesteps, FunCount, FunBatch
    global QMax, QMin, I0, QInitValues, QFinalValues
 
    ret = 0 # Return zero by default
 
    # Project the unconstrained vector of optimization parameters x into (QMin,QMax)
    q = (QMax + QMin) / 2. + (QMax - QMin) / 2. * np.cos(x)
 
    # Expand the set of control parameters with initial and final values
    q_ext_1 = np.hstack((QInitValues[0], q[:N], QFinalValues[0]))
    q_ext_2 = np.hstack((QInitValues[1], q[N:], QFinalValues[1]))
 
    # Compute principal components of inertia
    I_1 = 0.5 * I0 * (1 + q_ext_2*q_ext_2)
    I_2 = 0.5 * I0 * (1 + q_ext_1*q_ext_1)
    I_3 = 0.5 * I0 * (q_ext_1*q_ext_1 + q_ext_2*q_ext_2)
    nodal_values = np.hstack((I_1, I_2, I_3))
    print("Nodal values of principal moments of inertia:", nodal_values)
 
    # Compute cubic splines approximating the evolution of TOI
    xx, n_I1, n_I2, n_I3, xc, I1, dI1, I2, dI2, I3, dI3 = ComputeInertiaProfiles(nodal_values, ProgDuration, NTimesteps)
 
    arr = np.transpose(np.vstack((xc, I1, I2, I3, dI1, dI2, dI3)))
 
    # Expanded evolution incluidng the final piece with constant TOI
    I1_ext = np.hstack((I1, I1[len(I1)-1] * np.ones((int) (NTimesteps * (SymDuration - ProgDuration) / SymDuration))))
    I2_ext = np.hstack((I2, I2[len(I2)-1] * np.ones((int) (NTimesteps * (SymDuration - ProgDuration) / SymDuration))))
    I3_ext = np.hstack((I3, I3[len(I3)-1] * np.ones((int) (NTimesteps * (SymDuration - ProgDuration) / SymDuration))))
    T_I = np.arange(0, SymDuration, SymDuration / len(I1_ext)) # Timestamp
 
    # Save inertia control parameters of the simulation
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/inertia_profiles.txt', arr)
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/I_1.txt', I1_ext)
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/I_2.txt', I2_ext)
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/I_3.txt', I3_ext)
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/T_I.txt', T_I)
    np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/nodal_values.txt', nodal_values)
 
    # Run Mercury Driver file
    # (feed program duration, simulation duration, and initial angular velocity
    # as command line arguments)
    os.chdir(BuildDir + '/Drivers/Clump/ChangingTOI')
    command  = [BuildDir + '/Drivers/Clump/ChangingTOI/ChangingTOI']
    command.append('-p1')
    command.append(str(ProgDuration))
    command.append('-p2')
    command.append(str(SymDuration))
    command.append('-p3')
    command.append(str(BaseAngvel))
    print(command)
    subprocess.run(command)
 
    # Get the functional value
    f = open(BuildDir + '/Drivers/Clump/ChangingTOI/opt/functional.txt', "r")
    ret = float(f.read())
 
    # Save the functional value in the history batch
    if (FunCount == 0):
        FunBatch = ret
    else:
        FunBatch = np.append(FunBatch, ret)
        np.savetxt(BuildDir + '/Drivers/Clump/ChangingTOI/opt/fun_batch.txt', FunBatch)
 
    FunCount+=1 # Number of functional evaluations
    return ret
 
# Toy functional necessary for testing purposes

References ComputeInertiaProfiles(), Eigen::internal.print(), and compute_granudrum_aor.str.

ToyFun()

def MultiOpt.ToyFun

(

x

)

def ToyFun(x):
    global xc
    ret = 0
    for m in range(len(x)):
        ret = ret + (x[m] - xc[m])**2
    return ret
 
 

Variable Documentation

BaseAngvel

int MultiOpt.BaseAngvel = 5

BuildDir

string MultiOpt.BuildDir = "/home/iostanin/Desktop/MercuryGit/MercuryBuild"

delta

int MultiOpt.delta = np.pi/1000.

Referenced by Eigen::internal.bicgstabl(), Eigen::internal::TensorIntDivisor< int32_t, true >.calcMagic(), Eigen::internal::MatrixFunctionAtomic< MatrixType >.compute(), ChuteWithPeriodicInflow.computeInternalForces(), ScrewsymmetricIntersectionOfWalls.computeNormalRadialDeltaN(), Eigen::internal.convertArrayPointerBF16toF32Dup(), Eigen::internal.dogleg(), Eigen::internal.doubleword_div_fp(), Eigen::internal.erfc_double_large(), BaseInteraction.getOverlapVolume(), Eigen::internal.lmpar(), Eigen::internal.lmpar2(), Eigen::internal.loadBF16fromResult(), Eigen::LevenbergMarquardt< FunctorType_ >.minimizeOptimumStorageOneStep(), MandelbrotWidget.mouseMoveEvent(), HopperInsertionBoundary.placeParticle(), Eigen::internal::generic_fast_erfc< Scalar >.run(), UnstructuredImmersedEllipseProblem< ELEMENT >.set_boundary_velocity(), oomph::TetMeshBase.setup_boundary_coordinates(), oomph::XdaTetMesh< ELEMENT >.setup_boundary_coordinates(), ChargedBondedParticleUnitTest.setupInitialConditions(), oomph::TetMeshBase.snap_to_quadratic_surface(), Eigen::HybridNonLinearSolver< FunctorType, Scalar >.solveNumericalDiffOneStep(), Eigen::HybridNonLinearSolver< FunctorType, Scalar >.solveOneStep(), CGCoordinates.spaceEvenly(), storeBF16fromResult(), test_scalar_sugar_sub_div(), CFile.writeP4C(), and CFile.writeP4W().

FinalPhi

int MultiOpt.FinalPhi = delta

FinalTheta

int MultiOpt.FinalTheta = delta

FunBatch

bool MultiOpt.FunBatch = False

FunCount

int MultiOpt.FunCount = 0

I0

int MultiOpt.I0 = 10

InitialPhi

int MultiOpt.InitialPhi = np.pi/4

InitialTheta

MultiOpt.InitialTheta = np.arccos(1/3**0.5)

N

int MultiOpt.N = 5

NTimesteps

int MultiOpt.NTimesteps = ProgDuration / Timestep

ProgDuration

int MultiOpt.ProgDuration = 100

QFinalValues

MultiOpt.QFinalValues = np.array([1, 1])

QInitValues

MultiOpt.QInitValues = np.array([1, 1])

QMax

float MultiOpt.QMax = 1.5

QMin

float MultiOpt.QMin = 0.5

REAL_OPTIONS

dictionary MultiOpt.REAL_OPTIONS = {'xtol':1e-3, 'ftol':1e-3, 'maxiter':20}

SourceDir

string MultiOpt.SourceDir = "/home/iostanin/Desktop/MercuryGit/mercurydpm"

SymDuration

int MultiOpt.SymDuration = 105

TEST_MODE

bool MultiOpt.TEST_MODE = False

TEST_OPTIONS

dictionary MultiOpt.TEST_OPTIONS = {'xtol':1e-3, 'ftol':1e-4, 'maxiter':3}

Timestep

float MultiOpt.Timestep = 5.56833e-05 * 50

x0

MultiOpt.x0 = np.ones(2*N)

xc

MultiOpt.xc = np.random.rand(3*N)

Referenced by SCoupling< M, O >.computeSCouplingForcesFromTriangles(), and VolumeCoupling.computeWeightOnParticles().