double PronyModel(double t, matrix* beta, void *params)
{
    matrix *Ji, *taui;
    double *p;
    double J0, J, Jval, tauval;
    int n, i;
    p = (double*) params;
    J0 = *p;

    /* Pull out all the parameter values */
    n = nRows(beta)/2;
    Ji = CreateMatrix(n, 1);
    taui = CreateMatrix(n, 1);

    for(i=0; i<n; i++) {
        setval(Ji, val(beta, 2*i, 0), i, 0);
        setval(taui, val(beta, 2*i+1, 0), i, 0);
    }

    J = J0;
    for(i=0; i<n; i++) {
        Jval = val(Ji, i, 0)*val(Ji, i, 0);
        tauval = val(taui, i, 0)*val(taui, i, 0);
        J += Jval * (1-exp(-t/tauval));
    }
    DestroyMatrix(Ji);
    DestroyMatrix(taui);
    return J;
}

static unsigned long aout_addsym(char *name,taddr value,int bind,
                                 int info,int type,int desc,int be)
/* add a new symbol, return its symbol table index */
{
  struct SymbolNode **chain = &aoutsymlist.hashtab[hashcode(name)%SYMHTABSIZE];
  struct SymbolNode *sym;

  while (sym = *chain)
    chain = &sym->hashchain;
  /* new symbol table entry */
  *chain = sym = mycalloc(sizeof(struct SymbolNode));

  if (!name)
    name = emptystr;
  sym->name = name;
  sym->index = aoutsymlist.nextindex++;
  setval(be,sym->s.n_strx,4,aout_addstr(name));
  sym->s.n_type = type;
  /* GNU binutils don't use BIND_LOCAL/GLOBAL in a.out files! We do! */
  sym->s.n_other = ((bind&0xf)<<4) | (info&0xf);
  setval(be,sym->s.n_desc,2,desc);
  setval(be,sym->s.n_value,4,value);
  addtail(&aoutsymlist.l,&sym->n);
  return sym->index;
}

void _start(void)
{
    struct uart *const uart = get_base_uart(0);
    struct dma *const dma =((struct dma *)get_base_dmas()) + 0;

    /*开启uart0 的DMA received 功能*/
    //uart->UCON |= 0x2 << 2;
    uart->UCON = 9;

    /* 设置源 */
    dma->DISRC = (unsigned int)buf;
    dma->DISRCC = 0; /* 数据源在内存中，ahb */

    /* 设置目的地址 */
    dma->DIDST =(unsigned int)&uart->UTXH;
    dma->DIDSTC |= 3;  /* 目的地址uart连接在apb上*/

    setval(dma->DCON, 0x001, 3, 24); /* 接上uart0  */
    setval(dma->DCON, 0x1, 1, 23); /*  硬件DMA*/
    setval(dma->DCON, 0x1, 1, 22); /* 关闭重载 */
    setval(dma->DCON, sizeof(buf), 20, 0);

    /* 激活DMA  */
    dma->DMASKTRIG |= 0x2;

    return;
}

/* Of a 209-line interpreter, the actual interpreting function is only 17 lines... */
void run() {
    // run the program
    for (ip = 0; ip < n_of_lines; ip++) {
        switch (pgm[ip].instruction) {
        case LBL:
            break; // don't do anything, labels have already been stored.
        case CPY:
            setval(value(pgm[ip].op2),value(pgm[ip].op1));
            break;
        case INC:
            setval(value(pgm[ip].op1),memval(value(pgm[ip].op1))+1);
            break;
        case DEC:
            if (memval(value(pgm[ip].op1))) {
                setval(value(pgm[ip].op1), memval(value(pgm[ip].op1))-1);
            } else {
                ip = getlbl(pgm[ip].op2);
            };
            break;
        case OUT:
            putchar(value(pgm[ip].op1));
            break;
        case INP:
            setval(value(pgm[ip].op1),feof(stdin)?0:getchar());
            break;
        }
    }
}

