[Developers] Accessing final gradient information
Matthew Supernaw
matthew.supernaw at noaa.gov
Tue Oct 2 13:02:07 PDT 2012
All,
Is there anyway to access final gradient information from the function_minimizer class. Some people I support have interest in doing this from the report section.
Thanks.
Matthew Supernaw
Scientific Programmer
National Oceanic and Atmospheric Administration
National Marine Fisheries Service
Sustainable Fisheries Division
St. Petersburg, FL, 33701
Office 727-551-5606
Fax 727-824-5300
On Sep 17, 2012, at 3:00 PM, developers-request at admb-project.org wrote:
> Send Developers mailing list submissions to
> developers at admb-project.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.admb-project.org/mailman/listinfo/developers
> or, via email, send a message with subject or body 'help' to
> developers-request at admb-project.org
>
> You can reach the person managing the list at
> developers-owner at admb-project.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of Developers digest..."
>
>
> Today's Topics:
>
> 1. Fwd: [New post] Julia, Matlab, and C (John Sibert)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Mon, 17 Sep 2012 07:18:38 -1000
> From: John Sibert <sibert at hawaii.edu>
> To: ADMB Developers <developers at admb-project.org>
> Subject: [Developers] Fwd: [New post] Julia, Matlab, and C
> Message-ID: <50575B6E.8090606 at hawaii.edu>
> Content-Type: text/plain; charset=UTF-8; format=flowed
>
>
>
>
> -------- Original Message --------
> Subject: [New post] Julia, Matlab, and C
> Date: Mon, 17 Sep 2012 02:37:37 +0000
> From: Justin Domkes Weblog <comment-reply at wordpress.com>
> Reply-To: Justin Domkes Weblog
> <comment+_dvt62oavm4btqopdkeujz at comment.wordpress.com>
> To: sibert at hawaii.edu
>
>
>
> WordPress.com
> justindomke posted: "Julia is a new language in the same arena as Matlab
> or R. I've had failed attempts to quit the Matlab addiction in the past,
> making me generally quite conservative about new platforms. However,
> I've recently been particularly annoyed by Matlab's slow spee"
> Respond to this post by replying above this line
>
>
> New post on *Justin Domke's Weblog*
>
>
>
> <http://justindomke.wordpress.com/author/justindomke/>
>
>
> Julia, Matlab, and C
> <http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/>
>
> by justindomke <http://justindomke.wordpress.com/author/justindomke/>
>
> Julia <http://julialang.org/> is a new language in the same arena as
> Matlab or R. I've had failed attempts
> <http://justindomke.wordpress.com/2008/11/29/mandelbrot-in-scala/> to
> quit the Matlab addiction in the past, making me generally quite
> conservative about new platforms. However, I've recently been
> particularly annoyed by Matlab's slow speed, evil license manager
> errors, restrictions on parallel processes, C++ .mex file pain, etc.,
> and so I decided to check it out. It seems inevitable that Matlab will
> eventually displaced by /something/. The question is: is that something
> Julia?
>
> "We want a language that?s open source, with a liberal license. We
> want the speed of C with the dynamism of Ruby. We want a language
> that?s homoiconic, with true macros like Lisp, but with obvious,
> familiar mathematical notation like Matlab. We want something as
> usable for general programming as Python, as easy for statistics as
> R, as natural for string processing as Perl, as powerful for linear
> algebra as Matlab, as good at gluing programs together as the shell.
> Something that is dirt simple to learn, yet keeps the most serious
> hackers happy. We want it interactive and we want it compiled."
>
> Essentially, the goal seems to be a faster, freer Matlab that treats
> users like adults (macros!) and /doesn't require writing any .mex files
> in C++ or Fortan/. Sounds too good to be true? I decided to try it out.
> My comparisons are to Matlab (R2012a) and C (gcc 4.2 with -O2 and
> -fno-builtin to prevent compile-time computations) all on a recent
> MacBook Pro (2.3 GHz Intel Core i7 w/ 8GB RAM).
>
> Installation was trivial: I just grabbed a pre-compiled binary
> <https://github.com/JuliaLang/julia/downloads> and started it up.
>
> I deliberately used naive algorithms, since I am just testing raw speed.
> It should be a fair comparison, as long as the algorithm is constant.
> Please let me know about any bugs, though.
>
> Click here to skip to the results <#results>.
>
>
> First benchmark: Fibonnaci
>
> % Matlab
> function f=fib(n)
> if n <= 2
> f=1.0;
> else
> f=fib(n-1)+fib(n-2);
> end
> end
>
> // C
> double fib(int n){
> if(n<=2)
> return(1.0);
> else
> return(fib(n-2)+fib(n-1));
> }
>
> % julia
> function fib(n)
> if n <= 2
> 1.0
> else
> fib(n-1)+fib(n-2);
> end
> end
>
> Clarity is basically a tie. Running them for n=30 we get:
>
> time in matlab (fib): 14.344231
> time in c (fib): 0.005795
> time in julia (fib): 0.281221
>
>
> Second benchmark: Matrix Multiplication
>
> This is a test of the naive O(N^3) matrix multiplication algorithm.
>
> % matlab
> function C=mmult(A,B,C)
> [M,N] = size(A);
> for i=1:M
> for j=1:M
> for k=1:M
> C(i,j) = C(i,j) + A(i,k)*B(k,j);
> end
> end
> end
> end
>
> // C
> #define M 500
> void mmult(double A[M][M],double B[M][M], double C[M][M]){
> //double C[M][M];
> int i,j,k;
> for(i=0; i<M; i++)
> for(j=0; j<M; j++){
> C[i][j] = 0;
> for(k=0; k<M; k++)
> C[i][j] += A[i][k]*B[k][j];
> }
> }
>
> # julia
> function mmult(A,B)
> (M,N) = size(A);
> C = zeros(M,M);
> for i=1:M
> for j=1:M
> for k=1:M
> C[i,j] += A[i,k]*B[k,j];
> end
> end
> end
> C;
> end
>
> Here, I think that Matlab and Julia and a bit clearer, and Julia wins
> though the wonders of having "+=". The timing results on 500x500
> matrices are:
>
> time in matlab (matmult): 1.229571
> time in c (matmult): 0.168296
> time in julia (matmult): 0.505222
>
>
> Third Benchmark: numerical quadrature
>
> Here, we attempt to calculate the integral \int_{x=5}^{15} sin(x) dx by
> numerical quadrature, using a simple midpoint rule with computations at
> 10^7 points.
>
> % matlab
> function val=numquad(lb,ub,npoints)
> val = 0.0;
> for x=lb:(ub-lb)/npoints:ub
> val = val + sin(x)/npoints;
> end
> end
>
> // C
> double numquad(double lb,double ub,int npoints){
> double val = 0.0;
> int i;
> for(i=0; i<=npoints; i++){
> double x = lb + (ub-lb)*i/npoints;
> val += sin(x)/npoints;
> }
> return(val);
> }
>
> # julia
> function numquad(lb,ub,npoints)
> val = 0.0
> for x=lb:(ub-lb)/npoints:ub
> val += sin(x)/npoints
> end
> val
> end
>
> The timings are:
>
> time in matlab (numquad): 0.446151
> time in c (numquad): 0.201071
> time in julia (numquad): 0.273863
>
>
> Fourth Benchmark: Belief Propagation
>
> Finally, I decided to try a little algorithm similar to what I actually
> tend to implement for my research. Roughly speaking, Belief Propagation
> <http://en.wikipedia.org/wiki/Belief_propagation> is a repeated sequence
> of matrix multiplications, followed by normalization.
>
> % matlab
> function x=beliefprop(A,x,N)
> for i=1:N
> x = A*x;
> x = x/sum(x);
> end
> end
>
> // C
> void beliefprop(double A[25][25], double x[25], int N){
> int i, n, j;
> double x2[25];
> for(n=0; n<N; n++){
> for(i=0; i<25; i++){
> x2[i]=0;
> for(j=0; j<25; j++)
> x2[i] += A[i][j]*x[j];
> }
> for(i=0; i<25; i++)
> x[i]=x2[i];
> double mysum = 0;
> for(i=0; i<25; i++)
> mysum += x[i];
> for(i=0; i<25; i++)
> x[i] /= mysum;
> }
> return;
> }
>
> % julia
> function beliefprop(A,x,N)
> for i=1:N
> x = A*exp(x);
> x /= sum(x);
> end
> x
> end
>
> Here, I think we can agree that Matlab and Julia are clearer. (Please
> don't make fun of me for hardcoding the 25 dimensions in C.) Using a
> matrix package for C would probably improve clarity, but perhaps also
> slow things down. The results are:
>
> time in matlab (beliefprop): 0.627478
> time in c (beliefprop): 0.073564
> time in julia (beliefprop): 0.489665
>
>
> Fifth Benchmark: BP in log-space
>
> In practice, Belief Propagation is often implemented in log-space (to
> help avoid numerical under/over-flow.). To simulate an algorithm like
> this, I tried changing to propagation to take an exponent before
> multiplication, and a logarithm before storage.
>
> % matlab
> function x=beliefprop2(A,x,N)
> for i=1:N
> x = log(A*exp(x));
> x = x - log(sum(exp(x)));
> end
> end
>
> // C
> void beliefprop2(double A[25][25], double x[25], int N){
> int i, n, j;
> double x2[25];
> for(n=0; n<N; n++){
> for(i=0; i<25; i++){
> x2[i]=0;
> for(j=0; j<25; j++)
> x2[i] += A[i][j]*exp(x[j]);
> }
> for(i=0; i<25; i++)
> x[i]=log(x2[i]);
> double mysum = 0;
> for(i=0; i<25; i++)
> mysum += exp(x[i]);
> double mynorm = log(mysum);
> for(i=0; i<25; i++)
> x[i] -= mynorm;
> }
> return;
> }
>
> # julia
> function beliefprop2(A,x,N)
> for i=1:N
> x = log(A*exp(x));
> x -= log(sum(exp(x)));
> end
> x
> end
>
> Life is too short to write C code like that when not necessary. But how
> about the speed, you ask?
>
> time in matlab (beliefprop2): 0.662761
> time in c (beliefprop2): 0.646153
> time in julia (beliefprop2): 0.615113
>
>
> Sixth Benchmark: Markov Chain Monte Carlo
>
> Here, I implement a simple Metropolis algorithm. For no particular
> reason, I use the two-dimensional distribution:
>
> p(x) \propto \exp(\sin(5 x_1) - x_1^2 - x_2^2)
>
> <http://justindomke.files.wordpress.com/2012/09/dist.jpg>
>
> % matlab
> function mcmc(x,N)
> f = @(x) exp(sin(x(1)*5) - x(1)^2 - x(2)^2);
> p = f(x);
> for n=1:N
> x2 = x + .01*randn(size(x));
> p2 = f(x2);
> if rand < p2/p
> x = x2;
> p = p2;
> end
> end
> end
>
> // C
> double f(double *x){
> return exp(sin(x[0]*5) + x[1]*x[1]);
> }
> #define pi 3.141592653589793
> void mcmc(double *x,int N){
> double p = f(x);
> int n;
> double x2[2];
> for(n=0; n<N; n++){
> // run Box_Muller to get 2 normal random variables
> double U1 = ((double)rand())/RAND_MAX;
> double U2 = ((double)rand())/RAND_MAX;
> double R1 = sqrt(-2*log(U1))*cos(2*pi*U2);
> double R2 = sqrt(-2*log(U1))*sin(2*pi*U2);
> x2[0] = x[0] + .01*R1;
> x2[1] = x[1] + .01*R2;
> double p2 = f(x2);
> if(((double)rand())/RAND_MAX< p2/p){
> x[0] = x2[0];
> x[1] = x2[1];
> p = p2;
> }
> }
> }
>
> % julia
> function mcmc(x,N)
> f(x) = exp(sin(x[1]*5) - x[1]^2 - x[2]^2);
> p = f(x);
> for n=1:N
> x2 = x + .01*randn(size(x));
> p2 = f(x2);
> if rand() < p2/p
> x = x2;
> p = p2;
> end
> end
> end
>
> Again, I think that C is far less clear than Matlab or Julia. The
> timings are:
>
> time in matlab (mcmc): 7.747716
> time in c (mcmc): 0.646153
> time in julia (mcmc): 0.613223
>
>
> Table
>
> Matlab C Julia
> fib 14.344 0.005 0.281
> matmult 1.229 0.168 0.505
> numquad 0.446 0.201 0.273
> bp 0.627 0.073 0.489
> bp2 0.662 0.646 0.615
> mcmc 7.747 0.646 0.613
>
>
> Conclusions
>
> I'm sure all these programs can be sped up. In particular, I'd bet that
> an expert could optimize the C code to beat Julia on |bp2| and |mcmc|.
> These are a test of "how fast can Justin Domke make these programs", not
> the intrinsic capabilities of the languages. That said, Julia allows for
> optional type declarations. I did experiment with these but found
> absolutely no speed improvement. (Which is a good or a bad thing,
> depending on how you look at life.)
>
> Another surprise to me was how often Matlab's JIT managed a speed within
> a reasonable factor of C. (Except when it didn't...)
>
> The main thing that at Matlab programmer will miss in Julia is
> undoubtedly plotting. The Julia designers seem to understand the
> importance of this ("non-negotiable"). If Julia equalled Matlab's
> plotting facilities, Matlab would be in real trouble!
>
> Overall, I think that the killer features of freedom, kinda-sorta-C-like
> speed, and ease of use make Julia more likely as a Matlab-killer than
> other projects such as R, Sage, Octave, Scipy, etc. (Not to say that
> those projects have not succeeded in other ways!) Though Julia's
> designers also seem to be targeting current R users, my guess is that
> they will have more success with Matlab folks in the short term, since
> most Matlab functionality (other than plotting) already exists, while
> reproducing R's statistical libraries will be quite difficult. I also
> think that Julia would be very attractive to current users of languages
> like Lush <http://lush.sourceforge.net/>. Just to never write another
> .mex file, I'll very seriously consider Julia for new projects. Other
> benefits such as macros, better parallelism support are just bonuses. As
> Julia continues to develop, it will become yet more attractive.
>
> There was an interesting discussion on Lambda the Ultimate
> <http://lambda-the-ultimate.org/node/4452> about Julia back when it was
> announced
>
> *justindomke <http://justindomke.wordpress.com/author/justindomke/>* |
> September 17, 2012 at 2:37 am | Tags: c++
> <http://justindomke.wordpress.com/?tag=c>, efficiency
> <http://justindomke.wordpress.com/?tag=efficiency>, julia
> <http://justindomke.wordpress.com/?tag=julia>, matlab
> <http://justindomke.wordpress.com/?tag=matlab> | Categories:
> Uncategorized <http://justindomke.wordpress.com/?cat=1> | URL:
> http://wp.me/pgOXE-j8
>
> Comment
> <http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/#respond>
> See all comments
> <http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/#comments>
>
> Unsubscribe or change your email settings at Manage Subscriptions
> <https://subscribe.wordpress.com/?key=4220e3bc2f2ce5bd9d8f568e895ce9d4&email=sibert%40hawaii.edu>.
>
>
> *Trouble clicking?* Copy and paste this URL into your browser:
> http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/
>
>
> Thanks for flying with WordPress.com <http://wordpress.com>
>
>
>
>
>
>
>
>
> ------------------------------
>
> _______________________________________________
> Developers mailing list
> Developers at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/developers
>
>
> End of Developers Digest, Vol 43, Issue 5
> *****************************************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/developers/attachments/20121002/268a5aa9/attachment-0001.html>
More information about the Developers
mailing list