mfcalc.y

Linux/windows 环境下跨平台开发程序

	long npoints = 0;
	int length, res_mode = 0;
	unsigned int hn_x = HashValue((unsigned char *)"x");
	unsigned int hn_y = HashValue((unsigned char *)"y");

	if(x1 < x2) step = fabs(step);
	else step = -fabs(step);
	if(!(new_points = (lfPOINT*)calloc((iround(fabs(x2-x1)/fabs(step))+2), sizeof(lfPOINT))))
		return false;
	if(d) curr_data = d;
	push_parser();
	init_table();
	if(param) {
		length = strlen(param);
		if(!(buffer = (char*)malloc(length+2))){
			pop_parser();
			return false;
			}
		strcpy(buffer, param);	buffer[length++] = ';';
		buffer[length] = 0;	buff_pos = 0;
		do {
			yyparse();
			}while(buff_pos < length);
		free(buffer);		buffer = 0L;
		}		
	length = strlen(expr);
	buffer = expr;		sx = putsym(hn_x, VAR);
	for(x = x1; step > 0.0 ? x <= x2 : x >= x2; x += step) {
		if(sx){
			sx->SetValue(x);	buff_pos = 0;
			do {
				yyparse();
				}while(buff_pos < length);
			switch (res_mode) {
			case 1:
				y = sy->GetValue();	break;
			case 2:
				y = line_res.value;	break;
			default:
				if(sy = getsym(hn_y)) {
					y = sy->GetValue();	res_mode = 1;
					}
				else {
					y = line_res.value;	res_mode = 2;
					}
				break;
				}
			new_points[npoints].fx = (getsym(hn_x))->GetValue();
			new_points[npoints++].fy = y;
			}
		}
	*pts = new_points;	*npts = npoints;
	clear_table();
	pop_parser();
	if(curr_data) {
		curr_data->Command(CMD_CLEAR_ERROR, 0L, 0L);
		curr_data->Command(CMD_REDRAW, 0L, 0L);
		}
	return true;
}

bool do_func3D(DataObj *d, double x1, double x2, double xstep, double z1, double z2, double zstep, 
	char *expr, char *param)
{
	int length, nr, nc, r, c, res_mode=0;
	symrec *sx, *sz, *sy;
	double x, y, z;
	unsigned int hn_x = HashValue((unsigned char *)"x");
	unsigned int hn_y = HashValue((unsigned char *)"y");
	unsigned int hn_z = HashValue((unsigned char *)"z");

	if(!d || x2 <= x1 || z2 <= z1 || xstep <= 0.0 || zstep <= 0.0) return false;
	push_parser();
	init_table();
	if(param) {
		length = strlen(param);
		if(!(buffer = (char*)malloc(length+2))){
			pop_parser();
			return false;
			}
		strcpy(buffer, param);	buffer[length++] = ';';
		buffer[length] = 0;	buff_pos = 0;
		do {
			yyparse();
			}while(buff_pos < length);
		free(buffer);		buffer = 0L;
		}		
	length = strlen(expr);		buffer = expr;
	sx = putsym(hn_x, VAR);		sz = putsym(hn_z, VAR);
	nr = iround((z2-z1)/zstep)+1;	nc = iround((x2-x1)/xstep)+1;
	d->Init(nr, nc);
	for(r = 0, x = x1; r < nr; r++, x += xstep) {
		for(c = 0, z = z1; c < nc; c++, z+= zstep) {
			sx->SetValue(x);	sz->SetValue(z);	buff_pos = 0;
			do {
				yyparse();
				}while(buff_pos < length);
			switch (res_mode) {
			case 1:
				y = sy->GetValue();	break;
			case 2:
				y = line_res.value;	break;
			default:
				if(sy = getsym(hn_y)) {
					y = sy->GetValue();	res_mode = 1;
					}
				else {
					y = line_res.value;	res_mode = 2;
					}
				break;
				}
			d->SetValue(r, c, y);
			}
		} 
	clear_table();
	pop_parser();
	return true;
}

