gsl_reg.c

Go to the documentation of this file.

/* src/gsl_reg.c
 *
 * Copyright (C) 2011-2018 Dongdong Li
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundataion; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABLITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>

struct DATA_POINT {
  size_t n;
  EXPAW *expaw;
  EXPPH *expph;
  EXPSOL *expsol;
};

void
print_state (size_t iter, gsl_multifit_fdfsolver *s)
{
  int i;
  printf ("iter: %4lu |f(x)| = %g\n",
          iter, gsl_blas_dnrm2 (s->f));

  for (i = 0; i < fitted_param_num; i ++)
    {
      printf ("%0.10E\n", gsl_vector_get (s->x, i));
    }
}

void set_aqueous_species ()
{
  int species_num_index = 0;

  int i;
  for (i = 0; i < gas_phase_num; i ++)
    {
      species_num_index += phases.gas[i].species_num;
    }

  for (i = 0; i < aqueous_phase_num; i ++)
    {
      int j;

      for (j = 0; j < phases.aqueous[i].species_num; j++)
    {
      naq_eq[i][j] = n_eq[j + species_num_index];
    }

      species_num_index += phases.aqueous[i].species_num;
    }
}

int init_fitted_params ()
{
  int i, j, k;

  fitted_param = (FITPARAM *)malloc (sizeof (FITPARAM) * fitted_param_num);

  int index = 0;

  for (i = 0; i < aqueous_phase_num; i ++)
    {
      for (j = 0; j < phases.aqueous[i].ii_bin_param_num; j ++)
    {
      if (phases.aqueous[i].ii_bin_params[j].alpha > 0 && 
          phases.aqueous[i].ii_bin_params[j].alpha < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 0;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 0;
          fitted_param[index].coeff_index = 0;

          index ++;
        }
      if (phases.aqueous[i].ii_bin_params[j].alpha1 > 0 && 
          phases.aqueous[i].ii_bin_params[j].alpha1 < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 0;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 1;
          fitted_param[index].coeff_index = 0;

          index ++;
        }

      for (k = 0; k < 6; k ++)
        {
          if (phases.aqueous[i].ii_bin_params[j].B0ss[k] > 0 && 
          phases.aqueous[i].ii_bin_params[j].B0ss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 0;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 2;
          fitted_param[index].coeff_index = k;

          index ++;
        }
          if (phases.aqueous[i].ii_bin_params[j].B1ss[k] > 0 && 
          phases.aqueous[i].ii_bin_params[j].B1ss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 0;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 3;
          fitted_param[index].coeff_index = k;

          index ++;
        }
        }
    }

      for (j = 0; j < phases.aqueous[i].nn_bin_param_num; j ++)
    {
      for (k = 0; k < 6; k ++)
        {
          if (phases.aqueous[i].nn_bin_params[j].Wss[k] > 0 && 
          phases.aqueous[i].nn_bin_params[j].Wss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 1;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 0;
          fitted_param[index].coeff_index = k;

          index ++;
        }
          if (phases.aqueous[i].nn_bin_params[j].Uss[k] > 0 && 
          phases.aqueous[i].nn_bin_params[j].Uss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 1;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 1;
          fitted_param[index].coeff_index = k;

          index ++;
        }
        }
    }

      for (j = 0; j < phases.aqueous[i].nii_ter_param_num; j ++)
    {
      for (k = 0; k < 6; k ++)
        {
          if (phases.aqueous[i].nii_ter_params[j].Wsss[k] > 0 && 
          phases.aqueous[i].nii_ter_params[j].Wsss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 2;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 0;
          fitted_param[index].coeff_index = k;

          index ++;
        }
          if (phases.aqueous[i].nii_ter_params[j].Usss[k] > 0 && 
          phases.aqueous[i].nii_ter_params[j].Usss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 2;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 1;
          fitted_param[index].coeff_index = k;

          index ++;
        }

          if (phases.aqueous[i].nii_ter_params[j].Vsss[k] > 0 && 
          phases.aqueous[i].nii_ter_params[j].Vsss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 2;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 2;
          fitted_param[index].coeff_index = k;

          index ++;
        }
        }
    }

      for (j = 0; j < phases.aqueous[i].iii_ter_param_num; j ++)
    {
      for (k = 0; k < 6; k ++)
        {
          if (phases.aqueous[i].iii_ter_params[j].Wsss[k] > 0 && 
          phases.aqueous[i].iii_ter_params[j].Wsss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 3;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 0;
          fitted_param[index].coeff_index = k;

          index ++;
        }
        }
    }

      for (j = 0; j < phases.aqueous[i].niii_qua_param_num; j ++)
    {
      for (k = 0; k < 6; k ++)
        {
          if (phases.aqueous[i].niii_qua_params[j].Qssss[k] > 0 && 
          phases.aqueous[i].niii_qua_params[j].Qssss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 4;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 0;
          fitted_param[index].coeff_index = k;

          index ++;
        }

          if (phases.aqueous[i].niii_qua_params[j].Yssss[k] > 0 && 
          phases.aqueous[i].niii_qua_params[j].Yssss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 4;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 1;
          fitted_param[index].coeff_index = k;

          index ++;
        }
        }
    }

      for (j = 0; j < phases.aqueous[i].nnii_qua_param_num; j ++)
    {
      for (k = 0; k < 6; k ++)
        {
          if (phases.aqueous[i].nnii_qua_params[j].Yssss[k] > 0 && 
          phases.aqueous[i].nnii_qua_params[j].Yssss[k] < 1.0E-20)
        {
          fitted_param[index].value = 0.0;
          fitted_param[index].aq_index = i;
          fitted_param[index].param_type = 5;
          fitted_param[index].param_index = j;
          fitted_param[index].param_sub_type = 0;
          fitted_param[index].coeff_index = k;

          index ++;
        }
        }
    }
    }

  if (index != fitted_param_num)
    {
      printf ("Inconsistent Fitted Parameter Number!!!\n");
      exit (0);
    }
}

void renew_fitted_params (const gsl_vector *x)
{
  int i;

  for (i = 0; i < fitted_param_num; i ++)
    {
      fitted_param[i].value = gsl_vector_get(x, i);
    }

  for (i = 0; i < fitted_param_num; i ++)
    {
      switch (fitted_param[i].param_type)
    {
    case 0: // ii
      {
        if (fitted_param[i].param_sub_type == 0)
          {
        phases.aqueous[fitted_param[i].aq_index].ii_bin_params[fitted_param[i].param_index].alpha = fitted_param[i].value;
          }
        if (fitted_param[i].param_sub_type == 1)
          {
        phases.aqueous[fitted_param[i].aq_index].ii_bin_params[fitted_param[i].param_index].alpha1 = fitted_param[i].value;
          }
        if (fitted_param[i].param_sub_type == 2)
          {
        phases.aqueous[fitted_param[i].aq_index].ii_bin_params[fitted_param[i].param_index].B0ss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }
        if (fitted_param[i].param_sub_type == 3)
          {
        phases.aqueous[fitted_param[i].aq_index].ii_bin_params[fitted_param[i].param_index].B1ss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }

        break;
      }
    case 1: // nn
      {
        if (fitted_param[i].param_sub_type == 0)
          {
        phases.aqueous[fitted_param[i].aq_index].nn_bin_params[fitted_param[i].param_index].Wss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }
        if (fitted_param[i].param_sub_type == 1)
          {
        phases.aqueous[fitted_param[i].aq_index].nn_bin_params[fitted_param[i].param_index].Uss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }

        break;
      }
    case 2: // nii
      {
        if (fitted_param[i].param_sub_type == 0)
          {
        phases.aqueous[fitted_param[i].aq_index].nii_ter_params[fitted_param[i].param_index].Wsss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }
        if (fitted_param[i].param_sub_type == 1)
          {
        phases.aqueous[fitted_param[i].aq_index].nii_ter_params[fitted_param[i].param_index].Usss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }
        if (fitted_param[i].param_sub_type == 2)
          {
        phases.aqueous[fitted_param[i].aq_index].nii_ter_params[fitted_param[i].param_index].Vsss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }

        break;
      }
    case 3: // iii
      {
        if (fitted_param[i].param_sub_type == 0)
          {
        phases.aqueous[fitted_param[i].aq_index].iii_ter_params[fitted_param[i].param_index].Wsss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }

        break;
      }
    case 4: // niii
      {
        if (fitted_param[i].param_sub_type == 0)
          {
        phases.aqueous[fitted_param[i].aq_index].niii_qua_params[fitted_param[i].param_index].Qssss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }
        if (fitted_param[i].param_sub_type == 1)
          {
        phases.aqueous[fitted_param[i].aq_index].niii_qua_params[fitted_param[i].param_index].Yssss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }

        break;
      }
    case 5: // nnii
      {
        if (fitted_param[i].param_sub_type == 0)
          {
        phases.aqueous[fitted_param[i].aq_index].nnii_qua_params[fitted_param[i].param_index].Yssss[fitted_param[i].coeff_index] = fitted_param[i].value;
          }

        break;
      }
    }
    }
}

double objfun (const gsl_vector *x, void *data_point, int i)
{
  int exp_sum = 0;
  double OF = 0;

  renew_fitted_params (x);

  if (i < total_exp_aw_num)
    {
      EXPAW *expaw = ((struct DATA_POINT *)data_point)->expaw;

      int j, k;

      n_eq = (double *)malloc (sizeof(double) * total_species_num);
      for (k = 0; k < total_species_num; k ++)
    {
      n_eq[k] = 1.0e-10;
    }

      system_T = expaw[i].T;
      system_P = expaw[i].P;
      system_charge = 0;

      for (j = 0; j < total_comp_num; j ++)
    {
      total_components[j] = expaw[i].component[j];
    }

      // speciation equilibrium
      gem_ipopt ();
      if (!successed) {return 0.0;}

      set_aqueous_species ();

      for (k = 0; k < phases.aqueous[expaw[i].index].species_num; k ++)
    {
      if (phases.aqueous[expaw[i].index].aqueous_species[k].charge == 0 && 
          strcmp (phases.aqueous[expaw[i].index].aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
        {
          break;
        }
    }

      double aw = psc_a (phases.aqueous[expaw[i].index], k, naq_eq[expaw[i].index], expaw[i].T, expaw[i].P);

      OF = 1.0/expaw[i].sigma * (log(aw) - log(expaw[i].value));
    }

  exp_sum += total_exp_aw_num;
  if (i >= exp_sum && i < exp_sum + total_exp_ph_num)
    {
      EXPPH *expph = ((struct DATA_POINT *)data_point)->expph;

      int j, k;

      n_eq = (double *)malloc (sizeof(double) * total_species_num);
      for (k = 0; k < total_species_num; k ++)
    {
      n_eq[k] = 1.0e-10;
    }

      system_T = expph[i - exp_sum].T;
      system_P = expph[i - exp_sum].P;
      system_charge = 0;

      for (j = 0; j < total_comp_num; j ++)
    {
      total_components[j] = expph[i - exp_sum].component[j];
    }

      // speciation equilibrium
      gem_ipopt ();
      if (!successed) {return 0.0;}

      set_aqueous_species ();

      for (k = 0; k < phases.aqueous[expph[i - exp_sum].index].species_num; k ++)
    {
      if (phases.aqueous[expph[i - exp_sum].index].aqueous_species[k].charge != 0 && 
          strcmp (phases.aqueous[expph[i - exp_sum].index].aqueous_species[k].species_symbol, H_SYMBOL) == 0)
        {
          break;
        }
    }

      double ah = psc_a (phases.aqueous[expph[i - exp_sum].index], k, naq_eq[expph[i - exp_sum].index], expph[i - exp_sum].T, expph[i - exp_sum].P);

      OF = 1.0/expph[i].sigma * (-log(ah)/log(10) - expph[i - exp_sum].value);
    }

  exp_sum += total_exp_ph_num;
  if (i >= exp_sum && i < exp_sum + total_exp_sol_num)
    {
      EXPSOL *expsol = ((struct DATA_POINT *)data_point)->expsol;

      int j, k;

      n_eq = (double *)malloc (sizeof(double) * total_species_num);
      for (k = 0; k < total_species_num; k ++)
    {
      n_eq[k] = 1.0e-10;
    }

      system_T = expsol[i - exp_sum].T;
      system_P = expsol[i - exp_sum].P;
      system_charge = 0;

      for (j = 0; j < total_comp_num; j ++)
    {
      total_components[j] = expsol[i - exp_sum].component[j];
    }

      // speciation equilibrium
      gem_ipopt ();
      if (!successed) {return 0.0;}

      set_aqueous_species ();

      double *ai = (double *)malloc(sizeof (double) * phases.aqueous[expsol[i - exp_sum].index].species_num);

      for (k = 0; k < phases.aqueous[expsol[i - exp_sum].index].species_num; k ++)
    {
      ai[k] = psc_a (phases.aqueous[expsol[i - exp_sum].index], k, 
             naq_eq[expsol[i - exp_sum].index], 
             expsol[i - exp_sum].T, 
             expsol[i - exp_sum].P);

      OF += (log(ai[k])) * expsol[i - exp_sum].species[k];
    }

      OF = 1.0/expsol[i - exp_sum].sigma * (OF - expsol[i - exp_sum].lnk);
    }

  return OF;
}

int objfun_f (const gsl_vector *x, void *data_point, gsl_vector *f)
{
  EXPAW *expaw = ((struct DATA_POINT *)data_point)->expaw;
  EXPPH *expph = ((struct DATA_POINT *)data_point)->expph;
  EXPSOL *expsol = ((struct DATA_POINT *)data_point)->expsol;

  double Pi;
  int i, j, k;
  int exp_sum = 0;

  n_eq = (double *)malloc (sizeof(double) * total_species_num);
  for (k = 0; k < total_species_num; k ++)
    {
      n_eq[k] = 1.0e-10;
    }

  renew_fitted_params (x);

  for (i = 0; i < total_exp_aw_num; i ++)
    {
      Pi = 0;
      system_T = expaw[i].T;
      system_P = expaw[i].P;
      system_charge = 0;

      for (j = 0; j < total_comp_num; j ++)
    {
      total_components[j] = expaw[i].component[j];
    }

      // speciation equilibrium
      gem_ipopt ();
      //print_results ();

      set_aqueous_species ();

      for (k = 0; k < phases.aqueous[expaw[i].index].species_num; k ++)
    {
      if (phases.aqueous[expaw[i].index].aqueous_species[k].charge == 0 && 
          strcmp (phases.aqueous[expaw[i].index].aqueous_species[k].species_symbol, WATER_SYMBOL) == 0)
        {
          break;
        }
    }

      double aw = psc_a (phases.aqueous[expaw[i].index], k, naq_eq[expaw[i].index], expaw[i].T, expaw[i].P);

      Pi = 1.0/expaw[i].sigma * (log(aw) - log(expaw[i].value));

      if (!successed) {Pi = 0.0;}

      gsl_vector_set (f, i+exp_sum, Pi);
    }

  exp_sum += total_exp_aw_num;
  for (i = 0; i < total_exp_ph_num; i ++)
    {
      Pi = 0;
      system_T = expph[i].T;
      system_P = expph[i].P;
      system_charge = 0;

      for (j = 0; j < total_comp_num; j ++)
    {
      total_components[j] = expph[i].component[j];
    }

      // speciation equilibrium
      gem_ipopt ();
      //print_results ();

      set_aqueous_species ();

      for (k = 0; k < phases.aqueous[expph[i].index].species_num; k ++)
    {
      if (phases.aqueous[expph[i].index].aqueous_species[k].charge != 0 && 
          strcmp (phases.aqueous[expph[i].index].aqueous_species[k].species_symbol, H_SYMBOL) == 0)
        {
          break;
        }
    }

      double ah = psc_a (phases.aqueous[expph[i].index], k, naq_eq[expph[i].index], expph[i].T, expph[i].P);

      Pi = 1.0/expph[i].sigma * (-log(ah)/log(10) - expph[i].value);

      gsl_vector_set (f, i + exp_sum, Pi);
    }

  exp_sum += total_exp_ph_num;
  for (i = 0; i < total_exp_sol_num; i ++)
    {
      Pi = 0;
      system_T = expsol[i].T;
      system_P = expsol[i].P;
      system_charge = 0;

      for (j = 0; j < total_comp_num; j ++)
    {
      total_components[j] = expsol[i].component[j];
    }

      // speciation equilibrium
      gem_ipopt ();

      set_aqueous_species ();

      double *ai = (double *)malloc(sizeof (double) * phases.aqueous[expsol[i].index].species_num);

      for (k = 0; k < phases.aqueous[expsol[i].index].species_num; k ++)
    {
      ai[k] = psc_a (phases.aqueous[expsol[i].index], k, 
             naq_eq[expsol[i].index], 
             expsol[i].T, 
             expsol[i].P);

      Pi += (log(ai[k])) * expsol[i].species[k];
    }

      Pi = 1.0/expsol[i].sigma * (Pi - expsol[i].lnk);

      if (!successed) {Pi = 0.0;}

      gsl_vector_set (f, i + exp_sum, Pi);
    }

  return GSL_SUCCESS;
}

int objfun_df (const gsl_vector *x, void *data_point, gsl_matrix *J)
{
  double h = 1.0e-6;
  double F1, F2;
  int i, j;

  for (i = 0; i < total_exp_aw_num + total_exp_ph_num + total_exp_sol_num; i ++)
    {
      for (j = 0; j < fitted_param_num; j ++)
    {
      F1 = objfun (x, data_point, i);

      x->data[j] = x->data[j] + h;
      F2 = objfun (x, data_point, i);

      double Jij = (F2 - F1) / h;

      if (F2 == 0 || F1 == 0) {Jij = 0;}

      gsl_matrix_set (J, i, j, Jij);

      //printf ("%d\t%d\t%lf\n", i, j, Jij);
    }
    }

  return GSL_SUCCESS;
}

int objfun_fdf (const gsl_vector *x, void *data_point, gsl_vector *f, gsl_matrix *J)
{
  objfun_f (x, data_point, f);
  objfun_df (x, data_point, J);

  return GSL_SUCCESS;
}

void gsl_reg()
{
  init_fitted_params ();

  const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
  gsl_multifit_fdfsolver *s;

  int iter = 1;
  int status;
  size_t i;
  const size_t N = total_exp_aw_num + total_exp_ph_num + total_exp_sol_num;
  const size_t P = fitted_param_num;
  double chi, chi0, dof;

  total_components = (double *)malloc(sizeof (double) * total_comp_num);
  naq_eq = (double **)malloc(sizeof(double*) * aqueous_phase_num);

  for (i = 0; i < aqueous_phase_num; i ++)
    {
      naq_eq[i] = (double *)malloc(sizeof(double) * phases.aqueous[i].species_num);
    }

  struct DATA_POINT d = {N, expaw, expph, expsol};

  double *x_init = (double *)malloc(sizeof (double) * P);

  for(i = 0; i < fitted_param_num; i ++)
    {
      x_init[i] = 1.0;
    }

  gsl_vector_view x = gsl_vector_view_array (x_init, P);

  gsl_matrix *J = gsl_matrix_alloc (N, P);
  gsl_multifit_function_fdf f;

  f.f = &objfun_f;
  f.df = &objfun_df;
  f.fdf = &objfun_fdf;
  f.n = N;
  f.p = P;
  f.params = &d;

  s = gsl_multifit_fdfsolver_alloc (T, N, P);
  gsl_multifit_fdfsolver_set (s, &f, &x.vector);
  print_state (iter, s);

  do
    {
      iter ++;
      status = gsl_multifit_fdfsolver_iterate (s);

      printf ("status = %s\n", gsl_strerror (status));

      print_state (iter, s);

      if (status)
    break;

      status = gsl_multifit_test_delta (s->dx, s->x, 1.0e-4, 1.0e-4);
    }
  while (status == GSL_CONTINUE && iter < 1000);

  printf ("status = %s\n", gsl_strerror (status));

  chi = gsl_blas_dnrm2 (s->f); // 2 norm
  dof = N - P; // degree of freedom 

  printf ("OF = %g\n", sqrt(chi / dof));

  renew_fitted_params (s->x);
}