 load sched; load PLM_CACHE; load HUMAN_ZCONGR; l={'1.0 c/d','1.4 c/d','2.0 c/d','2.8 c/d','4.0 c/d'};
 descr=zeros(9600,5,20);for k=1:10 d=block_from_cache(sched1,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end;
 for k=11:20 d=block_from_cache(sched2,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end

*****   May 30, 2003   ******************************************
*****   Fitting PLM_Hebb1 to the empirical z-probability profiles

 Mparams.trigger=3.5;Nruns=6;o=zeros(9600,4,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));sst=PLM_summary(st);sst=PLM_summary(st);figure(1);plot_PLM_fit({sst.zP_congr',sst.zP_incon'});figure(2);plot_PLM_dynamics(sst);
 figure(3);M=0.6;for k=1:5;subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -M +M]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;end;b=10;c=Mparams.trigger;subplot(3,2,6);plot(1:9600/b,Wh(36,1:b:9600),'.',[0 0;1 1].*9600/b,[-1 1;-1 1].*c,'r-');axis([0 9600/b -12 +12]);set(gca,'xtick',bnd/b,'xticklabel',[]);grid on;title('Running Average');


%-- Parameter search
 Sparams=PLM_search_params
Sparams = 
    lsfun_name: 'PLM_paramsearch'
    model_name: 'PLM_Zcongr1'
        params: [1x1 struct]
     p2v_templ: {1x7 cell}
     v2p_templ: {1x7 cell}
        bounds: [2x7 double]
         optns: [1x1 struct]

 Sparams.p2v_templ'
ans = 
    'VAL = PARAMS.learn_rate *100 ;'
    'VAL = PARAMS.out_noise ;'
    'VAL = PARAMS.rep_noise ;'
    'VAL = PARAMS.W_init ;'
    'VAL = PARAMS.rep_gain ;'
    'VAL = PARAMS.criterion *10 ;'
    'VAL = PARAMS.trigger ;'


 x0=params2vector(Mparams,Sparams.p2v_templ)
x0 = [0.7000  0.2000  0.5000  0.2000  0.7000  0.7000  3.5000]

%%% 5 runs under sched1 + 5 runs under sched2
 for k=1:5 mdata1(k).stim=descr(:,:,k);mdata1(k).chev=make_chev(9600);end
 for k=1:5 mdata1(5+k).stim=descr(:,:,10+k);mdata1(5+k).chev=make_chev(9600);end
 for k=1:5 mdata2(k).stim=descr(:,:,k+5);mdata2(k).chev=make_chev(9600);end
 for k=1:5 mdata2(5+k).stim=descr(:,:,10+k+5);mdata2(5+k).chev=make_chev(9600);end

%%% Search for parameters using FMINCON
 [opt_par1,ssq1,opt_X1,det1]=paramsearch(mdata,Sparams)
ssq1 = 5.9034
opt_X1 = [0.1657  0.1620  0.0986  0.1739  0.4564  0.3335  3.2921]

 [opt_par2,ssq2,opt_X2,det2]=paramsearch(mdata2,Sparams)
ssq2 = 6.4786
opt_X2 = [0.3251  0.1737  0.2610  0.1790  0.4168  0.2953  3.8191]

%%% Document the performance of the model with these parameters
 Mparams=opt_par1;Nruns=20;o=zeros(9600,4,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));sst=PLM_summary(st);sst=PLM_summary(st);figure(1);plot_PLM_fit({sst.zP_congr',sst.zP_incon'});figure(2);plot_PLM_dynamics(sst);
 figure(3);M=0.6;for k=1:5;subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -M +M]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;end;b=10;c=Mparams.trigger;subplot(3,2,6);plot(1:9600/b,Wh(36,1:b:9600),'.',[0 0;1 1].*9600/b,[-1 1;-1 1].*c,'r-');axis([0 9600/b -12 +12]);set(gca,'xtick',bnd/b,'xticklabel',[]);grid on;title('Running Average');


%%%%%%%%  Include the d-prime profiles into the goodness-of-fit
% See PLM_SUMSQFIT2
 [opt_par2_1,ssq2_1,opt_X2_1,det]=paramsearch(mdata1,Sparams_dpr)
 [opt_par2_2,ssq2_2,opt_X2_2,det]=paramsearch(mdata2,Sparams_dpr)