anyResult *do_formula(DataObj *d, char *expr)
{
	int length;
	static anyResult ret, *pret = 0L;

	if(d) curr_data = d;
	ret.type = ET_ERROR;		ret.text = 0L;
	if(!expr || !expr[0]) return &ret;
	push_parser();		//make code reentrant
	init_table();		length = strlen(expr);
	if(!(buffer = (char*)malloc(length+2))){
		pop_parser();
		return &ret;
		}
	strcpy(buffer, expr);	if(buffer[length-1] != ';') buffer[length++] = ';';
	buffer[length] = 0;	buff_pos = 0;
	do {
		yyparse();
		}while(buff_pos < length);
	ret.type = ET_ERROR;		ret.text = 0L;
	if(curr_data && last_error) {
		if(!(strcmp(last_error, "parse error"))) curr_data->Command(CMD_ERROR, 0L, 0L);
		if(last_err_desc) pret = &line_res;
		else pret = &ret;
		}
	else pret = &line_res;
	last_error = last_err_desc = 0L;
	free(buffer);		buffer = 0L;
	clear_table();
	pop_parser();
	return pret;
}

bool MoveFormula(DataObj *d, char *of, char *nf, int dx, int dy, int r0, int c0)
{
	int length, tok, pos, i;
	char *res, desc1[2], desc2[2];

	if(d) curr_data = d;
	if(!curr_data || !of || !nf) return false;
	push_parser();		//make code reentrant
	init_table();		length = strlen(of);
	if(!(buffer = (char*)malloc(length+2))){
		pop_parser();
		return false;
		}
	strcpy(buffer, of);	buffer[length++] = ';';
	buffer[length] = 0;	buff_pos = pos = 0;
	res = (char *)calloc(length*2+10, sizeof(char));
	do {
		tok = yylex ();
		if(tok && tok < 256) {
			if(res[pos-1] == ' ') pos--;
			res[pos++] = (char)tok;
			}
		else switch(tok) {
			case NUM:
				pos += sprintf(res+pos, "%g", yylval.val);
				break;
			case FNCT:	case FUNC2:	case AFNCT:	case SFNCT:
				pos += sprintf(res+pos, "%s", curr_sym->name);
				break;
			case COLR:			case COLC:
				pos += sprintf(res+pos, ":");
				break;
			case PSEP:
				pos += sprintf(res+pos, ";");
				break;
			case CLVAL:
				pos += sprintf(res+pos, "$$");
				break;
			case CLAUSE:
				pos += sprintf(res+pos, " where ");
				break;
			case VAR:
				curr_sym->InitSS();
				if(curr_sym->col >= 0 && curr_sym->row >= 0) {
					desc1[0] = desc1[1] = desc2[0] = desc2[1] = 0;
					for(i=strlen(curr_sym->name)-1; i>0 && isdigit(curr_sym->name[i]); i--);
					if(curr_sym->name[0] == '$') desc1[0] = '$';
					if(curr_sym->name[i] == '$') desc2[0] = '$';
					pos += sprintf(res+pos, "%s%s%s%d", desc1, 
						Int2ColLabel(desc1[0] || curr_sym->col < c0 ? curr_sym->col : curr_sym->col+dx >=0 ?
						curr_sym->col+dx > c0 ? curr_sym->col+dx : c0 : 0, false),
						desc2, desc2[0] || curr_sym->row < r0 ? curr_sym->row+1 : curr_sym->row + dy >= 0 ? 
						curr_sym->row+dy > r0 ? curr_sym->row+1+dy : r0 : 1);
					}
				else pos += sprintf(res+pos, "%s ", curr_sym->name);
				break;
			case STR:
				pos += sprintf(res+pos, "\"%s\"", yylval.text && yylval.text[0] ? yylval.text : "");
				break;
			case SER:
				pos += sprintf(res+pos, "..");
				break;
			case INC:
				pos += sprintf(res+pos, "++");
				break;
			case DEC:
				pos += sprintf(res+pos, "--");
				break;
			case PI:
				pos += sprintf(res+pos, "pi");
				break;
			case E:
				pos += sprintf(res+pos, "e");
				break;
			case AND:
				pos += sprintf(res+pos, " && ");
				break;
			case OR:
				pos += sprintf(res+pos, " || ");
				break;
			case EQ:
				pos += sprintf(res+pos, " == ");
				break;
			case NE:
				pos += sprintf(res+pos, " != ");
				break;
			case GT:
				pos += sprintf(res+pos, ">");
				break;
			case GE:
				pos += sprintf(res+pos, ">=");
				break;
			case LT:
				pos += sprintf(res+pos, "<");
				break;
			case LE:
				pos += sprintf(res+pos, "<=");
				break;
			case IF: 
				pos += sprintf(res+pos, "if");
				break;
			case ELSE: 
				pos += sprintf(res+pos, "else");
				break;
			case BLOCK:
				pos += sprintf(res+pos, "{%s}", yylval.text && yylval.text[0] ? yylval.text : "");
				break;
			}
		}while(buff_pos < length);
	while((res[pos-1] == ';' || res[pos-1] == ' ') && pos > 0) { res[pos-1] = 0; pos--;} 
	strcpy(nf, res);	free(res);
	free(buffer);		buffer = 0L;
	clear_table();
	pop_parser();
	return true;
}

