Laboratoire RDP » Biophysics Team » Cell Division Rule

#include "compute_energy_and_forces.h"

#include "minimize_BFGS.h"

#include <nlopt.h>

#include <stdio.h>

#include <stdlib.h>

#include <cmath>

/* this function minimizes the mechanical energy */

void minimize_BFGS(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,

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

nlopt_opt opt;

opt=nlopt_create(NLOPT_LD_LBFGS,2*vertices_number);

nlopt_set_stopval(opt, -HUGE_VAL);

nlopt_set_maxeval(opt, -1.);

nlopt_set_maxtime(opt, -1.);

nlopt_set_xtol_rel(opt, -1.);

nlopt_set_xtol_abs1(opt, 1.e-6);

nlopt_set_ftol_rel(opt, -1.);

nlopt_set_ftol_abs(opt, -1.);

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;

                        int *wall_vertices; int wall_number; double *lengths; double *target_lengths; int *boundary; double *tissue_target; int anisotropy; double *relative_pressure;};

double min_energy; double tissue_target[2]={1.,1.};

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.wall_vertices=wall_vertices;

parameters.wall_number=wall_number;

parameters.lengths=lengths;

parameters.target_lengths=target_lengths;

parameters.boundary=boundary;

parameters.tissue_target=tissue_target;

parameters.anisotropy=0;

parameters.relative_pressure=relative_pressure;

void* pparameters=&parameters;

//double forces[2*vertices_number];

//double test=compute_energy_and_forces(0, vertices, forces, pparameters);

////for (int i=0; i<vertices_number;i++){ printf("%f %f %d %d\n",forces[2*i],forces[2*i+1],i,vertices_number);} printf("test\n");

//double test2;

//for (int i=0;i<2*vertices_number;i++){vertices[i]+=1.e-8;

//test2=compute_energy_and_forces(0, vertices, forces, pparameters);

//test2-=test;

//test2*=1.e8;

//vertices[i]-=1.e-8;

//printf("%f %f %f\n",test2,forces[i],(test2-forces[i])/forces[i]);

//}printf("\n");

int reso=nlopt_set_min_objective(opt, compute_energy_and_forces, pparameters);

nlopt_optimize(opt, vertices, &min_energy);

tissue_target[0]=0.; tissue_target[1]=0.;

int boundary_vertices_number=0.; double tissue_barycenter[2]={0.,0.};

for (int i=0;i<vertices_number_in_each_cell[cell_number];i++){

    tissue_barycenter[0]+=((double)boundary[i])*vertices[2*cells_vertices[i]];

    tissue_barycenter[1]+=((double)boundary[i])*vertices[2*cells_vertices[i]+1];

    boundary_vertices_number+=boundary[i];

}

tissue_barycenter[0]/=boundary_vertices_number;

tissue_barycenter[1]/=boundary_vertices_number;

for (int i=0;i<vertices_number_in_each_cell[cell_number];i++){

    tissue_target[0]+=((double)boundary[i])*(vertices[2*cells_vertices[i]]-tissue_barycenter[0])*(vertices[2*cells_vertices[i]]-tissue_barycenter[0]);

    tissue_target[1]+=((double)boundary[i])*(vertices[2*cells_vertices[i]+1]-tissue_barycenter[1])*(vertices[2*cells_vertices[i]+1]-tissue_barycenter[1]);

}

tissue_target[0]/=boundary_vertices_number;

tissue_target[1]/=boundary_vertices_number;

parameters.tissue_target=tissue_target;

parameters.anisotropy=1;

pparameters=&parameters;

reso=nlopt_set_min_objective(opt, compute_energy_and_forces, pparameters);

nlopt_optimize(opt, vertices, &min_energy);

nlopt_destroy(opt);

}