optimizing with non-linear constraints

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

optimizing with non-linear constraints

Rainer Krug-3
Hi

I have posted the following question on stackoverflow [1] but haven't
received a response yet. Has somebody here any idea, how to formulate
non-linear constraints using e.g. the packages alabama or nlopr? The
question is pasted below.

I am really stuck with this.

Thanks,

Rainer

,----
| I am successful in using constrOptim() for linear constraints (thanks to
| a comment in my earlier question), but I now have a problem with
| non-linear constraints.
|
| Consider the following example:
|
| I want to fit a function FUN(X1, X2, ... X12)
|
|  X1 = dep.a
|  X2 = dep.b
|  X3 = dep.c
|  ... for all z0. , na. and zjoint.
|
| and have the following constraint
| (link to formulas:
| http://stackoverflow.com/questions/32845169/how-to-optimize-with-non-linear-constraints-in-r
| )
| |
| With some re-formating and h > 0 and LAI > 0 I get for thew first one
|
| -X1 <= LAI^X2 / X3 < 1 - X1
|
| where LAI is simply a number which constant for each fit.
|
| I saw the packages alabama and nloptr but I don't manage to get my head
| around how I can specify these types of constraints?
|
| Also a web search did not yield any helpful information for me.
|
| Could somebody provide some info how I can convert these constraints
| into arguments for the fitting functions (e.g. hin, heq in the alabama
| package)?
`----


Footnotes:
[1]  http://stackoverflow.com/questions/32845169/how-to-optimize-with-non-linear-constraints-in-r

--
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :       +33 - (0)9 53 10 27 44
Cell:       +33 - (0)6 85 62 59 98
Fax :       +33 - (0)9 58 10 27 44

Fax (D):    +49 - (0)3 21 21 25 22 44

email:      [hidden email]

Skype:      RMkrug

PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

signature.asc (463 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: optimizing with non-linear constraints

Ravi Varadhan-2
Hi Rainer,
It is very simple to specify the constraints (linear or nonlinear) in "alabama" .  They are specified in a function called `hin', where the constraints are written such that they are positive.  Your two nonlinear constraints would be written as follows:

hin <- function(x, LAI) {
h <- rep(NA, 2)
h[1] <- LAI^x[2] / x[3] + x[1]
h[2] <- 1 - x[1] - LAI^x[2] / x[3]
h
}

Please take a look at the help page.  If it is still not clear, you can contact me offline.
Best,
Ravi

Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
Associate Professor,  Department of Oncology
Division of Biostatistics & Bionformatics
Sidney Kimmel Comprehensive Cancer Center
Johns Hopkins University
550 N. Broadway, Suite 1111-E
Baltimore, MD 21205
410-502-2619


        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optimizing with non-linear constraints

Rainer Krug-3
Ravi Varadhan <[hidden email]> writes:

> Hi Rainer,
> It is very simple to specify the constraints (linear or nonlinear) in
> "alabama" .  They are specified in a function called `hin', where the
> constraints are written such that they are positive.

OK - I somehow missed the part that, when the values x are valid, i.e. in
the range as defined by the conditions, the result of hin(x) that they
are all positive.

> Your two nonlinear constraints would be written as follows:
>
> hin <- function(x, LAI) {
> h <- rep(NA, 2)
> h[1] <- LAI^x[2] / x[3] + x[1]
> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
> h
> }

Makes perfect sense.

>
> Please take a look at the help page.  If it is still not clear, you can contact me offline.

Yup - I did. But I somehow missed the fact stated above.

I am using constrOptim() and constrOptim.nl() for a paper
and am compiling a separate document which explains how to get the
constraints for the two functions step by step - I will make it
available as a blog post and a pdf.

I might have further questions concerning the different fitting
functions and which ones are the most appropriate in my case.

Thanks a lot,

Rainer


> Best,
> Ravi
>
> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
> Associate Professor,  Department of Oncology
> Division of Biostatistics & Bionformatics
> Sidney Kimmel Comprehensive Cancer Center
> Johns Hopkins University
> 550 N. Broadway, Suite 1111-E
> Baltimore, MD 21205
> 410-502-2619
>
>
> [[alternative HTML version deleted]]
>
--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

signature.asc (463 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: optimizing with non-linear constraints

Ravi Varadhan-2
I would recommend that you use auglag() rather than constrOptim.nl() in the package "alabama."  It is a better algorithm, and it does not require feasible starting values.
Best,
Ravi  

-----Original Message-----
From: Rainer M Krug [mailto:[hidden email]]
Sent: Thursday, October 01, 2015 3:37 AM
To: Ravi Varadhan <[hidden email]>
Cc: '[hidden email]' <[hidden email]>
Subject: Re: optimizing with non-linear constraints

Ravi Varadhan <[hidden email]> writes:

> Hi Rainer,
> It is very simple to specify the constraints (linear or nonlinear) in
> "alabama" .  They are specified in a function called `hin', where the
> constraints are written such that they are positive.

OK - I somehow missed the part that, when the values x are valid, i.e. in the range as defined by the conditions, the result of hin(x) that they are all positive.

> Your two nonlinear constraints would be written as follows:
>
> hin <- function(x, LAI) {
> h <- rep(NA, 2)
> h[1] <- LAI^x[2] / x[3] + x[1]
> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
> h
> }

