> I'm trying to replicate a C++ code with R.Notes: (1) I'd recommend you make the code more modular. i.e. One function for initial data prep/modelling, one function for setting up and solving the QP, etc. This should be easier to debug. (However, you would probably have to do it to the C++ code first). (2) Your R code is not completely reproducible. i.e. AEJData (3) For the purposes of a reproducible example, your code can be simplified. i.e. Only one contributed R package should be attached. Regardless of (1) above, you should be able to identify at what point the C++ and R code becomes inconsistent. The simplest approach is to add print-based functions into both the C++ and R code, and print out state data, at each major step. Then all you need to do is compare the output for both.> Is there a better/another package for these types of problems?I'm not sure. After a quick search, this is the best I found: scam::scam scam::shape.constrained.smooth.terms
Thank you for giving me your time!
The problem is the quadratic optimization part. Something goes wrong along
the way. In C++ loops run from 0 and in R they run from 1, and I've tried
to take that into account. Still I'm having hard time figuring out the
mistake I make, cause I get a result from my R code. It's just not the same
that I get with the C++.
Here are the quadratic optimization parts for both codes.
C++
/* Set Up Quadratic Programing Problem */
Vector<double> hSmooth(J);
for(int j=0; j<J; j++) hSmooth(j) = -pow(kr.Phi(j),Rho);
Vector<double> Q(J);
for(int j=0; j<J; j++) Q(j) = exp(-Rho * (Beta * pow(Price(j),Eta + One) -
One) / (One + Eta));
SymmetricMatrix<double> H(J,Zero);
Vector<double> c(J,Zero);
Matrix<double> Aeq(0,J);
Vector<double> beq(0);
Matrix<double> Aneq(2*J-3,J,Zero);
Vector<double> bneq(2*J-3);
Vector<double> lb(J,-Inf);
Vector<double> ub(J,Inf);
for(int j=0; j<J; j++) H(j,j) = One;
for(int j=0; j<J; j++) c(j) = -hSmooth(j);
for(int j=1; j<J; j++)
{
Aneq(j-1,j-1) = -One;
Aneq(j-1,j) = One;
bneq[j-1] = Delta1;
}
for(int j=2; j<J; j++)
{
Aneq(J-1+j-2,j) = -One / (Q(j) - Q(j-1));
Aneq(J-1+j-2,j-1) = One / (Q(j) - Q(j-1)) + One / (Q(j-1) - Q(j-2));
Aneq(J-1+j-2,j-2) = -One / (Q(j-1) - Q(j-2));
bneq[J-1+j-2] = Delta2;
}
/* Solve Constrained Optimization Problem Using Quadratic Programming */
MinQuadProg qp(c,H,Aeq,beq,Aneq,bneq,lb,ub);
qp.PrintLevel = 0;
qp.Solve();
And R
J <- length(Price)
hs <- numeric(J)
for(j in 1:J){
hs[j] <-(-(gEst$KernelRegPartLin..Phi[j]^(-0.1)))
}
hs
Q <- rep(0,J)
for(j in 1:(length(Price))){
Q[j] <- exp((-0.1) * (Beta *Price[j]^(Eta + 1) - 1) / (1 + Eta))
}
Q
plot(Q)
Dmat <- matrix(0,nrow= J, ncol=J)
diag(Dmat) <- 1
dvec <- -hs
Aeq <- 0
beq <- 0
Amat <- matrix(0,J,2*J-3)
bvec <- rep(0,2*J-3)
for(j in 2:nrow(Amat)){
Amat[j-1,j-1] = -1
Amat[j,j-1] = 1
}
for(j in 3:nrow(Amat)){
Amat[j,J+j-3] = -1/(Q[j]-Q[j-1])
Amat[j-1,J+j-3] = 1/(Q[j]-Q[j-1])
Amat[j-2,J+j-3] = -1/(Q[j-1]-Q[j-2])
}
for(j in 2:nrow(bvec)) {
bvec[j-1] = Delta1
}
for(j in 3:nrow(bvec)) {
bvec[J-1+j-2] = Delta2
}
solution <- solve.QP(Dmat,dvec,Amat,bvec)
ke 23. syysk. 2020 klo 9.12 Abby Spurdle (spurdle.a at gmail.com) kirjoitti:
> > I'm trying to replicate a C++ code with R.
>
> Notes:
> (1) I'd recommend you make the code more modular.
> i.e. One function for initial data prep/modelling, one function for
> setting up and solving the QP, etc.
> This should be easier to debug.
> (However, you would probably have to do it to the C++ code first).
> (2) Your R code is not completely reproducible.
> i.e. AEJData
> (3) For the purposes of a reproducible example, your code can be
> simplified.
> i.e. Only one contributed R package should be attached.
>
> Regardless of (1) above, you should be able to identify at what point
> the C++ and R code becomes inconsistent.
> The simplest approach is to add print-based functions into both the
> C++ and R code, and print out state data, at each major step.
> Then all you need to do is compare the output for both.
>
> > Is there a better/another package for these types of problems?
>
> I'm not sure.
> After a quick search, this is the best I found:
>
> scam::scam
> scam::shape.constrained.smooth.terms
>
[[alternative HTML version deleted]]
Before going any further, I have to check, what is: MinQuadProg qp Also, if I'm following the C++ code correctly, H, is an identity matrix. This implies the input to the C++ solver, requires the QP in a different form to the R solver. In which case, the C++ inputs and the R inputs, should be different...? (A) It may be worthwhile comparing the solver output (for C++ and R) using a *much* simpler example. e.g. the examples from the quadprog package. (B) If they produce the same output (using a simple example), then that implies there's a difference in your inputs. So, you just need to work out which input values are different. Expanding on my previous post, just print them out. But check (A) above first. On Thu, Sep 24, 2020 at 11:15 PM Maija Sirkj?rvi <maija.sirkjarvi at gmail.com> wrote:> > Thank you for giving me your time! > > The problem is the quadratic optimization part. Something goes wrong along the way. In C++ loops run from 0 and in R they run from 1, and I've tried to take that into account. Still I'm having hard time figuring out the mistake I make, cause I get a result from my R code. It's just not the same that I get with the C++. > > Here are the quadratic optimization parts for both codes. > > C++ > > /* Set Up Quadratic Programing Problem */ > Vector<double> hSmooth(J); > for(int j=0; j<J; j++) hSmooth(j) = -pow(kr.Phi(j),Rho); > Vector<double> Q(J); > for(int j=0; j<J; j++) Q(j) = exp(-Rho * (Beta * pow(Price(j),Eta + One) - One) / (One + Eta)); > SymmetricMatrix<double> H(J,Zero); > Vector<double> c(J,Zero); > Matrix<double> Aeq(0,J); > Vector<double> beq(0); > Matrix<double> Aneq(2*J-3,J,Zero); > Vector<double> bneq(2*J-3); > Vector<double> lb(J,-Inf); > Vector<double> ub(J,Inf); > for(int j=0; j<J; j++) H(j,j) = One; > for(int j=0; j<J; j++) c(j) = -hSmooth(j); > > for(int j=1; j<J; j++) > { > Aneq(j-1,j-1) = -One; > Aneq(j-1,j) = One; > bneq[j-1] = Delta1; > } > for(int j=2; j<J; j++) > { > Aneq(J-1+j-2,j) = -One / (Q(j) - Q(j-1)); > Aneq(J-1+j-2,j-1) = One / (Q(j) - Q(j-1)) + One / (Q(j-1) - Q(j-2)); > Aneq(J-1+j-2,j-2) = -One / (Q(j-1) - Q(j-2)); > bneq[J-1+j-2] = Delta2; > } > > /* Solve Constrained Optimization Problem Using Quadratic Programming */ > MinQuadProg qp(c,H,Aeq,beq,Aneq,bneq,lb,ub); > qp.PrintLevel = 0; > qp.Solve(); > > And R > > J <- length(Price) > hs <- numeric(J) > for(j in 1:J){ > hs[j] <-(-(gEst$KernelRegPartLin..Phi[j]^(-0.1))) > } > hs > > Q <- rep(0,J) > for(j in 1:(length(Price))){ > Q[j] <- exp((-0.1) * (Beta *Price[j]^(Eta + 1) - 1) / (1 + Eta)) > } > Q > plot(Q) > Dmat <- matrix(0,nrow= J, ncol=J) > diag(Dmat) <- 1 > dvec <- -hs > Aeq <- 0 > beq <- 0 > Amat <- matrix(0,J,2*J-3) > bvec <- rep(0,2*J-3) > > for(j in 2:nrow(Amat)){ > Amat[j-1,j-1] = -1 > Amat[j,j-1] = 1 > } > > for(j in 3:nrow(Amat)){ > Amat[j,J+j-3] = -1/(Q[j]-Q[j-1]) > Amat[j-1,J+j-3] = 1/(Q[j]-Q[j-1]) > Amat[j-2,J+j-3] = -1/(Q[j-1]-Q[j-2]) > } > > for(j in 2:nrow(bvec)) { > bvec[j-1] = Delta1 > } > for(j in 3:nrow(bvec)) { > bvec[J-1+j-2] = Delta2 > } > > solution <- solve.QP(Dmat,dvec,Amat,bvec) > > > > > > ke 23. syysk. 2020 klo 9.12 Abby Spurdle (spurdle.a at gmail.com) kirjoitti: >> >> > I'm trying to replicate a C++ code with R. >> >> Notes: >> (1) I'd recommend you make the code more modular. >> i.e. One function for initial data prep/modelling, one function for >> setting up and solving the QP, etc. >> This should be easier to debug. >> (However, you would probably have to do it to the C++ code first). >> (2) Your R code is not completely reproducible. >> i.e. AEJData >> (3) For the purposes of a reproducible example, your code can be simplified. >> i.e. Only one contributed R package should be attached. >> >> Regardless of (1) above, you should be able to identify at what point >> the C++ and R code becomes inconsistent. >> The simplest approach is to add print-based functions into both the >> C++ and R code, and print out state data, at each major step. >> Then all you need to do is compare the output for both. >> >> > Is there a better/another package for these types of problems? >> >> I'm not sure. >> After a quick search, this is the best I found: >> >> scam::scam >> scam::shape.constrained.smooth.terms