text:
/**
 * This file is part of yasimSBML (http://www.labri.fr/perso/ghozlane/metaboflux/about/yasimSBML.php)
 * Copyright (C) 2010 Amine Ghozlane from LaBRI and University of Bordeaux 1
 *
 * yasimSBML is free software: you can redistribute it and/or modify
 * it under the terms of the Lesser GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * yasimSBML is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the Lesser GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * \file simulation.c
 * \brief Simulate a petri net
 * \author {Amine Ghozlane}
 * \version 1.0
 * \date 27 octobre 2009
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <time.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>
#include <sbml/SBMLTypes.h>
#include "especes.h"
#include "simulation.h"

/*TODO Verification des BoundaryConditions,Local parameter values, Global parameter values, Reaction, des equations */

/**
 * \fn void SBML_initEspeceAmounts(Model_t *mod, pEspeces molecules, int nbEspeces)
 * \author Amine Ghozlane
 * \brief  Alloc memory and initialize the struct Especes
 * \param  mod Model of the SBML file
 * \param  molecules Struct Especes
 * \param  nbEspeces Number of molecules
 */
void SBML_initEspeceAmounts(Model_t *mod, pEspeces molecules, int nbEspeces) {
  int i;
  Species_t *esp;

  /* Initialisation des quantites des especes*/
  for (i = 0; i < nbEspeces; i++) {
      esp = Model_getSpecies(mod, i);
      /*printf("espece: %s, compartiment : %s\n",Species_getId(esp),Species_getCompartment(esp));*/
      Especes_save(molecules, i, Species_getInitialAmount(esp),
          Species_getId(esp), Species_getCompartment(esp));
  }
}

/**
 * \fn void SBML_setReactions(Model_t *mod, pEspeces molecules, pScore result, double *reactions_ratio, int nbReactions, int nbEspeces)
 * \author Amine Ghozlane
 * \brief  Alloc memory and initialize the struct Especes
 * \param  mod Model of the SBML file
 * \param  molecules Struct Especes
 * \param  result Struct Score
 * \param  reactions_ratio List of computed reaction ratio
 * \param  nbReactions Number of reaction
 * \param  nbEspeces Number of molecules
 */
void SBML_setReactions(Model_t *mod, pEspeces molecules, int nbReactions, int nbEspeces) {
  int ref = 0, i, j;
  SpeciesReference_t *reactif;
  Species_t *especeId;
  Reaction_t *react;
  const char *kf;
  const ASTNode_t *km;
  KineticLaw_t *kl;

  /*fprintf(stdout, "Debut reaction\n");*/
  /* Recherche les reactions ou apparaissent chaque espece */
  for (i = 0; i < nbReactions; i++) {
      react = Model_getReaction(mod, i);
      for (j = 0; j < (int)Reaction_getNumReactants(react); j++) {
          reactif = Reaction_getReactant(react, j);
          especeId = Model_getSpeciesById(mod, SpeciesReference_getSpecies(
              reactif));
          ref = Especes_find(molecules, Species_getId(especeId), nbEspeces);
          /*printf("Reactif :Quantite de ref %d :%f \n", ref,Especes_getQuantite(molecules, ref));*/

          kl = Reaction_getKineticLaw(react);

          if (KineticLaw_isSetFormula(kl)) {
              kf = KineticLaw_getFormula(kl);
              /*printf("c du kf %s\n", kf);*/
              Especes_allocReactions(molecules, ref, react,
                  SBML_evalExpression(kf));
          } else {
              km = KineticLaw_getMath(kl);
              printf("C du kM... cas encore ignore\n");
              exit(1);
          }
      }
  }
  /*fprintf(stdout, "Fin reaction\n");*/
}

/**
 * \fn double SBML_evalExpression(const char *formule)
 * \author Amine Ghozlane
 * \brief  Get the reaction ratio define in the sbml
 * \param  formule Formule SBML
 * \return Return double value of the constraint
 */
double SBML_evalExpression(const char *formule) {
  return atof(formule);
}

/**
 * \fn int SBML_checkQuantite(Model_t *mod, Reaction_t *react, int nbEspeces, pEspeces molecules)
 * \author Amine Ghozlane
 * \brief  Determine the number of reaction for one molecule
 * \param  mod Model of the SBML file
 * \param  react Reaction id
 * \param  nbEspeces Number of molecules
 * \param  molecules Struct Especes
 * \return Number of reaction for one molecule
 */
int SBML_checkQuantite(Model_t *mod, Reaction_t *react, int nbEspeces, pEspeces molecules)
{
  double quantite = 0.0, minStep = 0.0, temp = 0.0;
  int ref = 0, i;
  SpeciesReference_t *reactif;
  Species_t *especeId;

  reactif = Reaction_getReactant(react, 0);
  especeId = Model_getSpeciesById(mod, SpeciesReference_getSpecies(reactif));
  ref = Especes_find(molecules, Species_getId(especeId), nbEspeces);

  if ((quantite = Especes_getQuantite(molecules, ref)) <= 0.0) {
      /*printf("je fais le return\n");*/
      return END;
  }
  /*printf("quantite : %d\n",quantite);*/
  minStep = quantite / SpeciesReference_getStoichiometry(reactif);
  /*printf("quantite : %d, minStep init : %d\n",quantite,minStep);*/
  /* Cas ou le nombre de pas est egal a 0 */
  if(minStep==0.0) return END;

  for (i = 1; i < (int)Reaction_getNumReactants(react); i++) {
      reactif = Reaction_getReactant(react, i);
      especeId = Model_getSpeciesById(mod, SpeciesReference_getSpecies(reactif));
      ref = Especes_find(molecules, Species_getId(especeId), nbEspeces);
      quantite= Especes_getQuantite(molecules, ref);
      /* La quantite est egale a 0 */
      if(quantite<=0.0) return END;
      temp = floor(Especes_getQuantite(molecules, ref)/SpeciesReference_getStoichiometry(reactif));
      /* Cas ou le nombre de pas est egal a 0 */
      if(temp==0.0) return END;
      /* Si le nouveau nombre est inferieur au precedent, on change la valeur de minStep */
      if (minStep > temp) minStep = temp;
  }

  return (int)minStep;
}

/**
 * \fn Reaction_t * SBML_reactChoice(pEspeces molecules, const gsl_rng * r, int ref)
 * \author Amine Ghozlane
 * \brief  Determine randomly the reaction to achieve for several nodes reactions
 * \param  molecules Struct Especes
 * \param  r Random number generator
 * \param  ref Number reference of one molecule
 * \return Id of the selected reaction
 */
Reaction_t * SBML_reactChoice(pEspeces molecules, const gsl_rng * r, int ref) {

  pReaction temp = NULL;
  pReaction Q = molecules[ref].system;
  double value = gsl_rng_uniform(r) * 100.0;
  double choice = 0.0;

  /*printf("Choix de la reaction :\n");
	 printf("value : %f\n", value);*/

  /* Choix de la reaction a realiser */
  if (value < Q->ratio)
    temp = Q;
  else {
      do {
          /*printf("reaction %s : ratio %f\n", Reaction_getId(Q->link), Q->ratio);*/
          choice += Q->ratio;
          /*printf("choice : %f\n", choice);*/
          Q = Q->suivant;
          temp = Q;
          /*printf("ratio suivant : %f\n",Q->ratio );*/
      } while (Q->suivant != NULL && value <= (choice + Q->suivant->ratio)
          && value > choice);
  }
  if (temp == NULL) {
      fprintf(stderr, "on a un probleme de ratio\n");
      exit(EXIT_FAILURE);
  }

  /*printf("Resultat :\n");
	 printf("value : %f, choice :%f\n", value, (choice + Q->ratio));
	 printf("reaction %s : ratio %f\n", Reaction_getId(temp->link), temp->ratio);*/
  /*if (Q->suivant == NULL)
	 printf("choice :%f", choice);*/
  /*printf("c bon\n");*/

  return (temp->link);

}

/**
 * \fn void SBML_reaction(Model_t *mod, pEspeces molecules, Reaction_t *react, int nbEspeces)
 * \author Amine Ghozlane
 * \brief  Simulation of a discrete transision
 * \param  mod Model of the SBML file
 * \param  molecules Struct Especes
 * \param  react Reaction id
 * \param  nbEspeces Number of molecules
 */
void SBML_reaction(Model_t *mod, pEspeces molecules, Reaction_t *react, int nbEspeces)
{
  SpeciesReference_t *reactif;
  Species_t *especeId;
  int i, ref = 0;

  /*boucle pour retirer des reactifs*/
  for (i = 0; i < (int)Reaction_getNumReactants(react); i++) {
      reactif = Reaction_getReactant(react, i);
      especeId = Model_getSpeciesById(mod, SpeciesReference_getSpecies(reactif));
      ref = Especes_find(molecules, Species_getId(especeId), nbEspeces);

      /*printf("Reactif : %s",Species_getId(especeId) );
		 printf("ref : %d\n",ref);
		 printf("Reactif :Quantite de ref %d :%d \n",ref,Espece_getQuantite(molecules,ref));*/
      Especes_setQuantite(molecules, ref,(Especes_getQuantite(molecules, ref)- SpeciesReference_getStoichiometry(reactif)));
      /*printf("Apres Reactif: Quantite de ref %d :%d \n",ref,Espece_getQuantite(molecules,ref));*/
  }

  /*boucle pour ajouter des produits */
  for (i = 0; i < (int)Reaction_getNumProducts(react); i++) {
      reactif = Reaction_getProduct(react, i);
      especeId = Model_getSpeciesById(mod, SpeciesReference_getSpecies(reactif));
      /*printf("Produit : %s",Species_getId(especeId) );*/
      ref = Especes_find(molecules, Species_getId(especeId), /*Species_getCompartment(especeId),*/ nbEspeces);
      /*printf("Produit : ref : %d\n",ref);
		 printf("Produit : Quantite de ref %d :%d \n",ref,Espece_getQuantite(molecules,ref));*/

      Especes_setQuantite(molecules, ref,(Especes_getQuantite(molecules, ref)+ SpeciesReference_getStoichiometry(reactif)));
      /*printf("Apres Produit : Quantite de ref %d :%d \n",ref,Espece_getQuantite(molecules,ref));*/
  }
}

/**
 * \fn void SBML_allocTest(pTestReaction T, int nbReactions)
 * \author Amine Ghozlane
 * \brief  Alloc memory and initialize the struct pTestReaction
 * \param  T Empty struct TestReaction
 * \param  nbReactions Number of reactions
 */
void SBML_allocTest(pTestReaction T, int nbReactions)
{
  int i;

  /* Initialisation des tableaux des reactions */
  T->tabReactions= (Reaction_t **) malloc(nbReactions * sizeof(Reaction_t *));
  T->minStepTab = (int*) malloc(nbReactions * sizeof(int));

  if (T->tabReactions == NULL)
    exit(EXIT_FAILURE);
  if (T->minStepTab == NULL)
    exit(EXIT_FAILURE);

  for (i = 0; i < nbReactions; i++) {
      T->tabReactions[i] = NULL;
      T->minStepTab[i] = 0;
  }
}

/**
 * \fn void SBML_freeTest(pTestReaction T)
 * \author Amine Ghozlane
 * \brief  Free memory of the struct TestReaction
 * \param T Struct TestReaction gives data on reaction
 */
void SBML_freeTest(pTestReaction T)
{
  if (T->tabReactions != NULL)
    free(T->tabReactions);
  if (T->minStepTab != NULL)
    free(T->minStepTab);
}

/**
 * \fn int SBML_EstimationReaction(Model_t *mod, pTestReaction T, pEspeces molecules, int ref, int nbEspeces)
 * \author Amine Ghozlane
 * \brief  Alloc memory and initialize the struct Especes
 * \param  mod Model of the SBML file
 * \param  T Struct TestReaction gives data on reaction
 * \param  molecules Struct Especes
 * \param  ref Number reference of one molecule
 * \param  nbEspeces Number of molecules
 * \return Estimated number of feasible step by reaction
 */
