Sync librpf
[openmx:openmx.git] / src / libifa-rpf.c
1 /*
2   Copyright 2012-2013 Joshua Nathaniel Pritikin and contributors
3
4   libifa-rpf is free software: you can redistribute it and/or modify
5   it under the terms of the GNU General Public License as published by
6   the Free Software Foundation, either version 3 of the License, or
7   (at your option) any later version.
8
9   This program is distributed in the hope that it will be useful,
10   but WITHOUT ANY WARRANTY; without even the implied warranty of
11   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12   GNU General Public License for more details.
13
14   You should have received a copy of the GNU General Public License
15   along with this program.  If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include <R.h>
19 #include <R_ext/BLAS.h>
20 #include <stdlib.h>
21 #include <math.h>
22 #include <string.h>
23 #include "libifa-rpf.h"
24
25 #ifndef M_LN2
26 #define M_LN2           0.693147180559945309417232121458        /* ln(2) */
27 #endif
28
29 #ifndef M_LN_SQRT_PI
30 #define M_LN_SQRT_PI    0.572364942924700087071713675677        /* log(sqrt(pi))
31                                                                    == log(pi)/2 */
32 #endif
33
34 static const double EXP_STABLE_DOMAIN = 35;
35
36 // This is far away from the item's difficulty so it is less
37 // interesting for estimation and the gradient becomes numerically
38 // unstable as athb < -25. For symmetry, it is necessary to clamp at
39 // 25 as well.
40 static const double GRADIENT_STABLE_DOMAIN = 25;
41
42 static void
43 irt_rpf_logprob_adapter(const double *spec,
44                         const double *restrict param, const double *restrict th,
45                         double *restrict out)
46 {
47   (*librpf_model[(int) spec[RPF_ISpecID]].prob)(spec, param, th, out);
48
49   int numOutcomes = spec[RPF_ISpecOutcomes];
50   for (int ox=0; ox < numOutcomes; ox++) {
51     out[ox] = log(out[ox]);
52   }
53 }
54
55 static double
56 dotprod(const double *v1, const double *v2, const int len)
57 {
58   double dprod = 0;
59   for (int dx=0; dx < len; dx++) {
60     dprod += v1[dx] * v2[dx];
61   }
62   return dprod;
63 }
64
65 static int
66 hessianIndex(int numParam, int row, int col)
67 {
68   return numParam + row*(row+1)/2 + col;
69 }
70
71 /**
72  * lognormal(x,sd) := (x*sd*sqrt(2*%pi))^-1 * exp(-(log(x)^2)/(2*sd^2))
73  * log(lognormal(x,sd)), logexpand=super;
74  */
75 static double lognormal_pdf(double aa, double sd)
76 {
77   double sd2 = sd*sd;
78   double loga = log(aa);
79   return -loga*loga/(2*sd2) - loga - log(sd) - M_LN_SQRT_PI - M_LN2/2;
80 }
81
82 /**
83  * diff(log(lognormal(a,sd))),a), logexpand=super;
84  */
85 static double lognormal_gradient(double aa, double sd)
86 {
87   double sd2 = sd*sd;
88   return -log(aa)/(aa * sd2) - 1/aa;
89 }
90
91 static double logit(double prob)
92 {
93   return log(prob/(1-prob));
94 }
95
96 // Prior for the lower asymptote parameter (Cai, Yang & Hansen, 2011, p. 246)
97
98 /**
99  * normal(x,mean,sd) := (sd*sqrt(2*%pi))^-1 * exp(-((x - mean)^2)/(2*sd^2))
100  * log(normal(log(x/(1-x)),mean,sd)), logexpand=super;
101  */
102 static double logitnormal_pdf(double cc, double mean, double sd)
103 {
104   double sd2 = sd * sd;
105   double offset = logit(cc) - mean;
106   return -(offset*offset)/(2*sd2) - log(sd) - M_LN_SQRT_PI - M_LN2/2;
107 }
108
109 /**
110  * factor(diff(log(normal(log(x/(1-x)),mean,sd)),x))
111  */
112 static double logitnormal_gradient(double cc, double mean, double sd)
113 {
114   double sd2 = sd * sd;
115   return (logit(cc) - mean) / (sd2 * (cc-1) * cc);
116 }
117
118 static int
119 irt_rpf_1dim_drm_numSpec(const double *spec)
120 { return RPF_ISpecCount; }
121
122 static int
123 irt_rpf_1dim_drm_numParam(const double *spec)
124 { return 4; }
125
126 static void
127 irt_rpf_1dim_drm_prob(const double *spec,
128                       const double *restrict param, const double *restrict th,
129                       double *restrict out)
130 {
131   double guessing = param[2];
132   double upper = param[3];
133   double athb = -param[0] * (th[0] - param[1]);
134   if (athb < -EXP_STABLE_DOMAIN) athb = -EXP_STABLE_DOMAIN;
135   else if (athb > EXP_STABLE_DOMAIN) athb = EXP_STABLE_DOMAIN;
136   double pp = guessing + (upper-guessing) / (1 + exp(athb));
137   out[0] = 1-pp;
138   out[1] = pp;
139 }
140
141 static void
142 set_deriv_nan(const double *spec, double *out)
143 {
144   int numParam = (*librpf_model[(int) spec[RPF_ISpecID]].numParam)(spec);
145
146   for (int px=0; px < numParam; px++) {
147     out[px] = nan("I");
148   }
149 }
150
151 static void
152 irt_rpf_1dim_drm_rescale(const double *spec, double *restrict param, const int *paramMask,
153                          const double *restrict mean, const double *restrict cov)
154 {
155   double thresh = param[1] * -param[0];
156   if (paramMask[0] >= 0) {
157     param[0] *= cov[0];
158   }
159   if (paramMask[1] >= 0) {
160     thresh += param[0] * mean[0];
161     param[1] = thresh / -param[0];
162   }
163 }
164
165 static int
166 irt_rpf_mdim_drm_numSpec(const double *spec)
167 { return RPF_ISpecCount; }
168
169 static int
170 irt_rpf_mdim_drm_numParam(const double *spec)
171 { return 3 + spec[RPF_ISpecDims]; }
172
173 static void
174 irt_rpf_mdim_drm_paramInfo(const double *spec, const int param,
175                            int *type, double *upper, double *lower)
176 {
177         int numDims = spec[RPF_ISpecDims];
178         *upper = nan("unset");
179         *lower = nan("unset");
180         if (param >= 0 && param < numDims) {
181                 *type = RPF_Slope;
182                 *lower = 1e-6;
183         } else if (param == numDims) {
184                 *type = RPF_Intercept;
185         } else if (param == numDims+1 || param == numDims+2) {
186                 *type = RPF_Bound;
187                 *lower = 1e-6;
188                 *upper = 1 - 1e-6;
189         }
190 }
191
192 static void
193 irt_rpf_mdim_drm_prob(const double *spec,
194                       const double *restrict param, const double *restrict th,
195                       double *restrict out)
196 {
197   int numDims = spec[RPF_ISpecDims];
198   double dprod = dotprod(param, th, numDims);
199   double diff = param[numDims];
200   double guessing = param[numDims+1];
201   double upper = param[numDims+2];
202   double athb = -(dprod + diff);
203   if (athb < -EXP_STABLE_DOMAIN) athb = -EXP_STABLE_DOMAIN;
204   else if (athb > EXP_STABLE_DOMAIN) athb = EXP_STABLE_DOMAIN;
205   double tmp = guessing + (upper-guessing) / (1 + exp(athb));
206   out[0] = 1-tmp;
207   out[1] = tmp;
208 }
209
210 static void
211 irt_rpf_mdim_drm_prob2(const double *spec,
212                        const double *restrict param, const double *restrict th,
213                        double *restrict out1, double *restrict out2)
214 {
215   int numDims = spec[RPF_ISpecDims];
216   double dprod = dotprod(param, th, numDims);
217   double diff = param[numDims];
218   double guessing = param[numDims+1];
219   double upper = param[numDims+2];
220   double athb = -(dprod + diff);
221   if (athb < -EXP_STABLE_DOMAIN) athb = -EXP_STABLE_DOMAIN;
222   else if (athb > EXP_STABLE_DOMAIN) athb = EXP_STABLE_DOMAIN;
223   double tmp = 1 / (1 + exp(athb));
224   out1[0] = 1-tmp;
225   out1[1] = tmp;
226   tmp = guessing + (upper-guessing) * tmp;
227   out2[0] = 1-tmp;
228   out2[1] = tmp;
229 }
230
231 static void
232 irt_rpf_mdim_drm_deriv1(const double *spec,
233                        const double *restrict param,
234                        const double *where,
235                        const double *weight, double *out)
236 {
237   int numDims = spec[RPF_ISpecDims];
238   double cc = param[numDims+1];
239   double upper = param[numDims+2];
240   double PQ[2];
241   double PQstar[2];
242   irt_rpf_mdim_drm_prob2(spec, param, where, PQstar, PQ);
243   double r1 = weight[0];
244   double r2 = weight[1];
245   double r1_P = r1/PQ[0];
246   double r1_P2 = r1/(PQ[0] * PQ[0]);
247   double r2_Q = r2/PQ[1];
248   double r2_Q2 = r2/(PQ[1] * PQ[1]);
249   for (int dx=0; dx < numDims; dx++) {
250     out[dx] += where[dx] * PQstar[0] * PQstar[1] * (upper-cc) * (r1_P - r2_Q);
251   }
252   out[numDims] += (upper-cc) * PQstar[0] * PQstar[1] * (r1_P - r2_Q);
253   out[numDims+1] += PQstar[0] * (r1_P - r2_Q);
254   out[numDims+2] += PQstar[1] * (r1_P - r2_Q);
255
256   int ox = numDims+2;
257   double ugD2 = (upper-cc);  // can factor out TODO
258   double ugD = (upper-cc);
259   double Pstar = PQstar[0];
260   double Pstar2 = Pstar * Pstar;
261   double Pstar3 = Pstar2 * Pstar;
262   double Qstar = PQstar[1];
263   double Qstar2 = Qstar * Qstar;
264
265   for(int ix=0; ix < numDims; ix++) {
266     for(int jx=0; jx <= ix; jx++) {
267       out[++ox] -= (r1_P * (ugD2 * where[ix] * where[jx] *
268                                    (Pstar - 3*Pstar2 + 2*Pstar3)) -
269                            r1_P2 * (ugD * where[ix] * (Pstar - Pstar2) *
270                                     (ugD * where[jx] * (Pstar - Pstar2))) +
271                            r2_Q * (ugD2 * where[ix] * where[jx] *
272                                    (-Pstar + 3*Pstar2 - 2*Pstar3)) -
273                            r2_Q2 * (ugD * where[ix] * (-Pstar + Pstar2) *
274                                     (ugD * where[jx] * (-Pstar + Pstar2))));
275     }
276   }
277   for(int ix=0; ix < numDims; ix++) {
278     out[++ox] -= (r1_P * (ugD2 * where[ix] * (Pstar - 3*Pstar2 + 2*Pstar3)) -
279                          r1_P2 * (ugD * where[ix] * (Pstar - Pstar2) *
280                                   (ugD * (Pstar - Pstar2))) +
281                          r2_Q * (ugD2 * where[ix] * (-Pstar + 3*Pstar2 - 2*Pstar3)) -
282                          r2_Q2 * (ugD * where[ix] * (-Pstar + Pstar2) *
283                                   (ugD * (-Pstar + Pstar2))));
284   }
285   double chunk1 = ugD * (Pstar - Pstar2);
286   out[++ox] -= (r1_P * (ugD2 * (Pstar - 3*Pstar2 + 2*Pstar3)) -
287                        r1_P2 * chunk1*chunk1 +
288                        r2_Q * (ugD2 * (-Pstar + 3*Pstar2 - 2*Pstar3)) -
289                        r2_Q2 * chunk1*chunk1);
290   for(int ix=0; ix < numDims; ix++) {
291     out[++ox] -= (r1_P * (where[ix] * (Pstar - Pstar2)) -
292                          r1_P2 * (ugD * where[ix] * (Pstar - Pstar2)) * Pstar +
293                          r2_Q * (where[ix] * (-Pstar + Pstar2)) +
294                          r2_Q2 * (ugD * where[ix] * (-Pstar + Pstar2) ) * Pstar);
295   }
296   out[++ox] -= (r1_P * (Pstar - Pstar2) -
297                        r1_P2 * (ugD * (Pstar - Pstar2)) * Pstar +
298                        r2_Q * (-Pstar + Pstar2) +
299                        r2_Q2 * (ugD * (-Pstar + Pstar2)) * Pstar);
300   out[++ox] += Pstar2 * (r1_P2 + r2_Q2);
301
302   for(int ix=0; ix < numDims; ix++) {
303     out[++ox] -= (r1_P * (where[ix] * (-Pstar + Pstar2)) -
304                          r1_P2 * (ugD * where[ix] * (Pstar - Pstar2)) * Qstar +
305                          r2_Q * (where[ix] * (Pstar - Pstar2)) -
306                          r2_Q2 * (ugD * where[ix] * (-Pstar + Pstar2) ) * (-1 + Pstar));
307   }
308
309   out[++ox] -= (r1_P * (-Pstar + Pstar2) -
310                        r1_P2 * (ugD * (Pstar - Pstar2)) * Qstar +
311                        r2_Q * (Pstar - Pstar2) -
312                        r2_Q2 * (ugD * (-Pstar + Pstar2)) * -Qstar);
313
314   out[++ox] -= (-r1_P2 * Pstar * Qstar + r2_Q2 * Pstar * (-1 + Pstar));
315   out[++ox] += Qstar2 * (r1_P2 + r2_Q2);
316 }
317
318 static void
319 irt_rpf_mdim_drm_deriv2(const double *spec,
320                         const double *restrict param,
321                         double *out)
322 {
323   int numDims = spec[RPF_ISpecDims];
324   const double *aa = param;
325   double cc = param[numDims+1];
326   double upper = param[numDims+2];
327
328   if (cc < 0 || cc >= 1 || upper <= 0 || upper > 1) {
329     set_deriv_nan(spec, out);
330     return;
331   }
332   for (int dx=0; dx < numDims; dx++) {
333     if (aa[dx] < 0) {
334       set_deriv_nan(spec, out);
335       return;
336     }
337   }
338   if (cc == 0) {
339     out[numDims+1] = nan("I");
340   }
341   if (upper == 1) {
342     out[numDims+2] = nan("I");
343   }
344 }
345
346 static void
347 irt_rpf_mdim_drm_rescale(const double *spec, double *restrict param, const int *paramMask,
348                          const double *restrict mean, const double *restrict cov)
349 {
350   int numDims = spec[RPF_ISpecDims];
351
352   double madj = dotprod(param, mean, numDims);
353
354   for (int d1=0; d1 < numDims; d1++) {
355     if (paramMask[d1] < 0) continue;
356     param[d1] = dotprod(param+d1, cov + d1 * numDims + d1, numDims-d1);
357   }
358
359   param[numDims] += madj;
360 }
361
362 static void
363 irt_rpf_mdim_drm_dTheta(const double *spec, const double *restrict param,
364                         const double *where, const double *dir,
365                         double *grad, double *hess)
366 {
367   int numDims = spec[RPF_ISpecDims];
368   double PQ[2];
369   double PQstar[2];
370   irt_rpf_mdim_drm_prob2(spec, param, where, PQstar, PQ);
371   double Pstar = PQstar[0];
372   double Qstar = PQstar[1];
373   const double *aa = param;
374   const double guess = param[numDims + 1];
375   const double upper = param[numDims + 2];
376   for (int ax=0; ax < numDims; ax++) {
377     double piece = dir[ax] * (upper-guess) * aa[ax] * (Pstar * Qstar);
378     grad[1] += piece;
379     grad[0] -= piece;
380     piece = dir[ax] * (2 * (upper - guess) * aa[ax]*aa[ax] * (Qstar * Qstar * Pstar) -
381                        (upper - guess) * aa[ax]*aa[ax] * (Pstar * Qstar));
382     hess[1] -= piece;
383     hess[0] += piece;
384   }
385 }
386
387 static void
388 irt_rpf_1dim_drm_dTheta(const double *spec, const double *restrict param,
389                         const double *where, const double *dir,
390                         double *grad, double *hess)
391 {
392   double nparam[4];
393   memcpy(nparam, param, sizeof(double) * 4);
394   nparam[1] = param[1] * -param[0];
395   irt_rpf_mdim_drm_dTheta(spec, nparam, where, dir, grad, hess);
396 }
397
398 static int
399 irt_rpf_mdim_grm_numSpec(const double *spec)
400 { return RPF_ISpecCount; }
401
402 static int
403 irt_rpf_mdim_grm_numParam(const double *spec)
404 { return spec[RPF_ISpecOutcomes] + spec[RPF_ISpecDims] - 1; }
405
406 static void
407 irt_rpf_mdim_grm_paramInfo(const double *spec, const int param,
408                            int *type, double *upper, double *lower)
409 {
410         int numDims = spec[RPF_ISpecDims];
411         const int numOutcomes = spec[RPF_ISpecOutcomes];
412         *upper = nan("unset");
413         *lower = nan("unset");
414         if (param >= 0 && param < numDims) {
415                 *type = RPF_Slope;
416                 *lower = 1e-6;
417         } else {
418                 *type = RPF_Intercept;
419         }
420 }
421
422 static void
423 irt_rpf_mdim_grm_prob(const double *spec,
424                       const double *restrict param, const double *restrict th,
425                       double *restrict out)
426 {
427   const int numDims = spec[RPF_ISpecDims];
428   const int numOutcomes = spec[RPF_ISpecOutcomes];
429   const double *slope = param;
430   const double dprod = dotprod(slope, th, numDims);
431   const double *kat = param + (int) spec[RPF_ISpecDims];
432
433   double athb = -(dprod + kat[0]);
434   if (athb < -EXP_STABLE_DOMAIN) athb = -EXP_STABLE_DOMAIN;
435   else if (athb > EXP_STABLE_DOMAIN) athb = EXP_STABLE_DOMAIN;
436   double tmp = 1 / (1 + exp(athb));
437   out[0] = 1-tmp;
438   out[1] = tmp;
439
440   for (int kx=2; kx < numOutcomes; kx++) {
441     double athb = -(dprod + kat[kx-1]);
442     if (athb < -EXP_STABLE_DOMAIN) athb = -EXP_STABLE_DOMAIN;
443     else if (athb > EXP_STABLE_DOMAIN) athb = EXP_STABLE_DOMAIN;
444     double tmp = 1 / (1 + exp(athb));
445     out[kx-1] -= tmp;
446     out[kx] = tmp;
447   }
448
449   // look for crazy stuff
450   for (int kx=0; kx < numOutcomes; kx++) {
451     if (out[kx] <= 0) {
452       int bigk = -1;
453       double big = 0;
454       for (int bx=0; bx < numOutcomes; bx++) {
455         if (out[bx] > big) {
456           bigk = bx;
457           big = out[bx];
458         }
459       }
460       for (int fx=0; fx < numOutcomes; fx++) {
461               if (out[fx] < 0) {
462                       set_deriv_nan(spec, out);
463                       return;
464               }
465         if (out[fx] == 0) {
466           double small = 1 / (1 + exp(EXP_STABLE_DOMAIN));
467           out[bigk] -= small;
468           out[fx] += small;
469         }
470       }
471       return;
472     }
473   }
474 }
475
476 static void
477 irt_rpf_mdim_grm_rawprob(const double *spec,
478                          const double *restrict param, const double *restrict th,
479                          double *restrict out)
480 {
481   int numDims = spec[RPF_ISpecDims];
482   const int numOutcomes = spec[RPF_ISpecOutcomes];
483   const double dprod = dotprod(param, th, numDims);
484   const double *kat = param + (int) spec[RPF_ISpecDims];
485
486   out[0] = 1;
487   for (int kx=0; kx < numOutcomes-1; kx++) {
488     double athb = -(dprod + kat[kx]);
489     if (athb < -EXP_STABLE_DOMAIN) athb = -EXP_STABLE_DOMAIN;
490     else if (athb > EXP_STABLE_DOMAIN) athb = EXP_STABLE_DOMAIN;
491     double tmp = 1 / (1 + exp(athb));
492     out[kx+1] = tmp;
493   }
494   out[numOutcomes] = 0;
495 }
496
497 static void
498 irt_rpf_mdim_grm_deriv1(const double *spec,
499                         const double *restrict param,
500                         const double *where,
501                         const double *weight, double *out)
502 {
503   int nfact = spec[RPF_ISpecDims];
504   int outcomes = spec[RPF_ISpecOutcomes];
505   int nzeta = spec[RPF_ISpecOutcomes] - 1;
506   double P[nzeta+2];
507   double PQfull[nzeta+2];
508   irt_rpf_mdim_grm_rawprob(spec, param, where, P);
509   PQfull[0] = 0;
510   PQfull[outcomes] = 0;
511   for (int kx=1; kx <= nzeta; kx++) PQfull[kx] = P[kx] * (1-P[kx]);
512   for (int jx = 0; jx <= nzeta; jx++) {
513     double Pk_1 = P[jx];
514     double Pk = P[jx + 1];
515     double PQ_1 = PQfull[jx];
516     double PQ = PQfull[jx + 1];
517     double Pk_1Pk = Pk_1 - Pk;
518     if (Pk_1Pk < 1e-10) Pk_1Pk = 1e-10;
519     double dif1 = weight[jx] / Pk_1Pk;
520     double dif1sq = weight[jx] / (Pk_1Pk * Pk_1Pk);
521     if(jx < nzeta) {
522       double Pk_p1 = P[jx + 2];
523       double PQ_p1 = PQfull[jx + 2];
524       double Pk_Pkp1 = Pk - Pk_p1;
525       if(Pk_Pkp1 < 1e-10) Pk_Pkp1 = 1e-10;
526       double dif2 = weight[jx+1] / Pk_Pkp1;
527       double dif2sq = weight[jx+1] / (Pk_Pkp1 * Pk_Pkp1);
528       out[nfact + jx] += PQ * (dif1 - dif2);
529
530       int d2base = (nfact + nzeta) + (nfact+jx) * (nfact+jx+1)/2;
531       out[d2base + nfact + jx] -= (-1 * PQ * PQ * (dif1sq + dif2sq) -
532                                           (dif1 - dif2) * (Pk * (1.0 - Pk) * (1.0 - 2.0*Pk)));
533       if (jx < (nzeta - 1)) {
534         int d2base1 = (nfact + nzeta) + (nfact+jx+1) * (nfact+jx+2)/2;
535         out[d2base1 + nfact + jx] -= dif2sq * PQ_p1 * PQ;
536       }
537       double tmp1 = (-1.0) * dif2sq * PQ * (PQ - PQ_p1);
538       double tmp2 = dif1sq * PQ * (PQ_1 - PQ);
539       double tmp3 = (dif1 - dif2) * (Pk * (1.0 - Pk) * (1.0 - 2.0*Pk));
540       for(int kx = 0; kx < nfact; kx++){
541         out[d2base + kx] -= (tmp1 + tmp2 - tmp3) * where[kx];
542       }
543     }
544     for(int kx = 0; kx < nfact; kx++) {
545       out[kx] -= dif1 * (PQ_1 - PQ) * where[kx];
546     }
547
548     double temp[nfact];
549     for(int ix = 0; ix < nfact; ix++)
550       temp[ix] = PQ_1 * where[ix] - PQ * where[ix];
551
552     int d2x = nfact + nzeta;
553     for(int i = 0; i < nfact; i++) {
554       for(int j = 0; j <= i; j++) {
555         double outer = where[i]*where[j];
556         out[d2x++] -= (-1 * dif1sq * temp[i] * temp[j] +
557                               (dif1 * (Pk_1 * (1.0 - Pk_1) * (1.0 - 2.0 * Pk_1) * 
558                                        outer - Pk * (1.0 - Pk) * (1.0 - 2.0 * Pk) * outer)));
559       }
560     }
561   }
562 }
563
564 static void
565 irt_rpf_mdim_grm_deriv2(const double *spec,
566                         const double *restrict param,
567                         double *out)
568 {
569   int nfact = spec[RPF_ISpecDims];
570   int nzeta = spec[RPF_ISpecOutcomes] - 1;
571   const double *aa = param;
572   for (int dx=0; dx < nfact; dx++) {
573     if (aa[dx] < 0) {
574       set_deriv_nan(spec, out);
575       return;
576     }
577   }
578   for (int zx=0; zx < nzeta-1; zx++) {
579     if (param[nfact+zx] < param[nfact+zx+1]) {
580       set_deriv_nan(spec, out);
581       return;
582     }
583   }
584 }
585
586 static void
587 irt_rpf_mdim_grm_dTheta(const double *spec, const double *restrict param,
588                         const double *where, const double *dir,
589                         double *grad, double *hess)
590 {
591   int numDims = spec[RPF_ISpecDims];
592   int outcomes = spec[RPF_ISpecOutcomes];
593   const double *aa = param;
594   double P[outcomes+1];
595   irt_rpf_mdim_grm_rawprob(spec, param, where, P);
596   for (int jx=0; jx < numDims; jx++) {
597     for (int ix=0; ix < outcomes; ix++) {
598       double w1 = P[ix] * (1-P[ix]) * aa[jx];
599       double w2 = P[ix+1] * (1-P[ix+1]) * aa[jx];
600       grad[ix] += dir[jx] * (w1 - w2);
601       hess[ix] += dir[jx] * (aa[jx]*aa[jx] * (2 * P[ix] * (1 - P[ix])*(1 - P[ix]) -
602                                               P[ix] * (1 - P[ix]) -
603                                               2 * P[ix+1] * (1 - P[ix+1])*(1 - P[ix+1]) +
604                                               P[ix+1] * (1 - P[ix+1])));
605     }
606   }
607 }
608
609 static void
610 irt_rpf_mdim_grm_rescale(const double *spec, double *restrict param, const int *paramMask,
611                          const double *restrict mean, const double *restrict cov)
612 {
613   int numDims = spec[RPF_ISpecDims];
614   int nzeta = spec[RPF_ISpecOutcomes] - 1;
615
616   double madj = dotprod(param, mean, numDims);
617
618   for (int d1=0; d1 < numDims; d1++) {
619     if (paramMask[d1] < 0) continue;
620     param[d1] = dotprod(param+d1, cov + d1 * numDims + d1, numDims-d1);
621   }
622
623   for (int tx=0; tx < nzeta; tx++) {
624     int px = numDims + tx;
625     if (paramMask[px] >= 0) param[px] += madj;
626   }
627 }
628
629 static int
630 irt_rpf_nominal_numSpec(const double *spec)
631 {
632   int outcomes = spec[RPF_ISpecOutcomes];
633   int Tlen = (outcomes - 1) * (outcomes - 1);
634   return RPF_ISpecCount + 4 * Tlen;
635 }
636
637 static int
638 irt_rpf_nominal_numParam(const double *spec)
639 { return spec[RPF_ISpecDims] + 2 * (spec[RPF_ISpecOutcomes]-1); }
640
641
642 static void
643 irt_rpf_nominal_paramInfo(const double *spec, const int param,
644                           int *type, double *upper, double *lower)
645 {
646         int numDims = spec[RPF_ISpecDims];
647         const int numOutcomes = spec[RPF_ISpecOutcomes];
648         *upper = nan("unset");
649         *lower = nan("unset");
650         if (param >= 0 && param < numDims) {
651                 *type = RPF_Slope;
652                 *lower = 1e-6;
653         } else if (param < numDims + numOutcomes - 1) {
654                 *type = RPF_Slope;
655         } else {
656                 *type = RPF_Intercept;
657         }
658 }
659
660 static void
661 _nominal_rawprob(const double *spec,
662                  const double *restrict param, const double *restrict th,
663                  double discr, double *ak, double *num)
664 {
665   int numDims = spec[RPF_ISpecDims];
666   int numOutcomes = spec[RPF_ISpecOutcomes];
667   const double *alpha = param + numDims;
668   const double *gamma = param + numDims + numOutcomes - 1;
669   const double *Ta = spec + RPF_ISpecCount;
670   const double *Tc = spec + RPF_ISpecCount + (numOutcomes-1) * (numOutcomes-1);
671
672   for (int kx=0; kx < numOutcomes; kx++) {
673     ak[kx] = 0;
674     double ck = 0;
675     if (kx) {
676       for (int tx=0; tx < numOutcomes-1; tx++) {
677         int Tcell = tx * (numOutcomes-1) + kx-1;
678         ak[kx] += Ta[Tcell] * alpha[tx];
679         ck += Tc[Tcell] * gamma[tx];
680       }
681     }
682
683     double z = discr * ak[kx] + ck;
684     if (z < -EXP_STABLE_DOMAIN) z = -EXP_STABLE_DOMAIN;
685     else if (z > EXP_STABLE_DOMAIN) z = EXP_STABLE_DOMAIN;
686     num[kx] = z;
687   }
688 }
689
690
691 static void
692 irt_rpf_nominal_prob(const double *spec,
693                      const double *restrict param, const double *restrict th,
694                      double *restrict out)
695 {
696   int numOutcomes = spec[RPF_ISpecOutcomes];
697   int numDims = spec[RPF_ISpecDims];
698   double num[numOutcomes];
699   double ak[numOutcomes];
700   double discr = dotprod(param, th, numDims);
701   _nominal_rawprob(spec, param, th, discr, ak, num);
702   double den = 0;
703
704   for (int kx=0; kx < numOutcomes; kx++) {
705     num[kx] = exp(num[kx]);
706     den += num[kx];
707   }
708   for (int kx=0; kx < numOutcomes; kx++) {
709     out[kx] = num[kx]/den;
710   }
711 }
712
713 static void
714 irt_rpf_nominal_logprob(const double *spec,
715                         const double *restrict param, const double *restrict th,
716                         double *restrict out)
717 {
718   int numOutcomes = spec[RPF_ISpecOutcomes];
719   int numDims = spec[RPF_ISpecDims];
720   double num[numOutcomes];
721   double ak[numOutcomes];
722   double discr = dotprod(param, th, numDims);
723   _nominal_rawprob(spec, param, th, discr, ak, num);
724   double den = 0;
725
726   for (int kx=0; kx < numOutcomes; kx++) {
727     den += exp(num[kx]);
728   }
729   den = log(den);
730
731   for (int kx=0; kx < numOutcomes; kx++) {
732     out[kx] = num[kx] - den;
733   }
734 }
735
736 static double makeOffterm(const double *dat, const double p, const double aTheta,
737                           const int ncat, const int cat)
738 {
739   double ret = 0;
740   for (int CAT = 0; CAT < ncat; CAT++) {
741     if (CAT == cat) continue;
742     ret += dat[CAT] * p * aTheta;
743   }
744   return(ret);
745 }
746
747 static double makeOffterm2(const double *dat, const double p1, const double p2, 
748                            const double aTheta, const int ncat, const int cat)
749 {
750   double ret = 0;
751   for (int CAT = 0; CAT < ncat; CAT++) {
752     if (CAT == cat) continue;
753     ret += dat[CAT] * p1 * p2 * aTheta;
754   }
755   return(ret);
756 }
757
758 static void
759 irt_rpf_nominal_deriv1(const double *spec,
760                        const double *restrict param,
761                        const double *where,
762                        const double *weight, double *out)
763 {
764   int nfact = spec[RPF_ISpecDims];
765   int ncat = spec[RPF_ISpecOutcomes];
766   double aTheta = dotprod(param, where, nfact);
767   double aTheta2 = aTheta * aTheta;
768
769   double num[ncat];
770   double ak[ncat];
771   _nominal_rawprob(spec, param, where, aTheta, ak, num);
772
773   double P[ncat];
774   double P2[ncat];
775   double P3[ncat];
776   double ak2[ncat];
777   double dat_num[ncat];
778   double numsum = 0;
779   double numakD = 0;
780   double numak2D2 = 0;
781   double numakDTheta_numsum[nfact];
782
783   for (int kx=0; kx < ncat; kx++) {
784     num[kx] = exp(num[kx]);
785     ak2[kx] = ak[kx] * ak[kx];
786     dat_num[kx] = weight[kx]/num[kx];
787     numsum += num[kx];
788     numakD += num[kx] * ak[kx];
789     numak2D2 += num[kx] * ak2[kx];
790   }
791   double numsum2 = numsum * numsum;
792
793   for (int kx=0; kx < ncat; kx++) {
794     P[kx] = num[kx]/numsum;
795     P2[kx] = P[kx] * P[kx];
796     P3[kx] = P2[kx] * P[kx];
797   }
798
799   double sumNumak = dotprod(num, ak, ncat);
800   for (int fx=0; fx < nfact; fx++) {
801     numakDTheta_numsum[fx] = sumNumak * where[fx] / numsum;
802   }
803
804   for (int jx = 0; jx < nfact; jx++) {
805     double tmpvec = 0;
806     for(int i = 0; i < ncat; i++) {
807       tmpvec += dat_num[i] * (ak[i] * where[jx] * P[i] -
808                               P[i] * numakDTheta_numsum[jx]) * numsum;
809     }
810     out[jx] -= tmpvec;
811   }
812   for(int i = 1; i < ncat; i++) {
813     double offterm = makeOffterm(weight, P[i], aTheta, ncat, i);
814     double offterm2 = makeOffterm(weight, P[i], 1, ncat, i);
815     double tmpvec = dat_num[i] * (aTheta * P[i] - P2[i] * aTheta) * numsum - offterm;
816     double tmpvec2 = dat_num[i] * (P[i] - P2[i]) * numsum - offterm2;
817     out[nfact + i - 1] -= tmpvec;
818     out[nfact + ncat + i - 2] -= tmpvec2;
819   }
820
821   int hessbase = nfact + (ncat-1)*2;
822   int d2ind = 0;
823   //a's
824   for (int j = 0; j < nfact; j++) {
825     for (int k = 0; k <= j; k++) {
826       double tmpvec = 0;
827       for (int i = 0; i < ncat; i++) {
828         tmpvec += dat_num[i] * (ak2[i] * where[j] * where[k] * P[i] -
829                                 ak[i] * where[j] * P[i] * numakDTheta_numsum[k] -
830                                 ak[i] * where[k] * P[i] * numakDTheta_numsum[j] + 
831                                 2 * P[i] * numakD * where[j] * numakD * where[k] / numsum2 -
832                                 P[i] * numak2D2 * where[j] * where[k] / numsum) * numsum - 
833           dat_num[i] * (ak[i] * where[j] * P[i] - P[i] * numakDTheta_numsum[j]) *
834           numsum * ak[i] * where[k] +
835           dat_num[i] * (ak[i] * where[j] * P[i] - P[i] * numakDTheta_numsum[j]) *
836           numakD * where[k];
837       }
838       out[hessbase + d2ind++] -= tmpvec;
839     }
840   }
841   //a's with ak and d
842   for(int k = 1; k < ncat; k++){
843     int akrow = hessbase + (nfact+k)*(nfact+k-1)/2;
844     int dkrow = hessbase + (nfact+ncat+k-1)*(nfact+ncat+k-2)/2;
845     for(int j = 0; j < nfact; j++){
846       double tmpvec = 0;
847       double tmpvec2 = 0;
848       for(int i = 0; i < ncat; i++){
849         if(i == k){
850           tmpvec += dat_num[i] * (ak[i]*where[j] * aTheta*P[i] -
851                                      aTheta*P[i]*numakDTheta_numsum[j] +
852                                      where[j]*P[i] - 2*ak[i]*where[j]*aTheta*P2[i] +
853                                      2*aTheta*P2[i]*numakDTheta_numsum[j] -
854                                      where[j]*P2[i])*numsum -
855             dat_num[i]*(aTheta*P[i] - aTheta*P2[i])*numsum*ak[i]*where[j] +
856             dat_num[i]*(aTheta*P[i] - aTheta*P2[i])*(numakD*where[j]);
857           tmpvec2 += dat_num[i]*(ak[i]*where[j]*P[i] -
858                                       2*ak[i]*where[j]*P2[i] -
859                                       P[i]*numakDTheta_numsum[j] +
860                                       2*P2[i]*numakDTheta_numsum[j])*numsum -
861             dat_num[i]*(P[i] - P2[i])*numsum*ak[i]*where[j] +
862             dat_num[i]*(P[i] - P2[i])*(numakD*where[j]);
863         } else {
864           tmpvec += -weight[i]*ak[k]*aTheta*where[j]*P[k] +
865             weight[i]*P[k]*aTheta*numakDTheta_numsum[j] -
866             weight[i]*P[k]*where[j];
867           tmpvec2 += -weight[i]*ak[k]*where[j]*P[k] +
868             weight[i]*P[k]*numakDTheta_numsum[j];
869         }
870       }
871       out[akrow + j] -= tmpvec;
872       out[dkrow + j] -= tmpvec2;
873     }
874   }
875   //ak's and d's
876   for(int j = 1; j < ncat; j++){
877     int akrow = hessbase + (nfact+j)*(nfact+j-1)/2;
878     int dkrow = hessbase + (nfact+ncat+j-1)*(nfact+ncat+j-2)/2;
879
880     double tmpvec = makeOffterm(weight, P2[j], aTheta2, ncat, j);
881     double tmpvec2 = makeOffterm(weight, P[j], aTheta2, ncat, j);
882     double offterm = tmpvec - tmpvec2;
883     tmpvec = makeOffterm(weight, P2[j], 1, ncat, j);
884     tmpvec2 = makeOffterm(weight, P[j], 1, ncat, j);
885     double offterm2 = tmpvec - tmpvec2;
886
887     out[akrow + nfact + j - 1] -=
888       (dat_num[j]*(aTheta2*P[j] - 3*aTheta2*P2[j] +
889                           2*aTheta2*P3[j])*numsum - weight[j]/num[j] *
890               (aTheta*P[j] - aTheta*P2[j])*numsum*aTheta + weight[j] *
891               (aTheta*P[j] - aTheta*P2[j])*aTheta + offterm);
892
893     out[dkrow + nfact + ncat + j - 2] -=
894       (dat_num[j]*(P[j] - 3*P2[j] + 2*P3[j])*numsum - weight[j]/num[j] *
895               (P[j] - P2[j])*numsum + weight[j] *
896               (P[j] - P2[j]) + offterm2);
897
898     for(int i = 1; i < ncat; i++) {
899       if(j > i) {
900         offterm = makeOffterm2(weight, P[j], P[i], aTheta2, ncat, i);
901         offterm2 = makeOffterm2(weight, P[j], P[i], 1, ncat, i);
902         tmpvec = dat_num[i] * (-aTheta2*P[i]*P[j] + 2*P2[i] *aTheta2*P[j])*numsum + 
903           dat_num[i] * (aTheta*P[i] - P2[i] * aTheta)*aTheta*num[j]+offterm;
904         tmpvec2 = dat_num[i] * (-P[i]*P[j] + 2*P2[i] *P[j]) * numsum +
905           dat_num[i] * (P[i] - P2[i]) * num[j] + offterm2;
906         out[akrow + nfact + i - 1] -= tmpvec;
907         out[dkrow + nfact + ncat + i - 2] -= tmpvec2;
908       }
909       if (abs(j-i) == 0) {
910         tmpvec = makeOffterm(weight, P2[i], aTheta, ncat, i);
911         tmpvec2 = makeOffterm(weight, P[i], aTheta, ncat, i);
912         offterm = tmpvec - tmpvec2;
913         tmpvec = dat_num[i]*(aTheta*P[i] - 3*aTheta*P2[i] +
914                              2*aTheta*P3[i]) * numsum - dat_num[i] *
915           (aTheta*P[i] - aTheta*P2[i])*numsum + weight[i] *
916           (P[i] - P2[i])*aTheta + offterm;
917         out[dkrow + nfact + i - 1] -= tmpvec;
918       } else {
919         offterm = makeOffterm2(weight, P[j], P[i], aTheta, ncat, i);
920         tmpvec = dat_num[i] * (-aTheta*P[i]*P[j] + 2*P2[i] *aTheta*P[j]) * numsum + 
921           dat_num[i] * (P[i] - P2[i]) * aTheta * num[j] + offterm;
922         out[dkrow + nfact + i - 1] -= tmpvec;
923       }
924     }
925   }
926 }
927
928 static void
929 pda(const double *ar, int rows, int cols) {   // column major order
930         for (int rx=0; rx < rows; rx++) {
931                 for (int cx=0; cx < cols; cx++) {
932                         Rprintf("%.6g, ", ar[cx * rows + rx]);
933                 }
934                 Rprintf("\n");
935         }
936 }
937
938 static void
939 irt_rpf_nominal_deriv2(const double *spec,
940                        const double *restrict param,
941                        double *out)
942 {
943   int nfact = spec[RPF_ISpecDims];
944   int nzeta = spec[RPF_ISpecOutcomes] - 1;
945   const double *aa = param;
946
947   for (int dx=0; dx < nfact; dx++) {
948     if (aa[dx] < 0) {
949       set_deriv_nan(spec, out);
950       return;
951     }
952   }
953
954   const double *Ta = spec + RPF_ISpecCount;
955   const double *Tc = spec + RPF_ISpecCount + nzeta * nzeta;
956   const int numParam = irt_rpf_nominal_numParam(spec);
957   double rawOut[numParam];
958   memcpy(rawOut, out, sizeof(double) * numParam);
959
960   // gradient
961   for (int tx=0; tx < nzeta; tx++) {
962     double ak1=0;
963     double ck1=0;
964     for (int kx=0; kx < nzeta; kx++) {
965       int Tcell = tx * nzeta + kx;
966       ak1 += rawOut[nfact + kx] * Ta[Tcell];
967       ck1 += rawOut[nfact + nzeta + kx] * Tc[Tcell];
968     }
969     out[nfact + tx] = ak1;
970     out[nfact + nzeta + tx] = ck1;
971   }
972
973   // don't need to transform the main a parameters TODO
974   double *dmat = Realloc(NULL, 3 * numParam * numParam, double);
975   const int hsize = hessianIndex(0, numParam-1, numParam-1);
976   {
977     int row=0;
978     int col=0;
979     for (int dx=0; dx <= hsize; dx++) {
980       dmat[numParam * col + row] = out[numParam + dx];
981       if (row == col) {
982         col=0; ++row;
983       } else {
984         dmat[numParam * row + col] = out[numParam + dx];
985         ++col;
986       }
987     }
988   }
989
990   double *tmat = dmat + numParam * numParam;
991   for (int dx=0; dx < numParam * numParam; dx++) tmat[dx] = 0;
992   for (int dx=0; dx < nfact; dx++) {
993     tmat[dx * numParam + dx] = 1;
994   }
995   for (int rx=0; rx < nzeta; rx++) {
996     for (int cx=0; cx < nzeta; cx++) {
997       tmat[(cx + nfact)*numParam + nfact + rx] = Ta[rx * nzeta + cx];
998       tmat[(cx + nfact + nzeta)*numParam + nfact + nzeta + rx] = Tc[rx * nzeta + cx];
999     }
1000   }
1001
1002   double *dest = dmat + 2 * numParam * numParam;
1003
1004   // It is probably possible to do this more efficiently than dgemm
1005   // since we know that we only care about the lower triangle.
1006   // I'm not sure whether this is worth optimizing. TODO
1007
1008   char normal = 'n';
1009   char transpose = 't';
1010   double one = 1;
1011   double zero = 0;
1012   F77_CALL(dgemm)(&normal, &normal, &numParam, &numParam, &numParam,
1013                   &one, tmat, &numParam, dmat, &numParam, &zero, dest, &numParam);
1014   F77_CALL(dgemm)(&normal, &transpose, &numParam, &numParam, &numParam,
1015                   &one, dest, &numParam, tmat, &numParam, &zero, dmat, &numParam);
1016
1017   {
1018     int row=0;
1019     int col=0;
1020     for (int dx=0; dx <= hsize; dx++) {
1021       out[numParam + dx] = dmat[numParam * col + row];
1022       if (row == col) {
1023         col=0; ++row;
1024       } else {
1025         ++col;
1026       }
1027     }
1028   }
1029
1030   Free(dmat);
1031 }
1032
1033 static void
1034 irt_rpf_mdim_nrm_dTheta(const double *spec, const double *param,
1035                         const double *where, const double *dir,
1036                         double *grad, double *hess)
1037 {
1038   int numDims = spec[RPF_ISpecDims];
1039   int outcomes = spec[RPF_ISpecOutcomes];
1040   const double *aa = param;
1041   double num[outcomes];
1042   double ak[outcomes];
1043   double discr = dotprod(param, where, numDims);
1044   _nominal_rawprob(spec, param, where, discr, ak, num);
1045
1046   double den = 0;
1047   for (int kx=0; kx < outcomes; kx++) {
1048     num[kx] = exp(num[kx]);
1049     den += num[kx];
1050   }
1051
1052   double P[outcomes];
1053   for (int kx=0; kx < outcomes; kx++) {
1054     P[kx] = num[kx]/den;
1055   }
1056
1057   for(int jx=0; jx < numDims; jx++) {
1058     double jak[outcomes];
1059     double jak2[outcomes];
1060     for (int ax=0; ax < outcomes; ax++) {
1061       jak[ax] = ak[ax] * aa[jx];
1062       jak2[ax] = jak[ax] * jak[ax];
1063     }
1064     double numjak = dotprod(num, jak, outcomes);
1065     double numjakden2 = numjak / den;
1066     numjakden2 *= numjakden2;
1067     double numjak2den = dotprod(num, jak2, outcomes) / den;
1068
1069     for(int ix=0; ix < outcomes; ix++) {
1070       grad[ix] += dir[jx] * (ak[ix] * aa[jx] * P[ix] - P[ix] * numjak / den);
1071       hess[ix] += dir[jx] * (ak[ix]*ak[ix] * aa[jx]*aa[jx] * P[ix] -
1072                              2 * ak[ix] * aa[jx] * P[ix] * numjak / den +
1073                              2 * P[ix] * numjakden2 - P[ix] * numjak2den);
1074     }
1075   }
1076 }
1077
1078 static void
1079 irt_rpf_mdim_nrm_rescale(const double *spec, double *restrict param, const int *paramMask,
1080                          const double *restrict mean, const double *restrict cov)
1081 {
1082   int numDims = spec[RPF_ISpecDims];
1083   int nzeta = spec[RPF_ISpecOutcomes] - 1;
1084   double *alpha = param + numDims;
1085   double *gamma = param + numDims + nzeta;
1086   const double *Ta  = spec + RPF_ISpecCount;
1087   const double *Tc  = spec + RPF_ISpecCount + nzeta * nzeta;
1088   const double *iTc = spec + RPF_ISpecCount + 3 * nzeta * nzeta;
1089
1090   double madj = dotprod(param, mean, numDims);
1091
1092   for (int d1=0; d1 < numDims; d1++) {
1093     if (paramMask[d1] < 0) continue;
1094     param[d1] = dotprod(param+d1, cov + d1 * numDims + d1, numDims-d1);
1095   }
1096
1097   double ak[nzeta];
1098   double ck[nzeta];
1099   memset(ak, 0, sizeof(double)*nzeta);
1100   memset(ck, 0, sizeof(double)*nzeta);
1101
1102   for (int kx=0; kx < nzeta; kx++) {
1103     for (int tx=0; tx < nzeta; tx++) {
1104       int Tcell = tx * nzeta + kx;
1105       ak[kx] += Ta[Tcell] * alpha[tx];
1106       ck[kx] += Tc[Tcell] * gamma[tx];
1107     }
1108   }
1109
1110   for (int kx=0; kx < nzeta; kx++) {
1111     ck[kx] += madj * ak[kx];
1112   }
1113
1114   for (int kx=0; kx < nzeta; kx++) {
1115     int px = numDims + nzeta + kx;
1116     if (paramMask[px] < 0) continue;
1117
1118     param[px] = 0;
1119
1120     for (int tx=0; tx < nzeta; tx++) {
1121       int Tcell = tx * nzeta + kx;
1122       param[px] += iTc[Tcell] * ck[tx];
1123     }
1124   }
1125 }
1126
1127 static void noop() {}
1128 static void notimplemented() { error("Not implemented"); }
1129
1130 const struct rpf librpf_model[] = {
1131   { "drm1-",
1132     irt_rpf_1dim_drm_numSpec,
1133     irt_rpf_1dim_drm_numParam,
1134     irt_rpf_mdim_drm_paramInfo,
1135     irt_rpf_1dim_drm_prob,
1136     irt_rpf_logprob_adapter,
1137     notimplemented,
1138     notimplemented,
1139     notimplemented,
1140     irt_rpf_1dim_drm_rescale,
1141   },
1142   { "drm",
1143     irt_rpf_mdim_drm_numSpec,
1144     irt_rpf_mdim_drm_numParam,
1145     irt_rpf_mdim_drm_paramInfo,
1146     irt_rpf_mdim_drm_prob,
1147     irt_rpf_logprob_adapter,
1148     irt_rpf_mdim_drm_deriv1,
1149     irt_rpf_mdim_drm_deriv2,
1150     irt_rpf_mdim_drm_dTheta,
1151     irt_rpf_mdim_drm_rescale,
1152   },
1153   { "grm",
1154     irt_rpf_mdim_grm_numSpec,
1155     irt_rpf_mdim_grm_numParam,
1156     irt_rpf_mdim_grm_paramInfo,
1157     irt_rpf_mdim_grm_prob,
1158     irt_rpf_logprob_adapter,
1159     irt_rpf_mdim_grm_deriv1,
1160     irt_rpf_mdim_grm_deriv2,
1161     irt_rpf_mdim_grm_dTheta,
1162     irt_rpf_mdim_grm_rescale,
1163   },
1164   { "nominal",
1165     irt_rpf_nominal_numSpec,
1166     irt_rpf_nominal_numParam,
1167     irt_rpf_nominal_paramInfo,
1168     irt_rpf_nominal_prob,
1169     irt_rpf_nominal_logprob,
1170     irt_rpf_nominal_deriv1,
1171     irt_rpf_nominal_deriv2,
1172     irt_rpf_mdim_nrm_dTheta,
1173     irt_rpf_mdim_nrm_rescale,
1174   }
1175 };
1176
1177 const int librpf_numModels = (sizeof(librpf_model) / sizeof(struct rpf));