--- /dev/null
+all:
+ make -C src
+ make -C test
+
+clean:
+ make -C src clean
+ make -C test clean
--- /dev/null
+- Composantes fortement connexes
+
+- Precision arbitraire
+ http://www.mpfr.org
+
+- Calculer les derivées par les differences
+ pour la moyenne
+ plus efficace que calculer la derivée directement
+
+- Réutiliser les termes dans la dérivation des équations (DAGuefier)
+ Ajouter un compteur de références dans les équations
+ ou
+ Forme normale (éviter les produits entre expressions compliquées)
+
+- Analyse de singularité
+ dépend des CFC
+
+- Schröder
+ quand la singularité est en racine
+ marche bien que pile sur la singularité, sinon il faut basculer à Newton simple quand on est proche
+
+- Remplacer la dichotomie par la resulution de l'equation |Y - (dH/dY)(x)| = 0
+
+- Vérifier le resultat des appels à malloc
+
+- Générer directement le code du jacobien sans créer l'AST correspondant
+ Réutiliser le code des sous-expressions communes
--- /dev/null
+CFLAGS=-Wall -fPIC -g
+LDFLAGS=-lm -Wl,-rpath .
+
+all: libcombstruct.a libcombstruct.so libcombstruct_compile.so \
+ libcombstruct_blas.so libcombstruct_compile_blas.so
+
+%_compile.o: %.c
+ $(CC) -c $(CFLAGS) -DCOMPILE -o $@ $<
+
+common.o: common.c common.h
+
+diff.o: diff.c diff.h equation.h combstruct.h system.h common.h
+
+equation.o: equation.c equation.h combstruct.h
+
+eval.o: eval.c eval.h combstruct.h common.h equation.h
+
+iter.o: iter.c combstruct.h eval.h system.h jacob.h diff.h equation.h
+
+iter_blas.o: iter.c combstruct.h eval.h system.h jacob.h diff.h equation.h
+ $(CC) -c $(CFLAGS) -DBLAS -o $@ $<
+
+iter_compile.o: iter.c combstruct.h eval.h system.h jacob.h diff.h \
+ equation.h compile.h
+
+iter_compile_blas.o: iter.c combstruct.h eval.h system.h jacob.h diff.h \
+ equation.h compile.h
+ $(CC) -c $(CFLAGS) -DCOMPILE -DBLAS -o $@ $<
+
+jacob.o: jacob.c jacob.h combstruct.h equation.h system.h diff.h
+
+print.o: print.c print.h combstruct.h jacob.h common.h equation.h \
+ system.h
+
+system.o: system.c system.h combstruct.h equation.h
+
+system_compile.o: system.c system.h combstruct.h equation.h compile.h
+
+compile.o: compile.c compile.h
+ $(CC) -c $(CFLAGS) -DCOMPILE -o $@ $<
+
+libcombstruct.a: system.o iter.o jacob.o diff.o eval.o print.o \
+ equation.o common.o
+ ar rcs $@ $^
+
+libcombstruct.so: system.o iter.o jacob.o diff.o eval.o print.o \
+ equation.o common.o
+ $(CC) $(LDFLAGS) -shared -o $@ $^
+
+libcombstruct_compile.so: compile.o system_compile.o iter_compile.o \
+ jacob.o diff.o print.o equation.o common.o
+ $(CC) $(LDFLAGS) -ldl -shared -o $@ $^
+
+libcombstruct_blas.so: system.o iter_blas.o jacob.o diff.o eval.o print.o \
+ equation.o common.o
+ $(CC) $(LDFLAGS) -lblas -shared -o $@ $^
+
+libcombstruct_compile_blas.so: compile.o system_compile.o iter_compile_blas.o \
+ jacob.o diff.o print.o equation.o common.o
+ $(CC) $(LDFLAGS) -lblas -ldl -shared -o $@ $^
+
+clean:
+ rm -f *.o *.so
--- /dev/null
+#ifndef _COMBSTRUCT_H
+# define _COMBSTRUCT_H
+
+/* Equation */
+typedef struct comb_eq *comb_eq;
+
+comb_eq nil();
+
+comb_eq epsilon();
+
+comb_eq atom();
+
+comb_eq sum(comb_eq, comb_eq);
+
+comb_eq prod(comb_eq, comb_eq);
+
+comb_eq coef(double, comb_eq);
+
+comb_eq seq(comb_eq);
+
+comb_eq ref(int);
+
+void free_eq(comb_eq);
+
+/* System */
+typedef struct comb_sys *comb_sys;
+
+comb_sys sys_new(int);
+
+void sys_add_eq(comb_sys, int, comb_eq);
+
+int sys_len(comb_sys);
+
+double sys_val(comb_sys, int);
+
+void free_sys(comb_sys);
+
+/* Iter */
+// TODO: retourner le max des 1/(1 - y') ?
+void eval_sys(comb_sys, double, double, int *);
+
+void eval_point_sys(comb_sys, double, double, int *);
+
+double sing_sys(comb_sys, double, double);
+
+// last argument must be twice the size of the system
+double mean_sys(comb_sys, double, double, double);
+
+# ifndef NOPRINT
+/* Print */
+void print_sys(comb_sys);
+# endif
+#endif
--- /dev/null
+#include "common.h"
+
+#ifdef NOPRINT
+void error(char *s)
+{
+}
+#endif
--- /dev/null
+#ifndef _COMMON_H
+# define _COMMON_H
+void my_error(char *);
+#endif
--- /dev/null
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include <lightning.h>
+#include "compile.h"
+#include "common.h"
+#include "equation.h"
+// TODO: rendre dynamique la taille du buffer
+#define BUFFER_SIZE (1<<28)
+#define STACK_SIZE 1024
+
+double zero(double *v, double z, int *err)
+{
+ return 0.;
+}
+
+static int stack_ptr;
+
+// TODO: utiliser les registres libres avant d'aller en mémoire
+
+void stack_push(int reg)
+{
+ jit_stxi_d (stack_ptr, JIT_FP, reg);
+ stack_ptr += sizeof (double);
+}
+
+void stack_pop(int reg)
+{
+ stack_ptr -= sizeof (double);
+ jit_ldxi_d (reg, JIT_FP, stack_ptr);
+}
+
+void compile_expr(comb_eq c)
+{
+ switch (c->type) {
+ case Epsilon:
+ jit_movi_d(JIT_FPR0, 1.0);
+ break;
+ case Atom:
+ jit_movr_d(JIT_FPR0, JIT_FPR2);
+ break;
+ case Ref:
+ // JIT_FPR0 = v[c->ref-1]
+ jit_ldxi_d(JIT_FPR0, JIT_R0, sizeof(double) * (c->ref-1));
+ break;
+ case Sum:
+ compile_expr(c->operand1);
+ stack_push(JIT_FPR0);
+ compile_expr(c->operand2);
+ stack_pop(JIT_FPR1);
+ jit_addr_d(JIT_FPR0, JIT_FPR1, JIT_FPR0);
+ break;
+ case Prod:
+ compile_expr(c->operand1);
+ stack_push(JIT_FPR0);
+ compile_expr(c->operand2);
+ stack_pop(JIT_FPR1);
+ jit_mulr_d(JIT_FPR0, JIT_FPR1, JIT_FPR0);
+ break;
+ case Coef:
+ compile_expr(c->operand1);
+ jit_movi_d(JIT_FPR1, c->coef);
+ jit_mulr_d(JIT_FPR0, JIT_FPR1, JIT_FPR0);
+ break;
+ case Seq:
+ {
+ jit_insn *ref;
+
+ compile_expr(c->operand1);
+ // IF JIT_FPR0 >= 1 THEN *err = 1
+ jit_floorr_d_i(JIT_R2, JIT_FPR0);
+ ref = jit_blti_i(jit_forward(), JIT_R2, 1);
+ jit_str_i(JIT_R1, JIT_R2);
+ jit_patch(ref);
+ jit_movi_d(JIT_FPR1, 1.0);
+ jit_subr_d(JIT_FPR0, JIT_FPR1, JIT_FPR0);
+ jit_divr_d(JIT_FPR0, JIT_FPR1, JIT_FPR0);
+ break;
+ }
+ default:
+ my_error("ERROR: eval_eq unknown equation type\n");
+ }
+}
+
+f_eval compile_eq(comb_eq c)
+{
+ f_eval fn;
+ int in;
+
+ if (c->type == Nil)
+ return zero;
+ fn = (f_eval) (jit_get_ip().dptr);
+ jit_leaf(3);
+ in = jit_arg_p();
+ // TODO: verifier que ca ne coute pas trop cher
+ stack_ptr = jit_allocai(STACK_SIZE * sizeof (double));
+
+ jit_getarg_p(JIT_R0, in);
+ in = jit_arg_d();
+ jit_getarg_d(JIT_FPR2, in);
+ in = jit_arg_p();
+ jit_getarg_p(JIT_R1, in);
+ compile_expr(c);
+ jit_movr_d(JIT_FPRET, JIT_FPR0);
+ jit_ret();
+ return fn;
+}
+
+void compile_sys(comb_sys s)
+{
+ int i, j;
+ int l;
+ static jit_insn codeBuffer[BUFFER_SIZE];
+ f_eval *fp;
+
+ if (s->f)
+ return;
+ jit_set_ip(codeBuffer);
+ l = s->len;
+ s->f = malloc(l*sizeof(f_eval));
+ s->fj = malloc(l*sizeof(f_eval *));
+ fp = malloc(l*l*sizeof(f_eval));
+ for (i = 0; i < l; i++) {
+ s->f[s->refs[i]] = compile_eq(s->eqs[i]);
+ s->fj[i] = &fp[s->refs[i]*l];
+ for (j = 0; j < l; j++) {
+ s->fj[s->refs[i]][s->refs[j]] = compile_eq(s->jacob[i][j]);
+ }
+ }
+ jit_flush_code(codeBuffer, jit_get_ip().ptr);
+}
--- /dev/null
+#ifndef _COMPILE_H
+# define _COMPILE_H
+# include "combstruct.h"
+# include "system.h"
+
+void compile_sys(comb_sys);
+
+#endif
--- /dev/null
+#include <stdlib.h>
+#include "diff.h"
+#include "common.h"
+#include "equation.h"
+#include "system.h"
+
+comb_eq copy_eq(comb_eq c)
+{
+ comb_eq res;
+
+ switch (c->type) {
+ case Epsilon:
+ res = epsilon();
+ break;
+ case Atom:
+ res = atom();
+ break;
+ case Ref:
+ res = ref(c->ref);
+ break;
+ case Sum:
+ res = sum(copy_eq(c->operand1), copy_eq(c->operand2));
+ break;
+ case Prod:
+ res = prod(copy_eq(c->operand1), copy_eq(c->operand2));
+ break;
+ case Coef:
+ res = coef(c->coef, copy_eq(c->operand1));
+ break;
+ case Seq:
+ res = seq(copy_eq(c->operand1));
+ break;
+ default:
+ my_error("ERROR: copy_eq unknown equation type\n");
+ return NULL;
+ }
+ return res;
+}
+
+comb_eq diff_eq_aux(comb_eq c, int d, int l)
+{
+ comb_eq res;
+
+ switch (c->type) {
+ case Nil:
+ res = NULL;
+ break;
+ case Epsilon:
+ res = NULL;
+ break;
+ case Atom:
+ if (d == -1)
+ res = epsilon();
+ else
+ res = NULL;
+ break;
+ case Ref:
+ if (d == -1) {
+ res = ref(c->ref + l);
+ } else if (d == c->ref) {
+ res = epsilon();
+ } else {
+ res = NULL;
+ }
+ break;
+ case Sum:
+ {
+ comb_eq dop1, dop2;
+
+ dop1 = diff_eq_aux(c->operand1, d, l);
+ dop2 = diff_eq_aux(c->operand2, d, l);
+ if (dop1 == NULL)
+ res = dop2;
+ else if (dop2 == NULL)
+ res = dop1;
+ else
+ res = sum(dop1, dop2);
+ }
+ break;
+ case Prod:
+ {
+ comb_eq dop1, dop2;
+
+ dop1 = diff_eq_aux(c->operand1, d, l);
+ dop2 = diff_eq_aux(c->operand2, d, l);
+ if (dop1 == NULL && dop2 == NULL)
+ res = NULL;
+ else if (dop1 == NULL)
+ res = prod(copy_eq(c->operand1),dop2);
+ else if (dop2 == NULL)
+ res = prod(dop1, copy_eq(c->operand2));
+ else
+ res = sum(prod(dop1,copy_eq(c->operand2)), prod(copy_eq(c->operand1),dop2));
+ }
+ break;
+ case Coef:
+ {
+ comb_eq dop;
+
+ dop = diff_eq_aux(c->operand1, d, l);
+ if (dop == NULL)
+ res = NULL;
+ else
+ res = coef(c->coef, dop);
+ }
+ break;
+ case Seq:
+ {
+ comb_eq dop;
+
+ dop = diff_eq_aux(c->operand1, d, l);
+ if (dop == NULL)
+ res = NULL;
+ else
+ res = prod(copy_eq(c), prod(dop, copy_eq(c)));
+ }
+ break;
+ default:
+ my_error("ERROR: diff_eq unknown equation type\n");
+ return NULL;
+ }
+ return res;
+}
+
+comb_eq diff_eq_ref(comb_eq c, int r)
+{
+ comb_eq res;
+
+ res = diff_eq_aux(c, r, 0);
+ if (res == NULL)
+ return nil();
+ return res;
+}
+
+comb_eq diff_eq_z(comb_eq c, int l)
+{
+ comb_eq res;
+
+ res = diff_eq_aux(c, -1, l);
+ if (res == NULL)
+ return nil();
+ return res;
+}
+
+comb_sys diff_sys_z(comb_sys s)
+{
+ comb_sys res;
+ int i, l;
+
+ l = s->len;
+ res = sys_new(2*l);
+ for (i = 0; i < l; i++)
+ sys_add_eq(res, i, copy_eq(s->eqs[i]));
+ for (i = 0; i < l; i++)
+ sys_add_eq(res, l+i, diff_eq_z(s->eqs[i], l));
+ return res;
+}
--- /dev/null
+#ifndef _DIFF_H
+# define _DIFF_H
+# include "equation.h"
+# include "system.h"
+
+comb_eq diff_eq_ref(comb_eq, int);
+
+comb_sys diff_sys_z(comb_sys);
+
+#endif
--- /dev/null
+#include <stdlib.h>
+#include "equation.h"
+
+comb_eq eq_init(constr type, comb_eq operand1, comb_eq operand2, int ref, double coef)
+{
+ comb_eq res;
+
+ res = malloc(sizeof(struct comb_eq));
+ res->type = type;
+ res->operand1 = operand1;
+ res->operand2 = operand2;
+ res->ref = ref;
+ res->coef = coef;
+ return res;
+}
+
+comb_eq nil()
+{
+ return eq_init(Nil, NULL, NULL, -1, 1.0);
+}
+
+comb_eq epsilon()
+{
+ return eq_init(Epsilon, NULL, NULL, -1, 1.0);
+}
+
+comb_eq atom()
+{
+ return eq_init(Atom, NULL, NULL, -1, 1.0);
+}
+
+comb_eq sum(comb_eq operand1, comb_eq operand2)
+{
+ return eq_init(Sum, operand1, operand2, -1, 1.0);
+}
+
+comb_eq prod(comb_eq operand1, comb_eq operand2)
+{
+ return eq_init(Prod, operand1, operand2, -1, 1.0);
+}
+
+comb_eq coef(double c, comb_eq operand)
+{
+ return eq_init(Coef, operand, NULL, -1, c);
+}
+
+comb_eq seq(comb_eq operand)
+{
+ return eq_init(Seq, operand, NULL, -1, 1.0);
+}
+
+comb_eq ref(int r)
+{
+ return eq_init(Ref, NULL, NULL, r, 1.0);
+}
+
+void free_eq(comb_eq c)
+{
+ if (!c)
+ return;
+ free_eq(c->operand1);
+ free_eq(c->operand2);
+ free(c);
+}
--- /dev/null
+#ifndef _EQUATION_H
+# define _EQUATION_H
+# include "combstruct.h"
+
+typedef enum { Epsilon = 0x0BEEF0,
+ Atom,
+ Sum,
+ Prod,
+ Coef,
+ Seq,
+ Ref,
+ Nil } constr;
+
+struct comb_eq {
+ constr type;
+ struct comb_eq *operand1;
+ struct comb_eq *operand2;
+ int ref;
+ double coef;
+};
+
+#endif
--- /dev/null
+#include "eval.h"
+#include "common.h"
+#include "equation.h"
+#include "system.h"
+
+static comb_sys system;
+
+/* TODO: ajouter mutex ou autre */
+void set_sys(comb_sys s)
+{
+ system = s;
+}
+
+double eval_eq(comb_eq c, double *v, double z, int *err)
+{
+ switch (c->type) {
+ case Nil:
+ return 0.;
+ case Epsilon:
+ return 1.;
+ case Atom:
+ return z;
+ case Ref:
+ return v[system->refs[c->ref-1]];
+ case Sum:
+ return eval_eq(c->operand1, v, z, err) + eval_eq(c->operand2, v, z, err);
+ case Prod:
+ return eval_eq(c->operand1, v, z, err) * eval_eq(c->operand2, v, z, err);
+ case Coef:
+ return c->coef * eval_eq(c->operand1, v, z, err);
+ case Seq:
+ {
+ double t;
+
+ t = eval_eq(c->operand1, v, z, err);
+ if (t >= 1.)
+ *err = 1;
+ return 1./(1. - t);
+ }
+ default:
+ my_error("ERROR: eval_eq unknown equation type\n");
+ return 0.;
+ }
+}
--- /dev/null
+#ifndef _EVAL_H
+# define _EVAL_H
+# include "combstruct.h"
+
+void set_sys(comb_sys);
+
+double eval_eq(comb_eq, double *, double, int *);
+#endif
--- /dev/null
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#ifdef BLAS
+# include <cblas.h>
+#endif
+#ifndef COMPILE
+# include "eval.h"
+#endif
+#include "system.h"
+#include "jacob.h"
+#include "diff.h"
+
+void multm(double **a, double **b, double **c, int n)
+{
+#ifdef BLAS
+ cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1, &b[0][0], n, &c[0][0], n, 0, &a[0][0], n);
+#else
+ int i, j, k;
+
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++) {
+ a[i][j] = 0;
+ for (k = 0; k < n; k++)
+ a[i][j] += b[i][k]*c[k][j];
+ }
+#endif
+}
+
+void iter_newt(comb_sys s, double *y1, double **dy, double **u1, double z,
+ double *y2, double **u2, int *err)
+{
+ int i, j, k, l;
+ double *tmp1;
+ double **tmp2;
+
+ l = s->len;
+ /* u2 will be filled later on */
+ tmp2 = u2;
+ // tmp2 = A*U - U + 1
+ multm(tmp2, dy, u1, l);
+ for (i = 0; i < l; i++) {
+ for (j = 0; j < l; j++) {
+ tmp2[i][j] += - u1[i][j] + (i==j?1:0);
+ }
+ }
+ // dy = U*tmp2
+ multm(dy, u1, tmp2, l);
+ // u2 = u1 + dy
+ for (i = 0; i < l; i++) {
+ for (j = 0; j < l; j++) {
+ u2[i][j] = u1[i][j] + dy[i][j];
+ }
+ }
+ /* we don't use the values of u1 in the rest of the function */
+ tmp1 = u1[0];
+ // tmp1 = H(y1) - y1
+ for (i = 0; i < l; i++) {
+#ifdef COMPILE
+ tmp1[i] = (s->f[i])(y1, z, err) - y1[i];
+#else
+ set_sys(s);
+ tmp1[s->refs[i]] = eval_eq(s->eqs[i], y1, z, err) - y1[s->refs[i]];
+#endif
+ }
+ // y2 = u2*tmp1 + y1
+ for (i = 0; i < l; i++) {
+ y2[i] = y1[i];
+ for (k = 0; k < l; k++) {
+ y2[i] += u2[i][k]*tmp1[k];
+ }
+ }
+ // dy = H(y2)
+ for (i = 0; i < l; i++) {
+ for (j = 0; j < l; j++) {
+#ifdef COMPILE
+ dy[i][j] = (s->fj[i][j])(y2, z, err);
+#else
+ set_sys(s);
+ dy[i][j] = eval_eq(s->jacob[i][j], y2, z, err);
+#endif
+ }
+ }
+}
+
+void eval_sys_noalloc(comb_sys s, double z, double prec, double *y, double *w,
+ double **dy, double **u, double **v, int *err)
+{
+ // y va contenir le resultat final, w sert de stockage temporaire
+ int i, j, l;
+ double e;
+ double *y1, *y2;
+ double **u1, **u2;
+ void *tmp;
+
+ l = s->len;
+ y2 = y; y1 = w;
+ u2 = u; u1 = v;
+ // w = 0; y2 = 0
+ for (i = 0; i < l; i++) {
+ w[i] = 0.;
+ y2[i] = 0.;
+ }
+ l = s->len;
+ // dy = H(0); u1 = I
+ for (i = 0; i < l; i++) {
+ for (j = 0; j < l; j++) {
+#ifdef COMPILE
+ dy[i][j] = (s->fj[i][j])(y2, z, err);
+#else
+ set_sys(s);
+ dy[i][j] = eval_eq(s->jacob[i][j], y2, z, err);
+#endif
+ u1[i][j] = (i==j?1.:0.);
+ u2[i][j] = 0.;
+ }
+ }
+ do {
+ e = 0.;
+ iter_newt(s, y1, dy, u1, z, y2, u2, err);
+ if (*err)
+ break;
+ for (i = 0; i < l; i++) {
+ if (e < y2[i]-y1[i])
+ e = y2[i]-y1[i];
+ }
+ tmp = y1; y1 = y2; y2 = tmp;
+ tmp = u1; u1 = u2; u2 = tmp;
+ } while (e > prec);
+ if (y != y2)
+ for (i = 0; i < l; i++)
+ y[i] = y2[i];
+}
+
+void init_eval_sys(comb_sys s, double **w, double ***dy, double ***u, double ***v)
+{
+ int i, l;
+ double *p;
+
+ l = s->len;
+ jacob_sys(s);
+#ifdef COMPILE
+ compile_sys(s);
+#endif
+ *w = malloc(l*sizeof(double));
+ *dy = malloc(l*sizeof(double *));
+ p = malloc(l*l*sizeof(double));
+ for (i = 0; i < l; i++)
+ (*dy)[i] = &p[i*l];
+ *u = malloc(l*sizeof(double *));
+ p = malloc(l*l*sizeof(double));
+ for (i = 0; i < l; i++)
+ (*u)[i] = &p[i*l];
+ *v = malloc(l*sizeof(double *));
+ p = malloc(l*l*sizeof(double));
+ for (i = 0; i < l; i++)
+ (*v)[i] = &p[i*l];
+}
+
+void end_eval_sys(comb_sys s, double *w, double **dy, double **u, double **v)
+{
+ free(w);
+ free(dy[0]); free(u[0]); free(v[0]);
+ free(dy); free(u); free(v);
+}
+
+void eval_sys(comb_sys s, double z, double prec, int *err)
+{
+ double *w, **dy, **u, **v, *y;
+
+ y = s->y;
+ init_eval_sys(s, &w, &dy, &u, &v);
+ eval_sys_noalloc(s, z, prec, y, w, dy, u, v, err);
+ end_eval_sys(s, w, dy, u, v);
+}
+
+double sing_sys(comb_sys s, double prec1, double prec2)
+{
+ double zmin = 0.;
+ double zmax = 1.;
+ double z;
+ int i;
+ double *w, **dy, **u, **v, *y;
+ int err;
+
+ y = s->y;
+ init_eval_sys(s, &w, &dy, &u, &v);
+ do {
+ z = (zmin + zmax)/2.;
+ err = 0;
+ eval_sys_noalloc(s, z, prec2, y, w, dy, u, v, &err);
+ if (err)
+ zmax = z;
+ for (i = 0; z != zmax && i < sys_len(s); i++)
+ if (y[i] < 0 || y[i] > 1./prec1)
+ zmax = z;
+ if (zmax != z)
+ zmin = z;
+ } while (zmax-zmin > prec1);
+ eval_sys_noalloc(s, zmin, prec2, y, w, dy, u, v, &err);
+ end_eval_sys(s, w, dy, u, v);
+ return zmin;
+}
+
+/* TODO: mettre à jour */
+void eval_point_sys(comb_sys s, double z, double prec, int *err)
+{
+ comb_sys ss;
+ double *w, **dy, **u, **v, *y;
+
+ y = s->y;
+ ss = diff_sys_z(s);
+ init_eval_sys(ss, &w, &dy, &u, &v);
+ eval_sys_noalloc(ss, z, prec, y, w, dy, u, v, err);
+ end_eval_sys(ss, w, dy, u, v);
+ free_sys(ss);
+}
+
+/* TODO: mettre à jour */
+double mean_sys(comb_sys s, double mean, double prec1, double prec2)
+{
+ double zmin = 0.;
+ double zmax = 1.;
+ double *y;
+ double z, m;
+ int i;
+ comb_sys ss;
+ double *w, **dy, **u, **v;
+ int err;
+
+ y = s->y;
+ ss = diff_sys_z(s);
+ init_eval_sys(ss, &w, &dy, &u, &v);
+ do {
+ z = (zmin + zmax)/2.;
+ err = 0;
+ eval_sys_noalloc(ss, z, prec2, y, w, dy, u, v, &err);
+ m = z*y[s->len]/y[0];
+ if (err)
+ zmax = z;
+ for (i = 0; z != zmax && i < sys_len(s); i++)
+ if (y[i] < 0 || y[i] > 1./prec1 || m > mean)
+ zmax = z;
+ if (zmax != z)
+ zmin = z;
+ } while (zmax-zmin > prec1);
+ eval_sys_noalloc(ss, zmin, prec2, y, w, dy, u, v, &err);
+ end_eval_sys(ss, w, dy, u, v);
+ free_sys(ss);
+ return z;
+}
--- /dev/null
+#include <stdlib.h>
+#include "jacob.h"
+#include "equation.h"
+#include "system.h"
+#include "diff.h"
+
+void jacob_sys(comb_sys s)
+{
+ int i, j, l;
+ comb_eq *p;
+
+ if (s->jacob != NULL)
+ return;
+ l = s->len;
+ s->jacob = malloc(l*sizeof(comb_eq *));
+ p = malloc(l*l*sizeof(comb_eq));
+ for (i = 0; i < l; i++) {
+ s->jacob[s->refs[i]] = &p[s->refs[i]*l];
+ for (j = 0; j < l; j++) {
+ s->jacob[s->refs[i]][s->refs[j]] = diff_eq_ref(s->eqs[i], j+1);
+ }
+ }
+}
--- /dev/null
+#ifndef _JACOB_H
+# define _JACOB_H
+# include "combstruct.h"
+
+void jacob_sys(comb_sys);
+
+#endif
--- /dev/null
+#include <stdio.h>
+#include "print.h"
+#include "equation.h"
+#include "system.h"
+#include "jacob.h"
+
+void print_eq(comb_eq c)
+{
+ switch (c->type) {
+ case Nil:
+ printf("0");
+ break;
+ case Epsilon:
+ printf("Epsilon");
+ break;
+ case Atom:
+ printf("Z");
+ break;
+ case Ref:
+ printf("Y%d", c->ref);
+ break;
+ case Sum:
+ printf("Sum(");
+ print_eq(c->operand1);
+ printf(",");
+ print_eq(c->operand2);
+ printf(")");
+ break;
+ case Prod:
+ printf("Prod(");
+ print_eq(c->operand1);
+ printf(",");
+ print_eq(c->operand2);
+ printf(")");
+ break;
+ case Coef:
+ printf("%lf*", c->coef);
+ print_eq(c->operand1);
+ break;
+ case Seq:
+ printf("Seq(");
+ print_eq(c->operand1);
+ printf(")");
+ break;
+ default:
+ my_error("ERROR: print_eq unknown equation type\n");
+ }
+}
+
+void print_jacob(comb_sys s)
+{
+ int i, j;
+
+ for (i = 0; i < s->len; i++) {
+ for (j = 0; j < s->len; j++) {
+ printf("dY%d/dY%d = ", i+1, j+1);
+ print_eq(s->jacob[s->refs[i]][s->refs[j]]);
+ puts("");
+ }
+ }
+}
+
+void print_sys(comb_sys s)
+{
+ int i;
+
+ for (i = 0; i < s->len; i++) {
+ printf("Y%d = ", i+1);
+ print_eq(s->eqs[i]);
+ puts("");
+ }
+}
+
+void my_error(char *s)
+{
+ fprintf(stderr, s);
+}
--- /dev/null
+#ifndef _PRINT_H
+# define _PRINT_H
+# include "combstruct.h"
+# include "jacob.h"
+
+# include "common.h"
+
+void print_eq(comb_eq);
+
+void print_jacob(comb_sys);
+
+#endif
--- /dev/null
+#include <stdlib.h>
+#include "system.h"
+#include "equation.h"
+
+#define SYS_DEF_SIZE 128
+
+comb_sys sys_new(int size)
+{
+ comb_sys res;
+
+ res = malloc(sizeof(struct comb_sys));
+ res->len = 0;
+ res->eqs = malloc(size*sizeof(comb_eq));
+ res->refs = malloc(size*sizeof(int));
+ res->y = malloc(size*sizeof(double));
+ res->jacob = NULL;
+#ifdef COMPILE
+ res->f = NULL;
+ res->fj = NULL;
+#endif
+ return res;
+}
+
+void sys_add_eq(comb_sys s, int r, comb_eq c)
+{
+ s->eqs[s->len] = c;
+ s->refs[r-1] = s->len;
+ s->len++;
+}
+
+int sys_len(comb_sys s)
+{
+ return s->len;
+}
+
+double sys_val(comb_sys s, int i)
+{
+ if (i > s->len)
+ return 0.0;
+ return s->y[s->refs[i-1]];
+}
+
+void free_sys(comb_sys s)
+{
+ int i, j;
+
+ for (i = 0; i < s->len; i++)
+ free_eq(s->eqs[i]);
+ free(s->eqs);
+ free(s->refs);
+ free(s->y);
+ for (i = 0; i < s->len; i++)
+ for (j = 0; j < s->len; j++)
+ free_eq(s->jacob[i][j]);
+ free(s->jacob[0]);
+ free(s->jacob);
+#ifdef COMPILE
+ free(s->f);
+ free(s->fj[0]);
+ free(s->fj);
+#endif
+ free(s);
+}
--- /dev/null
+#ifndef _SYSTEM_H
+# define _SYSTEM_H
+# include "combstruct.h"
+# ifdef COMPILE
+
+typedef double (*f_eval)(double *, double, int *);
+
+void compile_sys(comb_sys);
+# endif
+
+struct comb_sys {
+ int len;
+ comb_eq *eqs;
+ int *refs;
+ double *y;
+ comb_eq **jacob;
+# ifdef COMPILE
+ f_eval *f;
+ f_eval **fj;
+# endif
+};
+
+#endif
--- /dev/null
+CFLAGS=-Wall -g
+LDFLAGS=-L../src -lm -Wl,-rpath ../src
+CC=gcc
+
+all: test
+
+%: %.c ../src/combstruct.h
+ $(CC) $(CFLAGS) $(LDFLAGS) -lcombstruct -o $@ $<
+
+%_compile_blas: %.c ../src/combstruct.h
+ $(CC) $(CFLAGS) $(LDFLAGS) -lcombstruct_compile_blas -o $@ $<
+
+%_blas: %.c ../src/combstruct.h
+ $(CC) $(CFLAGS) $(LDFLAGS) -lcombstruct_blas -o $@ $<
+
+%_compile: %.c ../src/combstruct.h
+ $(CC) $(CFLAGS) $(LDFLAGS) -lcombstruct_compile -o $@ $<
+
+clean:
+ rm -f test *.o
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include "../src/combstruct.h"
+
+int main(int argc, char **argv)
+{
+ comb_sys s;
+ double *y;
+ double z = 0.5;
+ int err = 0;
+
+ s = empty_sys();
+ // Y = 1/(1 - Z) + 1/(1 - 4*Z)
+ add_sys(s, sum(seq(atom()), seq(prod(sum(sum(epsilon(),epsilon()),sum(epsilon(),epsilon())), atom()))));
+ y = malloc(2*sys_len(s)*sizeof(double));
+ z = sing_sys(s, 10e-12, 10e-5, y);
+ printf("z:%.20lf y:%.20lf\n", z, y[0]);
+ eval_point_sys(s, 0.5, 10e-5, y, &err);
+ printf("z:0.5 y:%.20lf y':%.20lf err:%d\n", y[0], y[1], err);
+ z = mean_sys(s, 500, 10e-12, 10e-5, y);
+ printf("z:%.20lf y:%.20lf\n", z, y[0]);
+ free(y);
+ free_sys(s);
+ return 0;
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include "../src/combstruct.h"
+
+int main(int argc, char **argv)
+{
+ comb_sys s;
+ double z = 4./27.;
+
+ s = sys_new(13);
+ // 1: package_uml
+ sys_add_eq(s, 1, coef(0.01577, prod(atom(), seq(ref(2)))));
+ // 2: packageableElement
+ sys_add_eq(s, 2, sum(ref(1), sum(ref(3), ref(6))));
+ // 3: class_uml
+ sys_add_eq(s, 3, prod(atom(), prod(seq(ref(5)), prod(seq(ref(12)), seq(ref(4)))))),
+ // 4: generalization
+ sys_add_eq(s, 4, atom());
+ // 5: property
+ sys_add_eq(s, 5, prod(atom(), coef(3.0, sum(epsilon(), ref(7)))));
+ // 6: association
+ sys_add_eq(s, 6, atom());
+ // 7: valueSpecification
+ sys_add_eq(s, 7, sum(ref(8), sum(ref(9), sum(ref(10), ref(11)))));
+ // 8: literalBoolean
+ sys_add_eq(s, 8, atom());
+ // 9: literalNull
+ sys_add_eq(s, 9, atom());
+ // 10: literalInteger
+ sys_add_eq(s, 10, atom());
+ // 11: literalString
+ sys_add_eq(s, 11, atom());
+ // 12: operation
+ sys_add_eq(s, 12, prod(atom(), coef(2.0, seq(ref(13)))));
+ // 13: parameter
+ sys_add_eq(s, 13, prod(atom(), coef(3.0, sum(epsilon(), ref(7)))));
+ z = sing_sys(s, 10e-12, 10e-5);
+ printf("=== SING ===\n");
+ printf("#define X %.20lf\n", z);
+ printf("#define PACKAGE_UML %.10lf\n", sys_val(s, 1));
+ printf("#define PACKAGEABLE_ELEMENT %.10lf\n", sys_val(s, 2));
+ printf("#define CLASS_UML %.10lf\n", sys_val(s, 3));
+ printf("#define GENERALIZATION %.10lf\n", sys_val(s, 4));
+ printf("#define PROPERTY %.10lf\n", sys_val(s, 5));
+ printf("#define ASSOCIATION %.10lf\n", sys_val(s, 6));
+ printf("#define VALUE_SPECIFICATION %.10lf\n", sys_val(s, 7));
+ printf("#define LITERAL_BOOLEAN %.10lf\n", sys_val(s, 8));
+ printf("#define LITERAL_NULL %.10lf\n", sys_val(s, 9));
+ printf("#define LITERAL_INTEGER %.10lf\n", sys_val(s, 10));
+ printf("#define LITERAL_STRING %.10lf\n", sys_val(s, 11));
+ printf("#define OPERATION %.10lf\n", sys_val(s, 12));
+ printf("#define PARAMETER %.10lf\n", sys_val(s, 13));
+/*
+ y = realloc(y, sys_len(s)*2*sizeof(double));
+ z = mean_sys(s, 30000, 10e-12, 10e-5, y);
+ printf("=== MEAN ===\n");
+ printf("#define X %.20lf\n", z);
+ printf("#define PACKAGE_UML %.10lf\n", y[0]);
+ printf("#define PACKAGEABLE_ELEMENT %.10lf\n", y[1]);
+ printf("#define CLASS_UML %.10lf\n", y[2]);
+ printf("#define GENERALIZATION %.10lf\n", y[3]);
+ printf("#define PROPERTY %.10lf\n", y[4]);
+ printf("#define ASSOCIATION %.10lf\n", y[5]);
+ printf("#define VALUE_SPECIFICATION %.10lf\n", y[6]);
+ printf("#define LITERAL_BOOLEAN %.10lf\n", y[7]);
+ printf("#define LITERAL_NULL %.10lf\n", y[8]);
+ printf("#define LITERAL_INTEGER %.10lf\n", y[9]);
+ printf("#define LITERAL_STRING %.10lf\n", y[10]);
+ printf("#define OPERATION %.10lf\n", y[11]);
+ printf("#define PARAMETER %.10lf\n", y[12]);
+ eval_point_sys(s, z, 10e-3, y, &err);
+ printf("Moyennes:\n");
+ printf("PACKAGE_UML %.10lf\n", z*y[13]/y[0]);
+ printf("PACKAGEABLE_ELEMENT %.10lf\n", z*y[14]/y[1]);
+ printf("CLASS_UML %.10lf\n", z*y[15]/y[2]);
+ printf("GENERALIZATION %.10lf\n", z*y[16]/y[3]);
+ printf("PROPERTY %.10lf\n", z*y[17]/y[4]);
+ printf("ASSOCIATION %.10lf\n", z*y[18]/y[5]);
+ printf("VALUE_SPECIFICATION %.10lf\n", z*y[19]/y[6]);
+ printf("LITERAL_BOOLEAN %.10lf\n", z*y[20]/y[7]);
+ printf("LITERAL_NULL %.10lf\n", z*y[21]/y[8]);
+ printf("LITERAL_INTEGER %.10lf\n", z*y[22]/y[9]);
+ printf("LITERAL_STRING %.10lf\n", z*y[23]/y[10]);
+ printf("OPERATION %.10lf\n", z*y[24]/y[11]);
+ printf("PARAMETER %.10lf\n", z*y[25]/y[12]); */
+ free_sys(s);
+ return 0;
+}
projects / boltzoc / commitdiff