static char *txt_formula;	//function to fit
static double **parval;		//pointers to parameter values
static void fcurve(double x, double z, double **a, double *y, double dyda[], int ma)
{
	int i, length;
	double tmp, y1, y2;
	symrec *symx, *sy=0L;
	unsigned int hn_x = HashValue((unsigned char *)"x");
	unsigned int hn_y = HashValue((unsigned char *)"y");

	if(!(symx = getsym(hn_x))) symx = putsym(hn_x, VAR);
	//swap parameters to requested set
	if(a != parval) for(i = 0; i < ma; i++) {
		tmp = *parval[i];	*parval[i]  = *a[i];	*a[i] = tmp;
		}
	//calc result
	symx->SetValue(x);	buffer = txt_formula;
	buff_pos = 0;		length = strlen(txt_formula);
	do {	yyparse();	}while(buff_pos < length);
	if(sy = getsym(hn_y)) *y = sy->GetValue();
	else *y = line_res.value;
	if(*y == HUGE_VAL || *y == -HUGE_VAL) {
		for(i = 0, *y = 0.0; i < ma; dyda[i++] = 0.0);
		return;
		}
	//partial derivatives for each parameter by numerical differentiation
	for(i = 0; i < ma; i++) {
		if(*parval[i] != 0.0) {
			tmp = *parval[i];
			*parval[i] = tmp*.995;
			buff_pos = 0;
			do {	yyparse();	}while(buff_pos < length);
			y1 = sy ? sy->GetValue() : line_res.value;
			*parval[i] = tmp*1.005;
			buff_pos = 0;
			do {	yyparse();	}while(buff_pos < length);
			y2 = sy ? sy->GetValue() : line_res.value;
			*parval[i] = tmp;
			dyda[i] = (y2-y1)*100.0/tmp;
			}
		else dyda[i] = 0.0;
		}
	//swap parameters back to original
	if(a != parval) for(i = 0; i < ma; i++) {
		tmp = *parval[i];	*parval[i]  = *a[i];	*a[i] = tmp;
		}
}