matrix* CreateElementMatrix(struct fe *p, Elem2D *elem, matrix *guess)
{
    basis *b;
    b = p->b;

    int nvars = p->nvars;

    int i, j;
    double value = 0;
    matrix *m;

    m = CreateMatrix(b->n*nvars, b->n*nvars);

    for(i=0; i<b->n*nvars; i+=nvars) {
        for(j=0; j<b->n*nvars; j+=nvars) {
            /* dRx/dx */
            value = quad2d3generic(p, guess, elem, &ElemJdRxdx, i/2, j/2);
            setval(m, value, i, j);

            /* dRy/dx */
            value = quad2d3generic(p, guess, elem, &ElemJdRydx, i/2, j/2);
            setval(m, value, i+1, j);

            /* dRx/dy */
            value = quad2d3generic(p, guess, elem, &ElemJdRxdy, i/2, j/2);
            setval(m, value, i, j+1);

            /* dRy/dy */
            value = quad2d3generic(p, guess, elem, &ElemJdRydy, i/2, j/2);
            setval(m, value, i+1, j+1);
        }
    }

    return m;
}

void Filter::sliderValueChanged (Slider* sliderThatWasMoved)
{
    //[UsersliderValueChanged_Pre]
    //[/UsersliderValueChanged_Pre]

    if (sliderThatWasMoved == lower_cutoffslider)
    {
        //[UserSliderCode_lower_cutoffslider] -- add your slider handling code here..
      setval(filter,lower_cutoff);
        //[/UserSliderCode_lower_cutoffslider]
    }
    else if (sliderThatWasMoved == upper_cutoffslider)
    {
        //[UserSliderCode_upper_cutoffslider] -- add your slider handling code here..
      setval(filter,upper_cutoff);
        //[/UserSliderCode_upper_cutoffslider]
    }
    else if (sliderThatWasMoved == sharpnessslider)
    {
        //[UserSliderCode_sharpnessslider] -- add your slider handling code here..
      setval(filter,sharpness);
        //[/UserSliderCode_sharpnessslider]
    }

    //[UsersliderValueChanged_Post]
    //[/UsersliderValueChanged_Post]
}

int main(int argc, char *argv[])
{
    int i, j, k;
    double Ti, Mj, percent;
    vector *T, *M;
    matrix *t, *Jij, *betaij, *output, *ttmp;

    /*
    if(argc < 3) {
        printf("Usage:\n"
               "fitcreep: <file> <t1> <t2> ... <tn>\n"
               "<file>: Filename containing the creep function data.\n"
               "<t1>: First retardation time\n"
               "<t2>: Second retardation time\n"
               "...\n"
               "<tn>: Nth retardation time.\n");
        exit(0);
    }
    */

    T = linspaceV(333, 363, 10);
    M = linspaceV(.05, .4, 10);

    ttmp = linspace(1e-3, 1e3, 1000);
    t = mtxtrn(ttmp);
    DestroyMatrix(ttmp);

    output = CreateMatrix(len(T)*len(M), 2+5);

    for(i=0; i<len(T); i++) {
        Ti = valV(T, i);
        for(j=0; j<len(M); j++) {
            Mj = valV(M, j);
            Jij = makedata(t, Ti, Mj);
            betaij = fitdata(t, Jij);

            setval(output, Ti, i*len(M)+j, 0);
            setval(output, Mj, i*len(M)+j, 1);
            setval(output, val(Jij, 0, 0), i*len(M)+j, 2);
            for(k=0; k<nRows(betaij); k++)
                setval(output, pow(val(betaij, k, 0), 2), i*len(T)+j, k+3);
            DestroyMatrix(Jij);
            DestroyMatrix(betaij);

            /* Print the percent done */
            percent = (1.*i*len(M)+j)/(len(M)*len(T))*100.;
            printf("%3.2f %%\r", percent);
            fflush(stdout);
        }
    }
    
    DestroyMatrix(t);
    DestroyVector(T);
    DestroyVector(M);
    mtxprntfilehdr(output, "output.csv", "T,M,J0,J1,tau1,J2,tau2\n");
    DestroyMatrix(output);
    return 0;
}

/* Create the load vector */
matrix* CreateElementLoad(double Pe, double h) {
    matrix *f;

    f = CreateMatrix(2, 1);

    setval(f, 0, 0, 0);
    setval(f, 0, 1, 0);

    return f;
}

matrix* testload(double Pe, double h) {
    matrix *f;

    f = CreateMatrix(2, 1);

    setval(f, h/2, 0, 0);
    setval(f, h/2, 1, 0);

    return f;
}

static void addRel64(int be,taddr o,taddr a,taddr i)
{
  struct RelocNode *rn = mymalloc(sizeof(struct RelocNode));

  setval(be,rn->r.r_offset,8,o);
#if RELA
  setval(be,rn->r.r_addend,8,a);
#endif
  setval(be,rn->r.r_info,8,i);
  addtail(&relalist,&(rn->n));
}

