Laboratoire RDP » Biophysics Team » Cell Division Rule

#include "compute_energy_and_forces.h"

#include "minimize_steepest_descent.h"

#include <nlopt.h>

#include <stdio.h>

#include <stdlib.h>

#include <cmath>

/* this function minimizes the mechanical energy. It is an alternative to the BFGS algorithm, just for debugging */

void minimize_steepest_descent(int vertices_number, int cell_number, double *vertices, int *cells_vertices, int *vertices_number_in_each_cell, double *barycenters, double *shapes, double *targets, double *areas, double *cell_areas, double *target_areas,

                   int *wall_vertices, int wall_number, double *lengths, double *target_lengths){

double step_size=1.e-2;

double precision=1.e-2;

struct parameters_list {int vertices_number; int cell_number; int *cells_vertices; int *vertices_number_in_each_cell; double *barycenters; double *shapes; double *targets; double *areas; double *cell_areas; double *target_areas;

                        int *wall_vertices; int wall_number; double *lengths; double *target_lengths;};

struct parameters_list parameters;

parameters.vertices_number=vertices_number;

parameters.cell_number=cell_number;

parameters.cells_vertices=cells_vertices;

parameters.vertices_number_in_each_cell=vertices_number_in_each_cell;

parameters.barycenters=barycenters;

parameters.shapes=shapes;

parameters.targets=targets;

parameters.areas=areas;

parameters.cell_areas=cell_areas;

parameters.target_areas=target_areas;

parameters.wall_vertices=wall_vertices;

parameters.wall_number=wall_number;

parameters.lengths=lengths;

parameters.target_lengths=target_lengths;

void* pparameters=&parameters;

double forces[2*vertices_number];

double energy=compute_energy_and_forces(0, vertices, forces, pparameters);

double max_step=0.; for (int i=0;i<2*vertices_number;i++){ if (fabs(forces[i])>max_step){ max_step=fabs(forces[i]);}}

int step_number=0;

while (max_step>precision){

    step_number++;

    for (int i=0;i<2*vertices_number;i++){ vertices[i]-=step_size*forces[i];}

    energy=compute_energy_and_forces(0, vertices, forces, pparameters);

    max_step=0;

    for (int i=0;i<2*vertices_number;i++){ if (fabs(forces[i])>max_step){ max_step=fabs(forces[i]);}}

}

printf("opt %f %d\n",energy,step_number);

}