Makes perfect sense.

>
> Please take a look at the help page.  If it is still not clear, you can contact me offline.

Yup - I did. But I somehow missed the fact stated above.

I am using constrOptim() and constrOptim.nl() for a paper and am compiling a separate document which explains how to get the constraints for the two functions step by step - I will make it available as a blog post and a pdf.

I might have further questions concerning the different fitting functions and which ones are the most appropriate in my case.

Thanks a lot,

Rainer


> Best,
> Ravi
>
> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
> Associate Professor,  Department of Oncology Division of Biostatistics
> & Bionformatics Sidney Kimmel Comprehensive Cancer Center Johns
> Hopkins University
> 550 N. Broadway, Suite 1111-E
> Baltimore, MD 21205
> 410-502-2619
>
>
> [[alternative HTML version deleted]]
>

--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: optimizing with non-linear constraints

Rainer M Krug-6


Envoyé de mon iPhone

> Le 1 oct. 2015 à 15:17, Ravi Varadhan <[hidden email]> a écrit :
>
> I would recommend that you use auglag() rather than constrOptim.nl() in the package "alabama."  It is a better algorithm, and it does not require feasible starting values.

Thanks - that was one question I wanted to ask later.

I will do so,

Rainer

> Best,
> Ravi  
>
> -----Original Message-----
> From: Rainer M Krug [mailto:[hidden email]]
> Sent: Thursday, October 01, 2015 3:37 AM
> To: Ravi Varadhan <[hidden email]>
> Cc: '[hidden email]' <[hidden email]>
> Subject: Re: optimizing with non-linear constraints
>
> Ravi Varadhan <[hidden email]> writes:
>
>> Hi Rainer,
>> It is very simple to specify the constraints (linear or nonlinear) in
>> "alabama" .  They are specified in a function called `hin', where the
>> constraints are written such that they are positive.
>
> OK - I somehow missed the part that, when the values x are valid, i.e. in the range as defined by the conditions, the result of hin(x) that they are all positive.
>
>> Your two nonlinear constraints would be written as follows:
>>
>> hin <- function(x, LAI) {
>> h <- rep(NA, 2)
>> h[1] <- LAI^x[2] / x[3] + x[1]
>> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
>> h
>> }
>
> Makes perfect sense.
>
>>
>> Please take a look at the help page.  If it is still not clear, you can contact me offline.
>
> Yup - I did. But I somehow missed the fact stated above.
>
> I am using constrOptim() and constrOptim.nl() for a paper and am compiling a separate document which explains how to get the constraints for the two functions step by step - I will make it available as a blog post and a pdf.
>
> I might have further questions concerning the different fitting functions and which ones are the most appropriate in my case.
>
> Thanks a lot,
>
> Rainer
>
>
>> Best,
>> Ravi
>>
>> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
>> Associate Professor,  Department of Oncology Division of Biostatistics
>> & Bionformatics Sidney Kimmel Comprehensive Cancer Center Johns
>> Hopkins University
>> 550 N. Broadway, Suite 1111-E
>> Baltimore, MD 21205
>> 410-502-2619
>>
>>
>>    [[alternative HTML version deleted]]
>>
>
> --
> Rainer M. Krug
> email: Rainer<at>krugs<dot>de
> PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Bug in auglag?

Rainer Krug-3
In reply to this post by Ravi Varadhan-2
Hi Ravi,

I would like come back to your offer. I have a problem which possibly is
caused by a bug or by something I don't understand:

My function to be minimised is executed even when an element in hin() is
negative.

My hin looks as follow:

--8<---------------cut here---------------start------------->8---
hinMahat <- function(x, hauteur, na, zjoint, y, LAI, ...) {
    if (x[1] < 0) {
        cat(names(list(...)), "\n")
        cat(..., "\n")
        cat(x, "|", hauteur, LAI, y, "\n")
    }

    h <- rep(NA, 8)
    if (!missing(na)) {
        x <- c(na, x )
    }
    if (!missing(y)) {
        x <- c(x, y)
    }
    if (!missing(zjoint)) {
        x <- c(x[1], zjoint, x[2])
    }
   
    ##
    dep <- hauteur * (0.05 + LAI^0.02 / 2) + (x[3] - 1)/20
    h[1] <- dep
    h[2] <- hauteur - dep
    ## if (h[2]==0) {
    ##     h[2] <- -1
    ## }
    ##
    z0 <- hauteur * (0.23 + LAI^0.25 / 10) + (x[3] - 1)/67
    h[3] <- z0
    ## if (h[3]==0) {
    ##     h[3] <- -1
    ## }
    h[4] <- hauteur - z0
    ##
    h[5] <- x[1]
    ##
    h[6] <- x[2]
    h[7] <- hauteur - x[2]
    ##
    h[8] <- hauteur - dep - z0
    if (any(h<=0)) {
        cat(h, "\n")
        cat("\n")
    }
    return(h)
}
--8<---------------cut here---------------end--------------->8---

the x contains up to three elements: c(na=, zjoint=, y=) and I fit these
three, unless one or two are specified explicitely.

The values going into hin are:

,----
| ... (z  u ua za z0sol )
| 3 11 17 23 29 37 0.315 0.422 0.458 0.556 1.567 1.747 1.747 37 0.001
|
| x(na, zjoint): -8.875735 24.51316
| hauteur: 28
| na:      8.1
| y:       3
|
| the resulting hin() is:
| 16.09815 11.90185 11.19352 16.80648 -8.875735 24.51316 3.486843 0.708335
`----