ssq2_1 = 6.5104
opt_X2_1 = [0.1756  0.2018  0.1441  0.1767  0.7999  0.4736  3.3350]
ssq2_2 = [6.7752
opt_X2_2 = [0.2088  0.2138  0.0786  0.1670  0.7661  0.4560  3.2987]

%%% Document the performance of the model with the new parameters
 descr=zeros(9600,5,100);for k=1:50 d=block_from_cache(sched1,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end;
 for k=51:100 d=block_from_cache(sched2,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end;

 Mparams=opt_par2_1;Nruns=100;o=zeros(9600,4,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));sst=PLM_summary(st);sst=PLM_summary(st);figure(1);plot_PLM_fit({sst.zP_congr',sst.zP_incon'});figure(2);plot_PLM_dynamics(sst);
 figure(3);M=0.6;for k=1:5;subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -M +M]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;end;b=10;c=Mparams.trigger;subplot(3,2,6);plot(1:9600/b,Wh(36,1:b:9600),'.',[0 0;1 1].*9600/b,[-1 1;-1 1].*c,'r-');axis([0 9600/b -12 +12]);set(gca,'xtick',bnd/b,'xticklabel',[]);grid on;title('Running Average');
sst2_1 =
       N_runs: 100
        accZP: [0.8494 0.9137 0.9529 1.2487 0.4864 -0.1277]
         accA: [0.5001 0.5145 0.5221 -0.0161 -0.1943 -0.3381]
         accI: [0.4996 0.5150 0.5220 -0.0147 -0.1942 -0.3368]
sst2_2 =                    % time = 56,000 sec = 15.6 hours
       N_runs: 100
        accZP: [0.8597 0.9102 0.9380 1.2301 0.4781 -0.1220]
         accA: [0.4990 0.5112 0.5177 8.3403e-004 -0.1810 -0.3272]
         accI: [0.4989 0.5129 0.5185 7.9533e-004 -0.1817 -0.3269]


********************************************************************
% June 2, 2003 -- Fix TRIGGER=3.3 and REP_NOISE=0.10
% Run six searches with different CHEVs

 load sched; load PLM_CACHE; load HUMAN_ZCONGR;
 descr=zeros(9600,5,60);for k=1:30 d=block_from_cache(sched1,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end;
 for k=31:60 d=block_from_cache(sched2,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end;foo=squeeze(descr(1,:,:));openvar foo

 for n=1:6 for k=1:5 Mdata(n,k).stim=descr(:,:,k+5*n-5);Mdata(n,k).chev=make_chev(9600);Mdata(n,k+5).stim=descr(:,:,k+5*n+25);Mdata(n,k+5).chev=make_chev(9600);end;end;clear n

 Sparams=PLM_search_params =
    lsfun_name: 'PLM_paramsearch'
    sumsq_name: 'PLM_sumsqfit2'
    model_name: 'PLM_Zcongr1'
        params: [1x1 struct]
     p2v_templ: {1x5 cell}
     v2p_templ: {1x5 cell}
        bounds: [2x5 double]
         optns: [1x1 struct]

 Sparams.params =
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.8000
     rep_noise: 0.1000
     out_noise: 0.2000
      W_minmax: [-1 1]
        W_init: 0.1700
       W0_seed: [35x1 double]
    learn_rate: 0.0020
    runav_rate: 0.9800
     criterion: 0.0450
       trigger: 3.3000
      blk_size: 300

 Sparams.p2v_templ' =
    'VAL = PARAMS.learn_rate *100 ;'
    'VAL = PARAMS.out_noise ;'
    'VAL = PARAMS.W_init ;'
    'VAL = PARAMS.rep_gain ;'
    'VAL = PARAMS.criterion *10 ;'

 for k=1:6 disp(['Data set ',int2str(k),':  ',datestr(now)]);[Sres(k).opt_par,Sres(k).ssq,Sres(k).opt_X,Sres(k).det]=paramsearch(Mdata(k,:),Sparams);[Sres(k).opt_X,Sres(k).ssq],save par6search;end
opt_X = [0.1786 0.2098 0.1670 0.7792 0.4163] ; sumsq2 = 6.3195  % time 4.7 hrs
opt_X = [0.1899 0.2051 0.1658 0.7683 0.4526] ; sumsq2 = 6.2303  % time 12 hrs
% terminated after the second run

 Mparams=Sres(1).opt_par,Nruns=60;o=zeros(9600,4,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));sst=PLM_summary(st);sst=PLM_summary(st);figure(1);plot_PLM_fit({sst.zP_congr',sst.zP_incon'});figure(2);plot_PLM_dynamics(sst);
 figure(3);M=0.6;for k=1:5;subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -M +M]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;end;c=Mparams.trigger;subplot(3,2,6);h=plot(1:9600,Wh(36,1:9600),'.',[0 0;1 1].*9600,[-1 1;-1 1].*c,'r-');set(h(1),'MarkerSize',3);axis([0 9600 -12 +12]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;title('Running Average');

sst1.accZP = [0.8563 0.9313 0.9504 1.2018 0.4461 -0.1426]
E1 = [1.91 1.85  0.97]    % as reported by plot_PLM_fit
sst2.accZP = [0.8653 0.9126 0.9498 1.2415 0.4800 -0.1176]
E2 = [1.97 1.68  0.98]    % as reported by plot_PLM_fit

% Setting TRIGGER=0 because the RUNAV scatterplot reveals that it isn't engaged anyway.
sst1.accZP = [0.8626 0.9386 0.9654 1.2341 0.4622 -0.1363]
E1 = [1.85 2.06  0.99]
sst2.accZP = [0.8635 0.9253 0.9497 1.2830 0.4887 -0.1133]
E2 = [1.88 2.19  1.05]


********************************************************************
% June 4, 2003 -- SSQ2 = SSQ_zP + 2*SSQ_dpr
% Run four more searches with different CHEVs

 Sparams.optns.MaxFunEvals=120;
 for n=1:6 for k=1:5 Mdata(k,n).stim=descr(:,:,k+5*n-5);Mdata(k,n).chev=make_chev(9600);Mdata(n,k+5).stim=descr(:,:,k+5*n+25);Mdata(n,k+5).chev=make_chev(9600);end;end;clear n
 for k=1:6 disp(['Data set ',int2str(k),':  ',datestr(now)]);[Sres(k).opt_par,Sres(k).ssq,Sres(k).opt_X,Sres(k).det]=paramsearch(Mdata(k,:),Sparams);[Sres(k).opt_X,Sres(k).ssq],save par6search;end

   learn_rate  out_noise  W_init  rep_gain  criterion     sumsq2
----------------------------------------------------------------
-   0.001786    0.2098    0.1670    0.7792    0.04163    (6.3195) %...+1*SSQ_dpr
-   0.001899    0.2051    0.1658    0.7683    0.04526    (6.2303) %...+1*SSQ_dpr
1   0.001716    0.2095    0.1628    0.8242    0.04373     7.6846
2   0.002504    0.2168    0.1703    0.7059    0.04310     7.6699
3   0.001728    0.1957    0.1662    0.7847    0.04382     7.5638
4   0.001907    0.2073    0.1651    0.7797    0.04638     7.8081
----------------------------------------------------------------
    0.001923    0.2074    0.1662    0.7737    0.04399

 mean_X=mean(cat(1,Sres_trig33(:).opt_X))
mean_X = [0.1923  0.2074  0.1662  0.7737  0.4399]
 Mparams=vector2params(mean_X,Sparams.params,Sparams.v2p_templ)
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.7737
     rep_noise: 0.1000
     out_noise: 0.2074
      W_minmax: [-1 1]
        W_init: 0.1662
       W0_seed: [35x1 double]
    learn_rate: 0.0019
    runav_rate: 0.9800
     criterion: 0.0440
       trigger: 3.3000
      blk_size: 300

 Nruns=60;o=zeros(9600,4,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));sst=PLM_summary(st);figure(1);plot_PLM_fit({sst.zP_congr',sst.zP_incon'});figure(2);plot_PLM_dynamics(sst);
 figure(3);M=0.6;for k=1:5;subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -M +M]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;end;c=Mparams.trigger;subplot(3,2,6);h=plot(1:9600,Wh(36,1:9600),'.',[0 0;1 1].*9600,[-1 1;-1 1].*c,'r-');set(h(1),'MarkerSize',2);axis([0 9600 -12 +12]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;title('Running Average');
sst = 
       N_runs: 60
        accZP: [0.8541 0.9279 0.9464 1.2295 0.4679 -0.1266]
         accI: [0.4927 0.5068 0.5129 -0.0053 -0.1842 -0.3257]

 Hstats=Hebb_stats(Mparams);print_Hebb_stats(Hstats);

Condition   dprime   bias  : mnYcon  mnYinc  :   sdY   SDcon  SDinc 
-------------------------------------------------------------------
Wsame
  Cp=.23    1.506  -0.047  :  0.111  -0.205  :  0.210  0.025  0.035
  Cp=.15    1.083   0.008  :  0.122  -0.106  :  0.210  0.025  0.042
  Cp=.10    0.659   0.058  :  0.128  -0.011  :  0.210  0.025  0.040
  Mean  +   1.083   0.006  +  0.120  -0.107  +  0.210  0.025  0.039
  Diff     -0.847          :  0.017   0.194  :

Wother
  Cp=.23    1.198   0.254  :  0.380   0.128  :  0.210  0.025  0.044
  Cp=.15    0.796   0.291  :  0.375   0.207  :  0.211  0.028  0.045
  Cp=.10    0.455   0.319  :  0.367   0.271  :  0.210  0.030  0.041
  Mean  +   0.816   0.288  +  0.374   0.202  +  0.211  0.027  0.043
  Diff     -0.743          : -0.013   0.143  :

Wmix
  Cp=.23    1.352   0.103  :  0.245  -0.039  :  0.210  0.025  0.039
  Cp=.15    0.940   0.149  :  0.248   0.051  :  0.210  0.026  0.043
  Cp=.10    0.557   0.189  :  0.247   0.130  :  0.210  0.027  0.040
  Mean  +   0.950   0.147  +  0.247   0.047  +  0.210  0.026  0.041
  Diff     -0.795          :  0.002   0.169  :

W0step
  Cp=.23    1.810   0.321  :  0.517   0.125  :  0.216  0.051  0.070
  Cp=.15    1.174   0.410  :  0.537   0.283  :  0.217  0.052  0.071
  Cp=.10    0.666   0.474  :  0.546   0.402  :  0.216  0.053  0.067
  Mean  +   1.217   0.402  +  0.533   0.270  +  0.216  0.052  0.069
  Diff     -1.145          :  0.029   0.277  :

W0lin
  Cp=.23    1.334   0.253  :  0.395   0.111  :  0.213  0.044  0.054
  Cp=.15    0.837   0.322  :  0.412   0.233  :  0.213  0.046  0.055
  Cp=.10    0.461   0.371  :  0.420   0.322  :  0.213  0.046  0.053
  Mean  +   0.878   0.315  +  0.409   0.222  +  0.213  0.045  0.054
  Diff     -0.873          :  0.025   0.211  :
-------------------------------------------------------------------
Condition   dprime   bias  : mnYcon  mnYinc  :   sdY   SDcon  SDinc 


%%%% Curve fitting
 dpr=mean(sst.dprime,2);edpr=mean(HUMAN_ZCONGR)'.*2; blk1=[1,10:17,26:31]';blk2=[2:9,18:25,32]';

%% Three time variables: Clock time T (in blocks), time since last switch Ts, and time in each context Tc
 T=[1:32]';Ts=[Inf,1:8,1:8,1:8,1:6,1]';Tc=[1,1:8,2:9,9:16,10:15,17]';

%% Switch-cost hypothesis: general curve + switch costs
 par1=lsqcurvefit('PLE_exponent',[2 1 .1 -.3 .5],[T Ts],dpr),par1=lsqcurvefit('PLE_exponent',par1,[T Ts],dpr)
par1 = [1.7837  0.4780  11.6629  0.0884  1.6192]
par1:  F_asympt=1.7837, g=0.4780, F_init=0.9311, tau=11.6629
       swcost=0.0884, delta_dpr=0.1577, tau_s=1.6192

 mdpr1=PLE_exponent(par1);resid1=dpr-mdpr1;mrse1=std(resid1,1),cc=corrcoef(mdpr1,dpr);cc(2)^2
mrse1 = 0.0368 ; R2 = 0.9760

 cc=corrcoef(mdpr1,edpr);cc(2)^2    %  R2 = 0.8642   Sw.cost vs. empir.data
 cc=corrcoef(dpr,edpr);cc(2)^2      %  R2 = 0.8602   Model vs. empir

 plot(blk1,edpr(blk1),'ro',blk2,dpr(blk2),'bs',1:32,dpr,'kx',2:9,mdpr1(2:9),'k.-');
 axis([0 33 0 2.5]);set(gca,'xtick',[0:2:32],'xticklabel','0||4||8||12||16||20||24||28||32');
 xlabel('Block number');ylabel('d''');legend('Data, context A','Data, context B','Model prediction','Model fit',4);
 hold on;plot(1,mdpr1(1),'k.',10:17,mdpr1(10:17),'k.-',18:25,mdpr1(18:25),'k.-',26:31,mdpr1(26:31),'k.-',32,mdpr1(32),'k.');hold off;
% Set page size=[1 3.25 6.5 4.5]; Marker size for data=8; Line width=2 and 1.


%%%%%%%%%%  By difficulty level  %%%%%%%%%%%%%%%%%%%

 dpr3=sst.dprime;d=HUMAN_ZCONGR';edpr3=[d(:,1)+d(:,2),d(:,3)+d(:,4),d(:,5)+d(:,6)];clear d
 [par3,rmse3] = PLE_curvefit('PLE_exponent',dpr3,par1);[par3,par3(:,1).*(1-par3(:,2)),par3(:,1).*par3(:,4),rmse3]

=====  Model
Cp    F_asym    g      tau     swcost  tau_s    F_init  d_dpr    RMSE
-----------------------------------------------------------------------
.23   2.5420  0.4523  11.4255  0.0852  1.2486   1.3923  0.2166  0.0526
.15   1.7517  0.4888  12.1670  0.0995  1.0611   0.8955  0.1743  0.0427
.10   1.0901  0.5229  13.0626  0.0930  3.7201   0.5201  0.1014  0.0435

===== Experimental data, see DPRIME_FITS.TXT
 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.3941  0.4343  1.3543    9.9001    0.2032  0.4865  1.1881    0.1280
.15   1.8690  0.4964  0.9412    8.8682    0.1595  0.2980  1.4397    0.1073
.10   1.2684  0.5339  0.5912   13.4278    0.1690  0.2143  0.5999    0.1072

%%%%%%  Three levels fitted together with common tau, s, and tau_s 
%% 7 parameters total:  params=[D1 D2 D3 , g , tau s tau_s]
%
 foo = [2.5 2.0 1.5 , 0.5 , 10 .2 1] ;
 [par7,rmse7] = PLE_curvefit('PLE_exponent',dpr3(:),foo),[par7,rmse7] = PLE_curvefit('PLE_exponent',dpr3(:),par7)
par7 = [2.5819  1.7341  1.0207   0.4696  11.9176  0.0881  1.3861]
rmse7 =  0.0505

=====  Model
 Cp   F_asym    g     F_init    tau         s      s*D    tau_s     RMSE
-------------------------------------------------------------------------
agg   1.7837  0.4780  0.9311   11.6629    0.0884  0.1577  1.6192   0.0368

.23   2.5819  0.4696  1.3694   11.9176    0.0881  0.2275  1.3861   0.0549
.15   1.7341  0.4696  0.9198   11.9176    0.0881  0.1528  1.3861   0.0452
.10   1.0207  0.4696  0.5414   11.9176    0.0881  0.0899  1.3861   0.0509
-------------------------------------------------------------------------


===== Experimental data, see DPRIME_FITS.TXT and Table 1 in the paper
 Cp   F_asym    g     F_init    tau         s      s*D    tau_s     RMSE
-------------------------------------------------------------------------
agg   1.8415  0.4749  0.9670   10.3764    0.1799  0.3313  1.1393   0.0825

.23   2.4145  0.4660  1.2894   10.1551    0.1837  0.4435  1.1784   0.1308
.15   1.8892  0.4660  1.0089   10.1551    0.1837  0.3470  1.1784   0.1102
.10   1.1974  0.4660  0.6395   10.1551    0.1837  0.2200  1.1784   0.1124
-------------------------------------------------------------------------

 mdpr7=reshape(PLE_exponent(par7),32,3);foo=(dpr3-mdpr7).^2;sqrt(mean(foo))
RMSE = [0.2118  0.1772  0.1877]

 plot(1:32,edpr3(:,1),'b^',1:32,edpr3(:,2),'ro',1:32,edpr3(:,3),'gs',1:32,dpr3(:,1),'bx',2:9,mdpr7(2:9,:),'k.-',1:32,dpr3(:,2),'rx',1:32,dpr3(:,3),'gx',10:17,mdpr7(10:17,:),'k.-',18:25,mdpr7(18:25,:),'k.-',26:31,mdpr7(26:31,:),'k.-',[1 32]',mdpr7([1 32],:),'k.');
 legend('Data, Cp=0.245','Data, Cp=0.160','Data, Cp=0.106','Model prediction','Model fits',4);
 axis([0 33 0 3]);set(gca,'xtick',[0:4:32]);xlabel('Block number');ylabel('d''');
% Set page size=[1 3.25 6.5 4.5]; Line width=1; Open symbols.


********************************************************************
% June 6, 2003 -- Fix TRIGGER=0 to see whether it can be eliminated altogether
% Run five searches with different CHEVs

 Sparams.params.trigger=0;
 for k=1:5 disp(['Data set ',int2str(k),':  ',datestr(now)]);[Sres(k).opt_par,Sres(k).ssq,Sres(k).opt_X,Sres(k).det]=paramsearch(Mdata(k,:),Sparams);[Sres(k).opt_X,Sres(k).ssq],save par6search;end

  learn_rate  out_noise  W_init  rep_gain  criterion     sumsq2
----------------------------------------------------------------
1  0.001318    0.1926    0.1521    0.8700    0.04452     7.4998
2  0.001692    0.2055    0.1987    0.7387    0.04743     8.1379
3  0.001478    0.1994    0.1630    0.7664    0.04253     7.6880
4  0.001595    0.2038    0.1603    0.8071    0.04164     7.3919
5  0.001483    0.1897    0.1579    0.7666    0.04262     7.9744
----------------------------------------------------------------
mn 0.001513    0.1982    0.1664    0.7898    0.04375     7.7384
sd 0.000125    0.0061    0.0165    0.0457    0.00207     0.2810






***********************************************************************
************      ERROR-CORRECTING VERSION 

%%%%%%%%%%%   Finalize parameters and generate accurate fits
%
 clear all ; load PLM_CACHE ; load HUMAN_ZCONGR ; load sched ; l={'1.0 c/d','1.4 c/d','2.0 c/d','2.8 c/d','4.0 c/d'};
 descr=zeros(9600,5,100);for k=1:50 d=block_from_cache(sched1,PLM_CACHE);descr(:,:,k)=cat(1,d{:});d=block_from_cache(sched2,PLM_CACHE);descr(:,:,50+k)=cat(1,d{:});end

 Mparams=PLM_params
Mparams = 
       rep_size: [7 5 1]
        max_act: 0.5000
     fdbk_input: [-2 2]
       rep_gain: 1.5000
      rep_noise: 0
       out_gain: 1.5000
      out_noise: 0.4000
       W_minmax: [-1 1]
         W_init: 0.1600
        W0_seed: [36x1 double]
    switch_bias: 0.1900
      threshold: 0.0830
     learn_rate: 0.0050
     runav_rate: 0.0800
        trigger: 0.4000
       blk_size: 300

 N_runs=100;o=zeros(9600,4,N_runs);for k=1:N_runs [o(:,:,k),Wh]=PLModel1(descr(:,:,k),Mparams);end;st=PLM_stats(descr(:,:,1:N_runs),o(:,1,:),o(:,2,:),o(:,3,:),0);sst=PLM_summary(st);figure(1);plot_PLM_dynamics(sst);figure(2);plot_PLM_fit({sst.zP_congr',sst.zP_incon'})
sst =  N_runs: 100
        accZP: [0.9731 0.9881 0.9570   1.1277  0.5317 -0.0016]
         accA: [0.1946 0.1968 0.1944   0.0881  0.0061 -0.0701]
         accI: [0.6606 0.6689 0.6596   0.1980 -0.0439 -0.2697]
% SSError reported by PLOT_PLM_FIT:  Congruent E=2.53, Incongruent E=1.35; D-prime E=0.63

 figure(3);for k=1:5 subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -0.5 +0.5]);set(gca,'xtick',[0:1200:9600],'xticklabel',[],'ytick',[-5 -3 -1 0 1 3 5]/10);grid on;title(l{k});end;subplot(3,2,6);b=Wh(36,:)';rav=Wh(37,:)';th=sign(Wh(38,:)').*Mparams.trigger;hh=plot([b rav th],'.');set(hh,'MarkerSize',2);axis([0 9600 -.5 +.5]);set(gca,'xtick',[0:1200:9600],'xticklabel',[],'ytick',[-5 -3 -1 0 1 3 5]/10);

 Z100=zeros(6,32,100);for k=1:100 Z100(:,:,k)=PLM_Zcongr1(descr(:,:,k),Mparams);end
 Z=mean(Z100,3);plot_PLM_fit(Z)
% SSError reported by PLOT_PLM_FIT:  Congruent E=2.36, Incongruent E=1.29; D-prime E=0.52


*************************************************************************
*************************************************************************

%%%%%%%%%%%   Fit a reduced model to the three d-prime curves
% April 18, 2003
%
 clear all ; load PLM_CACHE ; load HUMAN_ZCONGR ; load sched ;
 descr=zeros(9600,5,20);for k=1:10 d=block_from_cache(sched1,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end
 for k=11:20 d=block_from_cache(sched2,PLM_CACHE);descr(:,:,k)=cat(1,d{:});end

 Sparams=PLM_search_params;Sparams.model_name='PLM_dprime3';
 Sparams.p2v_templ(3:4)=[];Sparams.v2p_templ(3:4)=[];Sparams.bounds(:,3:4)=[]
Sparams = 
    lsfun_name: 'PLM_paramsearch'
    model_name: 'PLM_dprime3'
        params: [1x1 struct]
     p2v_templ: {'VAL = PARAMS.learn_rate *100 ;'  'VAL = PARAMS.out_noise ;'}
     v2p_templ: {'PARAMS.learn_rate = VAL /100 ;'  'PARAMS.out_noise = VAL ;'}
        bounds: [0.1  0.1 ; 1.0 0.5]
         optns: [1x1 struct]

% Disable the running-average mechanism
 Sparams.params.trigger=9999;Sparams.params.switch_bias=0;Sparams.params.threshold=0;
Sparams.params = 
       rep_size: [7 5 1]
        max_act: 0.5000
     fdbk_input: [-2 2]
       rep_gain: 1.5000
      rep_noise: 0
       out_gain: 1.5000
      out_noise: 0.4000
       W_minmax: [-1 1]
         W_init: 0.1600
        W0_seed: [36x1 double]
    switch_bias: 0
      threshold: 0
     learn_rate: 0.0050
     runav_rate: 0.0800
        trigger: 9999
       blk_size: 300


% 6 optimizations, sum-of-squares calculated over 6 mdata sets within each optimization
 for k=1:6 idx=randperm(10);for l=1:3 mdata6(l,k).stim=descr(:,:,idx(l));mdata6(3+l,k).stim=descr(:,:,idx(5+l)+10);mdata6(l,k).chev=make_chev(9600);mdata6(3+l,k).chev=make_chev(9600);end;end
 opt_X=zeros(6,2);for k=1:6 disp(['Data set ',int2str(k),':  ',datestr(now)]);[opt_par(k),ssq(k),opt_X(k,:),det(k)]=paramsearch(mdata6(:,k),Sparams);opt_X,ssq,search_results(k).opt_params=opt_par;search_results(k).ssq=ssq;disp('Saving to foo.mat');save foo;end
 [det(:).time]/60
    3.12  3.26  3.30  3.54  4.33  4.61   % minutes on the GX400
 [det(:).N_evals] = [19  19  19  20  19  20]

 [ssq',opt_X]        % [ssq',opt_X./repmat([100 1],6,1)]
    sum_sq  learn_rate   o_noise
--------------------------------
    2.3698    0.004534    0.3884
    2.3676    0.004534    0.3883
    2.5137    0.004534    0.3883
    2.4279    0.004509    0.3895
    2.5055    0.004534    0.3884
    2.4562    0.004511    0.3895
--------------------------------
    2.4401    0.004526    0.3887


%%%%%%%%%%%   Finalize parameters and generate accurate fits with the reduced model
%
 clear all ; load PLM_CACHE ; load HUMAN_ZCONGR ; load sched ; l={'1.0 c/d','1.4 c/d','2.0 c/d','2.8 c/d','4.0 c/d'};
 descr=zeros(9600,5,100);for k=1:50 d=block_from_cache(sched1,PLM_CACHE);descr(:,:,k)=cat(1,d{:});d=block_from_cache(sched2,PLM_CACHE);descr(:,:,50+k)=cat(1,d{:});end

 Mparams=PLM_params;Mparams.trigger=9999;Mparams.switch_bias=0;Mparams.threshold=0;
 Mparams.learn_rate=0.00453;Mparams.out_noise=0.389
Mparams = 
       rep_size: [7 5 1]
        max_act: 0.5000
     fdbk_input: [-2 2]
       rep_gain: 1.5000
      rep_noise: 0
       out_gain: 1.5000
      out_noise: 0.3890
       W_minmax: [-1 1]
         W_init: 0.1600
        W0_seed: [36x1 double]
    switch_bias: 0
      threshold: 0
     learn_rate: 0.0045
     runav_rate: 0.0800
        trigger: 9999
       blk_size: 300

 N_runs=100;o=zeros(9600,4,N_runs);for k=1:N_runs [o(:,:,k),Wh]=PLModel1(descr(:,:,k),Mparams);end;st=PLM_stats(descr(:,:,1:N_runs),o(:,1,:),o(:,2,:),o(:,3,:),0);sst=PLM_summary(st);figure(1);plot_PLM_dynamics(sst);figure(2);plot_PLM_fit({sst.zP_congr',sst.zP_incon'})
sst =  N_runs: 100
        accZP: [0.9396 0.9712 0.9651 1.1245 0.5390 -0.0179]
         accA: [0.1231 0.1270 0.1280 0.1437 0.0717 -0.0032]
         accI: [0.4470 0.4609 0.4617 0.3540 0.1316 -0.0869]
% SSError reported by PLOT_PLM_FIT:  Congruent E=20.47, Incongruent E=17.95; D-prime E=0.59

 figure(3);for k=1:5 subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -0.5 +0.5]);set(gca,'xtick',[0:1200:9600],'xticklabel',[],'ytick',[-5 -3 -1 0 1 3 5]/10);grid on;title(l{k});end;subplot(3,2,6);b=Wh(36,:)';rav=Wh(37,:)';th=sign(Wh(38,:)').*Mparams.trigger;hh=plot([b rav th],'.');set(hh,'MarkerSize',2);axis([0 9600 -.5 +.5]);set(gca,'xtick',[0:1200:9600],'xticklabel',[],'ytick',[-5 -3 -1 0 1 3 5]/10);

 Z100=zeros(6,32,100);for k=1:100 Z100(:,:,k)=PLM_Zcongr1(descr(:,:,k),Mparams);end
 Z=mean(Z100,3);plot_PLM_fit(Z)
% SSError reported by PLOT_PLM_FIT:  Congruent E=19.22, Incongruent E=16.86; D-prime E=0.54

% Going back to default params
 Mparams.learn_rate=0.005;Mparams.out_noise=0.4
Mparams = 
       rep_size: [7 5 1]
        max_act: 0.5000
     fdbk_input: [-2 2]
       rep_gain: 1.5000
      rep_noise: 0
       out_gain: 1.5000
      out_noise: 0.4000
       W_minmax: [-1 1]
         W_init: 0.1600
        W0_seed: [36x1 double]
    switch_bias: 0
      threshold: 0
     learn_rate: 0.0050
     runav_rate: 0.0800
        trigger: 9999
       blk_size: 300

 Z100=zeros(6,32,100);for k=1:100 Z100(:,:,k)=PLM_Zcongr1(descr(:,:,k),Mparams);end
 Z=mean(Z100,3);plot_PLM_fit(Z)
% SSError reported by PLOT_PLM_FIT:  Congruent E=19.22, Incongruent E=16.86; D-prime E=0.54

 foo=Z([1 3 5 2 4 6],:)';foo=[foo (foo(:,1:3)+foo(:,4:6))./2];describe(foo)
 for k=1:32 fprintf('\n%2d    ',k);fprintf(' %6.3f',foo(k,:));end

% SSError between Z_FULL and Z_REDUC:  Congruent E=13.40, Incongruent E=13.87; D-prime E=0.10