int do_fitfunc(DataObj *d, char *rx, char *ry, char *rz, char **par, char *expr, double conv, int maxiter, double *chi_2)
{
	int length, i, j, k, l, ndata, nparam, r1, r2, r3, c1, c2, c3, *lista, itst, itst1;
	symrec *tab1, *tab2, *csr, **parsym;
	AccRange *arx, *ary, *arz;
	double *x, *y, *z, currx, curry, currz, alamda, chisq, ochisq;
	double **covar, **alpha;
	char tmp_txt[500];

	if(d) curr_data = d;
	if(chi_2) *chi_2 = 0.0;
	txt_formula = expr;
	if(!curr_data || !par || !expr || !rx || !ry) return 0;
	//process ranges and create arrays
	arx = ary = arz = 0L;	x = y = z = 0L;	parval = 0L;	parsym = 0L;
	if(!(arx = new AccRange(rx)))return 0;
	i = arx->CountItems()+1;
	if(!(ary = new AccRange(ry))){
		delete arx;	return 0;
		}
	if(rz && !(arz = new AccRange(rz))){
		delete ary;	delete arx;	return 0;
		}
	if(!(x = (double*)malloc(i * sizeof(double)))){
		if(arz) delete arz;
		delete ary;	delete arx;	return 0;
		}
	if(!(y = (double*)malloc(i * sizeof(double)))){
		if(arz) delete arz;
		free(x);	delete arx;	delete ary;	return 0;
		}
	if(rz && !(y = (double*)malloc(i * sizeof(double)))){
		if(arz) delete arz;
		free(y);	free(x);	delete arx;	delete ary;	return 0;
		}
	arx->GetFirst(&c1, &r1);	ary->GetFirst(&c2, &r2);
	if(rz) arz->GetFirst(&c3, &r3);
	for(ndata = j = 0; j < i; j++) {
		if(rz) {
			if(arx->GetNext(&c1, &r1) && ary->GetNext(&c2, & r2) && arz->GetNext(&c3, &r3) &&
				curr_data->GetValue(r1, c1, &currx) && curr_data->GetValue(r2, c2, &curry) &&
				curr_data->GetValue(r3, c3, &currz)) {
				x[ndata] = currx;	y[ndata] = curry;	z[ndata] = currz;	ndata++;
				}
			}
		else {
			if(arx->GetNext(&c1, &r1) && ary->GetNext(&c2, & r2) &&
				curr_data->GetValue(r1, c1, &currx) && curr_data->GetValue(r2, c2, &curry)) {
				x[ndata] = currx;	y[ndata] = curry;	ndata++;
				}
			}
		}
	//common initialization for parser tasks
	push_parser();		//make code reentrant
	init_table();		length = strlen(*par);
	//process parameters
	if(!(buffer = (char*)malloc(length+2))){
		clear_table();	pop_parser();
		if(arz) delete arz;
		free(y);	free(x);	delete arx;	delete ary;
		return 0;
		}
	strcpy(buffer, *par);	buffer[length++] = ';';
	buffer[length] = 0;	buff_pos = 0;
	tab1 = sym_table;
	do {
		yyparse();
		}while(buff_pos < length);
	tab2 = sym_table;	free(buffer);	buffer =0L;
	for(nparam = 0, csr=tab2; csr != tab1; nparam++, csr = csr->next);
	parsym = (symrec**)malloc((nparam+1)*sizeof(symrec*));
	parval = (double**)malloc((nparam+1)*sizeof(double*));
	for(i = 0, csr=tab2; csr != tab1 && i < nparam; i++, csr = csr->next){
		parsym[i] = csr;	parval[i] = &csr->var;
		}
	//do iteratations to optimize fit
	lista = (int*)malloc(sizeof(int)*nparam);
	for(i = 0; i< nparam; i++) lista[i] = i;
	covar = dmatrix(1, nparam, 1, nparam);
	alpha = dmatrix(1, nparam, 1, nparam);
	alamda = -1.0;		itst = 0;
	mrqmin(x, y, z, ndata, parval, nparam, lista, nparam, covar, alpha, &chisq, fcurve, &alamda);
	if(!Check_MRQerror()) {
		for(itst = itst1 = 0, ochisq = chisq; itst < maxiter && chisq > conv && ochisq >= chisq && itst1 < 9; itst++) {
			ochisq = chisq;
			mrqmin(x, y, z, ndata, parval, nparam, lista, nparam, covar, alpha, &chisq, fcurve, &alamda);
			if(ochisq == chisq) itst1++;
			else itst1 = 0;
			}
		alamda = 0.0;
		mrqmin(x, y, z, ndata, parval, nparam, lista, nparam, covar, alpha, &chisq, fcurve, &alamda);
		Check_MRQerror();
		}
	for(i = nparam-1, j = k = l = 0; i >= 0; l = 0, i--) {
		if(k > 20) {
			if(tmp_txt[j-1] == ' ') j--;
			if(tmp_txt[j-1] == ';') j--;
			l = sprintf(tmp_txt+j, "\n");
			j += l;		k = 0;
			}
		l += sprintf(tmp_txt+j, "%s%s=%g;", j && k ? " " : "", parsym[i]->name, parsym[i]->GetValue());
		j += l;			k += l;
		}
	free(*par);	*par = strdup(tmp_txt);
	if(chi_2) *chi_2 = chisq;
	//write back spreadsheet data if necessary
	buffer = *par;	length = strlen(buffer);
	do {
		yyparse();
		}while(buff_pos < length);
	buffer = 0L;
	free_dmatrix(alpha, 1, nparam, 1, nparam);
	free_dmatrix(covar, 1, nparam, 1, nparam);
	if(arz) delete arz;		if(z) free(z);
	free(y);	free(x);	delete arx;	delete ary;
	if(parval) free(parval);	if(parsym) free(parsym);
	clear_table();
	pop_parser();
	if(curr_data){
		curr_data->Command(CMD_CLEAR_ERROR, 0L, 0L);
		curr_data->Command(CMD_REDRAW, 0L, 0L);
		}
	return itst < maxiter ? itst+1 : maxiter;
}