/* Create the element matrix for the specified values of Pe and h */
matrix* CreateElementMatrix(double Pe, double h)
{
    matrix *elem;

    elem = CreateMatrix(2, 2);
    setval(elem, 1/h-Pe/2, 0, 0);
    setval(elem, -1/h+Pe/2, 0, 1);
    setval(elem, -1/h-Pe/2, 1, 0);
    setval(elem, 1/h+Pe/2, 1, 1);

    return elem;
}

/* Test functions to create the element and load matricies based on the equation
 * in Problem 1
 */
matrix* testelem(double Pe, double h)
{
    matrix *elem;

    elem = CreateMatrix(2, 2);
    setval(elem, -1/h+h/3, 0, 0);
    setval(elem, 1/h+h/6, 0, 1);
    setval(elem, 1/h+h/6, 1, 0);
    setval(elem, -1/h+h/3, 1, 1);

    return elem;
}

matrix* GetDeformedCoords(struct fe *p, matrix *Soln)
{
    matrix *Def;
    int i;
    Def = CreateMatrix(nRows(Soln)/2, 2);

    for(i=0; i<nRows(Soln)/2; i++) {
        setval(Def, val(Soln, 2*i, 0) + valV(GetNodeCoordinates(p->mesh, i), 0),  i, 0);
        setval(Def, val(Soln, 2*i+1, 0) + valV(GetNodeCoordinates(p->mesh, i), 1), i, 1);
    }

    return Def;
}

matrix* FormatDisplacements(struct fe *p, matrix *Soln)
{
    matrix *Def;
    int i;
    Def = CreateMatrix(nRows(Soln)/2, 2);

    for(i=0; i<nRows(Soln)/2; i++) {
        setval(Def, val(Soln, 2*i, 0),  i, 0);
        setval(Def, val(Soln, 2*i+1, 0), i, 1);
    }

    return Def;
}

/**
 * Fit the GAB parameters given water activity
 */
int main(int argc, char *argv[])
{
    matrix *data, *aw, *Xdb, *tmp0, *tmp1, *beta0, *beta;

    if(argc != 2) {
        puts("Usage:");
        puts("gab <aw.csv>");
    }
    //data = mtxloadcsv("Andrieu.csv", 0);
    data = mtxloadcsv(argv[1], 0);

    /* Get the water activity from column 1 and the moisture content from
     * column 6. */
    aw = ExtractColumn(data, 0);
    Xdb = ExtractColumn(data, 5);

    /* Stick the two matricies together and delete any rows that contain
     * empty values */
    tmp0 = AugmentMatrix(aw, Xdb);
    tmp1 = DeleteNaNRows(tmp0);
    DestroyMatrix(aw);
    DestroyMatrix(Xdb);

    /* Pull the two columns back apart */
    aw = ExtractColumn(tmp1, 0);
    Xdb = ExtractColumn(tmp1, 1);
    DestroyMatrix(tmp0);
    DestroyMatrix(tmp1);

    /* Set up the beta matrix with some initial guesses at the GAB constants.
     * The solver needs these to be pretty close to the actual values, or it
     * will fail to converge */
    beta0 = CreateOnesMatrix(3, 1);
    setval(beta0, 6, 0, 0);
    setval(beta0, .5, 1, 0);
    setval(beta0, .04, 2, 0);

    /* Attempt to fit the gab parameters to the supplied data */
    beta = fitnlm(&gab, aw, Xdb, beta0);

    /* Print out the fitted values */
    printf("C = %g\nk = %g\nXm = %g\n",
            val(beta, 0, 0),
            val(beta, 1, 0),
            val(beta, 2, 0));

    return 0;
}

struct symbol*
setsymbol(const char *symname, struct value *val)
{
	struct symbol *sym;

	/* search for a matching symbol name */
	sym = findsymbol(symname);

	/* if there isn't one, add a new to the end of the list */
	if(!sym){
		sym = getlastsymbol();
		sym->sym_next = xmalloc(sizeof(struct symbol));
		sym = sym->sym_next;

		/* set the name, ensure list integrity and move along as
		 * if nothing happened */
		sym->name = allocstring(symname);
		sym->sym_next = NULL;
	}
	/* if a symbol existed, clear the old value */
	else freeval(&sym->val);

	setval(&sym->val, val);
	return sym;
}