int SBML_EstimationReaction(Model_t *mod, pTestReaction T, pEspeces molecules, int ref, int nbEspeces)
{
  pReaction Q = molecules[ref].system;
  int i = 0, min = 0, curr = 0;

  /* Compte le nombre de reactions rattachees a une espece */
  while (Q != NULL) {
      T->tabReactions[i] = Q->link;
      T->minStepTab[i] = SBML_checkQuantite(mod, Q->link, nbEspeces,
          molecules);
      if (i == 0)
        min = T->minStepTab[i];
      curr = T->minStepTab[i];
      if (T->minStepTab[i] <= 0)
        return END;
      else if (curr < min)
        min = curr;
      Q = Q->suivant;
      i++;
  }

  return min;
}

/**
 * \fn int SBML_simulate(Model_t *mod, pEspeces molecules, const gsl_rng * r, pTestReaction T, char **banned, int nbBanned, int nbEspeces, int ref)
 * \author Amine Ghozlane
 * \brief  Simulate one step of petri net
 * \param  mod Model of the SBML file
 * \param  molecules Struct Especes
 * \param  r Random number generator
 * \param  T Struct TestReaction gives data on reaction
 * \param  banned List of banned compound
 * \param  nbBanned Number of banned compound
 * \param  nbEspeces Number of molecules
 * \param  ref Number reference of one molecule
 * \return Condition of stop/pursue
 */
int SBML_simulate(Model_t *mod, pEspeces molecules, const gsl_rng * r, pTestReaction T, int nbEspeces, int ref)
{
  Reaction_t *react = NULL;
  int minStep = 0, valid = 0/*, i=0*/;
  /*int nbReactions = Especes_getNbreactions(molecules, ref);*/
  int nbReactions = molecules[ref].nbReactions;

  /*printf("nbreactions : %d, ref : %d\n", nbReactions, ref);
	 printf("molecules : %s\n", molecules[ref].id);*/

  /* Probleme ATP, ADP, NADH, NAD+ */
  if (!strcmp(molecules[ref].id, "ADP") || !strcmp(molecules[ref].id, "ATP")
      || !strcmp(molecules[ref].id, "NADH") || !strcmp(molecules[ref].id,
          "NADplus") || !strcmp(molecules[ref].id, "NADPH") || !strcmp(
              molecules[ref].id, "NADPplus")|| !strcmp(molecules[ref].id, "NADplus_g")||
              !strcmp(molecules[ref].id, "NADH_g")|| !strcmp(molecules[ref].id, "ADP_g")||
              !strcmp(molecules[ref].id, "ATP_g")|| !strcmp(molecules[ref].id, "ADP_c")|| !strcmp(molecules[ref].id, "ATP_c")) {
      /*printf("END : molecules : %s\n", molecules[ref].id);*/
      return END;
  }

  /* Elimine les calculs sur ADP... */
  /*for(j=0;j<BANNI;j++){
	 if(!strcmp(tab[i],molecules[ref].id)){
	 printf("test : tab: %s molecules : %s\n",tab[i],molecules[ref].id );
	 return END;
	 }
	 }*/

  /* Variation selon le cas. */
  switch (nbReactions) {
  /* Cas ou aucune reaction n'est possible */
  case 0:
    /*printf("Cas N°1\n");*/
    return END;
    break;
    /* Cas ou une seule reaction est possible */
  case 1:
    /*printf("Cas N°2\n");*/
    react = molecules[ref].system->link;
    if ((minStep = SBML_checkQuantite(mod, react, nbEspeces, molecules))<= END)
      return END;
    /*printf("minstep = %d\n", minStep);*/

    while (minStep > 0) {
        SBML_reaction(mod, molecules, react, nbEspeces);
        minStep--;
    }
    break;
    /* Cas ou  plusieurs reactions sont possibles*/
  default:
    /*printf("Cas N°3\n");*/
    /* Allocation de la memoire au tableau des reactions */
    SBML_allocTest(T, nbReactions);
    valid = SBML_EstimationReaction(mod, T, molecules, ref, nbEspeces);
    if (valid <= END) {
        /*printf("on ne fait rien\n");*/
        return END;
    }
    while (Especes_getQuantite(molecules, ref) > 0.0 && valid > END) {
        react = SBML_reactChoice(molecules, r, ref);
        /*printf("molecule : %s, reaction %s\n", molecules[ref].id, Reaction_getId(react));*/
        /*printf("reaction : %s\n", Reaction_getId(react));*/
        /*printf("valid : %d, i : %d\n", valid, i);*/
        SBML_reaction(mod, molecules, react, nbEspeces);
        /*if(i>=valid) valid=SBML_EstimationReaction(mod, T, molecules, ref, nbEspeces);*/
        /*i++;*/
        valid--;
    }
    /*printf("Fin de simulation\n");*/

    /* Liberation de la memoire allouee au tableau des reactions */
    SBML_freeTest(T);
    break;
  }
  /*printf("On quitte la simulation\n");*/
  return PURSUE;
}

/**
 * \fn void SBML_compute_simulation(pScore result, Model_t *mod, double *reactions_ratio, gsl_rng * r, char **banned, int nbBanned)
 * \author Amine Ghozlane
 * \brief  Simulation of metabolic network
 * \param  result Struct Score
 * \param  mod Model of the SBML file
 * \param  reactions_ratio List of computed reaction ratio
 * \param  r Random number generator
 * \param  banned List of banned compound
 * \param  nbBanned Number of banned compound
 */
void compute_simulation(Model_t *mod)
{
  int i, nbReactions = 0, nbEspeces = 0, temp = 1, tempo = 0;
  pEspeces molecules;
  const gsl_rng_type *T;
  gsl_rng * r;
  pTestReaction TR = NULL;

  /* Verifications */
  /*TODO numGlobalParameters=Model_getNumParameters(mod);
	checkBoundaryConditions(mod);
	if(numGlobalParameters>0){
	    globalParVec=gsl_vector_alloc(numGlobalParameters);
	    setGlobalParVec();
	  }*/
  /* Allocation memoire et initialisation des generateurs de nombre aleatoire */
  gsl_rng_env_setup();
  if (!getenv("GSL_RNG_SEED"))
    gsl_rng_default_seed = time(NULL);
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);
  gsl_rng_default_seed = gsl_rng_uniform(r);

  /* Allocation memoire */
  TR = (pTestReaction) malloc(1 * sizeof(TestReaction));
  if (T == NULL)
    exit(EXIT_FAILURE);

  /* File information */
  nbReactions = Model_getNumReactions(mod);
  nbEspeces = Model_getNumSpecies(mod);
  /*printf("Nom du model : %s\n",Model_getId(mod));*/
  molecules = Especes_alloc(nbEspeces);
  /* Initialisation de la quantite des especes */
  SBML_initEspeceAmounts(mod, molecules, nbEspeces);
  /* Initialisation des reactions et des ratios*/
  SBML_setReactions(mod, molecules, nbReactions, nbEspeces);

  /*TODO checkRatio(mod,molecules);*/

  /* Test des donnees enregistrees */
  printf("\nEtat des especes :\n\n");
  Especes_print(molecules, nbEspeces);

  printf("\nDebut de simulation ...\n");
  /* Simulation des reactions */
  while (temp > END) {
      temp = 0;
      for (i = 0; i < nbEspeces; i++) {
          tempo = SBML_simulate(mod, molecules, r, TR, nbEspeces, i);
          /*if (tempo != END) {
              Especes_print_2(molecules, nbEspeces);
				 printf("\n");
          }*/
          temp += tempo;
      }
  }
  printf("Fin de simulation...\n");
  printf("\nEtat des especes :\n");
  Especes_print_2(molecules, nbEspeces);
  /* Liberation de la memoire des generateurs aleatoire */
  gsl_rng_free(r);

  /* Liberation de la memoire de la structure Especes */
  Especes_free(molecules, nbEspeces);
  if (TR != NULL)
    free(TR);
}

/**
 * \fn void SBML_score_add(pScore result, pScore result_temp, FILE *debugFile)
 * \author Amine Ghozlane
 * \brief  Add scores
 * \param  result Struct Score used for all the simulation
 * \param  result_temp Struct Score used at each simulation step
 * \param  debugFile File use for debug
 */
void SBML_score_add(pScore result, pScore result_temp, int nbEspeces)
{
  /* Addition des scores */
  int i;

  /* Copie des resultats */
  for(i=0;i<nbEspeces;i++){
      if(result->name[i]==NULL){
          result->name[i]=(char*)malloc(((int)strlen(result_temp->name[i])+1)*sizeof(char));
          assert(result->name[i]!=NULL);
          strcpy(result->name[i], result_temp->name[i]);
      }
      result->quantite[i]+=result_temp->quantite[i];
  }
}

/**
 * \fn void SBML_score_mean(pScore result, int n)
 * \author Amine Ghozlane
 * \brief  Mean quantities for score
 * \param  result Struct Score
 * \param  n Number of simulation step
 */
void SBML_score_mean(pScore result, int nbEspeces, int n)
{
  /* Moyenne des resultats */
  int i;
  for(i=0;i<nbEspeces;i++){
      result->quantite[i]/=n;
  }
}

/**
 * \fn pScore SBML_scoreAlloc(pListParameters a)
 * \author Amine Ghozlane
 * \brief  Allocation of the struct Score
 * \param  a Global parameters : struct ListParameters
 * \return Allocated struct Score
 */
pScore SBML_scoreAlloc(int nbEspeces)
{
  int i;
  /* Allocation de la structure de score */
  pScore result=NULL;
  result=(pScore)malloc(1*sizeof(Score));
  assert(result!=NULL);

  result->name=NULL;
  result->quantite=NULL;
  result->name=(char**)malloc(nbEspeces*sizeof(char*));
  assert(result->name!=NULL);
  result->quantite=(double*)malloc(nbEspeces*sizeof(double));
  assert(result->quantite!=NULL);
  for(i=0;i<nbEspeces;i++){
      result->name[i]=NULL;
      result->quantite[i]=0.0;
  }
  return result;
}

/**
 * \fn void SBML_scoreFree(pScore out)
 * \author Amine Ghozlane
 * \brief  Free the struct Score
 * \param  out Struct score
 */
void SBML_scoreFree(pScore out, int nbEspeces)
{
  int i;
  /* Libere la memoire alloue a la structure de score */
  for(i=0;i<nbEspeces;i++){
      if(out->name[i]!=NULL) free(out->name[i]);
  }
  if(out->name!=NULL) free(out->name);
  if(out->quantite!=NULL) free(out->quantite);
  if(out!=NULL) free(out);
}

/* Affiche le score moyen */

/**
 * \fn void SBML_scoreFree(pScore out)
 * \author Amine Ghozlane
 * \brief  Free the struct Score
 * \param  out Struct score
 * \param  nbEspeces
 */
void SBML_score_print(pScore result, int nbEspeces)
{
  int i;
  /* Print head*/
  /*for(i=0;i<nbEspeces;i++){
      printf("%s;",result->name[i]);
  }
  printf("\n");*/
  /* Print Value */
  for(i=0;i<nbEspeces;i++){
      /*printf("Molecule: %s - Biomass: %.0f\n",result->name[i],result->quantite[i]);*/
      printf("%.0f;",result->quantite[i]);
  }
  printf("\n");
}

/**
 * \fn void SBML_compute_simulation_mean(Model_t *mod, int nb_simulation)
 * \author Amine Ghozlane
 * \brief  X time simulation of metabolic network
 * \param  mod Model of the SBML file
 * \param  nb_simulation Number of simulation step
 */
