Fix Ox implementation

Remove .h from Ox, fix `lex` typo, and add samples for Ox.

samples/Ox/IJCEmet2009.oxh

@@ -0,0 +1,72 @@
+/** Replicate Imai, Jain and Ching Econometrica 2009 (incomplete).
+
+**/
+#include "IJCEmet2009.h"
+
+Kapital::Kapital(L,const N,const entrant,const exit,const KP){
+	StateVariable(L,N);
+	this.entrant = entrant;
+	this.exit = exit;
+	this.KP = KP;
+	actual = Kbar*vals/(N-1);
+	upper = log(actual~.Inf);
+	}
+
+Kapital::Transit(FeasA) {
+	decl ent =CV(entrant), stayout = FeasA[][exit.pos], tprob, sigu = CV(KP[SigU]);
+	if (!v && !ent) return { <0>, ones(stayout) };
+	tprob = ent ? probn( (upper-CV(KP[Kbe]))/sigu )
+	           : probn( (upper-(CV(KP[Kb0])+CV(KP[Kb2])*upper[v])) / sigu );
+	tprob = tprob[1:] - tprob[:N-1];
+	return { vals, tprob.*(1-stayout)+(1.0~zeros(1,N-1)).*stayout };
+	}
+
+FirmEntry::Run() {
+	Initialize();
+	GenerateSample();
+	BDP->BayesianDP();
+	}
+
+FirmEntry::Initialize() {
+	Rust::Initialize(Reachable,0);
+	sige = new StDeviations("sige",<0.3,0.3>,0);
+	entrant = new LaggedAction("entrant",d);
+	KP = new array[Kparams];
+		KP[Kbe] = new Positive("be",0.5);
+		KP[Kb0] = new Free("b0",0.0);
+		KP[Kb1] = new Determined("b1",0.0);
+		KP[Kb2] = new Positive("b2",0.4);
+		KP[SigU] = new Positive("sigu",0.4);
+	EndogenousStates(K = new Kapital("K",KN,entrant,d,KP),entrant);
+	SetDelta(new Probability("delta",0.85));
+	kcoef = new Positive("kcoef",0.1);
+	ecost = new Negative("ec",-0.4);
+	CreateSpaces();
+	}
+
+FirmEntry::GenerateSample() {
+	Volume = LOUD;
+	EM = new ValueIteration(0);
+//	EM -> Solve(0,0);
+	data = new DataSet(0,EM);
+	data->Simulate(DataN,DataT,0,FALSE);
+	data->Print("firmentry.xls");
+	BDP = new ImaiJainChing("FMH",data,EM,ecost,sige,kcoef,KP,delta);
+	}
+
+/** Capital stock can be positive only for incumbents.
+**/
+FirmEntry::Reachable()	{ return CV(entrant)*CV(K) ? 0 : new FirmEntry() ;	}
+
+/** The one period return.
+<DD>
+<pre>U = </pre>
+</DD>
+**/
+FirmEntry::Utility()  {
+	decl ent = CV(entrant),
+		 u =
+		     ent*CV(ecost)+(1-ent)*CV(kcoef)*AV(K)
+		   | 0.0;
+	return u;
+	}

samples/Ox/ParallelObjective.ox

