/* Steven Andrews, 9/98; updated 3/02	*/
/* Basis functions for spectral fitting */
/* See documentation called BasisFn doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <string.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <limits.h>
#include "string2.h"
#include "math2.h"
#include "Set.h"
#include "Spectra.h"
#include "BasisFn.h"
#include "BasisFns.h"
#include "Rn.h"
#include "Plot.h"
#include "random.h"
#include "RnLU.h"
#include "RnSort.h"
#include "Utility.h"
#include "SpectFit.h"
#include "DiskIO.h"
#include "VoidComp.h"

#define CHECK(A) {if(!(A)) goto failure;}

set basisfns=NULL;

basisptr BasisAlloc(int n) {
	basisptr b;
	int i;

	if(n<1) return NULL;
	CHECK(b=(basisptr) malloc(sizeof(struct basisfn)));
	b->name[0]='\0';
	b->proc[0]='\0';
	b->color[0]='g';
	b->color[1]='\0';
	b->desc[0]='\0';
	b->addr=NULL;
	b->model=NULL;
	b->n=n;
	b->param=b->pold=b->paramlo=b->paramhi=b->scratch=NULL;
	b->eqn=b->pname=NULL;
	b->freeze=NULL;
	b->isspec=0;
	b->spec=NULL;
	CHECK(b->param=allocV(n));
	CHECK(b->pold=allocV(n));
	CHECK(b->paramlo=allocV(n));
	CHECK(b->paramhi=allocV(n));
	CHECK(b->eqn=(char **) calloc(n,sizeof(char *)));
	for(i=0;i<n;i++) b->eqn[i]=NULL;
	CHECK(b->pname=(char **) calloc(n,sizeof(char *)));
	for(i=0;i<n;i++) b->pname[i]=NULL;
	CHECK(b->freeze=(int *) calloc(n,sizeof(int)));
	CHECK(b->scratch=allocV(n));
	for(i=0;i<n;i++)	{
		b->param[i]=b->pold[i]=0;
		b->paramlo[i]=-FLT_MAX;
		b->paramhi[i]=FLT_MAX;
		CHECK(b->eqn[i]=EmptyString());
		CHECK(b->pname[i]=EmptyString());
		b->freeze[i]=0;
		b->scratch[i]=0; }
	return b;

 failure:
 	BasisFree(b);
 	return NULL; }

void BasisFree(basisptr b) {
	int i;
	
	if(!b) return;
	if(b->param) freeV(b->param);
	if(b->pold) freeV(b->pold);
	if(b->paramlo) freeV(b->paramlo);
	if(b->paramhi) freeV(b->paramhi);
	if(b->eqn) {
		for(i=0;i<b->n;i++)	if(b->eqn[i]) free(b->eqn[i]);
		free(b->eqn); }
	if(b->pname) {
		for(i=0;i<b->n;i++)	if(b->pname[i]) free(b->pname[i]);
		free(b->pname); }
	if(b->freeze)	free(b->freeze);
	return; }

modelptr ModelAlloc(sptr s) {
	modelptr m;

	m=(modelptr) malloc(sizeof(struct modeltype));
	if(!m) return NULL;
	m->name[0]='\0';
	m->file[0]='\0';
	m->color[0]='r';
	m->color[1]='\0';
	m->spec=s;
	m->uncert=NULL;
	m->basis=NULL;
	m->np=0;
	m->sigma=0;
	if(s) {
		m->xmin=minV(s->x,s->n);
		m->xmax=maxV(s->x,s->n);
		m->dx=(m->xmax-m->xmin)/(s->n>20?s->n-1:20); }
	else {
		m->xmin=0;
		m->xmax=10;
		m->dx=0.5; }
	m->covar=NULL;
	m->basis=SetAlloc(&StringCmp,NULL);
	if(!m->basis) {ModelFree(m);return NULL;}
	return m; }

void ModelFree(modelptr m) {
	void *x,*k;
	struct scell *trace;

	if(!m) return;
	if(m->basis)
		for(trace=SetNext(NULL,&k,&x,m->basis);trace;trace=SetNext(trace,&k,&x,m->basis))
			BasisFree((basisptr)x);
	SetFree(m->basis,0,0);
	if(m->covar) freeM(m->covar);
	free(m);
	return; }

int CheckBasis(basisptr b,modelptr m) {
	void *x,*k;
	struct scell *trace;
	int i,er;

	er=0;
	if(b) {
		if(!b->isspec) er=0;
		else if(b->isspec==1&&b->spec&&!b->spec->cmplx) er=0;
		else if(!b->spec) er=16;
		else if(b->isspec==1&&b->spec->cmplx) er=8;
		else if(b->isspec==2&&!b->spec->cmplx) er=12;
		for(i=0;i<b->n&&!er;i++) er=(b->paramlo[i]<=b->paramhi[i])?0:45; }
	if(m)
		for(trace=SetNext(NULL,&k,&x,m->basis);trace&&!er;trace=SetNext(trace,&k,&x,m->basis)) {
			b=(basisptr)x;
			if(!b->isspec) er=0;
			else if(b->isspec==1&&b->spec&&!b->spec->cmplx) er=0;
			else if(!b->spec) er=16;
			else if(b->isspec==1&&b->spec->cmplx) er=8;
			else if(b->isspec==2&&!b->spec->cmplx) er=12;
			for(i=0;i<b->n&&!er;i++) er=(b->paramlo[i]<=b->paramhi[i])?0:45; }
	return er; }

int CheckMSpec(modelptr m) {
	if(!m||!m->spec) return 0;
	if(m->spec->cmplx) return 8;
	return 0; }

int CheckUncert(modelptr m) {
	if(!m||!m->uncert) return 0;
	if(!m->spec) return 16;
	if(m->spec->n!=m->uncert->n) return 11;
	if(!equalV(m->spec->x,m->uncert->x,m->spec->n)) return 11;
	if(m->spec->cmplx) return 8;
	if(minV(m->uncert->y,m->uncert->n)<=0) return 14;	
	return 0; }

int setupbasis() {
	set bset;

	bset=basisfns=SetAlloc(&StringCmp,NULL);
	if(!bset) return 1;
	return getbasis(bset); }

void PlotBasis(basisptr b) {
	float xa,xb,ya,yb,x,dx;

	if(!b->model) return;
	if(CheckBasis(b,NULL)) return;
	GetScales(&xa,&xb,&ya,&yb);
	if(b->model->xmin>xa)	xa=b->model->xmin;
	if(b->model->xmax<xb)	xb=b->model->xmax;
	dx=b->model->dx;
	SetColor(b->color[0]);
	PlotPt(xa,(b->addr)(xa,b->param,b->spec,NULL));
	for(x=xa;x<xb;x+=dx)
		PlotLine(x,(b->addr)(x,b->param,b->spec,NULL));
	PlotLine(xb,(b->addr)(xb,b->param,b->spec,NULL));
	SetColor('0');
	return; }

void PlotModel(modelptr m) {
	float xa,xb,ya,yb,x,dx;

	if(CheckBasis(NULL,m)) return;
	GetScales(&xa,&xb,&ya,&yb);
	if(m->xmin>xa)	xa=m->xmin;
	if(m->xmax<xb)	xb=m->xmax;
	dx=m->dx;
	SetColor(m->color[0]);
	PlotPt(xa,ModelValue(xa,m));
	for(x=xa;x<xb;x+=dx)
		PlotLine(x,ModelValue(x,m));
	PlotLine(xb,ModelValue(xb,m));	
	SetColor('0');
	return; }

void BasisRange(basisptr b,float *xa,float *xb,float *ya,float *yb,int fn) {
	float y,ylo,yhi,xlo,xhi,x;
	
	if(!b->model) return;
	if(CheckBasis(b,NULL)) return;
	xlo=b->model->xmin;
	if(!fn||fn<0&&*xa<xlo||fn>0&&*xa>xlo) *xa=xlo;
	xhi=b->model->xmax;
	if(!fn||fn<0&&*xb>xhi||fn>0&&*xb<xhi) *xb=xhi;
	if(xlo<*xa) xlo=*xa;
	if(xhi>*xb) xhi=*xb;
	ylo=yhi=(b->addr)(xlo,b->param,b->spec,NULL);
	for(x=xlo;x<xhi;x+=b->model->dx)	{
		y=(b->addr)(x,b->param,b->spec,NULL);
		if(y<ylo) ylo=y;
		else if(y>yhi) yhi=y; }
	if(!fn||fn<0&&*ya<ylo||fn>0&&*ya>ylo) *ya=ylo;
	if(!fn||fn<0&&*yb>yhi||fn>0&&*yb<yhi) *yb=yhi;
	return; }

void ModelRange(modelptr m,float *xa,float *xb,float *ya,float *yb,int fn) {
	float y,ylo,yhi,xlo,xhi,x;
	
	if(CheckBasis(NULL,m)) return;
	xlo=m->xmin;
	if(!fn||fn<0&&*xa<xlo||fn>0&&*xa>xlo) *xa=xlo;
	xhi=m->xmax;
	if(!fn||fn<0&&*xb>xhi||fn>0&&*xb<xhi) *xb=xhi;
	if(xlo<*xa) xlo=*xa;
	if(xhi>*xb) xhi=*xb;
	ylo=yhi=ModelValue(xlo,m);
	for(x=xlo;x<xhi;x+=m->dx)	{
		y=ModelValue(x,m);
		if(y<ylo) ylo=y;
		else if(y>yhi) yhi=y; }
	if(!fn||fn<0&&*ya<ylo||fn>0&&*ya>ylo) *ya=ylo;
	if(!fn||fn<0&&*yb>yhi||fn>0&&*yb<yhi) *yb=yhi;
	return; }

int CountFreParam(modelptr m,int **frzptr) {
	int i,j,nfre;
	void *k,*item;
	struct scell *trace;
	basisptr b;

	nfre=0;
	for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
		b=(basisptr)item;
		for(i=0;i<b->n;i++) if(!b->freeze[i]) nfre++; }
	if(frzptr) {
		*frzptr=(int *)calloc(m->np,sizeof(int));
		if(!*frzptr) return -1;
		j=0;
		for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
			b=(basisptr)item;
			for(i=0;i<b->n;i++) (*frzptr)[j++]=b->freeze[i]; }}
	return nfre; }

int FindBasisParam(basisptr *bptr,modelptr m,char* str) {
	int i;
	basisptr b;
	void *k,*item;
	struct scell *trace;

	b=*bptr;
	if(b)	{
		for(i=0;i<b->n;i++) if(!strcmp(b->pname[i],str)) return i;
		return -1; }
	if(!m) return -1;
	for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
		b=*bptr=(basisptr)item;
		for(i=0;i<b->n;i++) if(!strcmp(b->pname[i],str)) return i; }
	*bptr=NULL;
	return -1; }

int FindBasisParam2(basisptr *bptr,modelptr m,float* fltptr) {
	int i;
	basisptr b;
	void *k,*item;
	struct scell *trace;

	if(!m) return -1;
	for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
		b=*bptr=(basisptr)item;
		for(i=0;i<b->n;i++) if(&(b->param[i])==fltptr) return i; }
	*bptr=NULL;
	return -1; }

float ModelValue(float x,modelptr m) {
	float y=0;
	void *k,*item;
	struct scell *trace;
	basisptr b;

	for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
		b=(basisptr)item;
		y+=(b->addr)(x,b->param,b->spec,NULL); }
	return y; }

float ModelValueD(float x,modelptr m,float *d) {
	float y=0;
	int i,j;
	void *k,*item;
	struct scell *trace;
	basisptr b;
	
	j=0;
	for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
		b=(basisptr)item;
		y+=(b->addr)(x,b->param,b->spec,b->scratch);
		for(i=0;i<b->n;i++)
			if(!b->freeze[i]) d[j++]=b->scratch[i]; }
	return y; }

void TypeBasis(basisptr b) {
	int i;

	if(!b) {
		printf("NULL basis function.\n");
		return; }
	if(b->model) printf("%s is a basis function in model %s\n",b->name,b->model->name);
	else printf("%s is a basis function without a model\n",b->name);
	printf(" procedure: %s, color: %s, %i parameters\n",b->proc,b->color,b->n);
	printf(" description: %s\n",b->desc);
	for(i=0;i<b->n;i++)	{
		printf("\t%s=%f\t%s\n",b->pname[i],b->param[i],b->freeze[i]?"fix":"free");
		if(strlen(b->eqn[i])) printf("\t\teqn %sactive: %s\n",b->freeze[i]?"":"in",b->eqn[i]);
		if(b->paramlo[i]>-FLT_MAX&&b->paramhi[i]<FLT_MAX)
			printf("\t\t%f ² %s ² %f\n",b->paramlo[i],b->pname[i],b->paramhi[i]);
		else if(b->paramlo[i]>-FLT_MAX)
			printf("\t\t%f ² %s\n",b->paramlo[i],b->pname[i]);
		else if(b->paramhi[i]<FLT_MAX)
			printf("\t\t%s ² %f\n",b->pname[i],b->paramhi[i]); }
	if(b->isspec&&!b->spec)
		printf(" spectrum required, none assigned.\n");
	else if(b->isspec==1&&!b->spec->cmplx)
		printf(" spectrum: %s\n",b->spec);
	else if(b->isspec==1&&b->spec->cmplx)
		printf(" complex spectrum should be real: %s\n",b->spec);
	else if(b->isspec==2&&b->spec->cmplx)
		printf(" complex spectrum: %s\n",b->spec);
	else if(b->isspec==2&&!b->spec->cmplx)
		printf(" real spectrum should be complex: %s\n",b->spec);
	else if(b->isspec)
		printf(" spectrum: %s\n",b->spec);
	printf("\n"); }

void TypeAllBasis() {
	struct scell *trace;
	void *k,*x;
	
	if(!basisfns) {
		printf("Set of basis functions not loaded in.\n");
		return; }
	printf("%i basis functions loaded.\n",SetCard(basisfns));
	for(trace=SetNext(NULL,&k,&x,basisfns);trace;trace=SetNext(trace,&k,&x,basisfns))
		printf(" %s\t%s\n",(char *) k,((basisptr)x)->desc);
	return; }

void TypeModel(modelptr m) {
	struct scell *trace;
	void *k,*x;
	int n,nfre,i,j,npts;
	basisptr b;
	sptr s;

	if(!m) {
		printf("No model.\n");
		return; }
	if(m->spec) printf("Model: %s for spectrum %s\n",m->name,m->spec->name);
	else printf("Model: %s, no spectrum\n",m->name);
	if(strlen(m->file)) printf(" model file: %s\n",m->file);
	printf(" domain: %f to %f, color: %s\n",m->xmin,m->xmax,m->color);
	if(CheckUncert(m)) printf(" invalid uncertainties\n");
	else if(m->sigma>0&&m->uncert)
		printf(" sigma: %f, uncertainties in: %s\n",m->sigma,m->uncert->name);
	else if(m->sigma>0) printf(" sigma: %f\n",m->sigma);
	else if(m->uncert) printf(" uncertainties in: %s\n",m->uncert->name);
	if(m->spec) {
		s=m->spec;
		if(s->x[0]>=m->xmin) i=0;
		else {
			i=locateV(s->x,m->xmin,s->n);
			if(i<0) i=0;
			else if(s->x[i]<m->xmin) i++; }
		if(s->x[s->n-1]<=m->xmax) j=s->n-1;
		else j=locateV(s->x,m->xmax,s->n);
		npts=j-i+1;
		printf(" points: %i, params: %i with %i free\n",npts,m->np,CountFreParam(m,NULL));
		if(CheckBasis(NULL,m))
			printf(" a basis set is missing numerical data.\n");
		else if(m->uncert||m->sigma>0)
			printf(" chi^2: %e, rms error: %e\n",chisq(m,-1),rmserror(m));
		else printf(" rms error: %e\n",rmserror(m)); }
	else {
		npts=0;
		printf(" no spectrum\n",m->name); }
	n=m->np;
	j=0;
	for(trace=SetNext(NULL,&k,&x,m->basis);trace;trace=SetNext(trace,&k,&x,m->basis))	{
		b=(basisptr)x;
		printf(" %s\n",b->name);
		for(i=0;i<b->n;i++)	{
			if(m->covar&&!b->freeze[i])
				printf("  %s= %f ± %f\tfree\n",b->pname[i],b->param[i],sqrt(m->covar[n*j+j]));
			else
				printf("  %s= %f\t\t\t%s\n",b->pname[i],b->param[i],b->freeze[i]?"fix":"free");
			if(b->eqn[i]&&strlen(b->eqn[i]))
				printf("\t\t%sactive eqn:%s\n",b->freeze[i]?"":"in",b->eqn[i]);
			if(b->paramlo[i]>-FLT_MAX&&b->paramhi[i]<FLT_MAX)
				printf("   %f ² %s ² %f\n",b->paramlo[i],b->pname[i],b->paramhi[i]);
			else if(b->paramlo[i]>-FLT_MAX)
				printf("   %f ² %s\n",b->paramlo[i],b->pname[i]);
			else if(b->paramhi[i]<FLT_MAX)
				printf("   %s ² %f\n",b->pname[i],b->paramhi[i]);
			j++; }}
	printf("\n");
	return; }

basisptr BasisCopy(basisptr b) {
	basisptr a;
	int i;
	
	if(!b) return NULL;
	a=BasisAlloc(b->n);   // *****************
	if(!a) return NULL;
	strcpy(a->name,b->name);
	strcpy(a->proc,b->proc);
	strcpy(a->color,b->color);
	strcpy(a->desc,b->desc);
	a->model=b->model;
	a->addr=b->addr;
	for(i=0;i<a->n;i++)	{
		a->param[i]=b->param[i];
		a->pold[i]=b->pold[i];
		a->paramlo[i]=b->paramlo[i];
		a->paramhi[i]=b->paramhi[i];
		strcpy(a->pname[i],b->pname[i]);
		strcpy(a->eqn[i],b->eqn[i]);
		a->freeze[i]=b->freeze[i];
		a->scratch[i]=b->scratch[i]; }
	a->isspec=b->isspec;
	a->spec=b->spec;
	return a; }

modelptr ModelCopy(modelptr m) {
	modelptr a;
	basisptr b;
	struct scell *trace;
	void *key,*item;
	int ok;

	if(!m) return NULL;
	a=ModelAlloc(m->spec);
	if(!a) return NULL;
	a->uncert=m->uncert;
	strcpy(a->name,m->name);
	strcpy(a->file,m->file);
	strcpy(a->color,m->color);
	a->sigma=m->sigma;
	a->np=m->np;
	a->xmin=m->xmin;
	a->xmax=m->xmax;
	a->dx=m->dx;
	if(m->covar) {
		a->covar=allocM(a->np,a->np);
		if(!a->covar) {ModelFree(a);return NULL;}
		copyM(m->covar,a->covar,m->np,m->np); }
	for(trace=SetNext(NULL,&key,&item,m->basis);trace;trace=SetNext(trace,&key,&item,m->basis))	{
		b=BasisCopy((basisptr)item);
		if(b)	{
			b->model=a;
			ok=SetInsert((void *)b->name,(void *)b,a->basis);
			if(!ok) {ModelFree(a);return NULL;}}}
	return a; }

basisptr AddBasis(basisptr b,char *bproc,sptr s,modelptr m) {
	void *x;
	basisptr b2;
	int i,ok;

	if(b)	b=BasisCopy(b);
	else if(bproc) b=BasisCopy((basisptr)SetItem((void *)bproc,basisfns));
	if(!b) return NULL;
	if(!strchr(b->name,':'))
		sprintf(b->name+strlen(b->name),":");
	i=0;
	sprintf(strchr(b->name,':')+1,"%i",i++);
	for(x=SetItem((void *)b->name,m->basis);x;x=SetItem((void *)b->name,m->basis))
		sprintf(strchr(b->name,':')+1,"%i",i++);
	b->model=m;
	b->spec=s;
	if(!b->addr)	{
		b2=SetItem((void *)b->proc,basisfns);
		if(b2) b->addr=b2->addr;
		else {BasisFree(b);return NULL;}}
	ok=SetInsert((void *)b->name,(void *)b,m->basis);
	if(!ok) {BasisFree(b);return NULL;}
	m->np+=b->n;
	if(m->covar) {
		freeM(m->covar);
		m->covar=NULL; }
	return b; }

void RemoveBasis(basisptr b) {
	modelptr m;

	if(!b) return;
	m=b->model;
	if(!m) return;
	if(SetDelete((void*)b->name,(void*)b,b->model->basis,0,0)) {
		m->np-=b->n;
		if(m->covar) {
			freeM(m->covar);
			m->covar=NULL; }}
	b->model=NULL;
	return; }

int SaveModel(modelptr m) {
	FILE *fptr;
	char *name;
	int nmfre,er,i,j;
	char str[STRCHAR];
	struct scell *trace;
	void *k,*x;
	basisptr b;

	name=m->file;
	er=GetName(&name,&nmfre);
	if(er) {m->file[0]='\0';return er;}
	if(!strcmp(name,"stdout")) fptr=stdout;
	else {
		fptr=fopen(name,"r");
		if(fptr) {
			fclose(fptr);
			fprintf(stderr,"File exists. Do you want to overwrite it? ");
			scanf("%s",str);
			if(str[0]!='y'&&str[0]!='Y') {
				if(nmfre) free(name);
				m->file[0]='\0';
				return 2; }}
		fptr=fopen(name,"w"); }
	if(nmfre) free(name);
	if(!fptr) return 3;
	er=0;
	er=er||fprintf(fptr,"# Model saved by SpectFit 2.0.\n\n")<0;
	er=er||fprintf(fptr,"name: %s\n",m->name)<0;
	er=er||fprintf(fptr,"file: %s\n",m->file)<0;
	er=er||fprintf(fptr,"color: %s\n",m->color)<0;
	er=er||fprintf(fptr,"# data modeled: %s\n",m->spec?m->spec->name:"none")<0;
	er=er||fprintf(fptr,"# uncertainties: %s\n",m->uncert?m->uncert->name:"none")<0;
	er=er||fprintf(fptr,"sigma: %f\n",m->sigma)<0;
	er=er||fprintf(fptr,"xmin: %f\n",m->xmin)<0;
	er=er||fprintf(fptr,"xmax: %f\n",m->xmax)<0;
	er=er||fprintf(fptr,"dx: %f\n",m->dx)<0;
	if(CheckBasis(NULL,m))
		er=er||fprintf(fptr,"# a basis set is missing numerical data.\n")<0;
	else if(m->uncert||m->sigma>0)
		er=er||fprintf(fptr,"# chi^2: %e, rms error: %e\n",chisq(m,-1),rmserror(m))<0;
	else
		er=er||fprintf(fptr,"# rms error: %e\n",rmserror(m))<0;
	j=0;
	for(trace=SetNext(NULL,&k,&x,m->basis);trace;trace=SetNext(trace,&k,&x,m->basis))	{
		b=(basisptr)x;
		er=er||fprintf(fptr,"basis: %s\n",b->proc)<0;
		er=er||fprintf(fptr," name: %s\n",b->name)<0;
		er=er||fprintf(fptr," color: %s\n",b->color)<0;
		er=er||fprintf(fptr," desc: %s\n",b->desc)<0;
		if(b->spec)
			er=er||fprintf(fptr," spec: %s\n",b->spec->name)<0;
		for(i=0;i<b->n;i++)	{
			er=er||fprintf(fptr," param: %s\n",b->pname[i])<0;
			if(m->covar&&!b->freeze[i])
				er=er||fprintf(fptr,"  value: %f ± %f\n",b->param[i],sqrt(m->covar[m->np*j+j]))<0;
			else
				er=er||fprintf(fptr,"  value: %f\n",b->param[i])<0;
			j++;
			if(b->freeze[i])
				er=er||fprintf(fptr,"  fix\n")<0;
			if(strlen(b->eqn[i]))
				er=er||fprintf(fptr,"  eqn: %s\n",b->eqn[i])<0;
			if(b->paramlo[i]>-FLT_MAX)
				er=er||fprintf(fptr,"  min: %f\n",b->paramlo[i])<0;
			if(b->paramhi[i]<FLT_MAX)
				er=er||fprintf(fptr,"  max: %f\n",b->paramhi[i])<0; }
		er=er||fprintf(fptr," endbasis\n")<0; }
	if(er) er=3;
	fclose(fptr);
	return er; }

int LoadModel(modelptr *mptr,char *name,sptr s) {
	FILE *fptr;
	modelptr m;
	basisptr b;
	char line[STRCHAR],word[STRCHAR],word2[STRCHAR],*line2;
	int i,er,nmfre,itct,endflag;

	if(!mptr) return 6;
	m=*mptr=NULL;
	er=GetName(&name,&nmfre);
	if(er) return er;
	if(!strcmp(name,"stdin")) fptr=stdin;
	else fptr=fopen(name,"r");
	if(nmfre) free(name);
	if(!fptr) return 5;

	m=ModelAlloc(s);
	if(!m) {er=1;goto failure;}
	
	endflag=0;
	er=1000;
	while(!endflag&&fgets(line,STRCHAR,fptr)) {
		er++;
		if(strlen(line)<STRCHAR-1) line[strlen(line)-1]='\0';
		itct=sscanf(line,"%s",word);
		line2=strnword(line,2);
		if(line2) sscanf(line2,"%s",word2);
		if(itct<=0);
		else if(line[0]=='#');
		else if(!strcmp(word,"end"))
			endflag=1;
		else if(!line2) goto failure;
		else if(!strcmp(word,"name:"))
			strncpy(m->name,word2,STRCHAR);
		else if(!strcmp(word,"file:"))
			strncpy(m->file,word2,STRCHAR);
		else if(!strcmp(word,"color:"))
			strncpy(m->color,word2,STRCHAR);
		else if(!strcmp(word,"sigma:")) {
			CHECK(sscanf(line2,"%f",&m->sigma)==1); }
		else if(!strcmp(word,"xmin:")) {
			CHECK(sscanf(line2,"%f",&m->xmin)==1); }
		else if(!strcmp(word,"xmax:")) {
			CHECK(sscanf(line2,"%f",&m->xmax)==1); }
		else if(!strcmp(word,"dx:")) {
			CHECK(sscanf(line2,"%f",&m->dx)==1); }
		else if(!strcmp(word,"basis:")) {
			CHECK(b=AddBasis(NULL,word2,NULL,m));
			i=-1;
			while(!endflag&&fgets(line,STRCHAR,fptr)) {
				er++;
				if(strlen(line)<STRCHAR-1) line[strlen(line)-1]='\0';
				itct=sscanf(line,"%s",word);
				line2=strnword(line,2);
				if(line2) sscanf(line2,"%s",word2);
				if(itct<=0);
				else if(line[0]=='#');
				else if(!strcmp(word,"endbasis")) {
					for(i=0;i<b->n;i++) {
						CHECK(b->paramlo[i]<=b->paramhi[i]); }
					CHECK(!b->isspec||b->spec);
					endflag=1; }
				else if(!strcmp(word,"fix")&&i>-1)
					b->freeze[i]=1;
				else if(!strcmp(word,"free")&&i>-1)
					b->freeze[i]=0;
				else if(!line2) goto failure;
				else if(!strcmp(word,"name:")) {
					if(!SetItem((void *)word2,m->basis))
						strncpy(b->name,word2,STRCHAR); }
				else if(!strcmp(word,"color:"))
					strncpy(b->color,word2,STRCHAR);
				else if(!strcmp(word,"desc:"))
					strncpy(b->desc,line2,STRCHAR);
				else if(!strcmp(word,"spec:")) {
					CHECK(b->spec=Name2Spec(word2)); }
				else if(!strcmp(word,"param:")) {
					for(i=0;i<b->n&&strcmp(b->pname[i],word2);i++)
					CHECK(i<b->n); }
				else if(!strcmp(word,"value:")&&i>-1) {
					CHECK(sscanf(line2,"%f",&b->param[i])==1); }
				else if(!strcmp(word,"min:")&&i>-1) {
					CHECK(sscanf(line2,"%f",&b->paramlo[i])==1); }
				else if(!strcmp(word,"max:")&&i>-1) {
					CHECK(sscanf(line2,"%f",&b->paramhi[i])==1); }
				else if(!strcmp(word,"eqn:")&&i>-1)
					strncpy(b->eqn[i],line2,STRCHAR);
				else goto failure; } 
			endflag=0; }
		else goto failure; }
	if(fptr!=stdin) fclose(fptr);
	fptr=NULL;
	CHECK(m->xmin<m->xmax&&m->dx>0);
	*mptr=m;
	return 0;

 failure:
	if(fptr&&fptr!=stdin) fclose(fptr);
	ModelFree(m);
	*mptr=NULL;
	return er; }

int BasisMath(basisptr b,void *v,sptr *ansptr,float k,char *fn) {
	float dx,y,xlo,*param,*s3x,*s3y;
	int i,j,n,er;
	sptr s,s3,bspec;
	basisptr b2;

	if(!b||!ansptr||!fn||!b->model||!b->model->spec) return 6;
	if((er=CheckBasis(b,NULL))) return er;
	if(*ansptr) SpectFree(*ansptr);
	*ansptr=NULL;
	s3=SpectAlloc(b->name,NULL,fn,b->model->spec->xunit,b->model->spec->yunit);
	if(!s3) return 1;
	strncat(s3->name,"s",STRCHAR-strlen(s3->name));
	strncat(s3->desc,b->name,STRCHAR-strlen(s3->desc));
	strncat(s3->desc," ",STRCHAR-strlen(s3->desc));
	strncat(s3->desc,fn,STRCHAR-strlen(s3->desc));
	dx=b->model->dx;
	xlo=b->model->xmin;
	n=s3->n=(b->model->xmax-xlo)/dx+1;
	s3x=s3->x=allocV(n);
	s3y=s3->y=allocV(n);
	if(!s3->x||!s3->y) {SpectFree(s3);return 1;}
	for(i=0;i<n;i++) s3->x[i]=xlo+i*dx;
	param=b->param;
	bspec=b->spec;

	if(!strcmp(fn,"exp b"))
		for(i=0;i<n;i++) {
			y=exp((b->addr)(s3x[i],param,bspec,NULL));
			s3y[i]=y<=FLT_MAX?y:0; }
	else if(!strcmp(fn,"ln b"))
		for(i=0;i<n;i++) {
			y=(b->addr)(s3x[i],param,bspec,NULL);
			s3y[i]=y>0?log(y):0; }
	else if(!strcmp(fn,"10^b"))
		for(i=0;i<n;i++) {
			y=pow(10,(b->addr)(s3x[i],param,bspec,NULL));
			s3y[i]=y<=FLT_MAX?y:0; }
	else if(!strcmp(fn,"log b"))
		for(i=0;i<n;i++) {
			y=(b->addr)(s3x[i],param,bspec,NULL);
			s3y[i]=y>0?log10(y):0; }
	else if(!strcmp(fn,"sqrt b"))
		for(i=0;i<n;i++) {
			y=(b->addr)(s3x[i],param,bspec,NULL);
			s3y[i]=y>0?sqrt(y):0; }
	else if(!strcmp(fn,"b+s")||!strcmp(fn,"s+b"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3->y[i]=interpolate1(s->x,s->y,s->n,&j,s3x[i])+(b->addr)(s3x[i],param,bspec,NULL); }
	else if(!strcmp(fn,"b-s"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)-interpolate1(s->x,s->y,s->n,&j,s3x[i]); }
	else if(!strcmp(fn,"b*s")||!strcmp(fn,"s*b"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3y[i]=interpolate1(s->x,s->y,s->n,&j,s3x[i])*(b->addr)(s3x[i],param,bspec,NULL); }
	else if(!strcmp(fn,"b/s"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++) {
			y=interpolate1(s->x,s->y,s->n,&j,s3x[i]);
			s3y[i]=(y<0||y>0)?(b->addr)(s3x[i],param,bspec,NULL)/y:0; }}
	else if(!strcmp(fn,"s-b"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3y[i]=interpolate1(s->x,s->y,s->n,&j,s3x[i])-(b->addr)(s3x[i],param,bspec,NULL); }
	else if(!strcmp(fn,"s/b"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++) {
			y=(b->addr)(s3x[i],param,bspec,NULL);
			s3y[i]=(y<0||y>0)?interpolate1(s->x,s->y,s->n,&j,s3x[i])/y:0; }}
	else if(!strcmp(fn,"b+b"))	{
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)+(b2->addr)(s3x[i],b2->param,b2->spec,NULL); }
	else if(!strcmp(fn,"b-b")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)-(b2->addr)(s3x[i],b2->param,b2->spec,NULL); }
	else if(!strcmp(fn,"b*b"))	{
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)*(b2->addr)(s3x[i],b2->param,b2->spec,NULL); }
	else if(!strcmp(fn,"b/b")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++) {
			y=(b2->addr)(s3x[i],b2->param,b2->spec,NULL);
			s3y[i]=(y<0||y>0)?(b->addr)(s3x[i],b->param,b->spec,NULL)/y:0; }}
	else if(!strcmp(fn,"b+k")||!strcmp(fn,"k+b"))
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)+k;
	else if(!strcmp(fn,"b*k")||!strcmp(fn,"k*b"))
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)*k;
	else if(!strcmp(fn,"k-b"))
		for(i=0;i<n;i++)
			s3y[i]=k-(b->addr)(s3x[i],param,bspec,NULL);
	else if(!strcmp(fn,"k/b"))
		for(i=0;i<n;i++) {
			y=(b->addr)(s3x[i],param,bspec,NULL);
			s3y[i]=(y<0||y>0)?k/y:0; }
	else if(!strcmp(fn,"b-k"))
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)-k;
	else if(!strcmp(fn,"b/k")) {
		if(!(k<0||k>0)) {SpectFree(s3);return 9;}
		for(i=0;i<n;i++)
			s3y[i]=(b->addr)(s3x[i],param,bspec,NULL)/k; }
	else if(!strcmp(fn,"b^k"))
		if(k>0) for(i=0;i<n;i++)
				s3y[i]=pow((b->addr)(s3x[i],param,bspec,NULL),k);
		else if(k==floor(k)) for(i=0;i<n;i++) {
				y=(b->addr)(s3x[i],param,bspec,NULL);
				s3y[i]=(y<0||y>0)?pow(y,k):0; }
		else for(i=0;i<n;i++) {
				y=(b->addr)(s3x[i],param,bspec,NULL);
				s3y[i]=y>0?pow(y,k):0; }
	else {
		SpectFree(s3);
		return 7; }
	*ansptr=s3;
	return 0; }

int ModelMath(modelptr m,void *v,sptr *ansptr,float k,char *fn) {
	float dx,y,*s3x,*s3y;
	int i,j,n,er;
	basisptr b2;
	modelptr m2;
	sptr s3,s;

	if(!m||!ansptr||!fn) return 6;
	if((er=CheckBasis(NULL,m))) return er;
	if(*ansptr) SpectFree(*ansptr);
	*ansptr=NULL;
	if(m->spec) s3=SpectAlloc(m->name,NULL,fn,m->spec->xunit,m->spec->yunit);
	else s3=SpectAlloc(m->name,NULL,fn,"x","y");
	if(!s3) return 1;
	strncat(s3->name,"s",STRCHAR-strlen(s3->name));
	strncat(s3->desc,m->name,STRCHAR-strlen(s3->desc));
	strncat(s3->desc," ",STRCHAR-strlen(s3->desc));
	strncat(s3->desc,fn,STRCHAR-strlen(s3->desc));
	dx=m->dx;
	n=s3->n=(m->xmax-m->xmin)/dx+1;
	s3x=s3->x=allocV(n);
	s3y=s3->y=allocV(n);
	if(!s3x||!s3y) {SpectFree(s3);return 1;}
	for(i=0;i<n;i++) s3x[i]=m->xmin+i*dx;

	if(!strcmp(fn,"exp m"))
		for(i=0;i<n;i++) {
			y=exp(ModelValue(s3x[i],m));
			s3y[i]=y<=FLT_MAX?y:0; }
	else if(!strcmp(fn,"ln m"))
		for(i=0;i<n;i++) {
			y=ModelValue(s3x[i],m);
			s3y[i]=y>0?log(y):0; }
	else if(!strcmp(fn,"10^m"))
		for(i=0;i<n;i++) {
			y=pow(10,ModelValue(s3x[i],m));
			s3y[i]=y<=FLT_MAX?y:0; }
	else if(!strcmp(fn,"log m"))
		for(i=0;i<n;i++) {
			y=ModelValue(s3x[i],m);
			s3y[i]=y>0?log10(y):0; }
	else if(!strcmp(fn,"sqrt m"))
		for(i=0;i<n;i++) {
			y=ModelValue(s3x[i],m);
			s3y[i]=y>0?sqrt(y):0; }
	else if(!strcmp(fn,"m+s")||!strcmp(fn,"s+m"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3y[i]=interpolate1(s->x,s->y,s->n,&j,s3x[i])+ModelValue(s3x[i],m); }
	else if(!strcmp(fn,"m-s")) {
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)-interpolate1(s->x,s->y,s->n,&j,s3x[i]); }
	else if(!strcmp(fn,"m*s")||!strcmp(fn,"s*m"))	{
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3y[i]=interpolate1(s->x,s->y,s->n,&j,s3x[i])*ModelValue(s3x[i],m); }
	else if(!strcmp(fn,"m/s")) {
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++) {
			y=interpolate1(s->x,s->y,s->n,&j,s3x[i]);
			s3y[i]=(y<0||y>0)?ModelValue(s3x[i],m)/y:0; }}
	else if(!strcmp(fn,"s-m")) {
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++)
			s3y[i]=interpolate1(s->x,s->y,s->n,&j,s3x[i])-ModelValue(s3x[i],m); }
	else if(!strcmp(fn,"s/m")) {
		s=(sptr) v;
		if(!s) {SpectFree(s3);return 6;}
		if(s->cmplx) {SpectFree(s3);return 8;}
		j=-2;
		for(i=0;i<n;i++) {
			y=ModelValue(s3x[i],m);
			s3y[i]=(y<0||y>0)?interpolate1(s->x,s->y,s->n,&j,s3x[i])/y:0; }}
	else if(!strcmp(fn,"m+b")||!strcmp(fn,"b+m")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)+(b2->addr)(s3x[i],b2->param,b2->spec,NULL); }
	else if(!strcmp(fn,"m-b")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)-(b2->addr)(s3x[i],b2->param,b2->spec,NULL); }
	else if(!strcmp(fn,"m*b")||!strcmp(fn,"b*m")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)*(b2->addr)(s3x[i],b2->param,b2->spec,NULL); }
	else if(!strcmp(fn,"m/b")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++) {
			y=(b2->addr)(s3x[i],b2->param,b2->spec,NULL);
			s3y[i]=(y<0||y>0)?ModelValue(s3x[i],m)/y:0; }}
	else if(!strcmp(fn,"b-m")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=(b2->addr)(s3x[i],b2->param,b2->spec,NULL)-ModelValue(s3x[i],m); }
	else if(!strcmp(fn,"b/m")) {
		b2=(basisptr) v;
		if(!b2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++) {
			y=ModelValue(s3x[i],m);
			s3y[i]=(y<0||y>0)?(b2->addr)(s3x[i],b2->param,b2->spec,NULL)/y:0; }}
	else if(!strcmp(fn,"m+m")) {
		m2=(modelptr) v;
		if(!m2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)+ModelValue(s3x[i],m2); }
	else if(!strcmp(fn,"m-m")) {
		m2=(modelptr) v;
		if(!m2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)-ModelValue(s3x[i],m2); }
	else if(!strcmp(fn,"m*m")) {
		m2=(modelptr) v;
		if(!m2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)*ModelValue(s3x[i],m2); }
	else if(!strcmp(fn,"m/m")) {
		m2=(modelptr) v;
		if(!m2) {SpectFree(s3);return 6;}
		for(i=0;i<n;i++) {
			y=ModelValue(s3x[i],m2);
			s3y[i]=(y<0||y>0)?ModelValue(s3x[i],m)/y:0; }}
	else if(!strcmp(fn,"m+k")||!strcmp(fn,"k+m"))
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)+k;
	else if(!strcmp(fn,"m*k")||!strcmp(fn,"k*m"))
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)*k;
	else if(!strcmp(fn,"k-m"))
		for(i=0;i<n;i++)
			s3y[i]=k-ModelValue(s3x[i],m);
	else if(!strcmp(fn,"k/m"))
		for(i=0;i<n;i++) {
			y=ModelValue(s3x[i],m);
			s3y[i]=(y<0||y>0)?k/y:0; }
	else if(!strcmp(fn,"m-k"))
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)-k;
	else if(!strcmp(fn,"m/k")) {
		if(!(k<0||k>0)) {SpectFree(s3);return 9;}
		for(i=0;i<n;i++)
			s3y[i]=ModelValue(s3x[i],m)/k; }
	else if(!strcmp(fn,"m^k"))
		if(k>0) for(i=0;i<n;i++)
				s3y[i]=pow(ModelValue(s3x[i],m),k);
		else if(k==floor(k)) for(i=0;i<n;i++) {
				y=ModelValue(s3x[i],m);
				s3y[i]=(y<0||y>0)?pow(y,k):0; }
		else for(i=0;i<n;i++) {
				y=ModelValue(s3x[i],m);
				s3y[i]=y>0?pow(y,k):0; }
	else {
		SpectFree(s3);
		return 7; }
	*ansptr=s3;
	return 0; }

void endbasis() {
	struct scell *trace,*trace2;
	void *k,*x;

	for(trace=SetNext(NULL,&k,&x,basisfns);trace;trace=trace2) {
		trace2=SetNext(trace,&k,&x,basisfns);
		BasisFree((basisptr) x); }
	SetFree(basisfns,0,0);
	basisfns=NULL;
	return; }

float rmserror(modelptr m) {
	float diff,ch2,siginv,*sx,*sy,*uy;
	int lo,hi,er;
	int i;
	sptr s;

	if(!m) return -6;
	if(!m->spec) return -16;
	if(m->spec->cmplx) return -8;
	if((er=CheckBasis(NULL,m))) return -er;
	s=m->spec;
	sx=s->x;
	sy=s->y;
	if(sx[0]>=m->xmin) lo=0;
	else {
		lo=locateV(sx,m->xmin,s->n);
		if(lo<0) lo=0;
		else if(sx[lo]<m->xmin) lo++; }
	if(sx[s->n-1]<=m->xmax) hi=s->n-1;
	else hi=locateV(sx,m->xmax,s->n);
	ch2=0;
	for(i=lo;i<=hi;i++)	{
		diff=(sy[i]-ModelValue(sx[i],m));
		ch2+=diff*diff; }
	return sqrt(ch2/(hi-lo+1)); }

float chisq(modelptr m,float ch2max) {
	float diff,ch2,siginv,*sx,*sy,*uy;
	int lo,hi;
	int i;
	sptr s;

	if(ch2max<0) ch2max=FLT_MAX;
	s=m->spec;
	if(m->sigma>0) siginv=1.0/m->sigma;
	else siginv=1;
	sx=s->x;
	sy=s->y;
	if(m->uncert) uy=m->uncert->y;
	if(sx[0]>=m->xmin) lo=0;
	else {
		lo=locateV(sx,m->xmin,s->n);
		if(lo<0) lo=0;
		else if(sx[lo]<m->xmin) lo++; }
	if(sx[s->n-1]<=m->xmax) hi=s->n-1;
	else hi=locateV(sx,m->xmax,s->n);
	ch2=0;
	if(!m->uncert)
		for(i=lo;i<=hi&&ch2<ch2max;i++)	{
			diff=(sy[i]-ModelValue(sx[i],m))*siginv;
			ch2+=diff*diff; }
	else
		for(i=lo;i<=hi&&ch2<ch2max;i++)	{
			diff=(sy[i]-ModelValue(sx[i],m))*siginv/uy[i];
			ch2+=diff*diff; }
	return ch2<ch2max?ch2:ch2max; }

float ChiSqD(modelptr m,float ch2max,float *alpha,float *beta,int np) {
	float diff,ch2,siginv,*sx,*sy,*uy;
	float *deriv;
	int lo,hi;
	int i,k,k2;
	sptr s;

	if(ch2max<0) ch2max=FLT_MAX;
	s=m->spec;
	if(m->sigma>0) siginv=1.0/m->sigma;
	else siginv=1;
	sx=s->x;
	sy=s->y;
	if(m->uncert) uy=m->uncert->y;
	deriv=allocV(np);
	if(!deriv) return -1;
	if(sx[0]>=m->xmin) lo=0;
	else {
		lo=locateV(sx,m->xmin,s->n);
		if(lo<0) lo=0;
		else if(sx[lo]<m->xmin) lo++; }
	if(sx[s->n-1]<=m->xmax) hi=s->n-1;
	else hi=locateV(sx,m->xmax,s->n);
	setstdV(beta,np,0);
	setstdM(alpha,np,np,0);
	
	ch2=0;
	if(!m->uncert)
		for(i=lo;i<=hi&&ch2<ch2max;i++)	{
			diff=(sy[i]-ModelValueD(sx[i],m,deriv))*siginv;
			ch2+=diff*diff;
			multKV(siginv,deriv,deriv,np);
			for(k=0;k<np;k++)	{
				beta[k]+=diff*deriv[k];
				for(k2=0;k2<np;k2++) alpha[np*k+k2]+=deriv[k]*deriv[k2]; }}
	else if(m->uncert)
		for(i=lo;i<=hi&&ch2<ch2max;i++)	{
			diff=(sy[i]-ModelValueD(sx[i],m,deriv))*siginv/uy[i];
			ch2+=diff*diff;
			for(k=0;k<np;k++) deriv[k]*=siginv/uy[i];
			for(k=0;k<np;k++)	{
				beta[k]+=diff*deriv[k];
				for(k2=0;k2<np;k2++) alpha[np*k+k2]+=deriv[k]*deriv[k2]; }}
	freeV(deriv);
	for(k=0;k<np;k++) {
		if(!(beta[k]<FLT_MAX&&beta[k]>-FLT_MAX)) beta[k]=0;
		for(k2=0;k2<np;k2++)
			if(!(alpha[np*k+k2]<FLT_MAX&&alpha[np*k+k2]>-FLT_MAX)) alpha[np*k+k2]=0; }
	return ch2<ch2max?ch2:ch2max; }

int prefit(modelptr m,float ***pptr,float **ploptr,float **phiptr) {
	int nfre,ntot,i,jfre,er;
	float **p;
	void *k,*item;
	struct scell *trace;
	basisptr b;

	if(!m) return -6;
	if(!m->spec) return -16;
	if(m->spec->cmplx) return -8;
	if((er=CheckBasis(NULL,m))) return -er;
	if((er=CheckUncert(m))) return -er;
	p=*pptr=NULL;
	*ploptr=NULL;
	*phiptr=NULL;
	nfre=CountFreParam(m,NULL);
	if(!nfre) return 0;
	ntot=m->np;
	p=(float **) calloc(nfre,sizeof(float *));							/* pointers to parameters */
	*ploptr=allocV(nfre);
	*phiptr=allocV(nfre);
	if(!p||!*ploptr||!*phiptr) {
		if(p) free(p);
		if(*ploptr) freeV(*ploptr);
		if(*phiptr) freeV(*phiptr);
		return -1; }
	jfre=0;
	for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
		b=(basisptr)item;
		for(i=0;i<b->n;i++) {
			b->pold[i]=b->param[i];
			if(!b->freeze[i]) {
				if(b->param[i]<b->paramlo[i]) b->param[i]=b->paramlo[i];
				else if(b->param[i]>b->paramhi[i]) b->param[i]=b->paramhi[i];
				p[jfre]=&b->param[i];
				(*ploptr)[jfre]=b->paramlo[i];
				(*phiptr)[jfre]=b->paramhi[i];
				jfre++; }}}
	*pptr=p;
	return nfre; }