void SBML_compute_simulation_mean(Model_t *mod, int nb_simulation)
{
  /* Simulation du reseau metabolique  */
  int i, j,nbReactions = 0, nbEspeces = 0, temp = 1, tempo=0, test=0;
  pEspeces molecules=NULL;
  pTestReaction TR=NULL;
  const gsl_rng_type *T;
  gsl_rng * r;
  pScore result=NULL,result_temp=NULL;

  /* Allocation memoire et initialisation des generateurs de nombre aleatoire */
  gsl_rng_env_setup();
  if (!getenv("GSL_RNG_SEED"))
    gsl_rng_default_seed = time(NULL)+(double)getpid();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);
  /*gsl_rng_default_seed = time(NULL)+(double)getpid();*/

  /* Allocation memoire */
  TR=(pTestReaction)malloc(1*sizeof(TestReaction));
  assert(TR!=NULL);

  /* File information */
  nbReactions = (int)Model_getNumReactions(mod);
  nbEspeces = (int)Model_getNumSpecies(mod);
  molecules = Especes_alloc(nbEspeces);

  /* Initialisation de la structure de score temporaire */
  result=SBML_scoreAlloc(nbEspeces);
  result_temp=SBML_scoreAlloc(nbEspeces);

  /* Initialisation de la quantite des especes */
  SBML_initEspeceAmounts(mod, molecules, nbEspeces);
  /* Initialisation des reactions et des ratios*/
  SBML_setReactions(mod, molecules, nbReactions, nbEspeces);

  /*printf("Start the simulation\n");*/
  /* SIMULATION */
  for(j=0;j<nb_simulation;j++){
      /* Simulation des reactions */
      while (temp > END) {
          temp = 0;
          for (i = 0; i < nbEspeces; i++) {
              tempo = SBML_simulate(mod, molecules, r, TR, nbEspeces, i);
              temp +=tempo;
          }
          test+=1;
      }
      temp=1;
      tempo=0;
      /*printf("Nombre de tour %d\n",test);*/
      /*Score */
      /* Enregistre le score des especes */
      Especes_scoreSpecies(molecules, nbEspeces, result_temp->name, result_temp->quantite);
      SBML_score_add(result,result_temp,nbEspeces);
      SBML_initEspeceAmounts(mod, molecules, nbEspeces);
  }
  SBML_score_mean(result,nbEspeces,nb_simulation);

  /*printf("End of the simulation...\n");
  printf("Final state of the molecules :\n");*/
  SBML_score_print(result, nbEspeces);

  /* Liberation de la memoire des generateurs aleatoire */
  gsl_rng_free(r);

  /* Liberation de la memoire du score */
  SBML_scoreFree(result, nbEspeces);
  SBML_scoreFree(result_temp, nbEspeces);

  /* Liberation de la memoire de la structure Especes */
  Especes_free(molecules, nbEspeces);
  if(TR!=NULL) free(TR);
}

max_stars_repo_name: aghozlane/yasimSBML

max_issues_repo_name: aghozlane/yasimSBML

max_forks_repo_name: aghozlane/yasimSBML


text:
/**
 * @file lpc.h
 * @brief Speech recognition front end.
 * @author John McDonough, Matthias Woelfel, Kenichi Kumatani
 */

#include <gsl/gsl_vector.h>
#include "feature/feature.h"

#ifndef LPC_H
#define LPC_H

// ----- definition for class `BaseFeature' -----
//
class BaseFeature {
protected:
  BaseFeature(unsigned order, unsigned dim);
  virtual ~BaseFeature();

  void fftPower(float* power);
  //virtual void autoCorrelation(const float* X, float* LP, float* E, float warp ) = 0;

  unsigned					_order;
  unsigned					_dim;

private:
  unsigned					_log2Length;

  unsigned 					_npoints;
  unsigned					_npoints2;
  double* 					_temp;
};

// ----- definition for class `WarpFeature' -----
//
class WarpFeature : protected BaseFeature {
protected:
  WarpFeature(unsigned order, unsigned dim);
  virtual ~WarpFeature();

  void autoCorrelation(const float* X, float* LP, float* E, float warp );

private:
  float*            				_K;
  float*            				_R;
  float*            				_WX;
  float*            				_WXTEMP;
  gsl_matrix_float* 				_tmpA;
};

// ----- definition for class `BurgFeature' -----
//
class BurgFeature : protected BaseFeature {
protected:
  BurgFeature(unsigned order, unsigned dim);
  virtual ~BurgFeature();

  void autoCorrelation(const float* X, float* LP, float* E, float warp );

private:
  float*					_EF;
  float*					_EB;
  float*					_EFP;
  float*					_EBP;
  float*					_A_flip;
  float*					_K;
};

// ----- definition for class `MVDRFeature' -----
//
template <class AutoCovariance>
class MVDRFeature : public VectorFeatureStream, private AutoCovariance {
 public:
  MVDRFeature(const VectorFloatFeatureStreamPtr& src, unsigned order = 60, unsigned correlate = 0, float warp = 0.0, const String& nm = "MVDR");
  virtual ~MVDRFeature();

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { _src->reset(); VectorFeatureStream::reset(); }

private:
  VectorFloatFeatureStreamPtr			_src;
  unsigned					_temp_order;
  unsigned					_correlate;
  float						_warp;
  float*					_tmpR;
  float*					_tmpA;
  float*					_tmpPC;
  float*					_tmpPA;
  float*					_E;
};

// ----- methods for class `MVDRFeature' -----
//
template <class AutoCovariance>
MVDRFeature<AutoCovariance>::
MVDRFeature(const VectorFloatFeatureStreamPtr& src, unsigned order, unsigned correlate, float warp, const String& nm)
  : VectorFeatureStream(src->size()/2+1, nm), AutoCovariance(order, src->size()),
    _src(src), _temp_order(2*order+1),
    _correlate(correlate), _warp(warp),
    _tmpR(new float[order+1]), _tmpA(new float[order+1]),
    _tmpPC(new float[_temp_order]), _tmpPA(new float[AutoCovariance::_dim+1]),
    _E(new float[order+1])
{
  if (_correlate < 10) _correlate = AutoCovariance::_dim;
  if (AutoCovariance::_order >= AutoCovariance::_dim/2+1)
    throw jparameter_error("Order (%d) and dimension (%d) do not match.", AutoCovariance::_order, AutoCovariance::_dim/2+1);
}

template <class AutoCovariance>
MVDRFeature<AutoCovariance>::~MVDRFeature()
{
  delete[] _tmpR;  delete[] _tmpA;  delete[] _tmpPC;  delete[] _tmpPA;  delete[] _E;
}

template <class AutoCovariance>
const gsl_vector* MVDRFeature<AutoCovariance>::MVDRFeature::next(int frame_no)
{
  if (frame_no == frame_no_) return vector_;

  if (frame_no >= 0 && frame_no - 1 != frame_no_)
    throw jindex_error("Problem in Feature %s: %d != %d\n",
		       name().c_str(), frame_no - 1, frame_no_);

  const gsl_vector_float* block = _src->next(frame_no_ + 1);
  increment_();

  const float* X  = block->data;
  float* A  = _tmpA;
  float* PC = _tmpPC;
  float* PA = _tmpPA;

  AutoCovariance::autoCorrelation(X, _tmpA, _E, _warp);

  for (int i = 0; i <= AutoCovariance::_order; i++) {
    double temp = 0;
    for (int ii = 0; ii <= (AutoCovariance::_order-i); ii++)
      temp += (float)(int(AutoCovariance::_order)+1-i-2*ii)*A[ii]*A[ii+i];
    if (_E[0] > 0) {
      PC[AutoCovariance::_order+i] = -temp;
      // instead of -temp/E[order] in the paper by B. Musicus IEEE 1985
    } else {
      PC[AutoCovariance::_order+i] = 10000000;
    }
  }

  for (int i = 1; i <= AutoCovariance::_order; i++) PC[AutoCovariance::_order-i] = PC[AutoCovariance::_order+i];

  PA[0] = 0;
  for (unsigned i = 1; i <= _temp_order; i++)
    PA[i] = PC[i-1];  // [..-1] because of fft

  for (unsigned i = _temp_order+1; i <= AutoCovariance::_dim; i++)
    PA[i] = 0;

  /* ------------------------------------------------------------------------
     Calculate MVDR Power
     (also called Capon Max. Likelihood Method)
  */

  AutoCovariance::fftPower(PA);

  for (unsigned i = 0; i<= AutoCovariance::_dim/2; i++) {
    double temp = sqrt(PA[i]);
    if (temp > 0) {
      temp = _E[0]/temp;
    } else {
      temp = 10000000;
    }

    gsl_vector_set(vector_, i, temp);
  }

  return vector_;
}

typedef MVDRFeature<WarpFeature> WarpMVDRFeature;
typedef Inherit<WarpMVDRFeature, VectorFeatureStreamPtr> WarpMVDRFeaturePtr;

typedef MVDRFeature<BurgFeature> BurgMVDRFeature;
typedef Inherit<BurgMVDRFeature, VectorFeatureStreamPtr> BurgMVDRFeaturePtr;

// ----- definition for class `WarpedTwiceFeature' -----
//
class WarpedTwiceFeature : protected BaseFeature {
protected:
  WarpedTwiceFeature(unsigned order, unsigned dim );
  virtual ~WarpedTwiceFeature();

  void autoCorrelation(const float* X, float* LP, float* E, float warp, float rewarp=0.0 );

private:
  float*     _K;
  float*     _R;
  float*     _WX;
  float*     _WXTEMP;
  float**    _tmpA;
};

// ----- definition for class `WarpedTwiceMVDRFeature' -----
//
class WarpedTwiceMVDRFeature : public VectorFeatureStream, private WarpedTwiceFeature {
public:
  WarpedTwiceMVDRFeature(const VectorFloatFeatureStreamPtr& src, unsigned order = 60, unsigned correlate = 0, float warp = 0.0, bool warpFactorFixed=false, float sensibility = 0.1, const String& nm = "WTMVDR");
  virtual ~WarpedTwiceMVDRFeature();

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { _src->reset(); VectorFeatureStream::reset(); }

private:
  VectorFloatFeatureStreamPtr			_src;
  unsigned					_correlate;
  unsigned					_temp_order;
  float						_warp;
  float*					_tmpR;
  float*					_tmpA;
  float*					_tmpPC;
  float*					_tmpPA;
  float*					_E;
  float                                         _rewarp;
  bool                                          _warpFactorFixed;
  float                                         _sensibility;
};

typedef Inherit<WarpedTwiceMVDRFeature,VectorFeatureStreamPtr> WarpedTwiceMVDRFeaturePtr;

// ----- definition for class `LPCFeature' -----
//
template <class AutoCovariance>
class LPCFeature : public VectorFeatureStream, private AutoCovariance {
 public:
  LPCFeature(const VectorFloatFeatureStreamPtr& src, unsigned order, unsigned correlate, float warp, const String& nm = "LPC");
  virtual ~LPCFeature();

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { _src->reset(); VectorFeatureStream::reset(); }

 private:
  VectorFloatFeatureStreamPtr			_src;
  unsigned					_correlate;
  float						_warp;
  float*					_tmpA;
  float*					_tmpPA;
  float*					_E;
  const String					_lpmethod;
};

// ----- methods for class `LPCFeature' -----
//
template <class AutoCovariance>
LPCFeature<AutoCovariance>::
LPCFeature(const VectorFloatFeatureStreamPtr& src, unsigned order, unsigned correlate, float warp, const String& nm)
  : VectorFeatureStream(src->size()/2+1, nm),
    AutoCovariance(order, src->size()),
    _src(src),
    _correlate(correlate),
    _warp(warp),
    _tmpA(new float[order+1]),
    _tmpPA(new float[AutoCovariance::_dim+1]),
    _E(new float[order+1])
{
  if (_correlate < 10) _correlate = AutoCovariance::_dim;
  if (AutoCovariance::_order >= AutoCovariance::_dim/2+1)
    throw jparameter_error("Order (%d) and dimension (%d) do not match.", AutoCovariance::_order, AutoCovariance::_dim/2+1);
}