Which is negative in element 5 as x[2]=na is negative.

So I would expect that the function fn is not evaluated. But it is, and
raises an error:

,----
| Error in wpLELMahat(z = z, ua = ua, na = ifelse(missing(na), par[1], na),  :
|   na has to be larger or equal than zero!
`----

Is this a misunderstanding on my part, or is it an error in the function
auglag?


Below is the function which is doing the minimisation.

If I replace auglag() with constrOptim.nl(), the optimisation is working
as expected.

So I think this is a bug in auglag?

Let me know if you need further information.

Cheers,

Rainer

--8<---------------cut here---------------start------------->8---
fitAuglag.wpLEL.mahat.single <- function(
                                         z,
                                         u,
                                         LAI,
                                         initial = c(na=9, zjoint=0.2*2, y=3),
                                         na, zjoint, y,
                                         h      = 28,
                                         za     = 37,
                                         z0sol  = 0.001,
                                         hin,
                                         ...
                                         ) {
    if (missing(hin)) {
        hin <- hinMahat
    }

    wpLELMin <- function(par, na, zjoint, y, z, u, ua, hauteur, za, z0sol, LAI) {
        result <- NA
        try({
            p <- wpLELMahat(
                z      = z,
                ua     = ua,
                na     = ifelse(missing(na), par[1], na),
                zjoint = ifelse(missing(zjoint), par[2], zjoint),
                h      = hauteur,
                za     = za,
                z0sol  = z0sol,
                LAI    = LAI,
                y      = ifelse(missing(y), par[3], y)
            )
            result <- sum( ( (p$u - u)^2 ) / length(u) )
        },
        silent = FALSE
        )
        ## cat("From wpLELMin", par, "\n")
        return( result )
    }

    ua <- u[length(u)]
    result <- list()
    result$method <- "fitAuglag.wpLEL.mahat.single"
    result$initial <-  initial
    result$dot <- list(...)
    result$z <- z
    result$u <- u

    result$fit <- auglag(
        par = initial,
        fn    = wpLELMin,
        hin   = hin,
        na     = na,
        zjoint = zjoint,
        y      = y,
        ##
        z     = z,
        u     = u,
        ua    = ua,
        hauteur = h,
        za    = za,
        z0sol = z0sol,
        LAI   = LAI,
        ...
    )
    result$wp <- wpLELMahat(
        z      = z,
        ua     = ua,
        na     = ifelse ( missing(na), result$fit$par["na"], na),
        zjoint = ifelse ( missing(zjoint), result$fit$par["zjoint"], zjoint),
        h      = h,
        za     = za,
        z0sol  = z0sol,
        LAI    = LAI,
        y      = ifelse ( missing(y), result$fit$par["y"], y)
    )
   
    class(result) <- c(class(result), "wpLELFit")
    return(result)
}
#+end_src--8<---------------cut here---------------end--------------->8---



Ravi Varadhan <[hidden email]> writes:

> I would recommend that you use auglag() rather than constrOptim.nl()
> in the package "alabama."  It is a better algorithm, and it does not
> require feasible starting values.
> Best,
> Ravi  
>
> -----Original Message-----
> From: Rainer M Krug [mailto:[hidden email]]
> Sent: Thursday, October 01, 2015 3:37 AM
> To: Ravi Varadhan <[hidden email]>
> Cc: '[hidden email]' <[hidden email]>
> Subject: Re: optimizing with non-linear constraints
>
> Ravi Varadhan <[hidden email]> writes:
>
>> Hi Rainer,
>> It is very simple to specify the constraints (linear or nonlinear) in
>> "alabama" .  They are specified in a function called `hin', where the
>> constraints are written such that they are positive.
>
> OK - I somehow missed the part that, when the values x are valid,
>> i.e. in the range as defined by the conditions, the result of hin(x)
>> that they are all positive.
>
>> Your two nonlinear constraints would be written as follows:
>>
>> hin <- function(x, LAI) {
>> h <- rep(NA, 2)
>> h[1] <- LAI^x[2] / x[3] + x[1]
>> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
>> h
>> }
>
> Makes perfect sense.
>
>>
>> Please take a look at the help page.  If it is still not clear, you can contact me offline.
>
> Yup - I did. But I somehow missed the fact stated above.
>
> I am using constrOptim() and constrOptim.nl() for a paper and am
>> compiling a separate document which explains how to get the
>> constraints for the two functions step by step - I will make it
>> available as a blog post and a pdf.
>
> I might have further questions concerning the different fitting
>> functions and which ones are the most appropriate in my case.
>
> Thanks a lot,
>
> Rainer
>
>
>> Best,
>> Ravi
>>
>> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
>> Associate Professor,  Department of Oncology Division of Biostatistics
>> & Bionformatics Sidney Kimmel Comprehensive Cancer Center Johns
>> Hopkins University
>> 550 N. Broadway, Suite 1111-E
>> Baltimore, MD 21205
>> 410-502-2619
>>
>>
>> [[alternative HTML version deleted]]
>>
>
> --
> Rainer M. Krug
> email: Rainer<at>krugs<dot>de
> PGP: 0x0F52F982
>
--
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :       +33 - (0)9 53 10 27 44
Cell:       +33 - (0)6 85 62 59 98
Fax :       +33 - (0)9 58 10 27 44

Fax (D):    +49 - (0)3 21 21 25 22 44

email:      [hidden email]

Skype:      RMkrug

PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

signature.asc (463 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Bug in auglag?

Rainer Krug-3
Please ignore - list members - accidentally CCd.

Rainer


Rainer M Krug <[hidden email]> writes:

> Hi Ravi,
>
> I would like come back to your offer. I have a problem which possibly is
> caused by a bug or by something I don't understand:
>
> My function to be minimised is executed even when an element in hin() is
> negative.
>
> My hin looks as follow:
>
> hinMahat <- function(x, hauteur, na, zjoint, y, LAI, ...) {
>     if (x[1] < 0) {
>         cat(names(list(...)), "\n")
>         cat(..., "\n")
>         cat(x, "|", hauteur, LAI, y, "\n")
>     }
>
>     h <- rep(NA, 8)
>     if (!missing(na)) {
>         x <- c(na, x )
>     }
>     if (!missing(y)) {
>         x <- c(x, y)
>     }
>     if (!missing(zjoint)) {
>         x <- c(x[1], zjoint, x[2])
>     }
>    
>     ##
>     dep <- hauteur * (0.05 + LAI^0.02 / 2) + (x[3] - 1)/20
>     h[1] <- dep
>     h[2] <- hauteur - dep
>     ## if (h[2]==0) {
>     ##     h[2] <- -1
>     ## }
>     ##
>     z0 <- hauteur * (0.23 + LAI^0.25 / 10) + (x[3] - 1)/67
>     h[3] <- z0
>     ## if (h[3]==0) {
>     ##     h[3] <- -1
>     ## }
>     h[4] <- hauteur - z0
>     ##
>     h[5] <- x[1]
>     ##
>     h[6] <- x[2]
>     h[7] <- hauteur - x[2]
>     ##
>     h[8] <- hauteur - dep - z0
>     if (any(h<=0)) {
>         cat(h, "\n")
>         cat("\n")
>     }
>     return(h)
> }
>
> the x contains up to three elements: c(na=, zjoint=, y=) and I fit these
> three, unless one or two are specified explicitely.
>
> The values going into hin are:
>
> ,----
> | ... (z  u ua za z0sol )
> | 3 11 17 23 29 37 0.315 0.422 0.458 0.556 1.567 1.747 1.747 37 0.001
> |
> | x(na, zjoint): -8.875735 24.51316
> | hauteur: 28
> | na:      8.1
> | y:       3
> |
> | the resulting hin() is:
> | 16.09815 11.90185 11.19352 16.80648 -8.875735 24.51316 3.486843 0.708335
> `----
>
>
> Which is negative in element 5 as x[2]=na is negative.
>
> So I would expect that the function fn is not evaluated. But it is, and
> raises an error:
>
> ,----
> | Error in wpLELMahat(z = z, ua = ua, na = ifelse(missing(na), par[1], na),  :
> |   na has to be larger or equal than zero!
> `----
>
> Is this a misunderstanding on my part, or is it an error in the function
> auglag?
>
>
> Below is the function which is doing the minimisation.
>
> If I replace auglag() with constrOptim.nl(), the optimisation is working
> as expected.
>
> So I think this is a bug in auglag?
>
> Let me know if you need further information.
>
> Cheers,
>
> Rainer
>
> --8<---------------cut here---------------start------------->8---
> fitAuglag.wpLEL.mahat.single <- function(
>                                          z,
>                                          u,
>                                          LAI,
>                                          initial = c(na=9, zjoint=0.2*2, y=3),
>                                          na, zjoint, y,
>                                          h      = 28,
>                                          za     = 37,
>                                          z0sol  = 0.001,
>                                          hin,
>                                          ...
>                                          ) {
>     if (missing(hin)) {
>         hin <- hinMahat
>     }
>
>     wpLELMin <- function(par, na, zjoint, y, z, u, ua, hauteur, za, z0sol, LAI) {
>         result <- NA
>         try({
>             p <- wpLELMahat(
>                 z      = z,
>                 ua     = ua,
>                 na     = ifelse(missing(na), par[1], na),
>                 zjoint = ifelse(missing(zjoint), par[2], zjoint),
>                 h      = hauteur,
>                 za     = za,
>                 z0sol  = z0sol,
>                 LAI    = LAI,
>                 y      = ifelse(missing(y), par[3], y)
>             )
>             result <- sum( ( (p$u - u)^2 ) / length(u) )
>         },
>         silent = FALSE
>         )
>         ## cat("From wpLELMin", par, "\n")
>         return( result )
>     }
>
>     ua <- u[length(u)]
>     result <- list()
>     result$method <- "fitAuglag.wpLEL.mahat.single"
>     result$initial <-  initial
>     result$dot <- list(...)
>     result$z <- z
>     result$u <- u
>
>     result$fit <- auglag(
>         par = initial,
>         fn    = wpLELMin,
>         hin   = hin,
>         na     = na,
>         zjoint = zjoint,
>         y      = y,
>         ##
>         z     = z,
>         u     = u,
>         ua    = ua,
>         hauteur = h,
>         za    = za,
>         z0sol = z0sol,
>         LAI   = LAI,
>         ...
>     )
>     result$wp <- wpLELMahat(
>         z      = z,
>         ua     = ua,
>         na     = ifelse ( missing(na), result$fit$par["na"], na),
>         zjoint = ifelse ( missing(zjoint), result$fit$par["zjoint"], zjoint),
>         h      = h,
>         za     = za,
>         z0sol  = z0sol,
>         LAI    = LAI,
>         y      = ifelse ( missing(y), result$fit$par["y"], y)
>     )
>    
>     class(result) <- c(class(result), "wpLELFit")
>     return(result)
> }
> #+end_src--8<---------------cut here---------------end--------------->8---
>
>
>
> Ravi Varadhan <[hidden email]> writes:
>
>> I would recommend that you use auglag() rather than constrOptim.nl()
>> in the package "alabama."  It is a better algorithm, and it does not
>> require feasible starting values.
>> Best,
>> Ravi  
>>
>> -----Original Message-----
>> From: Rainer M Krug [mailto:[hidden email]]
>> Sent: Thursday, October 01, 2015 3:37 AM
>> To: Ravi Varadhan <[hidden email]>
>> Cc: '[hidden email]' <[hidden email]>
>> Subject: Re: optimizing with non-linear constraints
>>
>> Ravi Varadhan <[hidden email]> writes:
>>
>>> Hi Rainer,
>>> It is very simple to specify the constraints (linear or nonlinear) in
>>> "alabama" .  They are specified in a function called `hin', where the
>>> constraints are written such that they are positive.
>>
>> OK - I somehow missed the part that, when the values x are valid,
>>> i.e. in the range as defined by the conditions, the result of hin(x)
>>> that they are all positive.
>>
>>> Your two nonlinear constraints would be written as follows:
>>>
>>> hin <- function(x, LAI) {
>>> h <- rep(NA, 2)
>>> h[1] <- LAI^x[2] / x[3] + x[1]
>>> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
>>> h
>>> }
>>
>> Makes perfect sense.
>>
>>>
>>> Please take a look at the help page.  If it is still not clear, you can contact me offline.
>>
>> Yup - I did. But I somehow missed the fact stated above.
>>
>> I am using constrOptim() and constrOptim.nl() for a paper and am
>>> compiling a separate document which explains how to get the
>>> constraints for the two functions step by step - I will make it
>>> available as a blog post and a pdf.
>>
>> I might have further questions concerning the different fitting
>>> functions and which ones are the most appropriate in my case.
>>
>> Thanks a lot,
>>
>> Rainer
>>
>>
>>> Best,
>>> Ravi
>>>
>>> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
>>> Associate Professor,  Department of Oncology Division of Biostatistics
>>> & Bionformatics Sidney Kimmel Comprehensive Cancer Center Johns
>>> Hopkins University
>>> 550 N. Broadway, Suite 1111-E
>>> Baltimore, MD 21205
>>> 410-502-2619
>>>
>>>
>>> [[alternative HTML version deleted]]
>>>
>>
>> --
>> Rainer M. Krug
>> email: Rainer<at>krugs<dot>de
>> PGP: 0x0F52F982
>>
--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

signature.asc (463 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Bug in auglag?

Ravi Varadhan-2
In reply to this post by Rainer Krug-3
Dear Rainer,
This is NOT a bug in auglag.  I already mentioned that auglag() can work with infeasible starting values, which also implies that the function must be evaluable at infeasible values.  A simple solution to your problem would be to fix up your objective function such that it evaluates to `Inf' or some large value, when the parameter values are not in the constrained domain.  constrOptim.nl() is a barrier method so it forces the initial value and the subsequent iterates to be feasible.
Best,
Ravi
________________________________________
From: Rainer M Krug <[hidden email]>
Sent: Tuesday, October 6, 2015 9:20 AM
To: Ravi Varadhan
Cc: '[hidden email]'
Subject: Bug in auglag?

Hi Ravi,

I would like come back to your offer. I have a problem which possibly is
caused by a bug or by something I don't understand:

My function to be minimised is executed even when an element in hin() is
negative.

My hin looks as follow:

--8<---------------cut here---------------start------------->8---
hinMahat <- function(x, hauteur, na, zjoint, y, LAI, ...) {
    if (x[1] < 0) {
        cat(names(list(...)), "\n")
        cat(..., "\n")
        cat(x, "|", hauteur, LAI, y, "\n")
    }

    h <- rep(NA, 8)
    if (!missing(na)) {
        x <- c(na, x )
    }
    if (!missing(y)) {
        x <- c(x, y)
    }
    if (!missing(zjoint)) {
        x <- c(x[1], zjoint, x[2])
    }

    ##
    dep <- hauteur * (0.05 + LAI^0.02 / 2) + (x[3] - 1)/20
    h[1] <- dep
    h[2] <- hauteur - dep
    ## if (h[2]==0) {
    ##     h[2] <- -1
    ## }
    ##
    z0 <- hauteur * (0.23 + LAI^0.25 / 10) + (x[3] - 1)/67
    h[3] <- z0
    ## if (h[3]==0) {
    ##     h[3] <- -1
    ## }
    h[4] <- hauteur - z0
    ##
    h[5] <- x[1]
    ##
    h[6] <- x[2]
    h[7] <- hauteur - x[2]
    ##
    h[8] <- hauteur - dep - z0
    if (any(h<=0)) {
        cat(h, "\n")
        cat("\n")
    }
    return(h)
}
--8<---------------cut here---------------end--------------->8---

the x contains up to three elements: c(na=, zjoint=, y=) and I fit these
three, unless one or two are specified explicitely.

The values going into hin are:

,----
| ... (z  u ua za z0sol )
| 3 11 17 23 29 37 0.315 0.422 0.458 0.556 1.567 1.747 1.747 37 0.001
|
| x(na, zjoint): -8.875735 24.51316
| hauteur: 28
| na:      8.1
| y:       3
|
| the resulting hin() is:
| 16.09815 11.90185 11.19352 16.80648 -8.875735 24.51316 3.486843 0.708335
`----


Which is negative in element 5 as x[2]=na is negative.

So I would expect that the function fn is not evaluated. But it is, and
raises an error:

,----
| Error in wpLELMahat(z = z, ua = ua, na = ifelse(missing(na), par[1], na),  :
|   na has to be larger or equal than zero!
`----

Is this a misunderstanding on my part, or is it an error in the function
auglag?


Below is the function which is doing the minimisation.

If I replace auglag() with constrOptim.nl(), the optimisation is working
as expected.

So I think this is a bug in auglag?

Let me know if you need further information.

Cheers,

Rainer

--8<---------------cut here---------------start------------->8---
fitAuglag.wpLEL.mahat.single <- function(
                                         z,
                                         u,
                                         LAI,
                                         initial = c(na=9, zjoint=0.2*2, y=3),
                                         na, zjoint, y,
                                         h      = 28,
                                         za     = 37,
                                         z0sol  = 0.001,
                                         hin,
                                         ...
                                         ) {
    if (missing(hin)) {
        hin <- hinMahat
    }

    wpLELMin <- function(par, na, zjoint, y, z, u, ua, hauteur, za, z0sol, LAI) {
        result <- NA
        try({
            p <- wpLELMahat(
                z      = z,
                ua     = ua,
                na     = ifelse(missing(na), par[1], na),
                zjoint = ifelse(missing(zjoint), par[2], zjoint),
                h      = hauteur,
                za     = za,
                z0sol  = z0sol,
                LAI    = LAI,
                y      = ifelse(missing(y), par[3], y)
            )
            result <- sum( ( (p$u - u)^2 ) / length(u) )
        },
        silent = FALSE
        )
        ## cat("From wpLELMin", par, "\n")
        return( result )
    }

    ua <- u[length(u)]
    result <- list()
    result$method <- "fitAuglag.wpLEL.mahat.single"
    result$initial <-  initial
    result$dot <- list(...)
    result$z <- z
    result$u <- u

    result$fit <- auglag(
        par = initial,
        fn    = wpLELMin,
        hin   = hin,
        na     = na,
        zjoint = zjoint,
        y      = y,
        ##
        z     = z,
        u     = u,
        ua    = ua,
        hauteur = h,
        za    = za,
        z0sol = z0sol,
        LAI   = LAI,
        ...
    )
    result$wp <- wpLELMahat(
        z      = z,
        ua     = ua,
        na     = ifelse ( missing(na), result$fit$par["na"], na),
        zjoint = ifelse ( missing(zjoint), result$fit$par["zjoint"], zjoint),
        h      = h,
        za     = za,
        z0sol  = z0sol,
        LAI    = LAI,
        y      = ifelse ( missing(y), result$fit$par["y"], y)
    )

    class(result) <- c(class(result), "wpLELFit")
    return(result)
}
#+end_src--8<---------------cut here---------------end--------------->8---



Ravi Varadhan <[hidden email]> writes:

> I would recommend that you use auglag() rather than constrOptim.nl()
> in the package "alabama."  It is a better algorithm, and it does not
> require feasible starting values.
> Best,
> Ravi
>
> -----Original Message-----
> From: Rainer M Krug [mailto:[hidden email]]
> Sent: Thursday, October 01, 2015 3:37 AM
> To: Ravi Varadhan <[hidden email]>
> Cc: '[hidden email]' <[hidden email]>
> Subject: Re: optimizing with non-linear constraints
>
> Ravi Varadhan <[hidden email]> writes:
>
>> Hi Rainer,
>> It is very simple to specify the constraints (linear or nonlinear) in
>> "alabama" .  They are specified in a function called `hin', where the
>> constraints are written such that they are positive.
>
> OK - I somehow missed the part that, when the values x are valid,
>> i.e. in the range as defined by the conditions, the result of hin(x)
>> that they are all positive.
>
>> Your two nonlinear constraints would be written as follows:
>>
>> hin <- function(x, LAI) {
>> h <- rep(NA, 2)
>> h[1] <- LAI^x[2] / x[3] + x[1]
>> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
>> h
>> }
>
> Makes perfect sense.
>
>>
>> Please take a look at the help page.  If it is still not clear, you can contact me offline.
>
> Yup - I did. But I somehow missed the fact stated above.
>
> I am using constrOptim() and constrOptim.nl() for a paper and am
>> compiling a separate document which explains how to get the
>> constraints for the two functions step by step - I will make it
>> available as a blog post and a pdf.
>
> I might have further questions concerning the different fitting
>> functions and which ones are the most appropriate in my case.
>
> Thanks a lot,
>
> Rainer
>
>
>> Best,
>> Ravi
>>
>> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
>> Associate Professor,  Department of Oncology Division of Biostatistics
>> & Bionformatics Sidney Kimmel Comprehensive Cancer Center Johns
>> Hopkins University
>> 550 N. Broadway, Suite 1111-E
>> Baltimore, MD 21205
>> 410-502-2619
>>
>>
>>      [[alternative HTML version deleted]]
>>
>
> --
> Rainer M. Krug
> email: Rainer<at>krugs<dot>de
> PGP: 0x0F52F982
>

--
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :       +33 - (0)9 53 10 27 44
Cell:       +33 - (0)6 85 62 59 98
Fax :       +33 - (0)9 58 10 27 44

Fax (D):    +49 - (0)3 21 21 25 22 44

email:      [hidden email]

Skype:      RMkrug

PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Bug in auglag?

Rainer Krug-3
Thanks a lot Ravi for the clarification and it makes sense - bug in my
understanding  acknowledged.

I'll change the objective function as suggested.

Thanks,

Rainer

Ravi Varadhan <[hidden email]> writes:

> Dear Rainer,
> This is NOT a bug in auglag.  I already mentioned that auglag() can
> work with infeasible starting values, which also implies that the
> function must be evaluable at infeasible values.  A simple solution to
> your problem would be to fix up your objective function such that it
> evaluates to `Inf' or some large value, when the parameter values are
> not in the constrained domain.  constrOptim.nl() is a barrier method
> so it forces the initial value and the subsequent iterates to be
> feasible.
> Best,
> Ravi
> ________________________________________
> From: Rainer M Krug <[hidden email]>
> Sent: Tuesday, October 6, 2015 9:20 AM
> To: Ravi Varadhan
> Cc: '[hidden email]'
> Subject: Bug in auglag?
>
> Hi Ravi,
>
> I would like come back to your offer. I have a problem which possibly is
> caused by a bug or by something I don't understand:
>
> My function to be minimised is executed even when an element in hin() is
> negative.
>
> My hin looks as follow:
>
> hinMahat <- function(x, hauteur, na, zjoint, y, LAI, ...) {
>     if (x[1] < 0) {
>         cat(names(list(...)), "\n")
>         cat(..., "\n")
>         cat(x, "|", hauteur, LAI, y, "\n")
>     }
>
>     h <- rep(NA, 8)
>     if (!missing(na)) {
>         x <- c(na, x )
>     }
>     if (!missing(y)) {
>         x <- c(x, y)
>     }
>     if (!missing(zjoint)) {
>         x <- c(x[1], zjoint, x[2])
>     }
>
>     ##
>     dep <- hauteur * (0.05 + LAI^0.02 / 2) + (x[3] - 1)/20
>     h[1] <- dep
>     h[2] <- hauteur - dep
>     ## if (h[2]==0) {
>     ##     h[2] <- -1
>     ## }
>     ##
>     z0 <- hauteur * (0.23 + LAI^0.25 / 10) + (x[3] - 1)/67
>     h[3] <- z0
>     ## if (h[3]==0) {
>     ##     h[3] <- -1
>     ## }
>     h[4] <- hauteur - z0
>     ##
>     h[5] <- x[1]
>     ##
>     h[6] <- x[2]
>     h[7] <- hauteur - x[2]
>     ##
>     h[8] <- hauteur - dep - z0
>     if (any(h<=0)) {
>         cat(h, "\n")
>         cat("\n")
>     }
>     return(h)
> }
>
> the x contains up to three elements: c(na=, zjoint=, y=) and I fit these
> three, unless one or two are specified explicitely.
>
> The values going into hin are:
>
> ,----
> | ... (z  u ua za z0sol )
> | 3 11 17 23 29 37 0.315 0.422 0.458 0.556 1.567 1.747 1.747 37 0.001
> |
> | x(na, zjoint): -8.875735 24.51316
> | hauteur: 28
> | na:      8.1
> | y:       3
> |
> | the resulting hin() is:
> | 16.09815 11.90185 11.19352 16.80648 -8.875735 24.51316 3.486843 0.708335
> `----
>
>
> Which is negative in element 5 as x[2]=na is negative.
>
> So I would expect that the function fn is not evaluated. But it is, and
> raises an error:
>
> ,----
> | Error in wpLELMahat(z = z, ua = ua, na = ifelse(missing(na), par[1], na),  :
> |   na has to be larger or equal than zero!
> `----
>
> Is this a misunderstanding on my part, or is it an error in the function
> auglag?
>
>
> Below is the function which is doing the minimisation.
>
> If I replace auglag() with constrOptim.nl(), the optimisation is working
> as expected.
>
> So I think this is a bug in auglag?
>
> Let me know if you need further information.
>
> Cheers,
>
> Rainer
>
> --8<---------------cut here---------------start------------->8---
> fitAuglag.wpLEL.mahat.single <- function(
>                                          z,
>                                          u,
>                                          LAI,
>                                          initial = c(na=9, zjoint=0.2*2, y=3),
>                                          na, zjoint, y,
>                                          h      = 28,
>                                          za     = 37,
>                                          z0sol  = 0.001,
>                                          hin,
>                                          ...
>                                          ) {
>     if (missing(hin)) {
>         hin <- hinMahat
>     }
>
>     wpLELMin <- function(par, na, zjoint, y, z, u, ua, hauteur, za, z0sol, LAI) {
>         result <- NA
>         try({
>             p <- wpLELMahat(
>                 z      = z,
>                 ua     = ua,
>                 na     = ifelse(missing(na), par[1], na),
>                 zjoint = ifelse(missing(zjoint), par[2], zjoint),
>                 h      = hauteur,
>                 za     = za,
>                 z0sol  = z0sol,
>                 LAI    = LAI,
>                 y      = ifelse(missing(y), par[3], y)
>             )
>             result <- sum( ( (p$u - u)^2 ) / length(u) )
>         },
>         silent = FALSE
>         )
>         ## cat("From wpLELMin", par, "\n")
>         return( result )
>     }
>
>     ua <- u[length(u)]
>     result <- list()
>     result$method <- "fitAuglag.wpLEL.mahat.single"
>     result$initial <-  initial
>     result$dot <- list(...)
>     result$z <- z
>     result$u <- u
>
>     result$fit <- auglag(
>         par = initial,
>         fn    = wpLELMin,
>         hin   = hin,
>         na     = na,
>         zjoint = zjoint,
>         y      = y,
>         ##
>         z     = z,
>         u     = u,
>         ua    = ua,
>         hauteur = h,
>         za    = za,
>         z0sol = z0sol,
>         LAI   = LAI,
>         ...
>     )
>     result$wp <- wpLELMahat(
>         z      = z,
>         ua     = ua,
>         na     = ifelse ( missing(na), result$fit$par["na"], na),
>         zjoint = ifelse ( missing(zjoint), result$fit$par["zjoint"], zjoint),
>         h      = h,
>         za     = za,
>         z0sol  = z0sol,
>         LAI    = LAI,
>         y      = ifelse ( missing(y), result$fit$par["y"], y)
>     )
>
>     class(result) <- c(class(result), "wpLELFit")
>     return(result)
> }
> #+end_src--8<---------------cut here---------------end--------------->8---
>
>
>
> Ravi Varadhan <[hidden email]> writes:
>
>> I would recommend that you use auglag() rather than constrOptim.nl()
>> in the package "alabama."  It is a better algorithm, and it does not
>> require feasible starting values.
>> Best,
>> Ravi
>>
>> -----Original Message-----
>> From: Rainer M Krug [mailto:[hidden email]]
>> Sent: Thursday, October 01, 2015 3:37 AM
>> To: Ravi Varadhan <[hidden email]>
>> Cc: '[hidden email]' <[hidden email]>
>> Subject: Re: optimizing with non-linear constraints
>>
>> Ravi Varadhan <[hidden email]> writes:
>>
>>> Hi Rainer,
>>> It is very simple to specify the constraints (linear or nonlinear) in
>>> "alabama" .  They are specified in a function called `hin', where the
>>> constraints are written such that they are positive.
>>
>> OK - I somehow missed the part that, when the values x are valid,
>>> i.e. in the range as defined by the conditions, the result of hin(x)
>>> that they are all positive.
>>
>>> Your two nonlinear constraints would be written as follows:
>>>
>>> hin <- function(x, LAI) {
>>> h <- rep(NA, 2)
>>> h[1] <- LAI^x[2] / x[3] + x[1]
>>> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
>>> h
>>> }
>>
>> Makes perfect sense.
>>
>>>
>>> Please take a look at the help page.  If it is still not clear, you can contact me offline.
>>
>> Yup - I did. But I somehow missed the fact stated above.
>>
>> I am using constrOptim() and constrOptim.nl() for a paper and am
>>> compiling a separate document which explains how to get the
>>> constraints for the two functions step by step - I will make it
>>> available as a blog post and a pdf.
>>
>> I might have further questions concerning the different fitting
>>> functions and which ones are the most appropriate in my case.
>>
>> Thanks a lot,
>>
>> Rainer
>>
>>
>>> Best,
>>> Ravi
>>>
>>> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
>>> Associate Professor,  Department of Oncology Division of Biostatistics
>>> & Bionformatics Sidney Kimmel Comprehensive Cancer Center Johns
>>> Hopkins University
>>> 550 N. Broadway, Suite 1111-E
>>> Baltimore, MD 21205
>>> 410-502-2619
>>>
>>>
>>>      [[alternative HTML version deleted]]
>>>
>>
>> --
>> Rainer M. Krug
>> email: Rainer<at>krugs<dot>de
>> PGP: 0x0F52F982
>>
>
> --
> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany)
>
> Centre of Excellence for Invasion Biology
> Stellenbosch University
> South Africa
>
> Tel :       +33 - (0)9 53 10 27 44
> Cell:       +33 - (0)6 85 62 59 98
> Fax :       +33 - (0)9 58 10 27 44
>
> Fax (D):    +49 - (0)3 21 21 25 22 44
>
> email:      [hidden email]
>
> Skype:      RMkrug
>
> PGP: 0x0F52F982
>
--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

signature.asc (463 bytes) Download Attachment