int ModelCovar(modelptr m) {
	int rmssig,er,jfre,nfre,itot,jtot,ntot,*frz,*luindx;
	float det,*alpha,*alphalu,*covar,*ident,*beta;

	if(!m) return 6;
	if(!m->spec) return 16;
	if(m->spec->cmplx) return 8;
	if((er=CheckBasis(NULL,m))) return er;
	if((er=CheckUncert(m))) return er;
	nfre=CountFreParam(m,&frz);
	if(nfre<0) return -nfre;
	else if(!nfre) return 0;
	er=0;
	rmssig=0;
	alpha=allocM(nfre,nfre);
	alphalu=allocM(nfre,nfre);
	covar=allocM(nfre,nfre);
	ident=allocM(nfre,nfre);
	beta=allocV(nfre);
	luindx=NULL;
	if(!alpha||!alphalu||!covar||!ident||!beta)	{er=1;goto endfunction;}
	if(m->sigma<=0&&!m->uncert) {
		rmssig=1;
		m->sigma=rmserror(m); }
	ChiSqD(m,-1,alpha,beta,nfre);
	copyM(alpha,alphalu,nfre,nfre);
	det=LUdecomp(alphalu,nfre,&luindx);
	if(det==0) {er=9;goto endfunction;}
	setstdM(covar,nfre,nfre,1);
	LUsolveM(alphalu,luindx,covar,nfre,nfre);
	setstdM(ident,nfre,nfre,1);
	LUimprM(alpha,alphalu,luindx,covar,ident,nfre,nfre);

	ntot=m->np;
	if(!m->covar) m->covar=allocM(ntot,ntot);
	if(!m->covar)	{er=1;goto endfunction;}
	setstdM(m->covar,ntot,ntot,0);
	jfre=0;
	for(itot=0;itot<ntot;itot++)
		for(jtot=0;jtot<ntot;jtot++)
			if(!frz[itot]&&!frz[jtot]) m->covar[ntot*itot+jtot]=covar[jfre++];

 endfunction:
	if(rmssig) m->sigma=0;
	if(frz) free(frz);
	if(alpha) freeM(alpha);
	if(alphalu) freeM(alphalu);
	if(covar) freeM(covar);
	if(ident) freeM(ident);
	if(beta) freeV(beta);
	if(luindx) free(luindx);
	if(er&&m->covar) {
		freeM(m->covar);
		m->covar=NULL; }
	return er; }

void unfit(modelptr m) {
	int i;
	float f1;
	void *k,*item;
	struct scell *trace;
	basisptr b;

	for(trace=SetNext(NULL,&k,&item,m->basis);trace;trace=SetNext(trace,&k,&item,m->basis))	{
		b=(basisptr)item;
		for(i=0;i<b->n;i++) {
			f1=b->param[i];
			b->param[i]=b->pold[i];
			b->pold[i]=f1; }}
	if(m->covar) {
		freeM(m->covar);
		m->covar=NULL; }
	return; }

int RandomFit(modelptr m) {
	float ch2,f,f2;
	float **p,*plo,*phi,*step;
	int i,j,np,er;

	np=prefit(m,&p,&plo,&phi);
	if(np<0) return -np;
	else if(!np) return 0;
	step=allocV(np);
	if(!step) {free(p);freeV(plo);freeV(phi);return 1;}
	for(j=0;j<np;j++)	step[j]=fabs(*p[j]/50)+0.001;

	ch2=chisq(m,-1);
	for(i=0;i<np*30&&!inkey();i++) {
		j=intrand(np);
		f=*p[j];
		*p[j]+=binomrand(3,0,step[j]);
		if(*p[j]<plo[j]) *p[j]=plo[j];
		else if(*p[j]>phi[j]) *p[j]=phi[j];
		UpdateModel(m);
		if((f2=chisq(m,ch2))<ch2)	{
			i=0;
			step[j]*=1.2;
			ch2=f2; }
		else {
			*p[j]=f;
			UpdateModel(m);
			step[j]*=0.9; }}
	er=ModelCovar(m);
	freeV(step);
	free(p);
	freeV(plo);
	freeV(phi);
	return er; }

int LinearFit(modelptr m) {
	float *A,*Atrans,*alpha,*alphalu,*beta,*ans,*b;
	float **p,det,*plo,*phi,ch2,siginv,*uy;
	int i,j,lo,hi,np,n,er;
	int *luindx;
	sptr s;

	er=0;
	np=prefit(m,&p,&plo,&phi);
	if(np<0) return -np;
	else if(!np) return 0;
	s=m->spec;
	if(m->sigma>0) siginv=1.0/m->sigma;
	else siginv=1;
	if(m->uncert) uy=m->uncert->y;
	if(s->x[0]>=m->xmin) lo=0;
	else {
		lo=locateV(s->x,m->xmin,s->n);
		if(lo<0) lo=0;
		else if(s->x[lo]<m->xmin) lo++; }
	if(s->x[s->n-1]<=m->xmax) hi=s->n-1;
	else hi=locateV(s->x,m->xmax,s->n);
	n=hi-lo+1;
	A=allocM(n,np);
	Atrans=allocM(np,n);
	alpha=allocM(np,np);
	alphalu=allocM(np,np);
	beta=allocV(np);
	ans=allocV(np);
	b=allocV(n);
	luindx=NULL;
	if(!A||!Atrans||!alpha||!beta||!alphalu||!ans||!b) {er=1;goto endfunction; }

	ch2=chisq(m,-1);
	if(!m->uncert)
		for(i=lo;i<=hi;i++) {
			b[i-lo]=(s->y[i]-ModelValueD(s->x[i],m,ans))*siginv;
			for(j=0;j<np;j++) A[np*(i-lo)+j]=ans[j]*siginv; }
	else
		for(i=lo;i<=hi;i++) {
			b[i-lo]=(s->y[i]-ModelValueD(s->x[i],m,ans))*siginv/uy[i];
			for(j=0;j<np;j++) A[np*(i-lo)+j]=ans[j]*siginv/uy[i]; }
	transM(A,Atrans,n,np);
	dotMM(Atrans,A,alpha,np,n,np);
	dotMV(Atrans,b,beta,np,n);
	copyM(alpha,alphalu,np,np);
	det=LUdecomp(alphalu,np,&luindx);
	if(det)	{
		copyV(beta,ans,np);
		LUsolveV(alphalu,luindx,ans,np);
		LUimprV(alpha,alphalu,luindx,ans,beta,np);
		for(i=0;i<np;i++) {
			beta[i]=*p[i];
			ans[i]+=*p[i];
			if(ans[i]<plo[i]) *p[i]=plo[i];
			else if(ans[i]>phi[i]) *p[i]=phi[i];
			else *p[i]=ans[i]; }
		UpdateModel(m);
		if(!(chisq(m,ch2)<ch2))
			for(i=0;i<np;i++) *p[i]=beta[i];
		er=ModelCovar(m); }
	else er=9;

 endfunction:
	if(A) freeM(A);
	if(Atrans) freeM(Atrans);
	if(alpha) freeM(alpha);
	if(alphalu) freeM(alphalu);
	if(beta) freeV(beta);
	if(ans) freeV(ans);
	if(b) freeV(b);
	if(luindx) free(luindx);
	free(p);
	freeV(plo);
	freeV(phi);
	return er; }

int LMFit(modelptr m) {
	float *alpha,*beta,*alphab,*betab,*pwas,ch2,ch2b,lambda;
	float **p,det,*plo,*phi;
	int np,it,j,*luindx,er;

	er=0;
	np=prefit(m,&p,&plo,&phi);
	if(np<0) return -np;
	else if(!np) return 0;
	alpha=allocM(np,np);
	beta=allocV(np);
	alphab=allocM(np,np);
	betab=allocV(np);
	pwas=allocV(np);
	luindx=NULL;
	if(!alpha||!alphab||!beta||!betab||!pwas) {er=1;goto endfunction;}

	ch2=ChiSqD(m,-1,alpha,beta,np);
	copyM(alpha,alphab,np,np);
	copyV(beta,betab,np);
	lambda=0.001;
	for(it=0;it<np*10&&!inkey();it++) {
		for(j=0;j<np;j++)	alphab[np*j+j]*=(1.0+lambda);
		det=LUdecomp(alphab,np,&luindx);
		if(det==0) goto endfunction;
		LUsolveV(alphab,luindx,betab,np);
		free(luindx);
		luindx=NULL;
		for(j=0;j<np;j++) {
			pwas[j]=*p[j];
			if(betab[j]<FLT_MAX&&betab[j]>-FLT_MAX) {
				*p[j]+=betab[j];
				if(*p[j]<plo[j]) *p[j]=plo[j];
				else if(*p[j]>phi[j])	*p[j]=phi[j]; }}
		UpdateModel(m);
		if((ch2b=ChiSqD(m,ch2,alphab,betab,np))>=0&&ch2b<ch2) {
			lambda*=0.1;
			it=0;
			ch2=ch2b;
			copyM(alphab,alpha,np,np);
			copyV(betab,beta,np); }
		else {
			for(j=0;j<np;j++)	*p[j]=pwas[j];
			UpdateModel(m);
			lambda*=10;
			copyM(alpha,alphab,np,np);
			copyV(beta,betab,np); }}
	er=ModelCovar(m);

 endfunction:
	free(p);
	freeV(plo);
	freeV(phi);
	if(alpha) freeM(alpha);
	if(alphab) freeM(alphab);
	if(beta) freeV(beta);
	if(betab) freeV(betab);
	if(pwas) freeV(pwas);
	if(luindx) free(luindx);
	return er; }

int RandMultiFit(set s) {
	float ch2,f,f2;
	float **p,**pm,*plo,*phi,*plom,*phim,*step;
	int ntotal,np,npm,it,j,i;
	struct scell *trace;
	void *key,*item;
	modelptr m;

	if(!s) return 6;
	ntotal=0;
	for(trace=SetNext(NULL,&key,&item,s);trace;trace=SetNext(trace,&key,&item,s)) {
		m=(modelptr)item;
		if(!m) return 6;
		ntotal+=m->np; }
	if(!ntotal) return 0;
	p=(float **) calloc(ntotal,sizeof(float *));
	plo=allocV(ntotal);
	phi=allocV(ntotal);
	step=allocV(ntotal);
	if(!p||!plo||!phi||!step) {
		if(p) free(p);
		if(plo) freeV(plo);
		if(phi) freeV(phi);
		if(step) freeV(step);
		return 1; }
	np=0;
	for(trace=SetNext(NULL,&key,&item,s);trace;trace=SetNext(trace,&key,&item,s))	{
		m=(modelptr)item;
		npm=prefit(m,&pm,&plom,&phim);
		if(npm<0) {free(p);freeV(step);freeV(plo);freeV(phi);return -npm;}
		if(npm>0)	{
			for(i=0;i<npm;i++)	{
				p[np+i]=pm[i];
				plo[np+i]=plom[i];
				phi[np+i]=phim[i];
				step[np+i]=fabs(*pm[i]/100)+0.001; }
			free(pm);
			freeV(plom);
			freeV(phim);
			np+=npm; }
		if(m->covar) {
			freeM(m->covar);
			m->covar=NULL; }}

	ch2=0;
	for(trace=SetNext(NULL,&key,&item,s);trace;trace=SetNext(trace,&key,&item,s))
		ch2+=chisq((modelptr)item,-1);
	if(np)
		for(it=0;it<np*20&&!inkey();it++)	{
			j=intrand(np);
			f=*p[j];
			*p[j]+=binomrand(3,0,step[j]);
			if(*p[j]<plo[j]) *p[j]=plo[j];
			else if(*p[j]>phi[j]) *p[j]=phi[j];
			for(trace=SetNext(NULL,&key,&item,s);trace;trace=SetNext(trace,&key,&item,s))
				UpdateModel((modelptr)item);
			f2=0;
			for(trace=SetNext(NULL,&key,&item,s);trace&&ch2>f2;trace=SetNext(trace,&key,&item,s))
				f2+=chisq((modelptr)item,ch2-f2);
			if(f2<ch2) {
				it=0;
				step[j]*=1.2;
				ch2=f2; }
			else {
				*p[j]=f;
				for(trace=SetNext(NULL,&key,&item,s);trace;trace=SetNext(trace,&key,&item,s))
					UpdateModel((modelptr)item);
				step[j]*=0.9; }}
	free(p);
	freeV(step);
	freeV(plo);
	freeV(phi);
	return 0; }