first commit

contact_function.py

@@ -0,0 +1,167 @@
+from dataclasses import dataclass
+import numpy as np
+import time
+import ctypes
+
+
+clib = ctypes.cdll.LoadLibrary('./overlap_algorithm.so')
+clib.contact_function.argtypes = [
+    ctypes.c_double,
+    ctypes.c_double,
+    ctypes.c_double,
+    ctypes.c_double,
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+    ctypes.c_int,
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.intc, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.intc, ndim=1, flags='C_CONTIGUOUS'),
+]
+clib.contact_function.restype = ctypes.c_double
+
+clib.contact_function_multi.argtypes = [
+    ctypes.c_int,
+    ctypes.c_double,
+    ctypes.c_double,
+    ctypes.c_double,
+    ctypes.c_double,
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    ctypes.c_int,
+    np.ctypeslib.ndpointer(dtype=np.intc, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+]
+clib.contact_function_multi.restype = ctypes.c_void_p
+
+
+def minimizer_map(minimizer: str) -> int:
+    if minimizer == 'brent':
+        return 0
+    elif minimizer == 'brent_early':
+        return 1
+    elif minimizer == 'gss':
+        return 2
+    elif minimizer == 'gss_early':
+        return 3
+    else:
+        raise ValueError('Unknown minimizer: "' + minimizer + '"')
+
+
+@dataclass
+class Contact:
+    contact_f: float
+    min_eig_T: float
+    med_eig_T: float
+    min_vec_scalar: float
+    med_vec_scalar: float
+    feval_min: int
+    feval_med: int
+    branch: int
+
+
+def contact_function(a0, a1, b0, b1, coord0, orient0, coord1, orient1, minimizer='brent') -> Contact:
+    other_results = np.zeros(4, dtype=np.float64)
+    fevals = np.zeros(2, dtype=np.intc)
+    branch = np.zeros(1, dtype=np.intc)
+
+    contact_f = clib.contact_function(a0, a1, b0, b1, coord0, orient0, coord1, orient1,
+                                      minimizer_map(minimizer), other_results, fevals, branch)
+
+    return Contact(contact_f, other_results[0], other_results[1], other_results[2], other_results[3],
+                   fevals[0], fevals[1], branch[0])
+
+
+@dataclass
+class ContactMulti:
+    contact_f: np.ndarray
+    avg_evals_mineig: float
+    avg_evals_medeig: float
+    time: float
+
+
+def contact_function_multi(a0, a1, b0, b1, coord0, orient0, coord1, orient1, minimizer='brent') -> ContactMulti:
+    n = len(coord0)
+    results = np.zeros(n, dtype=np.float64)
+    all_evals = np.zeros(2, dtype=np.intc)
+
+    t0 = time.perf_counter()
+    clib.contact_function_multi(n, a0, a1, b0, b1, coord0, orient0, coord1, orient1,
+                                minimizer_map(minimizer), all_evals, results)
+    t1 = time.perf_counter()
+
+    return ContactMulti(results, all_evals[0] / n, all_evals[1] / n, t1 - t0)
+
+
+class EllipsePair:
+
+    def __init__(self, a0, a1, epsilon):
+        """
+        Class in which optimization of ellipse configurations on a sphere is performed.
+        :param a0: semi-major axis for the first elliptical cylinder
+        :param a1: semi-major axis for the second elliptical cylinder
+        :param epsilon: aspect ratio
+        """
+        self.a0 = max(a0, a1)  # defined so that always a0 > a1
+        self.a1 = min(a0, a1)
+        self.b0 = self.a0 / epsilon
+        self.b1 = self.a1 / epsilon
+        self.size_ratio = self.a0 / self.a1
+        self.epsilon = epsilon
+
+    def contact(self, coord0, orient0, coord1=np.array([0., 0., 1.]), orient1=np.array([1., 0., 0.]),
+                minimizer='brent') -> Contact:
+        return contact_function(self.a0, self.a1, self.b0, self.b1, coord0, orient0, coord1, orient1, minimizer)
+
+    def contact_multi_dist(self, distances, num_ang, axis=np.array([1, 0, 0]), minimizer='brent')\
+            -> (np.ndarray, np.ndarray, np.ndarray):
+
+        timings = np.zeros(len(distances))
+        avg_evals_mineig = np.zeros(len(distances))
+        avg_evals_medeig = np.zeros(len(distances))
+        for i, dist in enumerate(distances):
+            coords0, coords1, orients0, orients1 = generate_confgs(dist, num_ang, axis=axis)
+            contact_multi = contact_function_multi(self.a0, self.a1, self.b0, self.b1,
+                                                   coords0, orients0, coords1, orients1,
+                                                   minimizer=minimizer)
+            timings[i] = contact_multi.time
+            avg_evals_mineig[i] = contact_multi.avg_evals_mineig
+            avg_evals_medeig[i] = contact_multi.avg_evals_medeig
+
+        print(f'Total time: {np.sum(timings)}')
+
+        return timings, avg_evals_mineig, avg_evals_medeig
+
+
+def generate_confgs(dist, n, axis=np.array([1, 0, 0])):
+    """
+    At a given distance, generate all possible orientational configurations of two spherical ellipses,
+    with n steps in a half-circle rotation. This results in n x n configurations.
+    The first ellipse will always be centered at the north pole while the position of the second is determined by
+    parameter axis (and obviously distance).
+    """
+
+    angles = np.linspace(1e-7, np.pi - 1e-7, n)
+
+    coords0 = np.zeros((n, 3))
+    coords0[:, 2] = 1
+    orients0 = np.array([1., 0., 0.])[None, :] * np.cos(angles)[:, None] + \
+               np.array([0., 1., 0.])[None, :] * np.sin(angles)[:, None]
+
+    coords1 = coords0 * np.cos(dist) + \
+              np.cross(axis, coords0) * np.sin(dist) + \
+              axis[None, :] * np.einsum('m, nm -> n', axis, coords0)[:, None] * (1 - np.cos(dist))
+    orients1 = orients0 * np.cos(dist) + \
+               np.cross(axis, orients0) * np.sin(dist) + \
+               axis[None, :] * np.einsum('m, nm -> n', axis, orients0)[:, None] * (1 - np.cos(dist))
+
+    indices = np.indices((n, n))
+    coords0 = coords0[indices[0].flatten()]
+    orients0 = orients0[indices[0].flatten()]
+    coords1 = coords1[indices[1].flatten()]
+    orients1 = orients1[indices[1].flatten()]
+
+    return coords0, coords1, orients0, orients1

overlap_algorithm.c

@@ -0,0 +1,931 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#define R (sqrt(5)-1.)*0.5
+#define C (1.0-R)
+#define SHFT2(a,b,c) (a)=(b);(b)=(c);
+#define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
+
+
+//.........................COMPILE TO SHARED LIBRARY FOR USE WITH PYTHON.........................................
+
+// gcc -shared overlap_algorithm.c -o overlap_algorithm.so -lm -fPIC -O3
+
+
+//..................................LINE MINIMIZATION ALGORITHMS.................................................
+
+
+double gss_eigenvalue_minim(double a, double b, double eps, double t, double *matA, double *matB, int early_stop,
+                  double f ( double x, double *mat1, double *mat2 ), double *xmin, int *feval) {
+
+  double f1,f2,x0,x1,x2,x3,m,tol,t2;
+  int num_eval;
+
+  x0=a;
+  x3=b;
+
+  x1 = a + C * (b - a);
+  x2 = a + R * (b - a);
+
+  f1=f(x1, matA, matB);
+  f2=f(x2, matA, matB);
+  num_eval = 2;
+
+  for ( ; ; )
+  {
+    m = 0.5 * ( x1 + x2 ) ;
+    tol = eps * fabs ( m ) + t;
+    t2 = 2.0 * tol;
+/*
+  Check the stopping criterion.
+*/
+    if (fabs(x3-x0) < t2)
+    {
+      break;
+    }
+
+    if ( early_stop == 1 && (f1 < -1. || f2 < -1.)) {
+      break;
+    }
+
+    if (f2 < f1) {
+    SHFT3(x0,x1,x2,R*x1+C*x3);
+    SHFT2(f1,f2,f(x2, matA, matB));
+    num_eval++;
+    }
+    else {
+      SHFT3(x3,x2,x1,R*x2+C*x0);
+      SHFT2(f2,f1,f(x1, matA, matB));
+      num_eval++;
+    }
+  }
+
+  if (f1 < f2) {
+    *xmin=x1;
+    *feval = num_eval;
+    return f1;
+  }
+  else {
+    *xmin=x2;
+    *feval = num_eval;
+    return f2;
+  }
+}
+
+
+
+double brent_eigenvalue_minim ( double a, double b, double eps, double t,
+  double *matA, double *matB, int early_stop,
+  double f ( double x, double *mat1, double *mat2 ), double *x, int *feval)
+
+/*
+
+  Adapted version of LOCAL_MIN function for finding eigenvalue extrema.
+
+  Authors of original minimization code:
+
+    Orignal FORTRAN77 version by Richard Brent.
+    C version by John Burkardt.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license.
+
+  Parameters:
+
+    Input, double A, B, the endpoints of the interval.
+
+    Input, double EPS, a positive relative error tolerance.
+    EPS should be no smaller than twice the relative machine precision,
+    and preferably not much less than the square root of the relative
+    machine precision.
+
+    Input, double T, a positive absolute error tolerance.
+
+    Input, double F ( double x ,double *mat1, double *mat2 ), a user-supplied
+    function whose local minimum along x is being sought, with matrices mat1
+    and mat2 as parameters.
+
+    Output, double *X, the estimated value of an abscissa
+    for which F attains a local minimum value in [A,B].
+
+    Output, double *feval, the number of function evaluations.
+
+    Output, double LOCAL_MIN, the value F(X).
+*/
+{
+  double c;
+  double d;
+  double e;
+  double fu;
+  double fv;
+  double fw;
+  double fx;
+  double m;
+  double p;
+  double q;
+  double r;
+  double sa;
+  double sb;
+  double t2;
+  double tol;
+  double u;
+  double v;
+  double w;
+  int num_eval;
+/*
+  C is the square of the inverse of the golden ratio.
+*/
+  c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
+
+  sa = a;
+  sb = b;
+  *x = sa + c * ( b - a );
+  w = *x;
+  v = w;
+  e = 0.0;
+  fx = f ( *x, matA, matB );
+  fw = fx;
+  fv = fw;
+  num_eval = 1;
+
+  for ( ; ; )
+  {
+    m = 0.5 * ( sa + sb ) ;
+    tol = eps * fabs ( *x ) + t;
+    t2 = 2.0 * tol;
+/*
+  Check the stopping criterion.
+*/
+    if ( fabs ( *x - m ) <= t2 - 0.5 * ( sb - sa ) )
+    {
+      break;
+    }
+
+/*
+  If we are looking for the maximum of smallest eigenvalue and the
+  calculation goes over the value of 1 (under -1 because algorithm finds min),
+  we can stop as that means the ellipses certainly do not overlap
+*/
+    if ( early_stop == 1 && fx < -1.) {
+      break;
+    }
+/*
+  Fit a parabola.
+*/
+    r = 0.0;
+    q = r;
+    p = q;
+
+    if ( tol < fabs ( e ) )
+    {
+      r = ( *x - w ) * ( fx - fv );
+      q = ( *x - v ) * ( fx - fw );
+      p = ( *x - v ) * q - ( *x - w ) * r;
+      q = 2.0 * ( q - r );
+      if ( 0.0 < q )
+      {
+        p = - p;
+      }
+      q = fabs ( q );
+      r = e;
+      e = d;
+    }
+
+    if ( fabs ( p ) < fabs ( 0.5 * q * r ) &&
+         q * ( sa - *x ) < p &&
+         p < q * ( sb - *x ) )
+    {
+/*
+  Take the parabolic interpolation step.
+*/
+      d = p / q;
+      u = *x + d;
+/*
+  F must not be evaluated too close to A or B.
+*/
+      if ( ( u - sa ) < t2 || ( sb - u ) < t2 )
+      {
+        if ( *x < m )
+        {
+          d = tol;
+        }
+        else
+        {
+          d = - tol;
+        }
+      }
+    }
+/*
+  A golden-section step.
+*/
+    else
+    {
+      if ( *x < m )
+      {
+        e = sb - *x;
+      }
+      else
+      {
+        e = sa - *x;
+      }
+      d = c * e;
+    }
+/*
+  F must not be evaluated too close to X.
+*/
+    if ( tol <= fabs ( d ) )
+    {
+      u = *x + d;
+    }
+    else if ( 0.0 < d )
+    {
+      u = *x + tol;
+    }
+    else
+    {
+      u = *x - tol;
+    }
+
+    fu = f ( u, matA, matB );
+    num_eval++;
+/*
+  Update A, B, V, W, and X.
+*/
+    if ( fu <= fx )
+    {
+      if ( u < *x )
+      {
+        sb = *x;
+      }
+      else
+      {
+        sa = *x;
+      }
+      v = w;
+      fv = fw;
+      w = *x;
+      fw = fx;
+      *x = u;
+      fx = fu;
+    }
+    else
+    {
+      if ( u < *x )
+      {
+        sa = u;
+      }
+      else
+      {
+        sb = u;
+      }
+
+      if ( fu <= fw || w == *x )
+      {
+        v = w;
+        fv = fw;
+        w = u;
+        fw = fu;
+      }
+      else if ( fu <= fv || v == *x || v == w )
+      {
+        v = u;
+        fv = fu;
+      }
+    }
+  }
+  *feval = num_eval;
+  return fx;
+}
+
+
+double r8_epsilon ( )
+/*
+  Purpose:
+
+    R8_EPSILON returns the R8 round off unit.
+
+  Discussion:
+
+    R8_EPSILON is a number R which is a power of 2 with the property that,
+    to the precision of the computer's arithmetic,
+      1 < 1 + R
+    but
+      1 = ( 1 + R / 2 )
+
+  Parameters:
+
+    Output, double R8_EPSILON, the R8 round-off unit.
+*/
+{
+  const double value = 2.220446049250313E-016;
+
+  return value;
+}
+
+
+//.................................EIGENVALUE AND EIGENVECTOR CALCULATION.............................................
+
+
+double determinant(double* a){
+  double det = a[0] * (a[4]*a[8] - a[7]*a[5]) - a[3] * (a[1]*a[8] - a[7]*a[2]) + a[6] * (a[1]*a[5] - a[4]*a[2]);
+  return det;
+}
+
+
+void eigenvalues(double* matrix, double* eigs){
+
+  double p1 = matrix[1]*matrix[1] + matrix[2]*matrix[2] + matrix[5]*matrix[5];
+  double q, p2, p;
+  double aux[9];
+  double r, phi;
+  double pi = M_PI;
+
+  if (p1 == 0) {
+    eigs[0] = matrix[0];
+    eigs[1] = matrix[4];
+    eigs[2] = matrix[8];
+  }
+  else {
+    q = (matrix[0]+matrix[4]+matrix[8])/3;
+    p2 = (matrix[0]-q)*(matrix[0]-q) + (matrix[4]-q)*(matrix[4]-q) + (matrix[8]-q)*(matrix[8]-q) + 2*p1;
+    p = sqrt(p2 / 6);
+
+    for (int i=0; i<9; i++) {
+      aux[i] = matrix[i];
+      if (i % 4 == 0) {
+        aux[i] -= q;
+      }
+      aux[i] = (1 / p) * aux[i];
+    }
+
+    r = determinant(aux) / 2;
+
+    if (r <= -1){
+      phi = pi / 3;
+    }
+    else if (r >= 1){
+      phi = 0;
+    }
+    else {
+      phi = acos(r) / 3;
+    }
+
+    eigs[2] = q + 2 * p * cos(phi);
+    eigs[0] = q + 2 * p * cos(phi + (2*pi/3));
+    eigs[1] = 3 * q - eigs[0] - eigs[2];
+  }
+
+}
+
+
+void eigenvector(double* matrix, int which, double* vector){
+
+  double eigs[3];
+  double mat_lbd1[9];
+  double mat_lbd2[9];
+  double product[9];
+  double norm1, norm2, norm3, norm;
+  int idx;
+
+  eigenvalues(matrix, eigs);
+
+  for (int i=0; i<9; i++) {
+    mat_lbd1[i] = matrix[i];
+    mat_lbd2[i] = matrix[i];
+    if (i % 4 == 0) {
+      if (which == 1) {
+        mat_lbd1[i] -= eigs[1];
+        mat_lbd2[i] -= eigs[2];
+      }
+      else if (which == 2) {
+        mat_lbd1[i] -= eigs[0];
+        mat_lbd2[i] -= eigs[2];
+      }
+      else {
+        mat_lbd1[i] -= eigs[0];
+        mat_lbd2[i] -= eigs[1];
+      }
+    }
+  }
+
+  for (int i=0; i<3; i++) {
+    for (int k=0; k<3; k++) {
+      product[3*i+k] = mat_lbd1[3*i] * mat_lbd2[k] + mat_lbd1[3*i+1] * mat_lbd2[3+k] + mat_lbd1[3*i+2] * mat_lbd2[6+k];
+    }
+  }
+
+  norm1 = product[0]*product[0] + product[3]*product[3] + product[6]*product[6];
+  norm2 = product[1]*product[1] + product[4]*product[4] + product[7]*product[7];
+  norm3 = product[2]*product[2] + product[5]*product[5] + product[8]*product[8];
+
+  if (norm1 > norm2 && norm1 > norm3) {
+    idx = 0;
+    norm = norm1;
+  }
+  else if (norm2 > norm1 && norm2 > norm3) {
+    idx = 1;
+    norm = norm2;
+  }
+  else {
+    idx = 2;
+    norm = norm3;
+  }
+
+  norm = 1/sqrt(norm);
+
+  vector[0] = product[idx]*norm;
+  vector[1] = product[3+idx]*norm;
+  vector[2] = product[6+idx]*norm;
+}
+
+
+//.....................................EIGENVALUE EXTREMIZATION.................................................
+
+
+double neg_min_eig(double lbd, double* matA, double* matB){
+
+  double lcb[9];
+  double eigs[3];
+  double result;
+
+  for (int i=0; i<9; i++) {
+    lcb[i] = lbd * matA[i] + (1 - lbd) * matB[i];
+  }
+
+  eigenvalues(lcb, eigs);
+
+  result = -eigs[0];
+
+  return result;
+}
+
+
+double med_eig(double lbd, double* matA, double* matB){
+
+  double lcb[9];
+  double eigs[3];
+  double result;
+
+  for (int i=0; i<9; i++) {
+    lcb[i] = lbd * matA[i] + (1 - lbd) * matB[i];
+  }
+
+  eigenvalues(lcb, eigs);
+
+  result = eigs[1];
+
+
+  return result;
+}
+
+
+double max_min_eig_brent(double* matA, double* matB, int early_stop, double* t_point, int* feval){
+
+  double a = r8_epsilon();
+  double b = 1. - r8_epsilon();
+  double e = r8_epsilon(); // relative error
+  double t = 1.0E-7; // absolute error
+  double solution;
+  double initial = 0.;
+  double *y = &initial;
+
+
+  solution = -brent_eigenvalue_minim(a, b, e, t, matA, matB, early_stop, neg_min_eig, t_point, feval);
+
+
+  return solution;
+
+}
+
+
+double min_med_eig_brent(double* matA, double* matB, double* t_point, int* feval){
+
+  double a = r8_epsilon();
+  double b = 1. - r8_epsilon();
+  double e = r8_epsilon(); // relative error
+  double t = 1.0E-7; // absolute error
+  int early_stop = 0;
+  double solution;
+  double initial = 0.;
+  double *y = &initial;
+
+  solution = brent_eigenvalue_minim(a, b, e, t, matA, matB, early_stop, med_eig, t_point, feval);
+
+  return solution;
+
+}
+
+
+double max_min_eig_gss(double* matA, double* matB, int early_stop, double* t_point, int* feval){
+
+  double a = r8_epsilon();
+  double b = 1. - r8_epsilon();
+  double c = 0.5;
+  double m = 100.;
+  double mechep = r8_epsilon();
+  double e = r8_epsilon(); // relative error
+  double t = 1.0E-7; // absolute error
+  double solution;
+
+  solution = -gss_eigenvalue_minim(a, b, e, t, matA, matB, early_stop, neg_min_eig, t_point, feval);
+
+  return solution;
+
+}
+
+double min_med_eig_gss(double* matA, double* matB, double* t_point, int* feval){
+
+  double a = r8_epsilon();
+  double b = 1. - r8_epsilon();
+  double c = 0.5;
+  double m = 100.;
+  double mechep = r8_epsilon();
+  double e = r8_epsilon(); // relative error
+  double t = 1.0E-7; // absolute error
+  int early_stop = 0;
+  double solution;
+
+  solution = gss_eigenvalue_minim(a, b, e, t, matA, matB, early_stop, med_eig, t_point, feval);
+
+  return solution;
+
+}
+
+
+//.......................................CONTACT FUNCTION...................................................
+
+
+void cross_product(double* coord, double* orient, double* cross) {
+
+  cross[0] = coord[1]*orient[2] - coord[2]*orient[1];
+  cross[1] = coord[2]*orient[0] - coord[0]*orient[2];
+  cross[2] = coord[0]*orient[1] - coord[1]*orient[0];
+
+}
+
+
+double contact_function(double a0, double a1, double b0, double b1,
+  double* coord0, double* orient0, double* coord1, double* orient1, int minimizer,
+  double* other_results, int* fevals, int* branch) {
+
+  /*
+    Other results array: [min_eig_T, med_eig_T, min_vec_scalar, med_vec_scalar]
+    fevals array: [feval_min, feval_med]
+  */
+
+  double a_inter;
+  if (a0 < a1) {
+    a_inter = a0;
+    a0 = a1;
+    a1 = a_inter;
+  }
+
+  double cross0[3];
+  double cross1[3];
+  double lin_interpolation[9];
+  double qf0[9];
+  double qf1[9];
+  double min_eig;
+  double med_eig;
+  double eigvec_mineig[3];
+  double mineig_coord0=0.;
+  double mineig_coord1=0.;
+  double inv2a0=1/(a0 * a0);
+  double inv2a1=1/(a1 * a1);
+  double inv2b0=1/(b0 * b0);
+  double inv2b1=1/(b1 * b1);
+  double limit=inv2a0;
+
+  cross_product(coord0, orient0, cross0);
+  cross_product(coord1, orient1, cross1);
+
+  /* Trasformation of quadratic forms */
+  for (int i=0; i<3; i++){
+    for (int l=0; l<3; l++){
+      qf0[3*i+l] = orient0[i]*inv2a0*orient0[l] + cross0[i]*inv2b0*cross0[l];
+      qf1[3*i+l] = orient1[i]*inv2a1*orient1[l] + cross1[i]*inv2b1*cross1[l];
+    }
+  }
+
+  if (minimizer == 0) {
+    min_eig = max_min_eig_brent(qf0, qf1, 0, &other_results[0], &fevals[0]);
+  }
+  else if (minimizer == 1) {
+    min_eig = max_min_eig_brent(qf0, qf1, 1, &other_results[0], &fevals[0]);
+  }
+  else if (minimizer == 2) {
+    min_eig = max_min_eig_gss(qf0, qf1, 0, &other_results[0], &fevals[0]);
+  }
+  else if (minimizer == 3) {
+    min_eig = max_min_eig_gss(qf0, qf1, 1, &other_results[0], &fevals[0]);
+  }
+  else {
+    return -1.;
+  }
+
+
+  /* eigenvector corresponding to maximal 1st eigenvalue (1st intersection coordinate)*/
+  for (int i=0; i<9; i++) {
+    lin_interpolation[i] = other_results[0] * qf0[i] + (1 -   other_results[0]) * qf1[i];
+  }
+  eigenvector(lin_interpolation, 1, eigvec_mineig);
+
+  /* scalar product between eigenvectors and position vectors for both ellipses in the pair */
+  for (int i=0; i<3; i++) {
+    mineig_coord0 += eigvec_mineig[i]*coord0[i];
+    mineig_coord1 += eigvec_mineig[i]*coord1[i];
+  }
+
+  other_results[2] = mineig_coord0 * mineig_coord1;
+  other_results[3] = 0;
+
+
+  if (mineig_coord0 * mineig_coord1 > 0) {
+    other_results[1] = -1; // meaningless result for med_eig_T
+    fevals[1] = 0;
+    if (min_eig < limit) {
+      *branch = 0;
+      return min_eig;
+    }
+    else {
+      *branch = 2;
+      return limit;
+    }
+  }
+
+  else {
+
+    if (minimizer == 0 || minimizer == 1 || minimizer == 4 || minimizer == 5) {
+      med_eig = min_med_eig_brent(qf0, qf1, &other_results[1], &fevals[1]);
+    }
+    else if (minimizer == 2 || minimizer == 3) {
+      med_eig = min_med_eig_gss(qf0, qf1, &other_results[1], &fevals[1]);
+    }
+    else {
+      return -1.;
+    }
+
+    if (1.0E-4 < other_results[1] && other_results[1] < 1 - 1.0E-4 && med_eig < limit) {
+      *branch = 1;
+      return med_eig;
+    }
+    else {
+      *branch = 2;
+      other_results[1] = 1;
+      return limit;
+    }
+
+  }
+
+}
+
+
+void contact_function_multi(int n, double a0, double a1, double b0, double b1,
+  double* coords0, double* orients0, double* coords1, double* orients1, int minimizer,
+  int* all_evals, double* results) {
+
+  double other_results[4];
+  int branch=0;
+  int* branch_point = &branch;
+  int fevals[2];
+  fevals[0] = 0;
+  fevals[1] = 0;
+
+  for (int i=0; i<n; i++){
+    results[i] = contact_function(a0, a1, b0, b1, &coords0[3*i], &orients0[3*i], &coords1[3*i], &orients1[3*i],
+      minimizer, &other_results[0], fevals, &branch);
+    all_evals[0] += fevals[0];
+    all_evals[1] += fevals[1];
+    fevals[0] = 0;
+    fevals[1] = 0;
+  }
+
+}
+
+
+
+double contact_mono_packing(double inv2a, double inv2b,
+                            double* coord0, double* orient0,
+                            double* coord1, double* orient1,
+                            int try_second, int early_stop) {
+
+  /*
+    Simplified contact function calculation for monodispersed ellipse packing.
+  */
+
+  double cross0[3];
+  double cross1[3];
+  double lin_interpolation[9];
+  double qf0[9];
+  double qf1[9];
+  double min_eig;
+  double min_eig_T = 0.;
+  double* meT_point = &min_eig_T;
+  int fevals = 0;
+  int* fevals_point = &fevals;
+
+  cross_product(coord0, orient0, cross0);
+  cross_product(coord1, orient1, cross1);
+
+  /* Trasformation of quadratic forms */
+  for (int i=0; i<3; i++){
+    for (int l=0; l<3; l++){
+      qf0[3*i+l] = orient0[i]*inv2a*orient0[l] + cross0[i]*inv2b*cross0[l];
+      qf1[3*i+l] = orient1[i]*inv2a*orient1[l] + cross1[i]*inv2b*cross1[l];
+    }
+  }
+
+  min_eig = max_min_eig_brent(qf0, qf1, early_stop, meT_point, fevals_point);
+
+  if (try_second == 0 ){
+    return min_eig;
+  }
+
+  else {  // The case where the 2nd eigenvalue must also be calculated
+
+    double lin_interpolation[9];
+    double eigvec_mineig[3];
+    double mineig_coord0 = 0;
+    double mineig_coord1 = 0;
+
+    /* eigenvector corresponding to maximal 1st eigenvalue (1st intersection coordinate)*/
+    for (int i=0; i<9; i++) {
+      lin_interpolation[i] = min_eig_T * qf0[i] + (1 - min_eig_T) * qf1[i];
+    }
+    eigenvector(lin_interpolation, 1, eigvec_mineig);
+
+    /* scalar product between eigenvectors and position vectors for both ellipses in the pair */
+    for (int i=0; i<3; i++) {
+      mineig_coord0 += eigvec_mineig[i]*coord0[i];
+      mineig_coord1 += eigvec_mineig[i]*coord1[i];
+    }
+
+    if (mineig_coord0 * mineig_coord1 > 0) {
+      return min_eig;
+    }
+
+    else {
+      min_eig_T = 0.;
+      fevals = 0;
+      return min_med_eig_brent(qf0, qf1, meT_point, fevals_point);
+    }
+
+  }
+
+}
+
+
+double contact_mono_pack_multi(int n, double inv2a, double inv2b, double* coords, double* orients, double dist_lim) {
+
+  double energy=0.;
+  double cos_dist;
+  double contact_f;
+  int try_second=0;
+
+  if (inv2a < 2) {
+    try_second = 1;  // for large ellipses, 2nd eigenvalue might be the correct one
+  }
+
+  for (int i=0; i<n; i++) {
+    for (int j=i+1; j<n; j++) {
+      cos_dist = coords[3*i]*coords[3*j] + coords[3*i+1]*coords[3*j+1] + coords[3*i+2]*coords[3*j+2];
+      if (cos_dist > dist_lim) {
+        contact_f = contact_mono_packing(inv2a, inv2b, &coords[3*i], &orients[3*i], &coords[3*j], &orients[3*j],
+                                          try_second, 1);
+        if (contact_f < 1) {
+          energy += contact_f + 1 / contact_f - 2;
+        }
+      }
+    }
+  }
+
+  return energy;
+}
+
+
+double contact_bidisp_packing(double inv2a_0, double inv2b_0,
+                              double inv2a_1, double inv2b_1,
+                              double* coord0, double* orient0,
+                              double* coord1, double* orient1,
+                              int try_second, double omega) {
+
+  /*
+    Simplified contact function calculation for bidispersed ellipse packing.
+  */
+
+  double cross0[3];
+  double cross1[3];
+  double lin_interpolation[9];
+  double qf0[9];
+  double qf1[9];
+  double min_eig;
+  double med_eig;
+  double min_eig_T = 0.;
+  double* meT_point = &min_eig_T;
+  int fevals = 0;
+  int* fevals_point = &fevals;
+
+  cross_product(coord0, orient0, cross0);
+  cross_product(coord1, orient1, cross1);
+
+  /* Trasformation of quadratic forms */
+  for (int i=0; i<3; i++){
+    for (int l=0; l<3; l++){
+      qf0[3*i+l] = orient0[i]*inv2a_0*orient0[l] + cross0[i]*inv2b_0*cross0[l];
+      qf1[3*i+l] = orient1[i]*inv2a_1*orient1[l] + cross1[i]*inv2b_1*cross1[l];
+    }
+  }
+
+  min_eig = max_min_eig_brent(qf0, qf1, 1, meT_point, fevals_point);
+
+  if (try_second == 0 ){
+    if (min_eig < omega) {
+      return min_eig;
+      }
+    else {
+      return omega;
+      }
+  }
+
+  else {  // The case where the 2nd eigenvalue must also be calculated
+
+    double lin_interpolation[9];
+    double eigvec_mineig[3];
+    double mineig_coord0 = 0;
+    double mineig_coord1 = 0;
+
+    /* eigenvector corresponding to maximal 1st eigenvalue (1st intersection coordinate)*/
+    for (int i=0; i<9; i++) {
+      lin_interpolation[i] = min_eig_T * qf0[i] + (1 - min_eig_T) * qf1[i];
+    }
+    eigenvector(lin_interpolation, 1, eigvec_mineig);
+
+    /* scalar product between eigenvectors and position vectors for both ellipses in the pair */
+    for (int i=0; i<3; i++) {
+      mineig_coord0 += eigvec_mineig[i]*coord0[i];
+      mineig_coord1 += eigvec_mineig[i]*coord1[i];
+    }
+
+    if (mineig_coord0 * mineig_coord1 > 0) {
+      if (min_eig < omega) {
+        return min_eig;
+        }
+      else {
+        return omega;
+        }
+    }
+
+    else {
+      min_eig_T = 0.;
+      fevals = 0;
+      med_eig = min_med_eig_brent(qf0, qf1, meT_point, fevals_point);
+
+      if (med_eig < omega) {
+        return med_eig;
+        }
+      else {
+        return omega;
+        }
+
+    }
+  }
+}
+
+
+double contact_bidisp_pack_multi(int n, double* inv2a, double* inv2b,
+                                double* coords, double* orients, double dist_lim, double omega) {
+
+  double energy=0.;
+  double cos_dist;
+  double contact_f;
+  int try_second;
+
+  for (int i=0; i<n; i++) {
+    for (int j=i+1; j<n; j++) {
+
+      cos_dist = coords[3*i]*coords[3*j] + coords[3*i+1]*coords[3*j+1] + coords[3*i+2]*coords[3*j+2];
+
+      if (cos_dist > dist_lim) {
+
+        if (inv2a[i] * inv2a[j] < inv2a[i] + inv2a[j]) {  // derived from the condition alpha0 + alpha1 > pi / 2
+          try_second = 1;  // for large ellipses, 2nd eigenvalue might be the correct one
+        }
+        else {
+          try_second = 0;
+        }
+
+        contact_f = contact_bidisp_packing(inv2a[i], inv2b[i], inv2a[j], inv2b[j],
+                      &coords[3*i], &orients[3*i], &coords[3*j], &orients[3*j], try_second, omega);
+
+        if (contact_f < 1) {
+          energy += contact_f + 1 / contact_f - 2;
+        }
+
+      }
+    }
+  }
+
+  return energy;
+}

packing_simulation.py

@@ -0,0 +1,427 @@
+import os
+
+os.environ["OPENBLAS_NUM_THREADS"] = "1"  # mitigates the clash between numpy and multiprocessing  parallelization
+
+import numpy as np
+import ctypes
+import multiprocessing
+from functools import partial
+import time
+from scipy.optimize import minimize as minim, OptimizeResult
+from dataclasses import dataclass, field
+from abc import ABC, abstractmethod
+
+
+def ephi(phi):
+    u = np.zeros((len(phi), 3))
+    u[:, 0] = -np.sin(phi)
+    u[:, 1] = np.cos(phi)
+    return u
+
+
+def etheta(phi, theta):
+    u = np.zeros((len(phi), 3))
+    u[:, 0] = np.cos(phi) * np.cos(theta)
+    u[:, 1] = np.sin(phi) * np.cos(theta)
+    u[:, 2] = -np.sin(theta)
+    return u
+
+
+def sph_to_cart(sph):
+    """Transformation form spherical coordinates phi-theta on a sphere to cartesian coordinates."""
+    cart = np.zeros((len(sph), 3))
+    cart[:, 0] = np.cos(sph[:, 0]) * np.sin(sph[:, 1])
+    cart[:, 1] = np.sin(sph[:, 0]) * np.sin(sph[:, 1])
+    cart[:, 2] = np.cos(sph[:, 1])
+    return cart
+
+
+def vector_orientations(sph, xi):
+    """
+    Calculates vector_orientations of ellipse orientations given their positions in spherical coordinates
+    and orientation angles.
+    """
+    cosxi = np.cos(xi)[:, np.newaxis]
+    sinxi = np.sin(xi)[:, np.newaxis]
+    basis_theta = etheta(sph[:, 0], sph[:, 1])
+    basis_phi = ephi(sph[:, 0])
+    return cosxi * basis_theta + sinxi * basis_phi
+
+
+clib = ctypes.cdll.LoadLibrary('./overlap_algorithm.so')
+clib.contact_mono_pack_multi.argtypes = [
+    ctypes.c_int,
+    ctypes.c_double,
+    ctypes.c_double,
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    ctypes.c_double
+    ]
+clib.contact_mono_pack_multi.restype = ctypes.c_double
+
+clib.contact_bidisp_pack_multi.argtypes = [
+    ctypes.c_int,
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, flags='C_CONTIGUOUS'),
+    ctypes.c_double,
+    ctypes.c_double
+    ]
+clib.contact_bidisp_pack_multi.restype = ctypes.c_double
+
+
+def monodisp_config_energy(inv2a, inv2b, coords, orients, dist_lim) -> float:
+    if inv2a >= 1:
+        n = len(coords)
+        return clib.contact_mono_pack_multi(n, inv2a, inv2b, coords, orients, dist_lim)
+    raise ValueError("Semi-major axis should be smaller that pi/2, otherwise the ellipse becomes inverted.")
+
+
+def bidisp_config_energy(inv2a_array, inv2b_array, coords, orients, dist_lim, omega) -> float:
+    if inv2a_array[0] >= 1:
+        n = len(coords)
+        return clib.contact_bidisp_pack_multi(n, inv2a_array, inv2b_array, coords, orients, dist_lim, omega)
+    raise ValueError("Semi-major axis should be smaller that pi/2, otherwise the ellipse becomes inverted.")
+
+
+class EllipseMinim(ABC):
+
+    """Base class fir optimization of ellipse configurations on a sphere."""
+
+    @abstractmethod
+    def change_params(self, new_alpha: float) -> None:
+        pass
+
+    @abstractmethod
+    def evaluate_energy(self) -> float:
+        pass
+
+    @abstractmethod
+    def export_config(self, folder: str, idx: int) -> None:
+        pass
+
+    @abstractmethod
+    def minim(self, eps, gtol) -> OptimizeResult:
+        """
+        Method for finding optimal configuration with BFGS energy minimization method.
+        :param eps: Absolute step size used for numerical approximation of the jacobian in the BFGS method.
+        :param gtol: BFGS minimization is terminated after all gradient components are smaller than gtol (inf norm)
+        """
+        pass
+
+
+class EllipseMinimMonodisp(EllipseMinim):
+
+    def __init__(self, n, alpha, epsilon, initial: np.ndarray = None):
+        """
+        :param n: number of spherical ellipses
+        :param alpha: ellipse size (geodesic semi-major axis)
+        :param epsilon: ellipse aspect ratio
+        """
+        self.n = n
+        self.alpha = alpha
+        self.beta = alpha / epsilon
+        self.epsilon = epsilon
+
+        self.config_xi = initial
+        if self.config_xi is None:
+            np.random.seed()
+            self.config_xi = np.array([2 * np.pi * np.random.random(self.n),
+                                       np.arccos(2 * np.random.random(self.n) - 1),
+                                       2 * np.pi * np.random.random(self.n)]).flatten()
+
+        self.inv2a = 1 / np.sin(self.alpha) ** 2
+        self.inv2b = 1 / np.sin(self.beta) ** 2
+        self.dist_lim = np.cos(2 * self.alpha)
+
+    def change_params(self, new_alpha):
+        self.alpha = new_alpha
+        self.beta = new_alpha / self.epsilon
+        self.inv2a = 1 / np.sin(self.alpha) ** 2
+        self.inv2b = 1 / np.sin(self.beta) ** 2
+        self.dist_lim = np.cos(2 * self.alpha)
+
+    def evaluate_energy(self):
+        sph = self.config_xi[:2 * self.n].reshape((2, self.n)).T
+        positions = sph_to_cart(sph)
+        orientations = vector_orientations(sph, self.config_xi[2 * self.n:])
+        return monodisp_config_energy(self.inv2a, self.inv2b, positions, orientations, self.dist_lim)
+
+    def export_config(self, folder: str, idx: int):
+        sph = self.config_xi[:2 * self.n].reshape((2, self.n)).T
+        positions = sph_to_cart(sph)
+        orientations = vector_orientations(sph, self.config_xi[2 * self.n:])
+        np.savetxt(folder + f'/coord{idx}.dat', positions)
+        np.savetxt(folder + f'/orient{idx}.dat', orientations)
+
+    def minim(self, eps=1.4901161193847656e-08, gtol=1e-5):
+
+        def energy(xi):
+            sph = xi[:2 * self.n].reshape((2, self.n)).T
+            positions = sph_to_cart(sph)
+            orientations = vector_orientations(sph, xi[2 * self.n:])
+            return monodisp_config_energy(self.inv2a, self.inv2b, positions, orientations, self.dist_lim)
+
+        result = minim(energy, self.config_xi, method='BFGS', options={'gtol': gtol, 'eps': eps})
+        self.config_xi = result.x
+
+        return result
+
+
+class EllipseMinimBidisp(EllipseMinim):
+
+    def __init__(self, n, alpha1, alpha2, epsilon, num_frac=0.5, initial=None):
+        """
+        :param n: number of spherical ellipses
+        :param alpha1: geodesic semi-major axis for the 1st size of ellipses
+        :param alpha2: geodesic semi-major axis for the 2nd size of ellipses
+        :param epsilon: ellipse aspect ratio
+        :param num_frac: number fraction of larger ellipses
+        """
+
+        alpha1, alpha2 = max(alpha1, alpha2), min(alpha1, alpha2)
+        if epsilon < 1:
+            epsilon = 1 / epsilon
+
+        self.n = n
+        self.alpha1 = alpha1
+        self.alpha2 = alpha2
+        self.beta1 = alpha1 / epsilon
+        self.beta2 = alpha2 / epsilon
+        self.num_frac = num_frac
+        self.size_ratio = alpha1 / alpha2
+        self.epsilon = epsilon
+
+        self.config_xi = initial
+        if self.config_xi is None:
+            np.random.seed()
+            self.config_xi = np.array([2 * np.pi * np.random.random(self.n),
+                                       np.arccos(2 * np.random.random(self.n) - 1),
+                                       2 * np.pi * np.random.random(self.n)]).flatten()
+
+        self.num_large = int(n * num_frac)
+
+        self.inv2a = np.zeros(self.n)
+        self.inv2b = np.zeros(self.n)
+        self.inv2a[:self.num_large] = 1 / np.sin(self.alpha1) ** 2
+        self.inv2a[self.num_large:] = 1 / np.sin(self.alpha2) ** 2
+        self.inv2b[:self.num_large] = 1 / np.sin(self.beta1) ** 2
+        self.inv2b[self.num_large:] = 1 / np.sin(self.beta2) ** 2
+
+        self.dist_lim = np.cos(2 * self.alpha1)
+        self.omega = 1 / np.sin(self.alpha1) ** 2
+
+    def change_params(self, new_alpha1):
+        self.alpha1 = new_alpha1
+        self.alpha2 = new_alpha1 / self.size_ratio
+        self.beta1 = self.alpha1 / self.epsilon
+        self.beta2 = self.alpha2 / self.epsilon
+
+        self.inv2a[:self.num_large] = 1 / np.sin(self.alpha1) ** 2
+        self.inv2a[self.num_large:] = 1 / np.sin(self.alpha2) ** 2
+        self.inv2b[:self.num_large] = 1 / np.sin(self.beta1) ** 2
+        self.inv2b[self.num_large:] = 1 / np.sin(self.beta2) ** 2
+
+        self.dist_lim = np.cos(2 * self.alpha1)
+
+    def evaluate_energy(self):
+        sph = self.config_xi[:2 * self.n].reshape((2, self.n)).T
+        positions = sph_to_cart(sph)
+        orientations = vector_orientations(sph, self.config_xi[2 * self.n:])
+        return bidisp_config_energy(self.inv2a, self.inv2b, positions, orientations, self.dist_lim, self.omega)
+
+    def export_config(self, folder: str, idx: int):
+        sph = self.config_xi[:2 * self.n].reshape((2, self.n)).T
+        positions = sph_to_cart(sph)
+        orientations = vector_orientations(sph, self.config_xi[2 * self.n:])
+        np.savetxt(folder + f'/coord{idx}.dat', positions)
+        np.savetxt(folder + f'/orient{idx}.dat', orientations)
+
+    def minim(self, eps=1.4901161193847656e-08, gtol=1e-5):
+
+        def energy_eig(xi):
+            sph = xi[:2 * self.n].reshape((2, self.n)).T
+            positions = sph_to_cart(sph)
+            orientations = vector_orientations(sph, xi[2 * self.n:])
+            return bidisp_config_energy(self.inv2a, self.inv2b, positions, orientations, self.dist_lim, self.omega)
+
+        result = minim(energy_eig, self.config_xi, method='BFGS',  options={'gtol': gtol, 'eps': eps})
+        self.config_xi = result.x
+
+        return result
+
+
+def ellipse_minim_factory(n, alpha1, alpha2, epsilon, num_frac=0.5) -> EllipseMinim:
+    if alpha1 == alpha2:
+        return EllipseMinimMonodisp(n, alpha1, epsilon)
+    return EllipseMinimBidisp(n, alpha1, alpha2, epsilon, num_frac=num_frac)
+
+
+class PackingSimulation:
+
+    def __init__(self, n: int, epsilon: float, data_folder: str = None, size_ratio: float = 1., num_frac: float = 0.5):
+        self.n = n
+        self.epsilon = epsilon
+        self.num_frac = num_frac
+
+        self.size_ratio = size_ratio
+        if self.size_ratio < 1:
+            self.size_ratio = 1 / self.size_ratio
+
+        self.data_folder = data_folder
+        if self.data_folder is None:
+            self.data_folder = os.getcwd()
+
+        if size_ratio == 1:
+            self.folder = self.data_folder + f'/packing_results/monodisp/n{n}ratio{int(1000 * epsilon)}'
+        else:
+            self.folder = self.data_folder + \
+            f'/packing_results/bidisp{int(100 * size_ratio)}frac{int(100 * num_frac)}/n{n}ratio{int(1000 * epsilon)}'
+
+        try:
+            existing = np.loadtxt(self.folder + '/ellipse_sizes.dat')
+            self.existing_idx = len(existing)
+        except OSError:
+            self.existing_idx = 0
+        except TypeError:
+            self.existing_idx = 1
+
+    def run(self, num_configs, tol=1e-8, n_jobs=8):
+        os.makedirs(self.folder, exist_ok=True)
+
+        if self.existing_idx != 0:
+            print(f'{self.existing_idx} configurations for given parameters already exist.')
+
+        threads = multiprocessing.Pool(n_jobs)
+        ellipse_sizes = threads.map(partial(ellipse_packing, n=self.n, epsilon=self.epsilon, tol=tol,
+                                            size_ratio=self.size_ratio, num_frac=self.num_frac,
+                                            folder=self.folder, existing_idx=self.existing_idx),
+                                    range(num_configs))
+        threads.close()
+
+        with open(self.folder + '/ellipse_sizes.dat', 'ab') as f:
+            np.savetxt(f, ellipse_sizes)
+
+        self.existing_idx += num_configs
+
+
+@dataclass
+class PackingSimData:
+
+    minim_iter: list = field(init=False, default_factory=list)
+    minim_nfev: list = field(init=False, default_factory=list)
+    ellipse_scaling: list = field(init=False, default_factory=list)
+    scaling_step: list = field(init=False, default_factory=list)
+    energy: list = field(init=False, default_factory=list)
+
+    def append_all(self, minim_result: OptimizeResult, ellipse_scaling: float, scaling_step: float):
+        self.minim_iter.append(minim_result.nit)
+        self.minim_nfev.append(minim_result.nfev)
+        self.ellipse_scaling.append(ellipse_scaling)
+        self.scaling_step.append(scaling_step)
+        self.energy.append(minim_result.fun)
+
+    def export(self, folder: str, idx: int):
+        minim_iter = np.array(self.minim_iter)
+        minim_nfev = np.array(self.minim_nfev)
+        ellipse_scaling = np.array(self.ellipse_scaling)
+        scaling_step = np.array(self.scaling_step)
+        energy = np.array(self.energy)
+        np.savetxt(folder + f'/iter_data{idx}.dat',
+                   np.vstack((ellipse_scaling, scaling_step, minim_iter, minim_nfev, energy)).T)
+
+
+def ellipse_packing(idx, epsilon: float, n: int, tol: float, size_ratio: float,
+                    num_frac: float, folder: str, existing_idx: int = 0):
+
+    scaling = np.sqrt(2 * epsilon / n)  # we start with scaling for packing fraction cca phi<0.5
+    scaling_step0 = scaling / 100
+    scaling_step = scaling_step0
+
+    inflate_deflate = 1
+    counter = 0
+    n_iter = 0
+
+    # default values
+    eps = 1.4901161193847656e-08
+    gtol = 1e-5
+
+    packing_sim_data = PackingSimData()
+
+    # random initialization of ellipse positions and orientations
+    ellipse_minim = ellipse_minim_factory(n, scaling, scaling / size_ratio, epsilon, num_frac=num_frac)
+    minim_result = ellipse_minim.minim(eps=eps, gtol=gtol)
+    packing_sim_data.append_all(minim_result, scaling, scaling_step)
+
+    if minim_result.fun > 1e-4:
+        raise ValueError("Check initialization (packing fraction), the first minimization should get rid of overlaps.")
+
+    # first: lowering of eps and gtol as long as minimization does not reach the desired accuracy
+    while True:
+        if minim_result.fun < 1e-12 or eps < 2 * 1e-10:
+            break
+
+        eps = eps / 2
+        print(f'Lowering eps at first minim for idx {idx}, now {eps}')
+        minim_result = ellipse_minim.minim(eps=eps, gtol=gtol)
+        packing_sim_data.append_all(minim_result, scaling, scaling_step)
+
+        # lowering the desired accuracy of gradient calculation if lowering eps does not help
+        delta_en = packing_sim_data.energy[-2] - packing_sim_data.energy[-1]
+        if delta_en < packing_sim_data.energy[-2] / 10 and gtol > 2 * 1e-7:
+            gtol = gtol / np.sqrt(10)
+            print(f'Lowering gtol at first minim for idx {idx}, now {gtol}')
+
+    while True:
+        if n_iter == 300:
+            print('Exceeded max number of scaling iterations (n_iter_max=300).')
+            break
+
+        energy = ellipse_minim.evaluate_energy()
+        if tol < energy < 2 * tol:
+            break
+        elif energy < tol and inflate_deflate == -1:
+            scaling_step = scaling_step / 2
+            inflate_deflate = 1  # ellipse size is increased
+            counter = 0
+        elif energy > 2 * tol and inflate_deflate == 1:
+            scaling_step = scaling_step / 2
+            inflate_deflate = -1  # ellipse size is decreased
+            counter = 0
+        # if there is no change in scaling step for 8 minimizations, we increase it
+        elif counter == 8 and scaling_step < scaling_step0:
+            scaling_step = 2 * scaling_step
+            counter = 0
+
+        # additional adjustment of eps and gtol:
+        # we mandate strict energy accuracy during the first 10 iterations (there should still be no overlaps)
+        if n_iter <= 10:
+            if energy > 2 * 1e-12 and eps > 2 * 1e-10:
+                eps = eps / 2
+                print(f'Lowering eps at {n_iter + 1}. minim for idx {idx}, now {eps}')
+                minim_result = ellipse_minim.minim(eps=eps, gtol=gtol)
+                packing_sim_data.append_all(minim_result, scaling, scaling_step)
+
+                delta_en = packing_sim_data.energy[-2] - packing_sim_data.energy[-1]
+                if delta_en < packing_sim_data.energy[-2] / 10 and gtol > 2 * 1e-7:
+                    gtol = gtol / np.sqrt(10)
+                    print(f'Lowering gtol at {n_iter + 1}. minim minim for idx {idx}, now {gtol}')
+
+            if energy > tol:
+                raise ValueError('Check initialization, there should still be no overlaps for the first 10 iterations.')
+
+        # increase (or decrease) the size of ellipses
+        scaling += inflate_deflate * scaling_step
+        ellipse_minim.change_params(scaling)
+
+        minim_result = ellipse_minim.minim(eps=eps, gtol=gtol)
+        packing_sim_data.append_all(minim_result, scaling, scaling_step)
+        counter += 1
+        n_iter += 1
+
+    # save packing configuration and simulation data
+    ellipse_minim.export_config(folder, idx=idx + existing_idx)
+    packing_sim_data.export(folder, idx=idx + existing_idx)
+
+    return scaling