@@ -0,0 +1,63 @@
+/** Client and Server classes for parallel optimization using CFMPI.**/
+#include "ParallelObjective.h"
+
+/** Set up MPI Client-Server support for objective optimization.
+@param obj `Objective' to parallelize
+@param DONOTUSECLIENT TRUE (default): client node does no object evaluation<br>FALSE after putting servers to work Client node does one evaluation.
+**/
+ParallelObjective(obj,DONOTUSECLIENT) {
+	if (isclass(obj.p2p)) {oxwarning("P2P object already exists for "+obj.L+". Nothing changed"); return;}
+	obj.p2p = new P2P(DONOTUSECLIENT,new ObjClient(obj),new ObjServer(obj));
+	}
+
+ObjClient::ObjClient(obj) { this.obj = obj;    }
+
+ObjClient::Execute() {    }
+
+ObjServer::ObjServer(obj) {
+    this.obj = obj;
+    basetag = P2P::STOP_TAG+1;
+    iml = obj.NvfuncTerms;
+    Nparams = obj.nstruct;
+    }
+
+/** Wait on the objective client.
+**/
+ObjServer::Loop(nxtmsgsz) {
+    Nparams = nxtmsgsz;   //free param length is no greater than Nparams
+    if (Volume>QUIET) println("ObjServer server ",ID," Nparams ",Nparams);
+    Server::Loop(Nparams);
+    Recv(ANY_TAG);                      //receive the ending parameter vector
+    obj->Encode(Buffer[:Nparams-1]);   //encode it.
+    }
+
+/** Do the objective evaluation.
+Receive structural parameter vector and `Objective::Encode`() it.
+Call `Objective::vfunc`().
+@return Nparams (max. length of next expected message);
+**/
+ObjServer::Execute() {
+	obj->Decode(Buffer[:obj.nfree-1]);
+	Buffer = obj.cur.V[] = obj->vfunc();
+    if (Volume>QUIET) println("Server Executive: ",ID," vfunc[0]= ",Buffer[0]);
+	return obj.nstruct;
+	}
+
+CstrServer::CstrServer(obj) { ObjServer(obj);	}
+
+SepServer::SepServer(obj) { ObjServer(obj);	}
+
+CstrServer::Execute() {
+	obj->Encode(Buffer);
+	obj->Lagrangian(0);
+	return rows(Buffer = obj.cur->Vec());
+	}
+
+/** Separable objective evaluations.
+**/
+SepServer::Execute() {
+	obj.Kvar.v = imod(Tag-basetag,obj.K);
+	obj->Encode(Buffer,TRUE);
+	Buffer = obj.Kvar->PDF() * obj->vfunc();
+	return obj.NvfuncTerms;
+	}

samples/Ox/particle.oxo

@@ -0,0 +1,38 @@
+nldge::ParticleLogLikeli()
+{	decl it, ip,
+		 mss, mbas, ms, my, mx, vw, vwi, dws,
+		 mhi, mhdet, loglikeli, mData,
+		 vxm, vxs, mxm=<>, mxsu=<>, mxsl=<>,
+		 time, timeall, timeran=0, timelik=0, timefun=0, timeint=0, timeres=0;
+
+	mData = GetData(m_asY);
+	mhdet = sqrt((2*M_PI)^m_cY * determinant(m_mMSbE.^2));		// covariance determinant
+	mhi   = invert(m_mMSbE.^2);					// invert covariance of measurement shocks
+
+	ms 	  = m_vSss + zeros(m_cPar, m_cS);			// start particles
+	mx 	  = m_vXss + zeros(m_cPar, m_cX);			// steady state of state and policy
+
+	loglikeli = 0;							// init likelihood
+																								//timeall=timer();
+	for(it = 0; it < sizer(mData); it++)
+	{
+		mss = rann(m_cPar, m_cSS) * m_mSSbE;			// state noise
+		fg(&ms, ms, mx, mss);					// transition prior as proposal
+		mx = m_oApprox.FastInterpolate(ms); 			// interpolate
+		fy(&my, ms, mx, zeros(m_cPar, m_cMS));			// evaluate importance weights
+		my -= mData[it][];					// observation error
+
+		vw = exp(-0.5 * outer(my,mhi,'d')' )/mhdet;		// vw = exp(-0.5 * sumr(my*mhi .*my ) )/mhdet;
+
+		vw = vw .== .NaN .? 0 .: vw;				// no policy can happen for extrem particles
+		dws = sumc(vw);
+		if(dws==0) return -.Inf;				// or extremely wrong parameters
+		loglikeli += log(dws/m_cPar)	;			// loglikelihood contribution
+																										//timelik += (timer()-time)/100;
+																										//time=timer();
+		vwi = resample(vw/dws)-1;				// selection step in c++
+		ms = ms[vwi][];						// on normalized weights
+		mx = mx[vwi][];
+																	}
+	return loglikeli;
+}