/**
 * @brief Smash several vectors together as column vectors into a matrix.
 *
 * @param n Number of vectors to smash.
 * @returns An n by l matrix containing the smashed vectors
 */
matrix* CatColVector(int n, ...)
{
    int i, j;
    va_list args;
    vector** v;
    matrix* result;
    int lmax = 0; /* Maximum vector length */
    v = (vector**) calloc(sizeof(vector*), n);

    va_start(args, n); /* Start getting the vectors */
    for(i=0; i<n; i++) {
        v[i] = va_arg(args, vector*);
        if(lmax < len(v[i]))
            lmax = len(v[i]);
    }
    va_end(args); /* Over it! */

    /* Create a matrix large enough to hold all the elements of the longest
     * vector. */
    result = CreateMatrix(lmax, n);

    for(i=0; i<n; i++)
        for(j=0; j<len(v[i]); j++)
            setval(result, valV(v[i], j), j, i);

    /* Cleanup */
    free(v);

    return result;
}

int
qrecv(Lextok *n, int full)
{	int whichq = eval(n->lft)-1;

	if (whichq == -1)
	{	if (n->sym && !strcmp(n->sym->name, "STDIN"))
		{	Lextok *m;

			if (TstOnly) return 1;

			for (m = n->rgt; m; m = m->rgt)
			if (m->lft->ntyp != CONST && m->lft->ntyp != EVAL)
			{	int c = getchar();
				(void) setval(m->lft, c);
			} else
				fatal("invalid use of STDIN", (char *)0);

			whichq = 0;
			return 1;
		}
		printf("Error: receiving from an uninitialized chan %s\n",
			n->sym?n->sym->name:"");
		whichq = 0;
		return 0;
	}
	if (whichq < MAXQ && whichq >= 0 && ltab[whichq])
	{	ltab[whichq]->setat = depth;
		return a_rcv(ltab[whichq], n, full);
	}
	return 0;
}

static int 
countComponents(Agraph_t *g, int* max_degree, float *nontree_frac)
{
  int        nc = 0;
  int        sum_edges = 0;
  int        sum_nontree = 0;
  int        deg;
  int        n_edges;
  int        n_nodes;
  Agnode_t*  n;

  for (n = agfstnode(g); n; n = agnxtnode(n)) {
    if (!getval(n)) {
      nc++;
      n_edges = 0;
      n_nodes = label(n,0,&n_edges);
      sum_edges += n_edges;
      sum_nontree += (n_edges - n_nodes + 1);
    }
  }
  if (max_degree) {
    int maxd = 0;
    for (n = agfstnode(g); n; n = agnxtnode(n)) {
      deg = agdegree(n,TRUE,TRUE);
      if (maxd < deg) maxd = deg;
      setval(n,0);
    }
    *max_degree = maxd;
  }
  if (nontree_frac) {
    if (sum_edges > 0) *nontree_frac = (float)sum_nontree / (float)sum_edges;
    else *nontree_frac = 0.0;
  }
  return nc;
}

/*
 * Generic object newindex metamethod handler.
 * Lua Stack: ud, key, value
 */
int object_newindex_handler(lua_State *L)
{
  char sskey[16], *skey;
  if ((!lua_isstring(L, 2)) || (!lua_isstring(L, 3))){
    gr_object_t *ud = toobject(L, 1, NULL, STRICT);
    skey = agget(ud->p.p, ".attrib");
    if (!skey || (strlen(skey) == 0)){
      /* Let's create an attrib table on the fly if none exists */
      sprintf(sskey, "%p", ud->p.p);
      skey = agstrdup(sskey);
      agset(ud->p.p, ".attrib", skey);
      lua_pushstring(L, skey);             /* ud, key, value, skey */
      lua_newtable(L);                     /* ud, key, value, skey, stab */
      lua_rawset(L, LUA_REGISTRYINDEX);    /* ud, key, value, */
    }
    lua_pushstring(L, skey);          /* ud, key, value, skey */
    lua_rawget(L, LUA_REGISTRYINDEX); /* ud, key, value, stab */
    lua_pushvalue(L, 2);              /* ud, key, value, stab, key */
    lua_pushvalue(L, 3);              /* ud, key, value, stab, key, value */
    lua_rawset(L, -3);                /* ud, key, value, stab */
    lua_pop(L, -1);                   /* ud, key, value */
    return 0;
  }
  return setval(L);
}