template <class AutoCovariance>
LPCFeature<AutoCovariance>::~LPCFeature()
{
  delete[] _tmpA;  delete[] _tmpPA;  delete[] _E;
}

template <class AutoCovariance>
const gsl_vector* LPCFeature<AutoCovariance>::next(int frame_no)
{
  if (frame_no == frame_no_) return vector_;

  if (frame_no >= 0 && frame_no - 1 != frame_no_)
    throw jindex_error("Problem in Feature %s: %d != %d\n",
		       name().c_str(), frame_no - 1, frame_no_);

  const gsl_vector_float* block = _src->next(frame_no_ + 1);
  increment_();

  const float* X = block->data;
  float* A  = _tmpA;
  float* PA = _tmpPA;

  AutoCovariance::autoCorrelation(X, _tmpA, _E, _warp);

  PA[0] = 0;
  for (unsigned i = 1; i <= AutoCovariance::_order+1; i++)
    PA[i] = A[i-1];  // [..-1] because of fft

  for (unsigned i = AutoCovariance::_order+2; i <= AutoCovariance::_dim; i++)
    PA[i] = 0;

  AutoCovariance::fftPower(PA);

  for (unsigned i = 0; i <= AutoCovariance::_dim/2; i++) {
    double temp = PA[i];
    if (temp > 0) {
      temp = (2*_E[0])/(temp*AutoCovariance::_dim);
    } else {
      temp = 10000000;
    }
    gsl_vector_set(vector_, i, temp);
  }

  return vector_;
}

typedef LPCFeature<WarpFeature> WarpLPCFeature;
typedef Inherit<WarpLPCFeature, VectorFeatureStreamPtr> WarpLPCFeaturePtr;

typedef LPCFeature<BurgFeature> BurgLPCFeature;
typedef Inherit<BurgLPCFeature, VectorFeatureStreamPtr> BurgLPCFeaturePtr;

// ----- definition for class `SpectralSmoothing' -----
//
class SpectralSmoothing : public VectorFeatureStream {
 public:
  SpectralSmoothing(const VectorFeatureStreamPtr& adjustTo, const VectorFeatureStreamPtr& adjustFrom,
		    const String& nm = "Spectral Smoothing");
  virtual ~SpectralSmoothing();

  virtual const gsl_vector* next(int frame_no = -5);

  virtual void reset() { _adjustTo->reset(); _adjustFrom->reset(); VectorFeatureStream::reset(); }

 private:
  VectorFeatureStreamPtr			_adjustTo;
  VectorFeatureStreamPtr			_adjustFrom;
  float*					_R;
};

typedef Inherit<SpectralSmoothing, VectorFeatureStreamPtr> SpectralSmoothingPtr;

#endif

ext: h


text:
#include <gsl/gsl_rng.h>
#include "gsl-sprng.h"
#include <gsl/gsl_randist.h>
#include <mpi.h>

int main(int argc,char *argv[])
{
  int k,i,iters; double x,can,a,alpha; gsl_rng *r;
  FILE *s; char filename[15];
  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&k);
  if ((argc != 3)) {
    if (k == 0)
      fprintf(stderr,"Usage: %s <iters> <alpha>\n",argv[0]);
    MPI_Finalize(); return(EXIT_FAILURE);
  }
  iters=atoi(argv[1]); alpha=atof(argv[2]);
  r=gsl_rng_alloc(gsl_rng_sprng20);
  sprintf(filename,"chain-%03d.tab",k);
  s=fopen(filename,"w");
  if (s==NULL) {
    perror("Failed open");
    MPI_Finalize(); return(EXIT_FAILURE);
  }
  x = gsl_ran_flat(r,-20,20);
  fprintf(s,"Iter X\n");
  for (i=0;i<iters;i++) {
    can = x + gsl_ran_flat(r,-alpha,alpha);
    a = gsl_ran_ugaussian_pdf(can) / gsl_ran_ugaussian_pdf(x);
    if (gsl_rng_uniform(r) < a)
      x = can;
    fprintf(s,"%d %f\n",i,x);
  }
  fclose(s);
  MPI_Finalize(); return(EXIT_SUCCESS);
}

ext: c


text:
#include <pthread.h>
#include <math.h>
#include <fftw.h>
#include <stdio.h>

//C -> Very ugly code

//DEFINITION OF COMPLEX NUMBERS
typedef struct {
  double re, im;
} mycom;

//DEFINITION OF DATA ARRAY NEEDED FOR THE TRANSFORMATION (small space)
typedef struct {
  double* m;
  double* vi;
  double* vo;
  int*    nrop;
} transdata;

//DEFINITION OF DATA ARRAY NEEDED FOR THE TRANSFORMATION (large space)
typedef struct {
  double* m;
  double* vi;
  double* vo;
  long long int*    nrop;
} ltransdata;

//THREAD FUNCTION FOR FORWARD TRANSFORMATION (small space)
void transforeT(void* td){
  transdata* Td = (transdata *)td;
  double* Mat       = Td->m;
  double* Vec_in    = Td->vi;
  double* Vec_out   = Td->vo;
  int*       Nrop   = Td->nrop;
  int x, y;

  for(x = 0; x < 2* *Nrop; x++){
    Vec_out[x] = 0.;
  }

  for(x = 0; x < *Nrop; x++){
    for(y = 0; y < *Nrop; y++){
      Vec_out[2*x]    +=     Mat[(x* (*Nrop))+y]*Vec_in[2*y];
      Vec_out[2*x+1]  +=     Mat[(x* (*Nrop))+y]*Vec_in[2*y+1];
    }
  }

}

//THREAD FUNCTION FOR FORWARD TRANSFORMATION (large space)
void ltransforeT(void* td){
  ltransdata* Td = (ltransdata *)td;
  double* Mat       = Td->m;
  double* Vec_in    = Td->vi;
  double* Vec_out   = Td->vo;
  long long int*       Nrop   = Td->nrop;
  long long int x, y;

  for(x = 0; x < 2* *Nrop; x++){
    Vec_out[x] = 0.;
  }

  for(x = 0; x < *Nrop; x++){
    for(y = 0; y < *Nrop; y++){
      Vec_out[2*x]    +=     Mat[(x* (*Nrop))+y]*Vec_in[2*y];
      Vec_out[2*x+1]  +=     Mat[(x* (*Nrop))+y]*Vec_in[2*y+1];
    }
  }

}

//THREAD FUNCTION FOR BACKWARD TRANSFORMATION (small space)
void transbackT(void* td){
  transdata* Td = (transdata *)td;
  double* Mat       = Td->m;
  double* Vec_in    = Td->vi;
  double* Vec_out   = Td->vo;
  int*       Nrop   = Td->nrop;
  int x, y;

  for(x = 0; x < 2* *Nrop; x++){
    Vec_out[x] = 0.;
  }

  for(x = 0; x < *Nrop; x++){
    for(y = 0; y < *Nrop; y++){
      Vec_out[2*y]    +=     Mat[x* *Nrop+y]*Vec_in[2*x];
      Vec_out[2*y+1]  +=     Mat[x* *Nrop+y]*Vec_in[2*x+1];
    }
  }

}

//THREAD FUNCTION FOR BACKWARD TRANSFORMATION (large space)
void ltransbackT(void* td){
  ltransdata* Td = (ltransdata *)td;
  double* Mat       = Td->m;
  double* Vec_in    = Td->vi;
  double* Vec_out   = Td->vo;
  long long int*       Nrop   = Td->nrop;
  long long int x, y;

  for(x = 0; x < 2* *Nrop; x++){
    Vec_out[x] = 0.;
  }

  for(x = 0; x < *Nrop; x++){
    for(y = 0; y < *Nrop; y++){
      Vec_out[2*y]    +=     Mat[x* *Nrop+y]*Vec_in[2*x];
      Vec_out[2*y+1]  +=     Mat[x* *Nrop+y]*Vec_in[2*x+1];
    }
  }

}

//MAKE TRANSFORMATION MULTI THREADED (small space) DIR=1 is forward !!
void mk_thread_trans(int dir, double* mat1, double* mat2, mycom* vec_in1, mycom* vec_in2,
		     mycom* vec_out1, mycom* vec_out2, int* nrop1, int* nrop2){
  pthread_t thread1, thread2;
  int  iret1, iret2;
  transdata td1, td2;

  td1.m    = mat1;
  td1.vi   = (double* ) vec_in1;
  td1.vo   = (double* ) vec_out1;
  td1.nrop = nrop1;

  td2.m    = mat2;
  td2.vi   = (double* ) vec_in2;
  td2.vo   = (double* ) vec_out2;
  td2.nrop = nrop2;

  if(dir == 1){
    iret1 = pthread_create( &thread1, NULL,  (void* ) &transforeT , (void*) &td1);
    iret2 = pthread_create( &thread2, NULL,  (void* ) &transforeT , (void*) &td2);
  }else{
    iret1 = pthread_create( &thread1, NULL,  (void* ) &transbackT , (void*) &td1);
    iret2 = pthread_create( &thread2, NULL,  (void* ) &transbackT , (void*) &td2);
  }

  pthread_join( thread1, NULL);
  pthread_join( thread2, NULL);

}

//MAKE TRANSFORMATION MULTI THREADED (large space) DIR=1 is forward !!
void lmk_thread_trans(int dir, double* mat1, double* mat2, mycom* vec_in1, mycom* vec_in2,
		      mycom* vec_out1, mycom* vec_out2, long long int* nrop1, long long int* nrop2){
  pthread_t thread1, thread2;
  int  iret1, iret2;
  ltransdata td1, td2;

  td1.m    = mat1;
  td1.vi   = (double* ) vec_in1;
  td1.vo   = (double* ) vec_out1;
  td1.nrop = nrop1;

  td2.m    = mat2;
  td2.vi   = (double* ) vec_in2;
  td2.vo   = (double* ) vec_out2;
  td2.nrop = nrop2;

  if(dir == 1){
    iret1 = pthread_create( &thread1, NULL,  (void* ) &ltransforeT , (void*) &td1);
    iret2 = pthread_create( &thread2, NULL,  (void* ) &ltransforeT , (void*) &td2);
  }else{
    iret1 = pthread_create( &thread1, NULL,  (void* ) &ltransbackT , (void*) &td1);
    iret2 = pthread_create( &thread2, NULL,  (void* ) &ltransbackT , (void*) &td2);
  }

  pthread_join( thread1, NULL);
  pthread_join( thread2, NULL);

}

max_stars_repo_name: sklinkusch/scripts

max_issues_repo_name: sklinkusch/scripts

max_forks_repo_name: sklinkusch/scripts


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

/*
 *
 *
 *       Filename:  cp_gen.c
 *
 *    Description:  Bloch matrix for the completely positive generator
 *
 *        Version:  1.0
 *        Created:  28/04/2014 23:09:21
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <math.h>
#include "funcs.h"

/*
 *      FUNCTION
 *         Name:  cp_mat
 *  Description:  Creating the CP generator matrix in Bloch form, taking as
 *  		  arguments all the integrals, the physical parameters and
 *  		  the pointer to the matrix
 *
 */
