Hello,

I’ve been taking a closer look at the ins and outs of setting constraints and the output from my test files is not what I expect.  

I set up things as follows:

DataSet myData = ReadDataFile("/Users/rachel/HYPHY/test.seq");

DataSetFilter ElementFilter = CreateFilter(myData, 1, "0-71");
DataSetFilter ReferenceFilter = CreateFilter(myData, 1, "72-11366");

HarvestFrequencies(obsFreqsRef, ReferenceFilter, 1, 1, 1);
HarvestFrequencies(obsFreqsEle, ElementFilter, 1, 1, 1);


FMatrix = {{*,mu,mu,mu}{mu,*,mu,mu}{mu,mu,*,mu}{mu,mu,mu,*}};

Model FRef = (FMatrix, obsFreqsRef);
Tree myTreeRef = (((hg18, panTro2), rheMac2), (mm8, rn4), canFam2);

Model FEle = (FMatrix, obsFreqsEle);
Tree myTreeEle = (((hg18, panTro2), rheMac2), (mm8, rn4), canFam2);

I build and optimize the likelihood function:

LikelihoodFunction theLnLik = (ReferenceFilter, myTreeRef, ElementFilter, myTreeEle);
Optimize (params, theLnLik);


I now apply the constraints.  

global G=0;
ReplicateConstraint(“this1.?.mu:=G*this2.?.mu”, myTreeEle, myTreeRef);
Optimize (params, theLnLik);

I get (in part) the following output:

G=0.57508
myTreeRef.hg18.mu=0.00597841
myTreeEle.hg18.mu=G*myTreeRef.hg18.mu=0.00343806
myTreeEle.Node2.mu=G*myTreeRef.Node2.mu=0.00980094

I then add another constraint:

ClearConstraints(myTreeEle.hg18.mu);
global S:>1;
myTreeEle.hg18.mu := S*G*myTreeRef.hg18.mu;
Optimize (params, theLnLik);
:

I get (in part) the following output:

G=0.482008
S=32.4987
myTreeRef.hg18.mu=0.00548192
myTreeEle.hg18.mu=S*G*myTreeRef.hg18.mu=0.0858723
myTreeEle.Node2.mu=G*myTreeRef.Node2.mu=0.00822833

It seems that G and the branch lengths have been re-optimized.   I would like to hold G and the branch lengths constant (at the values obtained under the first set of constraints) and only optimize S.  Is there a way to do this and am I correct in my assessment of what is going on?  


Thanks very much in advance for your help.

Cheers,

Rachel

Dear Rachel,

Can you attach the batch file to the message because some of the characters were lost in copy/paste.

Cheers,
Sergei

Hi Sergei,

Sorry about that.  It looks ok on my browser.  

The batch file is attached to this message.  

Thanks

Rachel

---

Dear Rachel,

After


Code (]
ClearConstraints(myTreeEle.hg18.mu);
[/code):

add [code] G:=G__; /*constraint G to the current value; __ is the in-place evaluation construct */ ReplicateConstraint ("this1.?.?:=this2.?.?__", myTreeRef, myTreeRef); /* constrain branch lengths in myTreeRef to current values */

Cheers,
Sergei

That answers my question.  Thanks very much!

Rachel