#include "structs.h" #include #include #include int istrunc (char *objtype); int findpos (char *haystack, char *needle); int blankstr (char *instring); void conv2lower (char *par); void isort (unsigned long n, int arr[]); int par2num (char *parname); int read_cons (struct inpars input, struct cons *constr, struct fitpars *fpar) { int parconv (struct cons *constr, char *par, int comp, struct fitpars *fpar, int compnum); FILE *cfile; char instring[300], strptr, par[10]; char ops[] = "*/+-_"; float *f, *arr; char *p, optype; size_t len; struct fitpars *tfpar; char tok1[100], tok2[10], tok3[10], tok4[10], tok5[10], tok6[10]; int op_pos, com_pos, ncons, pos, comp, ncomps=0, comp1, comp2; float temp; constr->next = NULL; cfile = fopen (input.constraints, "r"); if (cfile == (FILE *) 0) return (1); tfpar = fpar; while (tfpar != NULL) { ncomps++; tfpar = tfpar->next; }; ncons = 0; while (fgets (instring, 300, cfile) != NULL) { conv2lower (instring); /**************************************************************\ * Don't be fooled when the op symbol comes after the comment * * field. * \**************************************************************/ com_pos = findpos (instring, "#"); op_pos = findpos (instring, ops); /**********************************************************\ * Create space to hold more constraint blocks as needed. * \**********************************************************/ if (ncons >= 1 && !blankstr(instring)) { constr->next = (struct cons *) malloc ((size_t)(sizeof (struct cons))); constr = constr->next; constr->next = NULL; }; /***************************************************************\ * Figure out whether the constraint is to relate 2 or more * * parameters or to deal with just a single one. * \***************************************************************/ if (com_pos != 1 && !blankstr(instring)) { /* Read in the first 6 tokens and parse them for an operator */ if ((p = strtok (instring, " \t")) != NULL) strcpy (tok1, p); else return (2); if ((p = strtok (NULL, " \t")) != NULL) strcpy (tok2, p); else return (2); if ((p = strtok (NULL, " \t")) != NULL) strcpy (tok3, p); else return (2); p = strtok (NULL, " \t"); if (p != NULL) strcpy (tok4, p); p = strtok (NULL, " \t"); if (p != NULL) strcpy (tok5, p); p = strtok (NULL, " \t"); if (p != NULL) strcpy (tok6, p); strcpy (constr->parname, tok2); /**************************************************************\ * Find operator in first token. If found, the pararameters * * are coupled together somehow. Dice it up to find out what * * is the relationship between the parameters. * \**************************************************************/ pos = findpos (tok1, ops); if (pos != 0) { constr->op = tok1[pos-1]; p = strtok (tok1, ops); comp1 = atoi (p); p = strtok (NULL, ops); comp2 = atoi (p); if (comp1 > ncomps || comp2 > ncomps) return (3); if (comp1 <= ncomps && comp2 <= ncomps) { constr->comp = (int *) calloc (2+OFFSET, (size_t)(sizeof(constr->comp))); constr->par = (int *) calloc (2+OFFSET, (size_t)(sizeof(constr->par))); constr->cval = (float *) calloc (2+OFFSET, (size_t)(sizeof(constr->cval))), constr->ncomp = 2; constr->comp[1] = comp1; parconv (constr, tok2, constr->comp[1], fpar, 1); constr->comp[2] = comp2; parconv (constr, tok2, constr->comp[2], fpar, 2); /**********************************\ * Soft, 2 parameter, constraints * \**********************************/ if (constr->op != '_') { constr->cval[1] = atof (tok3); constr->cval[2] = atof (tok4); } /**********************************************************\ * Hard, multi-parameter constraints can relate more than * * one parameter at a time. Figure out how many. Note * * that the first two already have storage spaces. * \**********************************************************/ if (constr->op == '_') { if (strncmp (tok3, "offset", 6) != 0 && strncmp (tok3, "ratio", 5) != 0) return (2); while ((p = strtok (NULL, "_ \t"))!= NULL) { comp = atoi(p); if (comp > ncomps) return (3); constr->ncomp++; constr->comp = (int *) realloc (constr->comp, (size_t)(sizeof (int) * (constr->ncomp + OFFSET))); constr->par = (int *) realloc (constr->par, (size_t)(sizeof (int) * (constr->ncomp + OFFSET))); constr->cval = (float *) realloc (constr->cval, (size_t)(sizeof (int) * (constr->ncomp + OFFSET))); constr->comp[constr->ncomp] = comp; parconv (constr, tok2, comp, fpar, constr->ncomp); }; /******************************************************\ * Use the numerically lowest component to lead. * * It doesn't matter which one is used at all, but * * algorithmically, this is the easiest to implement. * \******************************************************/ isort (constr->ncomp, constr->comp); isort (constr->ncomp, constr->par); if (strncmp (tok3, "offset", 6) == 0) constr->op = 'o'; else if (strncmp (tok3, "ratio", 5) == 0) constr->op = 'r'; else return (2); }; }; }; /***************************************************\ * No operator found. Parameters are not coupled. * * Soft constraints. * \***************************************************/ if (pos == 0) { comp = atoi (tok1); if (comp > 0 && comp <= ncomps) { constr->par = (int *) calloc (1 + OFFSET, (size_t)(sizeof(constr->par))); constr->cval = (float *) calloc (2 + OFFSET, (size_t)(sizeof(constr->cval))); constr->ncomp = 1; constr->op = 0; parconv (constr, tok2, comp, fpar, 1); constr->cval[1] = atof (tok3); if (findpos (tok4, "to") != 0) { constr->op = 't'; /* Constrain by min max */ constr->cval[2] = atof (tok5); } else constr->cval[2] = atof (tok4); } else return (2); }; /*****************************************************\ * Make sure the constraints are in the right order. * \*****************************************************/ if (constr->cval[1] > constr->cval[2]) { temp = constr->cval[1]; constr->cval[1] = constr->cval[2]; constr->cval[2] = temp; }; ncons++; }; }; return (0); } /****************************************************************************\ * Convert parameters from names like mag, re, etc. to parameter numbers * \****************************************************************************/ int parconv (struct cons *constr, char *par, int comp, struct fitpars *fpar, int compnum) { int parnumber, i, npars=0, quit; struct fitpars *ptr, *high; /* Convert parameter name to an actual number */ parnumber = par2num (par); /********************************************************************\ * Figure out what this parameter corresponds to in the linearlized * * array that mrqmin needs. * \********************************************************************/ ptr = fpar; npars = 0; for (i=1; i <= comp; i++) { quit = 0; /************************************************************\ * If this is a high order mode or the index is not at the * * component of interest.... * \************************************************************/ if ( i != comp || (strncmp (par, "f", 1) == 0 || (strncmp (par, "r", 1) == 0 && par[1] >= '0' && par[1] <= '9') || strncmp (par, "c", 1) == 0)) npars += ptr->npars; else if (i == comp) quit = 1; high = ptr->high; while (high != NULL && !quit) { /********************************************************\ * If at the component of interest and the parameter is * * a high order mode, go down the list.... * \********************************************************/ if (i == comp && ((strncmp (par, "f", 1) == 0 && strncmp(high->objtype, "four", 4)== 0) || ((strncmp (par, "r", 1) == 0 && par[1]>='0' && par[1]<='9') && (strncmp (high->objtype, "log", 3) == 0 || strncmp (high->objtype, "tahn", 4) == 0)) || (strncmp (par, "c", 1) == 0 && strncmp (high->objtype, "C0", 2) == 0))) quit = 1; else { /***********************************************************\ * If not at the parameter of interest, keep advancing. * * Keep in mind that the truncation parameters stored in * * the "fpar->high" structs are only pointers which are * * referencing the truncation components in the main * * object struct list. So "fpar.high.objtype=radial" pars * * should not get counted again. * \***********************************************************/ if (!istrunc (high->objtype)) npars += high->npars; high = high->next; }; }; ptr = ptr->next; }; constr->par[compnum] = npars+parnumber; return (parnumber); }