int cp_mat ( gsl_matrix* cp_mx , void* params )
{
	double integrals[12] ;
	if ( (integration(params,integrals)) != 0 )
		return -1;

	/* Setting the integrals */
	double rcc, rss, rsc, isc, rcs, ics, rc0 ;
	rcc = integrals[0] ;
	rss = integrals[2] ;
	rsc = integrals[4] ;
	isc = integrals[5] ;
	rcs = integrals[6] ;
	ics = integrals[7] ;
	rc0 = integrals[8] ;

	/* Copying the parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p ( pars, &o_c, &b, &O, &o_1 ) ;
	double D = sqrt(POW_2(o_1)-POW_2(O)) ;

	/* Ensuring that the matrix is zeroed */
	gsl_matrix_set_zero ( cp_mx ) ;

	/*  Building the matrix */
	unsigned int i ;
	for ( i = 0; i < 4; i++ )
		gsl_matrix_set ( cp_mx, 0 , i , 0 ) ;
	gsl_matrix_set ( cp_mx, 1, 0, 0 ) ;
	gsl_matrix_set ( cp_mx, 1, 1, 2*(O/o_1)*rss+(1+POW_2(O/o_1))*rcc+2*POW_2(D/o_1)*rc0 ) ;
	gsl_matrix_set ( cp_mx, 1, 2, o_1/4-(1+POW_2(O/o_1))*rcs+2*(O/o_1)*rsc ) ;
	gsl_matrix_set ( cp_mx, 1, 3, 0	) ;
	gsl_matrix_set ( cp_mx, 2, 0, 0	) ;
	gsl_matrix_set ( cp_mx, 2, 1, -(o_1/4-(1+POW_2(O/o_1))*rcs+2*(O/o_1)*rsc) ) ;
	gsl_matrix_set ( cp_mx, 2, 2, 2*(O/o_1)*rss+(1+POW_2(O/o_1))*rcc+2*POW_2(D/o_1)*rc0 ) ;
	gsl_matrix_set ( cp_mx, 2, 3, 0	) ;
	gsl_matrix_set ( cp_mx, 3, 0, (1+POW_2(O/o_1))*ics-2*(O/o_1)*isc ) ;
	gsl_matrix_set ( cp_mx, 3, 1, 0	) ;
	gsl_matrix_set ( cp_mx, 3, 2, 0	) ;
	gsl_matrix_set ( cp_mx, 3, 3, (1+POW_2(O/o_1))*rcc+2*(O/o_1)*rss )  ;

	/* Multiplies times -4 to obtain the Bloch matrix */
	gsl_matrix_scale ( cp_mx, -4 ) ;

	return 0;
}		/* -----  end of function cp_mat  ----- */


text:
/*
 * knnring_sequential.c
 *
 *  Created on: Nov 21, 2019
 *      Author: Lambis
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include "knnring.h"

knnresult kNN(double * X, double * Y, int n, int m, int d, int k) {

	knnresult knn;

    // Allocating memory for the knnresult.
	knn.nidx = (int *)malloc(m*k*sizeof(int));
	knn.ndist = (double *)malloc(m*k*sizeof(double));

	// Passing m and k values to knnresult.
	knn.m = m;
	knn.k = k;

    // Allocating memory for the distances array.
	double *D = (double *)malloc(m*n*sizeof(double));

	// Calculation of sum(X.^2,2).
	double *a = (double *)malloc(n*sizeof(double));
	for (int i=0; i<n; i++) {
		a[i] = cblas_dnrm2(d, &X[i*d], 1);
		a[i] = a[i]*a[i];
	}

    // Calculation of sum(Y.^2,2).
	double *b = (double *)malloc(m*sizeof(double));
	for (int i=0; i<m; i++) {
		b[i] = cblas_dnrm2(d, &Y[i*d], 1);
		b[i] = b[i]*b[i];
	}

    // Calculation of -2*X*Y.' multiplication.
    double *c = (double *)malloc(n*m*sizeof(double));
	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n, m, d, -2, X, d, Y, d, 0, c, m);

	// Adding sum(X.^2,2) sum(Y.^2,2).' to the c array.
    for (int i=0; i<n; i++)
		for(int j=0; j<m; j++)
			c[i*m+j] += a[i] + b[j];

    // Cleanup.
    free(a);
    free(b);

    // Applying elementwise square root to the c array.
	for (int i=0; i<n; i++)
		for(int j=0; j<m; j++)
			c[i*m+j] = sqrt(fabs(c[i*m+j]));

    // Creates a nxn identity array.
    double *iA = (double *)calloc(n*n, sizeof(double));
    for (int i=0; i<n; i++)
        iA[i*n+i] = 1;

    // Transposing c array and storing it in D.
	cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, n, 1, c, m, iA, n, 0, D, n);

    // More efficient than the multiplication above. This can only be used with openblas.
    //cblas_domatcopy(CblasRowMajor, CblasTrans, n, m, 1, c, m, D, n);

    // Cleanup.
    free(iA);
    free(c);

    // Setup of idx number for corpus points.
	int *idx = (int *)malloc(m*n*sizeof(int));
	for (int i=0; i<m; i++)
        for (int j=0; j<n; j++)
            idx[i*n+j]= j;

    // Sorting of all distances.
	for (int i=0; i<m; i++)
		quickSort(D, i*n, ((i+1)*n)-1, idx);

    // Passing the k nearest neighbors' indexes and distances to knnresult.
	for (int i=0; i<m; i++) {
		for (int j=0; j<k; j++) {
			knn.nidx[i*k+j] = idx[i*n+j];
			knn.ndist[i*k+j] = D[i*n+j];
		}
	}

	return knn;

}

// A utility function to swap two elements.
void swap(void* a, void* b, size_t s) {
    void* tmp = malloc(s);
    memcpy(tmp, a, s);
    memcpy(a, b, s);
    memcpy(b, tmp, s);
    free(tmp);
}

/* This function takes last element as pivot, places the pivot element at its correct position in sorted array
 * and places all smaller (smaller than pivot) to left of pivot and all greater elements to right of pivot */
int partition(double *arr, int low, int high, int *idx) {
    double pivot = arr[high];    		// Pivot
    int i = (low - 1);  			// Index of smaller element

    for (int j=low; j<=high-1; j++) {
        // If current element is smaller than the pivot
        if (arr[j] < pivot) {
            i++;    				// Increment index of smaller element
            swap(&arr[i], &arr[j], sizeof(double));
            swap(&idx[i], &idx[j], sizeof(int));
        }
    }
    swap(&arr[i + 1], &arr[high], sizeof(double));
    swap(&idx[i + 1], &idx[high], sizeof(int));
    return (i + 1);
}

// The main function that implements QuickSort arr[] --> Array to be sorted, low  --> Starting index, high  --> Ending index
void quickSort(double *arr, int low, int high, int *idx) {
    if (low < high) {
        // pi is partitioning index, arr[pi] is now at right place
        int pi = partition(arr, low, high, idx);

        // Separately sort elements before partition and after partition
        quickSort(arr, low, pi - 1, idx);
        quickSort(arr, pi + 1, high, idx);
    }
}


text:
/* specfunc/elljac.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 Foundation; 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
 * MERCHANTABILITY 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_elljac.h>

/* GJ: See [Thompson, Atlas for Computing Mathematical Functions] */

/* BJG 2005-07: New algorithm based on Algorithm 5 from Numerische
   Mathematik 7, 78-90 (1965) "Numerical Calculation of Elliptic
   Integrals and Elliptic Functions" R. Bulirsch.

   Minor tweak is to avoid division by zero when sin(x u_l) = 0 by
   computing reflected values sn(K-u) cn(K-u) dn(K-u) and using
   transformation from Abramowitz & Stegun table 16.8 column "K-u"*/

