text: //go to: cd /cluster/bioprogs/modeller9v5_RSR/modlib/modeller 
  //gcc -shared -Wall -fPIC `/cluster/user/armin/bin/modeller9.10/bin/mod9.10 --cflags --libs` `pkg-config --cflags glib-2.0` -I/usr/include/python2.4 cuser_form54.c cuser_form54_wrap.c -o _cuser_form54.so -lm
  // new restraints with phylogenetic weight - python version 2.4, see also cuser_form59.c for python version 2.6
  
  #include <glib.h>
  #include "modeller.h"
  #include "cuser_form54.h"
  #include <stdio.h>
  #include <math.h>
  #include <stdlib.h>
  
  static const float RT = 6.0;
  
  /* Decode parameters from Modeller parameter array */
  inline static void get_param(const float pcsr[7], float *w, float *mean, float *stdev, float *w1, float *mean1, float *stdev1, float *phyloW)
  {
    *w = pcsr[0] + 1E-5;  
    *mean =  pcsr[1]; 	
    *stdev = pcsr[2];
    *w1 = pcsr[3]+1E-5;	
    *mean1 =  pcsr[4]; 	
    *stdev1 = pcsr[5];
    *phyloW = pcsr[6];
    //*phyloW = 1.0;
    //check if sum w+w1 equals 1 
    float sum = *w + *w1;
    *w /= sum;
    *w1 /= sum;
  }
  
  float get_p(float stdev, float x, float delta)
  {  
    return exp(-0.5 * delta * delta / (stdev*stdev)) / x;
  } 
  
  ///////////////////////
  /* Evaluate the form */
  ///////////////////////
  static int myform_eval(void *data, const float *feat, const int *iftyp,
                         const int *modal, int n_feat, const float *pcsr,
                         int n_pcsr, gboolean deriv, float *fderv, float *val)
  {
    float w, mean, stdev, w1, mean1, stdev1, phyloW;
    get_param(pcsr, &w, &mean, &stdev, &w1, &mean1, &stdev1, &phyloW);
    if (feat[0]<=0) {fprintf (stderr, "Error in cuser_form.c, myform_eval: negative argument of log!\n");return 1;}
    float logdist = log(feat[0]);
    float delta = logdist - mean;
    float delta1 = logdist - mean1;
  
    float arg = delta1*delta1/(2.*stdev1*stdev1) - delta*delta/(2.*stdev*stdev);
    float argBar = delta1/(stdev1*stdev1*feat[0]) - delta/(stdev*stdev*feat[0]);
  
  
    float com = w*stdev1/stdev * exp(arg);
    //printf("com=%f, arg=%f, argBar=%f, w=%f, feat=%f, logdist=%f, mean=%f, mean1=%f, stdev=%f, stdev1=%f\n", com, arg, argBar, w, feat[0], logdist, mean, mean1, stdev, stdev1);
  
  
    *val = RT * (-log(com + 1 - w) ) * phyloW;		
    
    
    if (deriv) {
      fderv[0] = RT * ((-1) * (com * argBar) / (com + 1 - w)) * phyloW;
    }
  
    return 0; 
  }
  
  
  // 1 ///////////////////////////////////
       /* Evaluate the minimum violation*/
       ///////////////////////////////////
  static int myform_vmin(void *data, const float *feat, const int *iftyp,
                         const int *modal, int n_feat, const float *pcsr,
                         int n_pcsr, float *val)
  {
    float w, mean, stdev, w1, mean1, stdev1, phyloW;
    get_param(pcsr, &w, &mean, &stdev, &w1, &mean1, &stdev1, &phyloW);
  
    if (feat[0]<=0) {
      fprintf (stderr, "Error in cuser_form.c, negative argument of log!\n");
      return 1;
    }
  
    float prec = 1.0/(stdev*stdev);
    float prec1 = 1.0/(stdev1*stdev1);
  
    float xmin = exp((mean1*prec1 - mean*prec)/(prec1 - prec));
  
    *val = feat[0] - xmin;
  
    /* if (abs(logdist-mean) < abs(logdist-mean1)) { */
  /*     *val = logdist - mean; */
  /*   } */
  /*   else { */
  /*     *val = logdist - mean1; */
  /*   } */
   
    return 0;
  }
  
  
  // 2 /////////////////////////////////
       /* Evaluate the heavy violation*/
       /////////////////////////////////
  static int myform_vheavy(void *data, const float *feat, const int *iftyp,
                         const int *modal, int n_feat, const float *pcsr,
                         int n_pcsr, float *val)
  {
    float w, mean, stdev, w1, mean1, stdev1, phyloW;
    get_param(pcsr, &w, &mean, &stdev, &w1, &mean1, &stdev1, &phyloW);
  
  
    if (feat[0]<=0) {
      fprintf (stderr, "Error in cuser_form.c, negative argument of log!\n");
      return 1;
    }
  
    float prec = 1.0/(stdev*stdev);
    float prec1 = 1.0/(stdev1*stdev1);
  
    float xmin = exp((mean1*prec1 - mean*prec)/(prec1-prec));
  
    *val = feat[0] - xmin;
  
    //search global minimum 
  
    /* if (get_p(stdev, logdist, logdist-mean) < get_p(stdev1, logdist, logdist - mean1)) { */
  /*     *val = logdist - mean1; */
  /*   } */
  /*   else { */
  /*     *val = logdist - mean; */
  /*   } */
    //float max = w*get_p(stdev,x,x-mean) + w1*get_p(stdev1,x,x-mean1);
  
    return 0;
  }
  
  // 3 ////////////////////////////////////////////
       /* Evaluate the relative minimum violation*/
       ////////////////////////////////////////////
  static int myform_rvmin(void *data, const float *feat, const int *iftyp,
                          const int *modal, int n_feat, const float *pcsr,
                          int n_pcsr, float *val)
  {
    float w, mean, stdev, w1, mean1, stdev1, phyloW;
    get_param(pcsr, &w, &mean, &stdev, &w1, &mean1, &stdev1, &phyloW);
  
    if (feat[0]<=0) {
      fprintf (stderr, "Error in cuser_form.c, negative argument of log!\n");
      return 1;
    }
  
    float logdist = log(feat[0]);
  
    
    if (abs(logdist-mean) < abs(logdist-mean1)) {
      *val = (logdist - mean)/stdev;
    }
    else {
      *val = (logdist - mean1)/stdev1;
    }
    
    return 0;
  }
  
  // 4 //////////////////////////////////////////
       /* Evaluate the relative heavy violation*/
       //////////////////////////////////////////
  static int myform_rvheavy(void *data, const float *feat, const int *iftyp,
                          const int *modal, int n_feat, const float *pcsr,
                          int n_pcsr, float *val)
  {
    float w, mean, stdev, w1, mean1, stdev1, phyloW;
    get_param(pcsr, &w, &mean, &stdev, &w1, &mean1, &stdev1, &phyloW);
  
    if (feat[0]<=0) {
      fprintf (stderr, "Error in cuser_form.c, negative argument of log!\n");
      return 1;
    }
  
    float logdist = log(feat[0]); 
  
    if (get_p(stdev, logdist, logdist-mean) < get_p(stdev1, logdist, logdist - mean1)) {
      *val = (logdist - mean1)/stdev1;
    }
    else {
      *val = (logdist - mean)/stdev;
    }
   
    return 0;
  }
  
   //////////////////////////////////////
   /* Evaluate the minimum feature mean*/
   //////////////////////////////////////
  static int myform_minmean(void *data, const float *feat, const int *iftyp,
                            const int *modal, int n_feat, const float *pcsr,
                            int n_pcsr, float *val)
  {
    float w, mean, stdev, w1, mean1, stdev1, phyloW;
    get_param(pcsr, &w, &mean, &stdev, &w1, &mean1, &stdev1, &phyloW);
  
    if (feat[0]<=0) {
      fprintf (stderr, "Error in cuser_form.c, negative argument of log!\n");
      return 1;
    }
  
    //search nearest minimum
    float prec = 1/(stdev*stdev);
    float prec1 = 1/(stdev1*stdev1);
  
    float x = exp( (mean1*prec1 - mean*prec) / (prec1 - prec) );
   
    *val = x; 
  
    return 0;
  }
  
   ////////////////////////////////////
   /* Evaluate the heavy feature mean*/
   ////////////////////////////////////
  static int myform_heavymean(void *data, const float *feat, const int *iftyp,
                            const int *modal, int n_feat, const float *pcsr,
                            int n_pcsr, float *val)
  {
    float w, mean, stdev, w1, mean1, stdev1, phyloW;
    get_param(pcsr, &w, &mean, &stdev, &w1, &mean1, &stdev1, &phyloW);
  
    if (feat[0]<=0) {
      fprintf (stderr, "Error in cuser_form.c, negative argument of log!\n");
      return 1;
    }
  
    //search global minimum 
    float prec = 1/(stdev*stdev);
    float prec1 = 1/(stdev1*stdev1);
  
    float x = exp( (mean1*prec1 - mean*prec) / (prec1 - prec) );
  
    *val = x;
  
    return 0;
  }
  
  ////////////////////////////////////////////////////
  /* Create the new form, and return its identifier */
  ////////////////////////////////////////////////////
  int myform_create(void)
  {  
    return mod_user_form_new(myform_eval, NULL, myform_vmin, NULL, myform_vheavy,
                             NULL, myform_rvmin, NULL, myform_rvheavy, NULL,
                             myform_minmean, NULL, myform_heavymean, NULL);
  }
  