Skip to content

Commit

Permalink
fix: objective function is not none
Browse files Browse the repository at this point in the history
  • Loading branch information
rabii-chaarani committed May 22, 2024
1 parent 8c686be commit c2f143c
Show file tree
Hide file tree
Showing 12 changed files with 9,799 additions and 12,638 deletions.
2 changes: 1 addition & 1 deletion FoldOptLib/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
InputGeologicalKnowledge,
)
from FoldOptLib.utils import utils
from FoldOptLib.input import CheckInputData, InputDataProcessor
from FoldOptLib.input import CheckInputData, InputDataProcessor, InputData
from FoldOptLib.objective_functions import (
GeologicalKnowledgeFunctions,
VonMisesFisher,
Expand Down
8 changes: 6 additions & 2 deletions FoldOptLib/builders/builder.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,9 @@ def __init__(self, boundingbox: BoundingBox):
self.boundingbox = boundingbox

self.interpolator = LoopInterpolator(
self.boundingbox, dimensions=3, nelements=1000
self.boundingbox,
dimensions=self.boundingbox.dimensions,
nelements=1000
)

def set_constraints(self, constraints: InterpolationConstraints):
Expand All @@ -23,7 +25,9 @@ def evaluate_scalar_value(self, locations: numpy.ndarray) -> numpy.ndarray:
return self.interpolator.evaluate_scalar_value(locations)

def evaluate_gradient(self, locations: numpy.ndarray) -> numpy.ndarray:
return self.interpolator.evaluate_gradient(locations)
gradient = self.interpolator.evaluate_gradient(locations)
gradient /= numpy.linalg.norm(gradient, axis=1)[:, None]
return gradient

def min(self):
"""Calculate the min value of the fold frame
Expand Down
19,374 changes: 6,972 additions & 12,402 deletions FoldOptLib/examples/axial_surface_optimisation_part_1.ipynb

Large diffs are not rendered by default.

2,900 changes: 2,756 additions & 144 deletions FoldOptLib/examples/fourier_series_optimisation.ipynb

Large diffs are not rendered by default.

12 changes: 6 additions & 6 deletions FoldOptLib/examples/misc_functions.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
import numpy as np
import numpy


def sample_random_dataset(bounding_box, sample_size=2, seed=1806545646548564032015):
Expand All @@ -20,7 +20,7 @@ def sample_random_dataset(bounding_box, sample_size=2, seed=18065456465485640320
A sample_size x 3 array of random 3D coordinates within the specified bounding box.
"""
# Set the seed for the random number generator for reproducible results
rng = np.random.default_rng(seed)
rng = numpy.random.default_rng(seed)

# Extract the maximum x and y coordinates from the bounding box
xmax, ymax = bounding_box[1, 0], bounding_box[1, 1]
Expand All @@ -29,15 +29,15 @@ def sample_random_dataset(bounding_box, sample_size=2, seed=18065456465485640320
zmax = 0.0

# Generate 'sample_size' number of random x-coordinates
xn = rng.random.uniform(low=bounding_box[0, 0], high=xmax, size=sample_size)
xn = numpy.random.uniform(low=bounding_box[0, 0], high=xmax, size=sample_size)

# Generate 'sample_size' number of random y-coordinates
yn = rng.random(low=bounding_box[0, 1], high=ymax, size=sample_size)
yn = numpy.random.uniform(low=bounding_box[0, 1], high=ymax, size=sample_size)

# Create an array of z-coordinates, all set to 'zmax' (fixed z-coordinate for all points)
zn = np.tile(zmax, sample_size)
zn = numpy.tile(zmax, sample_size)

# Combine the x, y, and z coordinates into a single 2D array (shape: sample_size x 3)
random_xyz = np.array([xn, yn, zn]).T
random_xyz = numpy.array([xn, yn, zn]).T

return random_xyz
9 changes: 5 additions & 4 deletions FoldOptLib/fold_modelling/engine.py
Original file line number Diff line number Diff line change
Expand Up @@ -368,7 +368,7 @@ def fit_fourier_series(
)

# Check the type of knowledge and set the fit type flag
self.geological_knowledge.fittypeflag[fit_type] = True
fourier_optimiser.geological_knowledge.fittypeflag[fit_type] = True
opt = fourier_optimiser.optimise()

return opt.x
Expand Down Expand Up @@ -439,10 +439,11 @@ def get_predicted_foliation(self, axial_normal: numpy.ndarray) -> numpy.ndarray:
# Build the fold frame
self.build_fold_frame(axial_normal)

# Create and build fold event
self.create_and_build_fold_event()

# Calculate folded foliation vectors
predicted_foliation = self.calculate_folded_foliation_vectors()

# clean up memory
del self.axial_surface
gc.collect()

return predicted_foliation
8 changes: 4 additions & 4 deletions FoldOptLib/input/data_storage.py
Original file line number Diff line number Diff line change
Expand Up @@ -140,8 +140,8 @@ def set_axial_surface_field_constraints(self):
value_constraints = numpy.array(
[
[mean_x, mean_y, mean_z, 0.0, 1.0],
[min_x, min_y, min_z, -1.0, 1.0],
[max_x, max_y, max_z, 1.0, 1.0],
[min_x, min_y, min_z, -500., 1.0],
[max_x, max_y, max_z, 500.0, 1.0],
],
dtype=float,
)
Expand Down Expand Up @@ -171,8 +171,8 @@ def set_fold_axis_field_constraints(self):
value_constraints = numpy.array(
[
[mean_x, mean_y, mean_z, 0.0, 1.0],
[min_x, min_y, min_z, -1.0, 1.0],
[max_x, max_y, max_z, 1.0, 1.0],
[min_x, min_y, min_z, -500., 1.0],
[max_x, max_y, max_z, 500.0, 1.0],
],
dtype=float,
)
Expand Down
4 changes: 2 additions & 2 deletions FoldOptLib/objective_functions/geological_knowledge.py
Original file line number Diff line number Diff line change
Expand Up @@ -453,7 +453,7 @@ def setup_objective_functions_for_restricted_mode(
# Return the list of constraints
return constraints

def __call__(self, theta: numpy.ndarray) -> float:
def __call__(self, theta: numpy.ndarray):
"""
Calculate the total geological knowledge objective function value for all constraints by summing up the
objective function values for all constraints. This objective function represent only the
Expand Down Expand Up @@ -537,4 +537,4 @@ def __getitem__(self, knowledge_type: KnowledgeType):
Union[NormalDistribution, VonMisesFisherDistribution]
The knowledge constraints for the given knowledge type.
"""
return self.input_knowledge(knowledge_type)
return self.input_knowledge[knowledge_type]
2 changes: 1 addition & 1 deletion FoldOptLib/objective_functions/objective_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ def log_normal(

@beartype.beartype
@staticmethod
def vonmises(x: Union[list, numpy.ndarray]):
def vonmises(x: Union[int, float, list, numpy.ndarray]):
"""
Objective function for the axial surface.
This function calculates the loglikelihood of an axial surface using the VonMisesFisher distribution.
Expand Down
99 changes: 37 additions & 62 deletions FoldOptLib/optimisers/axial_surface_optimiser.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
import gc
from typing import Dict, Any, Callable
from typing import Dict, Any, Callable, Union
import numpy
from .fold_optimiser import BaseOptimiser
from ..objective_functions import GeologicalKnowledgeFunctions, ObjectiveFunction
Expand All @@ -19,40 +19,9 @@
import beartype


@beartype.beartype
def calculate_intersection_lineation(axial_surface, folded_foliation):
"""
Calculate the intersection lineation of the axial surface and the folded foliation.
Parameters:
axial_surface (np.ndarray): The normal vector of the axial surface.
folded_foliation (np.ndarray): The normal vector of the folded foliation.
Returns:
np.ndarray: The normalised intersection lineation vector.
"""
# Check if the inputs are numpy arrays
if not isinstance(axial_surface, numpy.ndarray):
raise TypeError("Axial surface vector must be a numpy array.")
if not isinstance(folded_foliation, numpy.ndarray):
raise TypeError("Folded foliation vector must be a numpy array.")

# Check if the inputs have the same shape
if axial_surface.shape != folded_foliation.shape:
raise ValueError(
"Axial surface and folded foliation arrays must have the same shape."
)

# Calculate cross product of the axial surface and folded foliation normal vectors
li = numpy.cross(axial_surface, folded_foliation)

# Normalise the intersection lineation vector
li /= numpy.linalg.norm(li, axis=1)[:, None]

return li


@beartype.beartype
class AxialSurfaceOptimiser(BaseOptimiser):
"""
Optimiser for Axial Surfaces.
Expand Down Expand Up @@ -91,14 +60,12 @@ def __init__(

super().__init__(method=method)
self.fold_engine = FoldModel(data, **kwargs)

self.data = data
# self.geological_knowledge = geological_knowledge
self.geological_knowledge = self.setup_geological_knowledge(
self.data[DataType.GEOLOGICAL_KNOWLEDGE]
)
self.optimisation_type = self.setup_optimisation_type()
self.gradient_data = self.data[["gx", "gy", "gz"]].to_numpy()
self.gradient_data = self.data[DataType.DATA][["gx", "gy", "gz"]].to_numpy()
self.objective_function = None
self.guess = None
self.bounds = None
Expand Down Expand Up @@ -134,7 +101,7 @@ def setup_optimisation_type(self):
@beartype.beartype
@staticmethod
def setup_geological_knowledge(
geological_knowledge: InputGeologicalKnowledge = None,
geological_knowledge: Union[InputGeologicalKnowledge, None] = None,
):
"""
Setup the geological knowledge.
Expand All @@ -157,19 +124,21 @@ def generate_bounds_and_initial_guess(self):
"""
self.bounds = [(0, 360), (0, 90)]

if self.geological_knowledge[KnowledgeType.AXIAL_SURFACE] is not None:
# Create a VonMisesFisher distribution with the given parameters
mu = self.geological_knowledge[KnowledgeType.AXIAL_SURFACE].mu
kappa = 5
vmf = VonMisesFisher(mu, kappa)
# Sample from the distribution
initial_guess = vmf.draw_samples(size=20, random_state=180)
initial_guess = normal_vector_to_strike_and_dip(initial_guess)
return initial_guess
if self.geological_knowledge is not None:

if self.geological_knowledge[KnowledgeType.AXIAL_SURFACE] is not None:
# Create a VonMisesFisher distribution with the given parameters
mu = self.geological_knowledge[KnowledgeType.AXIAL_SURFACE].mu
kappa = self.geological_knowledge[KnowledgeType.AXIAL_SURFACE].kappa
vmf = VonMisesFisher(mu, kappa)
# Sample from the distribution
initial_guess = vmf.draw_samples(size=20, random_state=180564327)
initial_guess = normal_vector_to_strike_and_dip(initial_guess)
return initial_guess

if self.geological_knowledge[KnowledgeType.AXIAL_SURFACE] is None:
# use the halton method to initialise the optimisation
return "halton"
if self.geological_knowledge is None:
# use the halton method to initialise the optimisation
return "halton"

@beartype.beartype
def get_predicted_foliation(self, unit_vector: numpy.ndarray):
Expand Down Expand Up @@ -222,7 +191,12 @@ def setup_optimisation_method(self):
)

# If no geological knowledge about the axial surface is available
elif self.geological_knowledge is None:
if self.geological_knowledge is None:
if self.optimisation_type is OptimisationType.ANGLE:
self.objective_function = self.build_optimisation_function(
ObjectiveFunction[ObjectiveType.ANGLE],
self.get_predicted_foliation,
)
# If the optimisation type is MLE, calculate the logpdf using the log normal distribution
if self.optimisation_type is OptimisationType.MLE:
self.objective_function = self.build_optimisation_function(
Expand All @@ -231,7 +205,7 @@ def setup_optimisation_method(self):
)

# If the optimisation type is VMF_MLE, calculate the logpdf using the Von Mises distribution
elif self.optimisation_type is OptimisationType.VM_MLE:
if self.optimisation_type is OptimisationType.VM_MLE:
self.objective_function = self.build_optimisation_function(
ObjectiveFunction[ObjectiveType.VON_MISES],
self.get_predicted_foliation,
Expand All @@ -244,6 +218,7 @@ def build_optimisation_function(
knowledge_function: GeologicalKnowledgeFunctions = None,
):
def optimisation_function(strike_dip):
print(strike_dip)
# Convert the strike-dip to a unit vector
unit_vector = strike_dip_to_vector(strike_dip[0], strike_dip[1])
# Normalize the unit vector
Expand All @@ -254,22 +229,23 @@ def optimisation_function(strike_dip):
angle_difference = ObjectiveFunction[ObjectiveType.ANGLE](
predicted_foliation, self.gradient_data
)
print("angle difference: ",angle_difference)

# If the optimisation type is angle, return the angle difference
if self.optimisation_type == OptimisationType.ANGLE:
return angle_difference
else:
if knowledge_function is None:
# Calculate the logpdf of the angle difference
logpdf = objective_function(angle_difference)
return logpdf

elif knowledge_function is not None:
if knowledge_function is not None:
# Calculate the logpdf of the angle difference
logpdf = objective_function(angle_difference) + knowledge_function(
unit_vector
)
return logpdf

elif knowledge_function is None:
# Calculate the logpdf of the angle difference
logpdf = objective_function(angle_difference)
return logpdf

# clean up memory
del predicted_foliation, unit_vector
Expand Down Expand Up @@ -307,12 +283,11 @@ def optimise(self):

self.setup_optimisation()

if self._solver is self.optimiser._solvers[SolverType.DIFFERENTIAL_EVOLUTION]:
return self._solver(self.objective_function, self._bounds, init=self._guess)
if self.solver is self._solvers[SolverType.DIFFERENTIAL_EVOLUTION]:
return self.solver(self.objective_function, self.bounds, init=self.guess)

elif (
self._solver is self.optimiser._solvers[SolverType.CONSTRAINED_TRUST_REGION]
):
return self._solver(self.objective_function, x0=self._guess)
if self.solver is self._solvers[SolverType.CONSTRAINED_TRUST_REGION]:

return self.solver(self.objective_function, x0=self.guess)

# TODO: ...add support for restricted optimisation mode...
9 changes: 4 additions & 5 deletions FoldOptLib/optimisers/fourier_optimiser.py
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,6 @@ def __init__(
)
# TODO: check how to initialise self.x = x in self.geological_knowledge
self.x = x
self.solver = None
self.objective_function = None
self.guess = None
self.bounds = None
Expand Down Expand Up @@ -214,12 +213,12 @@ def optimise(self):

self.setup_optimisation()

if self._solver is self.optimiser._solvers[SolverType.DIFFERENTIAL_EVOLUTION]:
return self._solver(self.objective_function, self._bounds, init=self._guess)
if self.solver is self._solvers[SolverType.DIFFERENTIAL_EVOLUTION]:
return self.solver(self.objective_function, self.bounds)

elif (
self._solver is self.optimiser._solvers[SolverType.CONSTRAINED_TRUST_REGION]
self.solver is self._solvers[SolverType.CONSTRAINED_TRUST_REGION]
):
return self._solver(self.objective_function, x0=self._guess)
return self.solver(self.objective_function, x0=self.guess)

# TODO: ...add support for restricted optimisation mode...
10 changes: 5 additions & 5 deletions FoldOptLib/solvers/solvers.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
from typing import Callable, Dict, Any, Tuple, Union, List
import numpy as np
import numpy
from ..utils.utils import *
from scipy.optimize import minimize, differential_evolution
from ..datatypes import SolverType
Expand Down Expand Up @@ -28,8 +28,8 @@ def __init__(self, solver="differential_evolution", **kwargs: Dict[str, Any]):
@staticmethod
def differential_evolution(
objective_function: Callable,
bounds: Union[Tuple, List],
init: str = "halton",
bounds: Union[Tuple, List, numpy.ndarray],
init: Union[str, numpy.ndarray] = "halton",
maxiter: int = 5000,
seed: int = 80,
polish: bool = True,
Expand Down Expand Up @@ -85,7 +85,7 @@ def differential_evolution(
@beartype.beartype
@staticmethod
def constrained_trust_region(
objective_function: Callable, x0: np.ndarray, constraints=None, **kwargs
objective_function: Callable, x0: numpy.ndarray, constraints=None, **kwargs
) -> Dict:
"""
Solves the optimisation problem using the trust region method.
Expand All @@ -95,7 +95,7 @@ def constrained_trust_region(
objective_function: Callable
The objective function to be optimised.
x0 : np.ndarray
x0 : numpy.ndarray
Initial guess of the parameters to be optimised.
Returns
Expand Down

0 comments on commit c2f143c

Please sign in to comment.