int
gsl_sf_elljac_e(double u, double m, double * sn, double * cn, double * dn)
{
  if(fabs(m) > 1.0) {
    *sn = 0.0;
    *cn = 0.0;
    *dn = 0.0;
    GSL_ERROR ("|m| > 1.0", GSL_EDOM);
  }
  else if(fabs(m) < 2.0*GSL_DBL_EPSILON) {
    *sn = sin(u);
    *cn = cos(u);
    *dn = 1.0;
    return GSL_SUCCESS;
  }
  else if(fabs(m - 1.0) < 2.0*GSL_DBL_EPSILON) {
    *sn = tanh(u);
    *cn = 1.0/cosh(u);
    *dn = *cn;
    return GSL_SUCCESS;
  }
  else {
    int status = GSL_SUCCESS;
    const int N = 16;
    double mu[16];
    double nu[16];
    double c[16];
    double d[16];
    double sin_umu, cos_umu, t, r;
    int n = 0;

    mu[0] = 1.0;
    nu[0] = sqrt(1.0 - m);

    while( fabs(mu[n] - nu[n]) > 4.0 * GSL_DBL_EPSILON * fabs(mu[n]+nu[n])) {
      mu[n+1] = 0.5 * (mu[n] + nu[n]);
      nu[n+1] = sqrt(mu[n] * nu[n]);
      ++n;
      if(n >= N - 1) {
        status = GSL_EMAXITER;
        break;
      }
    }

    sin_umu = sin(u * mu[n]);
    cos_umu = cos(u * mu[n]);

    /* Since sin(u*mu(n)) can be zero we switch to computing sn(K-u),
       cn(K-u), dn(K-u) when |sin| < |cos| */

    if (fabs(sin_umu) < fabs(cos_umu))
      {
        t = sin_umu / cos_umu;

        c[n] = mu[n] * t;
        d[n] = 1.0;

        while(n > 0) {
          n--;
          c[n] = d[n+1] * c[n+1];
          r = (c[n+1] * c[n+1]) / mu[n+1];
          d[n] = (r + nu[n]) / (r + mu[n]);
          }

        *dn = sqrt(1.0-m) / d[n];
        *cn = (*dn) * GSL_SIGN(cos_umu) / gsl_hypot(1.0, c[n]);
        *sn = (*cn) * c[n] /sqrt(1.0-m);
      }
    else
      {
        t = cos_umu / sin_umu;

        c[n] = mu[n] * t;
        d[n] = 1.0;

        while(n > 0) {
          --n;
          c[n] = d[n+1] * c[n+1];
          r = (c[n+1] * c[n+1]) / mu[n+1];
          d[n] = (r + nu[n]) / (r + mu[n]);
        }

        *dn = d[n];
        *sn = GSL_SIGN(sin_umu) / gsl_hypot(1.0, c[n]);
        *cn = c[n] * (*sn);
      }

    return status;
  }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef L_Katyusha_H
#define L_Katyusha_H

#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <sstream>
#include <math.h>

//This class implements the loopless variance reduced type methods with arbitrary sampling

/*
The optimization problem to solve is:

F(x):=\sum_{i=1}^n 1/n* phi_i(x)+ g(x)
*/

template<typename L, typename D>
class L_Katyusha
{

protected:

  // involved variables

  std::vector<D> proba_vector;

  std::vector<D>  tilde_proba_vector;

  std::vector<D> group_C;

  std::vector<D> index_group_C;

  std::vector<D> maxp_group_C;

  std::vector<D> isolated_I;

  std::vector<D> sump_group_C;
  std::vector<D> theta_S;   //theta_S in the paper

  std::vector<D> Li;  // the Lipchitz constant of phi_i

  std::vector<D> gk;     // the vector g^k in the paper

  std::vector<D> gradient_f_w;  // gradient of f at w
  std::vector<D> gradient_f_x;  // gradient of f at x
  std::vector<D> batch_delta_gradient; //  1/n(G(x^k)-G(w^k))theta_{S_k} I_{S_k} e in the paper

  L c;

  L noverc;

  // auxiliary variables

  L nb_of_iters_per_loop;

  L max_nb_loops;

  D max_p;

  L nsamples;
  L nfeatures;

  D running_time;

  L nb_loops;

  L print_every_N;

  vector<L> batch_i;
  vector<L> my_batch;

  L batch_size;

  L nb_groups;    //number of groups in group sampling

  string uniform;

  D Lf;

  D L2;

  D sumLi;

  D upper_bound;   // upper bound of F(x)-F^*

  D p; // the probability of changing w to x as in the paper
  L tau;  //number of threads on each node/computer

  D mu;

  D scaler;

  L current_nb_iters;

  L nb_iters;

  ofstream samp;

  D theta1;
  D theta2;
  D theta3;
  D eta;  //eta/L in the paper
  D oneovereta;

  D primal_value;

  std::vector<D> x;
  std::vector<D> w;
  std::vector<D> y;
  std::vector<D> z;
  std::vector<D> next_x;

public:

  gsl_rng * rng;
  virtual inline D value_of_phi_i(L) {return D(NULL);}
  virtual inline D value_of_g(){return D(NULL);}
  virtual inline void prox_of_g(D, vector<D> &, vector<D> &, vector<D> &){}   //prox_of_g(L,x, gr, y) computes y=argmin_u{L/2 ||u-(x-gr/L)||^2+g(u)}
  virtual inline void set_auxiliary_v(){}
  virtual inline void update_gradient(){}
  virtual inline void set_Li_Lf(){}
  virtual inline void compute_batch_delta_gradient(){}
  virtual inline void compute_full_gradient(vector<D> &, vector<D> &){}

  L_Katyusha()
  {

  }

  void update_x(){
    for(L i=0; i<nfeatures; i++){
      x[i]= theta1*z[i]+ theta2*w[i]+ theta3*y[i];
    }
  }

  void update_y(){
    for(L i=0; i<nfeatures; i++){
      y[i]=theta1*z[i]+ theta2*w[i]+ theta3*y[i];
    }
  }

    void update_gk(){
      batch_delta_gradient.clear();
      batch_delta_gradient.resize(nfeatures,0);
      compute_batch_delta_gradient();
      for(L i=0;i<nfeatures;i++){
        gk[i]=batch_delta_gradient[i]+gradient_f_w[i];
      }
    }

    void update_z(){
      prox_of_g(oneovereta,z,gk,z);
    }

    void update_w(){
      D yi=gsl_rng_uniform(rng);
      if(yi<=p){
        w=x;
        compute_full_gradient(w, gradient_f_w);
      }
    }

    void compute_upper_bound_of_optimality_gap(){
      compute_full_gradient(x, gradient_f_x);
      prox_of_g(Lf,x,gradient_f_x,next_x);
      D tmp=0;
      for(L i=0;i<nfeatures;i++){
        tmp+=(next_x[i]-x[i])*(next_x[i]-x[i]);
        x[i]=next_x[i];
      }
      upper_bound=tmp*Lf*Lf/mu;
   	}

    void set_rng()
    {
      gsl_rng_env_setup();
      const gsl_rng_type * T;
      T = gsl_rng_default;
      rng = gsl_rng_alloc(T);
      gsl_rng_set(rng,time(NULL));
      //gsl_rng_set(rng, 27432042);

    }

    void set_print_every_N(L i){print_every_N=i;}

    void set_optimal_probability()
    {
      tilde_proba_vector.clear();
      tilde_proba_vector.resize(nsamples,0);
      proba_vector.clear();
      proba_vector.resize(nsamples,0);
      D sum=0;
      for(L i=0; i<nsamples;i++)
      {
        tilde_proba_vector[i]=Li[i];
        sum+=Li[i];
      }
      max_p=0;
      for(L i=0; i<nsamples;i++)
      {
        tilde_proba_vector[i]=tilde_proba_vector[i]/sum;
        proba_vector[i]=1-pow(1-tilde_proba_vector[i],tau);
        if(max_p<tilde_proba_vector[i])
        {
          max_p=tilde_proba_vector[i];
        }
      }
      cout<<"sum="<<sum<<"; proba"<<tilde_proba_vector[0]<<endl;
    }

    void set_group_sampling_probability(){
      proba_vector.clear();
      proba_vector.resize(nsamples,0);
      std::vector<D> q(nsamples);
      D sumq=0;
      cout<<"start setting group sampling probablity"<<endl;
      for(L i=0; i<nsamples;i++)
      {
        q[i]=Li[i];
        sumq+=q[i];
      }
      D maxq=0;
      L nb=0;
      D tmp=0;
      for(L i=0; i<nsamples;i++)
      {
        q[i]=q[i]/sumq*tau;
        if(q[i]>maxq) maxq=q[i];
        if(q[i]>1) {
          nb++; //count the number of elements larger than 1
          tmp+=q[i]-1;
        }
      }
      cout<<"maxq="<<maxq<<endl;
      if(maxq<=1){
        for(L i=0; i<nsamples;i++)
        {
          proba_vector[i]=q[i];
        }
      }else{
        for(L i=0; i<nsamples;i++)
        {
          if(q[i]>1) q[i]=1;
          else {
            D deltaq=min(1-q[i],tmp);
            tmp=tmp-deltaq;
            q[i]+=deltaq;
          }
          proba_vector[i]=q[i];
        }
      }
    }

    void set_L2(L p_mod, L u){
      if(p_mod==0){ //batch sampling mode
        D tmp=0;
        D st;
        for(L i=0;i<nsamples;i++){
          st=Li[i]/tilde_proba_vector[i];
          if (st>tmp) tmp=st;
        }
        L2=tmp/nsamples/tau;
        cout<<"L2="<<L2<<"; "<<tilde_proba_vector[0]<<endl;
      }
      else{ //group sampling mode
        D tmp=0;
        D st;
        for(L i=0;i<nsamples;i++){
          st=Li[i]/proba_vector[i];
          if(isolated_I[i]==1) st=st-Li[i];
          if(st>tmp) tmp=st;
        }
        L2=tmp/nsamples;
      }
      cout<<"set_L2="<<L2<<endl;
    }

    void set_uniform_probability(L p_mod)
    {
      if(p_mod==0) {
        //batch sampling (stochastic process as defined in the paper)
        proba_vector.clear();
        tilde_proba_vector.clear();
        tilde_proba_vector.resize(nsamples,1./nsamples);
        proba_vector.resize(nsamples,1-pow(1-1.0/nsamples,tau));
        max_p=1.0/nsamples;
        cout<<"uniform="<<tilde_proba_vector[0]<<endl;
      }
      else{ //group sampling
        D pi=(tau*c+0.0)/nsamples;
        proba_vector.clear();
        proba_vector.resize(nsamples,pi);
      }
    }

    void sort_p(){
      vector<pair<D,L> >a;
      for (L i = 0 ;i < nsamples ; i++) {
        a.push_back(make_pair(proba_vector[i],i)); // k = value, i = original index
      }
      sort(a.begin(),a.end());
      group_C.clear();
      group_C.resize(nsamples);
      isolated_I.clear();
      isolated_I.resize(nsamples,0);
      D tmp=0;
      D maxpi=0;
      D previous_i=0;
      index_group_C.clear();
      index_group_C.push_back(0);
      maxp_group_C.clear();
      sump_group_C.clear();

      for (L i = 0 ;i < nsamples ; i++){
        group_C[i]=a[nsamples-1-i].second;
        D pi=a[nsamples-1-i].first;
        if(tmp+pi<=1){
          if(pi>maxpi) maxpi=pi;
          tmp+=pi;
        }
        else{
          index_group_C.push_back(i);
          maxp_group_C.push_back(maxpi);
          sump_group_C.push_back(tmp);
          tmp=pi;
          maxpi=pi;
          if(i-previous_i==1) isolated_I[group_C[previous_i]]=1;
          previous_i=i;
        }
      }
      index_group_C.push_back(nsamples);
      maxp_group_C.push_back(maxpi);
      sump_group_C.push_back(tmp);
      nb_groups=sump_group_C.size();
      if(previous_i==nsamples-1) isolated_I[group_C[previous_i]]=1;
      cout<<"size of group="<<nb_groups<<endl;
    }

    inline L sampling(L n)
    {
      L i=(floor)(gsl_rng_uniform(rng)*n);
      D y=gsl_rng_uniform(rng);
      while(y*max_p>tilde_proba_vector[i])
      {
        i=(floor)(gsl_rng_uniform(rng)*n);
        y=gsl_rng_uniform(rng);
      }
      return i;
    }

    void set_theta(){
      theta2=L2/2/max(Lf,L2);
      if(Lf<=L2/p){
        D tmp=mu/L2/p;
        cout<<"L2="<<L2<<endl;
        cout<<"tmp="<<tmp<<endl;
        if(tmp>=1){
          theta1=theta2;
        }
        else{
          theta1=sqrt(tmp)*theta2;
        }
      }else
      theta1=min(sqrt(mu/Lf),p/2);
      theta3=1-theta1-theta2;
      eta=1./(theta1*(Lf+2*max(L2,Lf)));
      oneovereta=theta1*(Lf+2*max(L2,Lf));
      cout<<"Lf="<<Lf<<"; L2="<<L2<<"; theta1="<<theta1<<"; theta2="<<theta2<<"; theta3="<<theta3<<endl;
      cout<<"eta="<<eta<<endl;
    }

    void set_p(D scal_p){
        p=scal_p*tau/(0.0+nsamples);
        cout<<"changing probablity: "<<p<<endl;
    }

    void initialize(vector<D> & x0, D val_mu, L max_nb, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
    {
      cout<<"start initializing loopless Katyusha"<<" u="<<u<<endl;
      c=nb_c;
      noverc=nsamples/c;
      if(nb_tau>noverc) perror("tau should be less than n over c");
      tau=nb_tau;
      nb_of_iters_per_loop=floor(nsamples/(c*(tau+0.0)));
      batch_i.clear();
      batch_i.resize(nsamples,0);

      gk.clear();
      gk.resize(nfeatures,0);
      upper_bound=std::numeric_limits<double>::max();

      /**setup parameters**/
      max_nb_loops=max_nb;
      mu=val_mu;
      cout<<"mu="<<mu<<endl;
      set_rng();

      set_Li_Lf();

      running_time=0;
      /**setup probability**/
      if(p_mod==0)   // batch sampling (stochastic process as defined in the paper)
      {
        if(u==0)   // uniform sampling
        {
          set_uniform_probability(p_mod);
          uniform="uniform";
        }
        else{
          set_optimal_probability();
          uniform="nonuniform";
        }
      }
      else{  //group sampling
        if(u==0){
          set_uniform_probability(p_mod);
          uniform="uniform";
        }else{
          set_group_sampling_probability();
          uniform="nonuniform";
        }
        sort_p();
        cout<<"sort p"<<endl;
      }

      set_L2(p_mod,u);
      scaler=scal_p;
      set_p(scal_p);
      set_thetaS(p_mod);
      cout<<"set thetaS"<<endl;

      x.clear();
      y.clear();
      z.clear();
      w.clear();
      x.resize(nfeatures);
      y.resize(nfeatures);
      z.resize(nfeatures);
      w.resize(nfeatures);
      next_x.clear();
      next_x.resize(nfeatures);

      gradient_f_w.clear();
      gradient_f_w.resize(nfeatures);
      gradient_f_x.clear();
      gradient_f_x.resize(nfeatures);

      for(L i=0;i<nfeatures;i++)
      {
        x[i]=x0[i];
        y[i]= x0[i];
        w[i]=x0[i];
        z[i]=x0[i];
      }

      set_theta();
      cout<<"set stepsizes"<<endl;

      current_nb_iters=0;
      compute_primal_value();
      cout<<"primal value="<<primal_value<<endl;
      cout<<"Initialization  loopless Katyusha is finished!"<<endl;
    }

void compute_primal_value()
{
  D res=0;
  for(L i=0;i<nsamples;i++)
  {
    res+=value_of_phi_i(i);
  }
  D res2=value_of_g();
  primal_value= res/nsamples+ res2;
}

void compute_and_record_result()
{
  if(nb_loops%print_every_N==0){
    compute_primal_value();
    cout<<setprecision(9)<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<";  "<<running_time<<" primal value="<<primal_value<<endl;
    samp<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<" "<<running_time<<" "<<primal_value<<endl;
  }
}

void compute_and_record_result2()
{
  if(nb_loops%print_every_N==0){
    compute_primal_value();
    compute_upper_bound_of_optimality_gap();
    cout<<setprecision(9)<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<";  "<<running_time<<" primal value="<<primal_value<<" upper bound of F(x)-F^*= "<< upper_bound<< endl;
    samp<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<" "<<upper_bound<<" "<<running_time<<" "<<primal_value<<endl;
  }

}

void batch_sampling()
{
  batch_size=0;
  L i;
  for(L it_t=0;it_t<tau;it_t++)
  {
    i=sampling(nsamples);
    batch_i[batch_size]=i;
    batch_size++;
  }
}

void set_thetaS(L p_mod){
  theta_S.clear();
  theta_S.resize(nsamples,0);
  if(p_mod==0){
    for(L i=0;i<nsamples;i++)
    theta_S[i]=1.0/(tau*tilde_proba_vector[i]);
  }
  else{
    for(L i=0;i<nsamples;i++)
    theta_S[i]=1.0/proba_vector[i];
  }
}

void group_sampling(){
  batch_size=0;
  std::vector<L> sampled_groups;
  for(L i=0;i<nb_groups;i++){
    D y=gsl_rng_uniform(rng);
    if(y<sump_group_C[i])
    {
      sampled_groups.push_back(i);
      batch_size++;
    }
  }
  for(L t=0;t<batch_size;t++){
    L group_i=sampled_groups[t];
    L s1=index_group_C[group_i];
    L s2=index_group_C[group_i+1];
    L i=s1+(floor)(gsl_rng_uniform(rng)*(s2-s1));
    D y=gsl_rng_uniform(rng);
    D maxpi=maxp_group_C[group_i];
    while(y*maxpi>proba_vector[group_C[i]])
    {
      i=s1+(floor)(gsl_rng_uniform(rng)*(s2-s1));
      y=gsl_rng_uniform(rng);
    }
    batch_i[t]=group_C[i];
  }
}

void loopless(vector<D> & x0, string filename, D val_mu, L max_nb, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
{
  initialize(x0, val_mu, max_nb,  nb_tau, nb_c,  u, p_mod, scal_p);
  string sampname="results/L_Katyusha"+filename+uniform;
  if(p_mod==0) sampname=sampname+"_batch";
  else sampname=sampname+"_group";
  string scal_str;
  stringstream  scal_convert;
  scal_convert<<scal_p;
  scal_str=scal_convert.str();
  sampname+=scal_str;
  cout<<"running Loopless Katyusha "<<" ; "<<sampname<<endl;
  samp.open(sampname.c_str());

  nb_loops=0;
  nb_iters=0;

  srand48(27432042);
  //srand(time(NULL));

  D start;
  while(nb_loops<max_nb_loops)
  {
    compute_and_record_result();

    nb_loops++;

    start = std::clock();

    for(L it=0;it<nb_of_iters_per_loop;it++)
    {
      if(p_mod==0)
      batch_sampling();
      else
      group_sampling();

      start = std::clock();
      update_x();
      update_gk();
      update_z();
      update_y();
      update_w();
      update_gradient();
      current_nb_iters++;
      nb_iters++;
      running_time+= ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
    }
  }
}

void loopless2(vector<D> & x0, string filename, D val_mu, L max_nb, D epsilon, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
{
  cout<<"max_nb="<<max_nb<<endl;
  initialize(x0, val_mu, max_nb,  nb_tau, nb_c,  u, p_mod, scal_p);
  string sampname="results/L_"+filename+uniform;
  if(p_mod==0) sampname=sampname+"_batch";
  else sampname=sampname+"_group";
  string scal_str;
  stringstream  scal_convert;
  scal_convert<<scal_p;
  scal_str=scal_convert.str();
  sampname+=scal_str;
  cout<<"running Loopless Katyusha 2"<<" ; "<<sampname<<endl;
  samp.open(sampname.c_str());

  nb_loops=0;
  nb_iters=0;

  cout<<setprecision(9)<<"initial: "<< "primal: "<<primal_value << endl;
  samp<<((0.0+nb_iters)*c*tau/(nsamples))<<" "<<running_time<<" "<<primal_value<<endl;
  //srand48(27432042);
  srand(time(NULL));

  D start;
  cout<<"nb_loops="<<nb_loops<<"; max_nb_loops="<<max_nb_loops<<"; upper bound="<<upper_bound<<"; epsilon="<<epsilon<<endl;
  while(nb_loops<max_nb_loops && upper_bound> epsilon)
  {
    compute_and_record_result2();
    start = std::clock();

    nb_loops++;

    //cout<<"before the loop time elapsed="<<duration<<endl;
    start = std::clock();

    //cout<<"nb_of_iters_per_loop="<<nb_of_iters_per_loop<<endl;
    for(L it=0;it<nb_of_iters_per_loop;it++)
    {
      if(p_mod==0)
      batch_sampling();
      else
      group_sampling();

      start = std::clock();
      update_x();
      update_gk();
      update_z();
      update_y();
      update_w();
      update_gradient();
      current_nb_iters++;
      //cout<< "test 1"<< endl;
      //cout<< "test 2"<< endl;
      nb_iters++;
      running_time+= ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
    }
    //cout<< "epoch "<< nb_loops<< " finished"<< endl;
    //cout<<"after the loop time elapsed="<<duration<<endl;
  }
}

};

#endif /* MIN_SMOOTH_CONVEX_H */


text:
/* -*- linux-c -*- */
/* triple.c

   Copyright (C) 2002-2004 John M. Fregeau

   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 Foundation; either version 2 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
   MERCHANTABILITY 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <getopt.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include "fewbody.h"
#include "triple.h"

/* print the usage */
void print_usage(FILE *stream)
{
	fprintf(stream, "USAGE:\n");
	fprintf(stream, "  triple [options...]\n");
	fprintf(stream, "\n");
	fprintf(stream, "OPTIONS:\n");
	fprintf(stream, "  -m --m000 <m000/MSUN>        : set mass of star 0 of inner binary of triple [%.6g]\n", FB_M000/FB_CONST_MSUN);
	fprintf(stream, "  -n --m001 <m001/MSUN>        : set mass of star 1 of inner binary of triple [%.6g]\n", FB_M001/FB_CONST_MSUN);
	fprintf(stream, "  -o --m01 <m01/MSUN>          : set mass of outer star of triple [%.6g]\n", FB_M01/FB_CONST_MSUN);
	fprintf(stream, "  -r --r000 <r000/RSUN>        : set radius of star 0 of inner binary of triple [%.6g]\n", FB_R000/FB_CONST_RSUN);
	fprintf(stream, "  -g --r001 <r001/RSUN>        : set radius of star 1 of inner binary of triple [%.6g]\n", FB_R001/FB_CONST_RSUN);
	fprintf(stream, "  -i --r01 <r01/RSUN>          : set radius of outer star of triple [%.6g]\n", FB_R01/FB_CONST_RSUN);
	fprintf(stream, "  -a --a00 <a00/AU>            : set inner semimajor axis of triple [%.6g]\n", FB_A00/FB_CONST_AU);
	fprintf(stream, "  -Q --a0 <a0/AU>              : set outer semimajor axis of triple [%.6g]\n", FB_A0/FB_CONST_AU);
	fprintf(stream, "  -e --e00 <e00>               : set inner eccentricity of triple [%.6g]\n", FB_E00);
	fprintf(stream, "  -F --e0 <e0>                 : set outer eccentricity of triple [%.6g]\n", FB_E0);
	fprintf(stream, "  -t --tstop <tstop/t_dyn>     : set stopping time [%.6g]\n", FB_TSTOP);
	fprintf(stream, "  -D --dt <dt/t_dyn>           : set approximate output dt [%.6g]\n", FB_DT);
	fprintf(stream, "  -c --tcpustop <tcpustop/sec> : set cpu stopping time [%.6g]\n", FB_TCPUSTOP);
	fprintf(stream, "  -A --absacc <absacc>         : set integrator's absolute accuracy [%.6g]\n", FB_ABSACC);
	fprintf(stream, "  -R --relacc <relacc>         : set integrator's relative accuracy [%.6g]\n", FB_RELACC);
	fprintf(stream, "  -N --ncount <ncount>         : set number of integration steps between calls\n");
	fprintf(stream, "                                 to fb_classify() [%d]\n", FB_NCOUNT);
	fprintf(stream, "  -z --tidaltol <tidaltol>     : set tidal tolerance [%.6g]\n", FB_TIDALTOL);
	fprintf(stream, "  -x --fexp <f_exp>            : set expansion factor of merger product [%.6g]\n", FB_FEXP);
	fprintf(stream, "  -k --ks                      : turn K-S regularization on or off [%d]\n", FB_KS);
	fprintf(stream, "  -s --seed                    : set random seed [%ld]\n", FB_SEED);
	fprintf(stream, "  -d --debug                   : turn on debugging\n");
	fprintf(stream, "  -V --version                 : print version info\n");
	fprintf(stream, "  -h --help                    : display this help text\n");
}

/* calculate the units used */
int calc_units(fb_obj_t *obj[1], fb_units_t *units)
{
	double m0, m00, m01, m000, m001, a00, a0;

	m0 = obj[0]->m;
	m00 = obj[0]->obj[0]->m;
	m01 = obj[0]->obj[1]->m;
	m000 = obj[0]->obj[0]->obj[0]->m;
	m001 = obj[0]->obj[0]->obj[1]->m;

	a0 = obj[0]->a;
	a00 = obj[0]->obj[0]->a;

	/* Unit of velocity is approximate relative orbital speed of inner binary,
	   unit of length is semimajor axis of inner binary;
	   therefore, unit of time is approximately 1 inner orbital period. */
	units->v = sqrt(FB_CONST_G*(m000+m001)/a00);
	units->l = a00;
	units->t = units->l / units->v;
	units->m = units->l * fb_sqr(units->v) / FB_CONST_G;
	units->E = units->m * fb_sqr(units->v);

	return(0);
}

/* the main attraction */
int main(int argc, char *argv[])
{
	int i, j;
	unsigned long int seed;
	double m000, m001, m01, r000, r001, r01, a00, a0, e00, e0;
	double Ei, Lint[3], Li[3], t;
	fb_hier_t hier;
	fb_input_t input;
	fb_ret_t retval;
	fb_units_t units;
	char string1[FB_MAX_STRING_LENGTH], string2[FB_MAX_STRING_LENGTH];
	gsl_rng *rng;
	const gsl_rng_type *rng_type=gsl_rng_mt19937;
	const char *short_opts = "m:n:o:r:g:i:a:Q:e:F:t:D:c:A:R:N:z:x:k:s:dVh";
	const struct option long_opts[] = {
		{"m000", required_argument, NULL, 'm'},
		{"m001", required_argument, NULL, 'n'},
		{"m01", required_argument, NULL, 'o'},
		{"r000", required_argument, NULL, 'r'},
		{"r001", required_argument, NULL, 'g'},
		{"r01", required_argument, NULL, 'i'},
		{"a00", required_argument, NULL, 'a'},
		{"a0", required_argument, NULL, 'Q'},
		{"e00", required_argument, NULL, 'e'},
		{"e0", required_argument, NULL, 'F'},
		{"tstop", required_argument, NULL, 't'},
		{"dt", required_argument, NULL, 'D'},
		{"tcpustop", required_argument, NULL, 'c'},
		{"absacc", required_argument, NULL, 'A'},
		{"relacc", required_argument, NULL, 'R'},
		{"ncount", required_argument, NULL, 'N'},
		{"tidaltol", required_argument, NULL, 'z'},
		{"fexp", required_argument, NULL, 'x'},
		{"ks", required_argument, NULL, 'k'},
		{"seed", required_argument, NULL, 's'},
		{"debug", no_argument, NULL, 'd'},
		{"version", no_argument, NULL, 'V'},
		{"help", no_argument, NULL, 'h'},
		{NULL, 0, NULL, 0}
	};

	/* set parameters to default values */
	m000 = FB_M000;
	m001 = FB_M001;
	m01 = FB_M01;
	r000 = FB_R000;
	r001 = FB_R001;
	r01 = FB_R01;
	a00 = FB_A00;
	a0 = FB_A0;
	e00 = FB_E00;
	e0 = FB_E0;
	input.ks = FB_KS;
	input.tstop = FB_TSTOP;
	input.Dflag = 0;
	input.dt = FB_DT;
	input.tcpustop = FB_TCPUSTOP;
	input.absacc = FB_ABSACC;
	input.relacc = FB_RELACC;
	input.ncount = FB_NCOUNT;
	input.tidaltol = FB_TIDALTOL;
	input.fexp = FB_FEXP;
	seed = FB_SEED;
	fb_debug = FB_DEBUG;

	while ((i = getopt_long(argc, argv, short_opts, long_opts, NULL)) != -1) {
		switch (i) {
		case 'm':
			m000 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'n':
			m001 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'o':
			m01 = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'r':
			r000 = atof(optarg) * FB_CONST_RSUN;
			break;
		case 'g':
			r001 = atof(optarg) * FB_CONST_RSUN;
			break;
		case 'i':
			r01 = atof(optarg) * FB_CONST_RSUN;
			break;
		case 'a':
			a00 = atof(optarg) * FB_CONST_AU;
			break;
		case 'Q':
			a0 = atof(optarg) * FB_CONST_AU;
			break;
		case 'e':
			e00 = atof(optarg);
			if (e00 >= 1.0) {
				fprintf(stderr, "e00 must be less than 1\n");
				return(1);
			}
			break;
		case 'F':
			e0 = atof(optarg);
			if (e0 >= 1.0) {
				fprintf(stderr, "e0 must be less than 1\n");
				return(1);
			}
			break;
		case 't':
			input.tstop = atof(optarg);
			break;
		case 'D':
			input.Dflag = 1;
			input.dt = atof(optarg);
			break;
		case 'c':
			input.tcpustop = atof(optarg);
			break;
		case 'A':
			input.absacc = atof(optarg);
			break;
		case 'R':
			input.relacc = atof(optarg);
			break;
		case 'N':
			input.ncount = atoi(optarg);
			break;
		case 'z':
			input.tidaltol = atof(optarg);
			break;
		case 'x':
			input.fexp = atof(optarg);
			break;
		case 'k':
			input.ks = atoi(optarg);
			break;
		case 's':
			seed = atol(optarg);
			break;
		case 'd':
			fb_debug = 1;
			break;
		case 'V':
			fb_print_version(stdout);
			return(0);
		case 'h':
			fb_print_version(stdout);
			fprintf(stdout, "\n");
			print_usage(stdout);
			return(0);
		default:
			break;
		}
	}

	/* check to make sure there was nothing crazy on the command line */
	if (optind < argc) {
		print_usage(stdout);
		return(1);
	}

	/* initialize a few things for integrator */
	t = 0.0;
	hier.nstarinit = 3;
	hier.nstar = 3;
	fb_malloc_hier(&hier);
	fb_init_hier(&hier);

	/* put stuff in log entry */
	snprintf(input.firstlogentry, FB_MAX_LOGENTRY_LENGTH, "  command line:");
	for (i=0; i<argc; i++) {
		snprintf(&(input.firstlogentry[strlen(input.firstlogentry)]),
			 FB_MAX_LOGENTRY_LENGTH-strlen(input.firstlogentry), " %s", argv[i]);
	}
	snprintf(&(input.firstlogentry[strlen(input.firstlogentry)]),
		 FB_MAX_LOGENTRY_LENGTH-strlen(input.firstlogentry), "\n");

	/* print out values of paramaters */
	fprintf(stderr, "PARAMETERS:\n");
	fprintf(stderr, "  ks=%d  seed=%ld\n", input.ks, seed);
	fprintf(stderr, "  a00=%.6g AU  e00=%.6g  m000=%.6g MSUN  m001=%.6g MSUN  r000=%.6g RSUN  r001=%.6g RSUN\n", \
		a00/FB_CONST_AU, e00, m000/FB_CONST_MSUN, m001/FB_CONST_MSUN, r000/FB_CONST_RSUN, r001/FB_CONST_RSUN);
	fprintf(stderr, "  a0=%.6g AU  e0=%.6g  m01=%.6g MSUN  r01=%.6g RSUN\n", \
		a0/FB_CONST_AU, e0, m01/FB_CONST_MSUN, r01/FB_CONST_RSUN);
	fprintf(stderr, "  tstop=%.6g  tcpustop=%.6g\n", \
		input.tstop, input.tcpustop);
	fprintf(stderr, "  tidaltol=%.6g  abs_acc=%.6g  rel_acc=%.6g  ncount=%d  fexp=%.6g\n\n", \
		input.tidaltol, input.absacc, input.relacc, input.ncount, input.fexp);

	/* initialize GSL rng */
	gsl_rng_env_setup();
	rng = gsl_rng_alloc(rng_type);
	gsl_rng_set(rng, seed);

	/* create hierarchies */
	hier.narr[2] = 1;
	hier.narr[3] = 1;
	/* inner binary of triple */
	hier.hier[hier.hi[2]+0].obj[0] = &(hier.hier[hier.hi[1]+0]);
	hier.hier[hier.hi[2]+0].obj[1] = &(hier.hier[hier.hi[1]+1]);
	hier.hier[hier.hi[2]+0].t = t;
	/* outer binary of triple */
	hier.hier[hier.hi[3]+0].obj[0] = &(hier.hier[hier.hi[2]+0]);
	hier.hier[hier.hi[3]+0].obj[1] = &(hier.hier[hier.hi[1]+2]);
	hier.hier[hier.hi[3]+0].t = t;

	/* give the objects some properties */
	for (j=0; j<hier.nstar; j++) {
		hier.hier[hier.hi[1]+j].ncoll = 1;
		hier.hier[hier.hi[1]+j].id[0] = j;
		snprintf(hier.hier[hier.hi[1]+j].idstring, FB_MAX_STRING_LENGTH, "%d", j);
		hier.hier[hier.hi[1]+j].n = 1;
		hier.hier[hier.hi[1]+j].obj[0] = NULL;
		hier.hier[hier.hi[1]+j].obj[1] = NULL;
		hier.hier[hier.hi[1]+j].Eint = 0.0;
		hier.hier[hier.hi[1]+j].Lint[0] = 0.0;
		hier.hier[hier.hi[1]+j].Lint[1] = 0.0;
		hier.hier[hier.hi[1]+j].Lint[2] = 0.0;
	}

	hier.hier[hier.hi[1]+0].R = r000;
	hier.hier[hier.hi[1]+1].R = r001;
	hier.hier[hier.hi[1]+2].R = r01;

	hier.hier[hier.hi[1]+0].m = m000;
	hier.hier[hier.hi[1]+1].m = m001;
	hier.hier[hier.hi[1]+2].m = m01;

	hier.hier[hier.hi[2]+0].m = m000 + m001;
	hier.hier[hier.hi[3]+0].m = m000 + m001 + m01;

	hier.hier[hier.hi[2]+0].a = a00;
	hier.hier[hier.hi[3]+0].a = a0;

	hier.hier[hier.hi[2]+0].e = e00;
	hier.hier[hier.hi[3]+0].e = e0;

	hier.nobj = 1;
	hier.obj[0] = &(hier.hier[hier.hi[3]+0]);
	hier.obj[1] = NULL;
	hier.obj[2] = NULL;

	/* get the units and normalize */
	calc_units(hier.obj, &units);
	fb_normalize(&hier, units);

	/* place triple at origin */
	for (j=0; j<3; j++) {
		hier.obj[0]->x[j] = 0.0;
		hier.obj[0]->v[j] = 0.0;
	}

	/* randomize binary orientations and downsync */
	fb_randorient(&(hier.hier[hier.hi[3]+0]), rng);
	fb_downsync(&(hier.hier[hier.hi[3]+0]), t);
	fb_randorient(&(hier.hier[hier.hi[2]+0]), rng);
	fb_downsync(&(hier.hier[hier.hi[2]+0]), t);

	fprintf(stderr, "UNITS:\n");
	fprintf(stderr, "  v=%.6g km/s  l=%.6g AU  t=t_dyn=%.6g yr\n", \
		units.v/1.0e5, units.l/FB_CONST_AU, units.t/FB_CONST_YR);
	fprintf(stderr, "  M=%.6g M_sun  E=%.6g erg\n\n", units.m/FB_CONST_MSUN, units.E);

	/* trickle down properties (not sure if this is actually needed here, but it doesn't harm anything) */
	fb_trickle(&hier, t);

	/* store the initial energy and angular momentum*/
	Ei = fb_petot(&(hier.hier[hier.hi[1]]), hier.nstar) + fb_ketot(&(hier.hier[hier.hi[1]]), hier.nstar) +
		fb_einttot(&(hier.hier[hier.hi[1]]), hier.nstar);
	fb_angmom(&(hier.hier[hier.hi[1]]), hier.nstar, Li);
	fb_angmomint(&(hier.hier[hier.hi[1]]), hier.nstar, Lint);
	for (j=0; j<3; j++) {
		Li[j] += Lint[j];
	}

	/* integrate along */
	fb_dprintf("calling fewbody()...\n");

	/* call fewbody! */
	retval = fewbody(input, &hier, &t);

	/* print information to screen */
	fprintf(stderr, "OUTCOME:\n");
	if (retval.retval == 1) {
		fprintf(stderr, "  encounter complete:  t=%.6g (%.6g yr)  %s  (%s)\n\n",
			t, t * units.t/FB_CONST_YR,
			fb_sprint_hier(hier, string1),
			fb_sprint_hier_hr(hier, string2));
	} else {
		fprintf(stderr, "  encounter NOT complete:  t=%.6g (%.6g yr)  %s  (%s)\n\n",
			t, t * units.t/FB_CONST_YR,
			fb_sprint_hier(hier, string1),
			fb_sprint_hier_hr(hier, string2));
	}

	fb_dprintf("there were %ld integration steps\n", retval.count);
	fb_dprintf("fb_classify() was called %ld times\n", retval.iclassify);

	fprintf(stderr, "FINAL:\n");
	fprintf(stderr, "  t_final=%.6g (%.6g yr)  t_cpu=%.6g s\n", \
		t, t*units.t/FB_CONST_YR, retval.tcpu);

	fprintf(stderr, "  L0=%.6g  DeltaL/L0=%.6g  DeltaL=%.6g\n", fb_mod(Li), retval.DeltaLfrac, retval.DeltaL);
	fprintf(stderr, "  E0=%.6g  DeltaE/E0=%.6g  DeltaE=%.6g\n", Ei, retval.DeltaEfrac, retval.DeltaE);
	fprintf(stderr, "  Rmin=%.6g (%.6g RSUN)  Rmin_i=%d  Rmin_j=%d\n", \
		retval.Rmin, retval.Rmin*units.l/FB_CONST_RSUN, retval.Rmin_i, retval.Rmin_j);
	fprintf(stderr, "  Nosc=%d (%s)\n", retval.Nosc, (retval.Nosc>=1?"resonance":"non-resonance"));

	/* free GSL stuff */
	gsl_rng_free(rng);

	/* free our own stuff */
	fb_free_hier(hier);

	/* done! */
	return(0);
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos