diff --git a/deb b/deb
index 4f2efc7763b33da14d6c63d5e37dbab5041f7928..aa07200151af66580d88ab2f0ad317b450286c78 100755
--- a/deb
+++ b/deb
@@ -1,22 +1,144 @@
-9,11d8
-< #if !defined(__base2z_h)
-< #define __base2z_h
-<
-13a11,12
-> #include "newmatap.h"
-> #include "newmatio.h"
-15c14,19
-< namespace MISCMATHS {
----
-> using namespace NEWMAT;
->
-> namespace UTILS{
->
-> #if !defined(__base2z_h)
-> #define __base2z_h
-38a43,44
->
-> #endif
-41,42d46
-<
-< #endif
+Log directory is: /usr/people/woolrich/fsl/nsrc/smmfast/logdir++++
+epibt =1000
+spatialdatafile =/usr/people/woolrich/matlab/smm/data/art_z_checker
+maskfile =/usr/people/woolrich/matlab/smm/data/art_mask_checker
+mrf_precision=0.25
+means=-0.187964 2.466296 -5.830409
+
+vars=0.781579 2.096998 1.603299
+
+props=0.559440 0.127918 0.312642
+
+it=1
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=1988.99
+smmfunc.evaluate_withmrfprec(latest)=-1702.47
+latest(nclasses*num_superthreshold+1)=71.6585
+mrf_precision=71.6585
+it=2
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-1702.47
+smmfunc.evaluate_withmrfprec(latest)=-1969.63
+latest(nclasses*num_superthreshold+1)=114.107
+mrf_precision=114.107
+it=3
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-1969.63
+smmfunc.evaluate_withmrfprec(latest)=-2069.99
+latest(nclasses*num_superthreshold+1)=135.829
+mrf_precision=135.829
+it=4
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2069.99
+smmfunc.evaluate_withmrfprec(latest)=-2119.04
+latest(nclasses*num_superthreshold+1)=147.931
+mrf_precision=147.931
+it=5
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2119.04
+smmfunc.evaluate_withmrfprec(latest)=-2168.2
+latest(nclasses*num_superthreshold+1)=161.158
+mrf_precision=161.158
+it=6
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2168.2
+smmfunc.evaluate_withmrfprec(latest)=-2202.21
+latest(nclasses*num_superthreshold+1)=170.851
+mrf_precision=170.851
+it=7
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2202.21
+smmfunc.evaluate_withmrfprec(latest)=-2257.59
+latest(nclasses*num_superthreshold+1)=187.959
+mrf_precision=187.959
+it=8
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2257.59
+smmfunc.evaluate_withmrfprec(latest)=-2291.83
+latest(nclasses*num_superthreshold+1)=199.568
+mrf_precision=199.568
+it=9
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2291.83
+smmfunc.evaluate_withmrfprec(latest)=-2320.1
+latest(nclasses*num_superthreshold+1)=209.342
+mrf_precision=209.342
+it=10
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2320.1
+smmfunc.evaluate_withmrfprec(latest)=-2339.78
+latest(nclasses*num_superthreshold+1)=216.569
+mrf_precision=216.569
+it=11
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2339.78
+smmfunc.evaluate_withmrfprec(latest)=-2362.63
+latest(nclasses*num_superthreshold+1)=225.299
+mrf_precision=225.299
+it=12
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2362.63
+smmfunc.evaluate_withmrfprec(latest)=-2374.07
+latest(nclasses*num_superthreshold+1)=229.721
+mrf_precision=229.721
+it=13
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2374.07
+smmfunc.evaluate_withmrfprec(latest)=-2385.38
+latest(nclasses*num_superthreshold+1)=234.241
+mrf_precision=234.241
+it=14
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2385.38
+smmfunc.evaluate_withmrfprec(latest)=-2397.25
+latest(nclasses*num_superthreshold+1)=239.055
+mrf_precision=239.055
+it=15
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2397.25
+smmfunc.evaluate_withmrfprec(latest)=-2412.04
+latest(nclasses*num_superthreshold+1)=245.244
+mrf_precision=245.244
+it=16
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2412.04
+smmfunc.evaluate_withmrfprec(latest)=-2422.29
+latest(nclasses*num_superthreshold+1)=249.555
+mrf_precision=249.555
+it=17
+sit=1
+"Updating distribution params"=Updating distribution params
+"Updating weights"=Updating weights
+smmfunc.evaluate_withmrfprec(latest)=-2422.29
diff --git a/minimize.cc b/minimize.cc
index a6dd5040943d5cb8b53b07b08d5d6d389eac3fa3..eb199a825514dbfd91f465c8d3cfe3a63d09a7e1 100644
--- a/minimize.cc
+++ b/minimize.cc
@@ -350,9 +350,7 @@ void fminsearch(ColumnVector& x, const EvalFunction& func){
}
-
-void scg(ColumnVector& x,const gEvalFunction& func){
- int niters=100; //default maybe make an options option
+void scg(ColumnVector& x,const gEvalFunction& func, float tol, float eps, int niters){
int fevals=0;
int gevals=0;
@@ -371,8 +369,7 @@ void scg(ColumnVector& x,const gEvalFunction& func){
float lambdamax = 1.0e15;
int j = 1;
float mu=0,kappa=0,sigma=0,gamma=0,alpha=0,delta=0,Delta,beta=0;
- float eps=1e-16;
- float tol=0.0000001;
+
// main loop..
while(j<niters){
diff --git a/minimize.h b/minimize.h
index c0c67b98343423ccf2a0e48bcad007d5ec37bb29..39de572199929d8f81f16c3b9b488f8de3ec3bde 100644
--- a/minimize.h
+++ b/minimize.h
@@ -81,7 +81,7 @@ ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func,float h,int
void fminsearch(ColumnVector& x, const EvalFunction& func);
-void scg(ColumnVector& x, const gEvalFunction& func);
+void scg(ColumnVector& x, const gEvalFunction& func, float tol = 0.0001, float eps=1e-16, int niters=100);
}
diff --git a/miscprob.h b/miscprob.h
index b6bd56047fea396b52e4d4e4763282605be9c88f..b50691da906d2fae6aa2cbb9204cf148bd2c0907 100644
--- a/miscprob.h
+++ b/miscprob.h
@@ -26,6 +26,7 @@ namespace MISCMATHS {
ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1,
const float mu = 0, const float sigma = 1);
+ // returns nsamps*nparams matrix:
ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp = 1);
diff --git a/sparsefn.cc b/sparsefn.cc
index 625d45c5b4ab1a5233cbb4f6190c84f93ae9120e..cf97bed1c8f791656bacac54b2f3021244489596 100644
--- a/sparsefn.cc
+++ b/sparsefn.cc
@@ -341,12 +341,12 @@ namespace MISCMATHS {
}
- float solvefortracex(const SparseMatrix& A, const SparseMatrix& b, SparseMatrix& x)
+ float solvefortracex(const SparseMatrix& A, const SparseMatrix& b, SparseMatrix& x, int nsamps, float tol)
{
Tracer_Plus trace("sparsefns::solvefortracex");
- int every = Max(1,A.Ncols()/100);
- //int every = 1;
+ int every = Max(1,A.Ncols()/nsamps);
+ // int every = 1;
OUT(every);
float tr = 0.0;
@@ -360,7 +360,7 @@ namespace MISCMATHS {
ColumnVector br = b.RowAsColumn(r);
ColumnVector xr = x.RowAsColumn(r);
- solveforx3(A,br,xr);
+ solveforx3(A,br,xr,tol);
for(int c = 1; c<=b.Ncols(); c++)
{
@@ -456,7 +456,7 @@ namespace MISCMATHS {
// write_ascii_matrix("rho",rho);
}
- void solveforx3(const SparseMatrix& A, const ColumnVector& b, ColumnVector& x)
+ void solveforx3(const SparseMatrix& A, const ColumnVector& b, ColumnVector& x, float tol)
{
Tracer_Plus trace("sparsefns::solveforx3");
@@ -478,7 +478,7 @@ namespace MISCMATHS {
ColumnVector w;
ColumnVector p = r;
- while(std::sqrt(rho(k))>0.2*norm2(b) && k < kmax)
+ while(std::sqrt(rho(k))>tol*norm2(b) && k < kmax)
{
k++;
//if(k>2)
diff --git a/sparsefn.h b/sparsefn.h
index 740a397ad8a6f24dbdaca48f74c6b51f0eb1105c..7c2309813632ef0f2c7923627f6506691d0c33ee 100644
--- a/sparsefn.h
+++ b/sparsefn.h
@@ -31,8 +31,8 @@ namespace MISCMATHS {
void solveforx(const SparseMatrix& A, const ColumnVector& b, ColumnVector& x);
void solveforx(const SparseMatrix& A, const ColumnVector& b, SparseMatrix& x);
void solveforx(const SparseMatrix& A, const SparseMatrix& b, SparseMatrix& x);
- float solvefortracex(const SparseMatrix& A, const SparseMatrix& b, SparseMatrix& x);
- void solveforx3(const SparseMatrix& A, const ColumnVector& b, ColumnVector& x);
+ float solvefortracex(const SparseMatrix& A, const SparseMatrix& b, SparseMatrix& x, int nsamps = 50, float tol = 0.001);
+ void solveforx3(const SparseMatrix& A, const ColumnVector& b, ColumnVector& x, float tol = 0.001);
void solve(const SparseMatrix& A, const Matrix& b, SparseMatrix& x);
void addto(SparseMatrix& A, const SparseMatrix& B, float S);
void symmetric_addto(SparseMatrix& A, const SparseMatrix& B, float S);