diff --git a/Beauchamp_Cell2020.m b/Beauchamp_Cell2020.m new file mode 100644 index 0000000..eb13d5c --- /dev/null +++ b/Beauchamp_Cell2020.m @@ -0,0 +1,176 @@ +% Beauchamp Cell 2020.m +% +% Beauchamp, M. S., Oswalt, D., Sun, P., Foster, B. L., Magnotti, J. F., Niketeghad, S., ... & Yoshor, D. (2020). +% Dynamic stimulation of visual cortex produces form vision in sighted and blind humans. Cell, 181(4), 774-783. + +clear all; +c.efthr = 0.05; % what mag of electric field goes through the model, just a speed thing +v.drawthr = 1; + +PARALLEL = 0; +if PARALLEL + numCores = feature('numcores'); + p = parpool(numCores); +end + +%% define cortical and visual space +c.cortexHeight = [-35,35]; % degrees top to bottom, degrees LR, +c.cortexLength = [-10, 80]; +c.pixpermm = 8; c.e.radius = 0.25; +c = p2p_c.define_cortex(c); + +v.visfieldHeight = [-40, 40]; v.visfieldWidth= [-40,40]; v.pixperdeg = 12; +v = p2p_c.define_visualmap(v); +[c, v] = p2p_c.generate_corticalmap(c, v); + +%% figure 4 +if 1% figure 4, YBN doing letters + % folder CharacterDiscriminationYBN + [all_v, xy] = Beauchamp_getDataFig4(); + letter(1).name = 'C'; letter(1).order = [5 1 3 8 12 20 18 ]; + letter(2).name = 'N'; letter(2).order = [3 8 12 20 11 5 9 17]; + letter(3).name = 'S'; letter(3).order =[1 3 8 12 11 10 9 13 17 18 19 20]; + letter(4).name = 'U'; letter(4).order = [8 12 20 19 14 9 5 ]; + n_electrodes = 20; + trl.amp = 1.5; trl.pw = 100*10^-6; trl.dur = 50*10^-3; trl.freq = 200; trl.lag = 50*10^-3; + lag = 100*10^-3; % time between each electrode + trl.order = -1; +else % 03-328, figure 6 + % folder CharacterDiscrimination03281 + trl.amp = 5; + trl.pw = 100*10^-6; + trl.dur = 100*10^-3; + trl.freq = 120; + lag=100*10^-3; % time between each electrode +end + +% do the temporal pulse trains +for tt = 1:12 + trl.lag = trl.lag +lag; + trl.simdur = 6; + tp = p2p_c.define_temporalparameters(); + trl = p2p_c.define_trial(tp,trl); + all_trl(tt) = p2p_c.convolve_model(tp, trl); % create response to the pulse +end + +% locations of the phosphenes, taken from Fig. 2 +figure(1); clf +for ee = 1:length(all_v.e) + polarplot(all_v.e(ee).ang.*pi/180, all_v.e(ee).ecc, 'ko'); hold on + text(all_v.e(ee).ang.*pi/180, all_v.e(ee).ecc, num2str(ee)); + ax = gca; ax.RLim = [0 13]; +end + +for ii = 1:n_electrodes + disp([' generating maps for letter ', num2str(ii), 'out of ', num2str(n_electrodes)]); + v.e = all_v.e(ii); + c = p2p_c.define_electrodes(c, v); + c = p2p_c.generate_ef(c); + v = p2p_c.generate_corticalelectricalresponse(c, v); % create receptive field map for each electrode + all_rfmaps(ii) = v; % generate the receptive field map for each electrode +end + +for ex = 1; %:length(letter) + figure(ex); clf + p2p_c.plotretgrid((img./max(img(:)))*256, v, gray(256), ex, ['';]); drawnow; hold on + % open the video + filename = ['Beauchamp_Cell2020', letter(ex).name]; vid = VideoWriter([filename, '.avi']); + vid.FrameRate = 30; open(vid); + for t = 0:(6*vid.FrameRate)-1 %for each video frame + tmpt = 1+ceil((t/30)*tp.tSamp); + img = zeros(size(squeeze(all_rfmaps(1).e.rfmap(:, :, 1)))); + for ii = 1:length(letter(ex).order) % add up the response of all of the electrodes at that moment in time + trl = all_trl(ii); % define the pulse train + v = all_rfmaps(letter(ex).order(ii)); % define the electrode + img = img + (mean(v.e.rfmap, 3)*trl.resp(tmpt)); + end + p2p_c.plotretgrid(img*500, v, gray(256), ex); drawnow + frame = getframe(gca); + writeVideo(vid,frame); + end + close(vid); +end + + + +% +% %% creating movies +% +% +% +% for t = 1:length(out(exp).val(1).resp_i) +% % sim +% img = zeros(size(squeeze(v.e(1).rfmap(:, :, 1)))); +% for ee=1:length(letter(exp).order) +% e = letter(exp).order(ee); +% img = img + max(v.e(e).rfmap,[], 3).*out(exp).val(end).resp_i(t); % just use the timing of an electrode simulated in the middle of the trial +% end p2p_c.plotretgrid(img*23, v, [], 1); +% +% t = text(-17, 15, 'simultaneous stimulation'); +% set(t, 'Color', [1 1 1]) +% frame = getframe(gca); +% writeVideo(vid,frame); +% end +% for t = 1:length(out(exp).val(1).resp_i) +% % seq +% img = zeros(size(squeeze(v.e(1).rfmap(:, :, 1)))); +% for ee=1:length(letter(exp).order) +% e = letter(exp).order(ee); +% img = img + max(v.e(e).rfmap,[], 3).*out(exp).val(ee).resp_i(t); % just use the timing of an electrode simulated in the middle of the trial +% end +% p2p_c.plotretgrid(img*23, v, [], 1); hold on +% t = text(-17, 15, 'sequential stimulation'); +% set(t, 'Color', [1 1 1]) +% frame = getframe(gca); +% writeVideo(vid,frame); +% end +% close(vid); +% end + + + + %% all the functions + function electrode = generate_resp(stim, tp, v) + + for ee=1:length(stim.order) + trl = []; trl.expname = 'Beauchamp_BioRxiv'; + trl.e = stim.order(ee); trl.lag = (ee-1)*100*10^-3;% 50ms blank period between each stimulation + trl = p2p_c.define_trial(tp,trl); + + + trl = p2p_c.generate_phosphene(v, tp, trl); + itpl= round(linspace(1,length(trl.resp), 30*trl.trialdur)); % interpolate to 30fps + resp_i = trl.resp(itpl); + electrode(ee).resp_i = resp_i; + end + end + + function [v, xy] = Beauchamp_getDataFig4() + xy = [5 8.678; 4.0805 7.842; 3.046 7.519; 2.3565 7.495; 6.4945 7.925; 5.2875 7.366; ... + 4.4255 7.049; 3.2185 6.6625; 7.931 6.653; 6.954 6.1595; 5.9195 5.951; 4.8275 5.741; ... + 8.9655 5.539; 7.701 5.3805; 6.8965 4.9505; 5.747 4.796; 10.23 3.916; 9.138 3.9935; ... + 7.701 3.7715; 6.5515 3.7315; 10.23 3.2265; 9.483 3.028; 8.3335 2.586; 6.954 2.309]; + + [theta,rho] = cart2pol(xy(:, 1),xy(:, 2)); + theta = theta * 180/pi; + + for ee = 1:size(xy, 1) + v.e(ee).ang = theta(ee); + v.e(ee).ecc = rho(ee); + end + end + + function [v, xy] = Beauchamp_getDataFig6() + xy = [13.37 2.11; 11.44 2.14; 9.77 2.13; 8.02 2.19; 5.35 2.21]; + + [theta,rho] = cart2pol(xy(:, 1),xy(:, 2)); + theta = theta * 180/pi; + + for ee = 1:size(xy, 1) + v.e(ee).ang = theta(ee); + v.e(ee).ecc = rho(ee); + end + end + + + diff --git a/Beauchamp_Cell2020C.avi b/Beauchamp_Cell2020C.avi new file mode 100644 index 0000000..3aa81ba Binary files /dev/null and b/Beauchamp_Cell2020C.avi differ diff --git a/Bosking_JNeuro_17.asv b/Bosking_JNeuro_17.asv deleted file mode 100644 index 84d365f..0000000 --- a/Bosking_JNeuro_17.asv +++ /dev/null @@ -1,164 +0,0 @@ -% Bosking_JNeuro_17.m -c.e.radius = 0.25; -c.efthr = 0.05; -n_e_sample = 10; % how many electrodes to simulate, if NaN simulates all of them -%% begin by defining the location and size of electrodes -all_v= Bosking_getData(42); - - -% add some electrodes near the fovea, since we couldn't sample those from the paper -morevals = linspace(.01, 2, 8); -for i=43:50 - all_v.e(i).ecc = morevals(i-42); - all_v.e(i).ang = (rand(1)*2*pi)-pi; -end - -ss = randperm(length(all_v.e), n_e_sample); % collect a random subsample of the electrodes - -%% define cortical and visual space -c.cortexHeight = [-35,35]; % degrees top to bottom, degrees LR, -c.cortexLength = [-10, 80]; -c.pixpermm = 6; -c = p2p_c.define_cortex(c); - -v.visfieldHeight = [-60,60]; v.visfieldWidth= [-60,60]; v.pixperdeg = 12; -v = p2p_c.define_visualmap(v); -[c, v] = p2p_c.generate_corticalmap(c, v); -% temporal parameters -tp = p2p_c.define_temporalparameters(); -ampList = [50 100 150 200 250 500,750,1000,1500,2000,3000,4000]; -trl.dur = 200*10^-3; % duration in s -trl.t = 0:tp.dt:trl.dur-tp.dt; -trl.pw = .1 * 10^-3; % pulse width in s -trl.order = 1; % 1 = cathodic first, -1 = anodic first -trl.freq = 200; % -1 for a single pulse, NaN if not using a temporal model - -% set up plots -figure(1); hold on -p2p_c.plotretgrid(0, v, gray(64), 1, ['']); - -for ii=1:n_e_sample - disp(fprintf('Electrode %d of %d',ii,n_e_sample)); -v.e = all_v.e(ss(ii)); -c.e.radius = 0.25; -c = p2p_c.define_electrodes(c, v); -c = p2p_c.generate_ef(c); - -ampList = [50 100 150 200 250 500,750,1000,1500,2000,3000,4000]; -durList = 200*10^-3*ones(size(ampList)); -pdList = exp(linspace(log(0.05), log(6.4), 40))/1000; -freq = standard_trl.freq*ones(size(pdList)); -dur = standard_trl.dur*ones(size(pdList)); - -Tsim_50 = table(pdList', freq', dur'); - -figure(1) % show electrode location in visual space -plot([v.e.ecc].*cos([v.e.ang]),[v.e.ecc].*sin([v.e.ang]),'ko','MarkerFaceColor','w'); - -v = p2p_c.generate_corticalelectricalresponse(c, v); - for tt=1:length(ampList) - trl.amp = ampList(tt); - trl = rmfield(trl, 'resp'); - trl = p2p_c.define_trial(tp,trl); - trl = p2p_c.generate_phosphene(v, tp, trl); - img = mean(trl.maxphos, 3); % average over both eye-dominant cells - p2p_c.plotretgrid((img./max(img(:)))*256, v, gray(256), 2, ['';]); - trl.sim_radius= mean([trl.ellipse(1).sigma_x trl.ellipse(1).sigma_y]); - trl.sim_diameter = 2 * trl.sim_radius; - trl.sim_brightness = max(trl.maxphos(:)); - trl.maxresp = max(trl.resp); - sim_sizes(ii,tt) = trl.sim_diameter; - disp(['amp = ', num2str(ampList(tt)), ' size ', num2str(sim_sizes(ii,tt))]); - drawnow; - end -end - -%% -% Plot normalized size as a function of Current (figure 3D) - -max_size = max(sim_sizes')'; -y = sim_sizes./repmat(max_size,1,length(ampList)); - -figure(2) -clf -hold on -plot(ampList(1:end)/1000,y(:,1:end),'k-','LineWidth',1); - -plot(ampList(1:end)/1000,y(:,1:end),'s','Color','k','MarkerFaceColor',.75*[1,1,1],'MarkerSize',18,'LineWidth',2); -xlabel('Current (mA)'); -ylabel('Normalized phosphene size'); -set(gca,'XLim',[0,4.100]); -set(gca,'XTick',[0:4]); -set(gca,'YLIm',[0,1.1]); - -set(gca,'FontSize',24); - - -%% -% Plot phosphene size as a function of phosphene eccentricity (for -% stimulation at 1000 microamps) - -figure(3) -clf -idx = find(ampList == 1000); - -hold on -for ee=1:length(v.e) - plot(v.e(ee).ecc, sim_sizes(ee,idx), 'ks', 'MarkerSize', 18,'MarkerFaceColor',[.75,.75,.75],'LineWidth',2); % Figure 4, Bosking 2017 -if ee>42 plot(v.e(ee).ecc, sim_sizes(ee,idx), 'ks', 'MarkerSize', 18,'MarkerFaceColor',[1,.5,.5],'LineWidth',2); % Figure 4, Bosking 2017 -end -end -xlabel('Eccentricity (deg)'); ylabel('Phosphene size (deg)'); -set(gca,'FontSize',20);set(gca, 'YLim', [0 6]) - - -%% -% plot example pulse train and response -dt = trl.t(2)-trl.t(1); -figure(4) -clf -plot(dt:dt:(dt*length(trl.pt)),trl.pt,'k-'); -set(gca,'YTick',[]); -widen -%heighten(2) -xlabel('Time (s)'); -set(gca,'FontSize',16); ylabel('Amplitude'); - -figure(5) -clf -plot(dt:dt:(dt*length(trl.resp)),trl.resp,'k-'); -set(gca,'YTick',[]); widen; -%heighten -xlabel('Time (s)'); set(gca,'FontSize',16); ylabel('Response'); - -%% - -function v = Bosking_getData(varargin) - -if nargin <1 % how many electrodes to simulate - n = 40; -else - n = varargin{1}; -end -% eccentricity and angle values from the Bosking paper -ecc = [2.186915888 2.973407843 2.373639079 2.374361692 3.485162347 3.161720782 3.903266211 2.654157433 2.842470373 ... - 3.071683206 3.535889777 4.046488101 4.091868195 4.414587147 4.833847191 4.974323153 5.392137971 3.956016957 4.467482416 4.931110897 5.58045091 4.098371712 ... - 4.608536468 5.722661143 6.047836979 2.88467097 5.307303208 5.355140187 6.143366413 5.172897196 6.056074766 6.610897004 7.446093073 7.769534637 8.557905386 ... - 3.21230369 8.461219771 9.30089604 9.441372001 10.50852683 8.608921861 9.028326428 8.193852972 6.940697562 14.13532132 13.99513441 14.22694865 14.37537335]; -ecc = repmat(ecc, 1, ceil(n/length(ecc))); -ecc= ecc(randperm(length(ecc))); % select a random subset of electrodes - -ang = [-2.103467982 -1.914786272 -2.00638409 -1.910671717 -1.987623682 -1.600120824 -1.629608464 -1.536051334 -1.626277634 -1.623534598 -1.424713451 ... - -1.423537864 -1.437399993 -1.452241779 -1.265911243 -1.282320478 -1.184061 -0.365460602 -0.195931163 -0.31682069 0.664892401 0.447114914 ... - 2.172582699 2.139960161 2.455605264 2.549358325 2.489207458 2.598977892 2.19330242 -0.76158443 -0.78054077 -0.921954087 -1.030548934 ... - -1.001649087 -1.096185873 -0.902801816 -0.980145642 -0.963344545 -1.029569278 -1.139143781 -1.578029586 -1.094422493 -0.904173334]; -ang = repmat(ang, 1, ceil(n/length(ang))); -ang = ang(randperm(length(ang)))*180/pi; - -for e = 1:n - v.e(e).ang = ang(e); - v.e(e).ecc = ecc(e); -end -end - - diff --git a/Bosking_JNeuro_17.m b/Bosking_JNeuro_17.m index a36ddfb..08b8589 100644 --- a/Bosking_JNeuro_17.m +++ b/Bosking_JNeuro_17.m @@ -1,105 +1,104 @@ % Bosking_JNeuro_17.m c.e.radius = 0.25; -c.efthr = 0.05; -n_e_sample = NaN; % how many electrodes to simulate, if NaN simulates all of them -%% begin by defining the location and size of electrodes -all_v= Bosking_getData(42); - - -% add some electrodes near the fovea, since we couldn't sample those from the paper -morevals = linspace(.01, 2, 8); -for i=43:50 - all_v.e(i).ecc = morevals(i-42); - all_v.e(i).ang = (rand(1)*2*pi)-pi; -end - -if ~isnan(n_e_sample) -ss = randperm(length(all_v.e), n_e_sample); % collect a random subsample of the electrodes -else - ss = 1:length(all_v.e); - n_e_sample = length(all_v.e); +c.efthr = 0.1; +v.drawthr = 5; +n_e_sample = 5; % how many electrodes to simulate, if NaN simulates all of them +% %% begin by defining the location and size of electrodes +% all_v= Bosking_getData(42); +% +% % add some electrodes near the fovea, since we couldn't sample those from the paper +% morevals = linspace(.01, 2, 8); +% for i=43:50 +% all_v.e(i).ecc = morevals(i-42); +% all_v.e(i).ang = (rand(1)*2*pi)-pi; +% end +% +% if ~isnan(n_e_sample); ss = randperm(length(all_v.e), n_e_sample); % collect a random subsample of the electrodes +% else ; ss = 1:length(all_v.e); n_e_sample = length(all_v.e); end + +ecc = linspace(0.01, 25, n_e_sample); +ang = (rand(size(ecc))*2*pi)-pi; +for i = 1:length(ecc) + all_v.e(i).ecc = ecc(i); + all_v.e(i).ang = 0; end + ss = 1:length(all_v.e); %% define cortical and visual space c.cortexHeight = [-35,35]; % degrees top to bottom, degrees LR, -c.cortexLength = [-10, 80]; -c.pixpermm = 6; +c.cortexLength = [-10, 80]; +c.pixpermm = 8; c = p2p_c.define_cortex(c); - -v.visfieldHeight = [-60,60]; v.visfieldWidth= [-60,60]; v.pixperdeg = 12; +vLim = max([all_v.e(:).ecc])*2; +v.visfieldHeight = [-vLim, vLim]; v.visfieldWidth= [-vLim,vLim]; v.pixperdeg = 12; v = p2p_c.define_visualmap(v); [c, v] = p2p_c.generate_corticalmap(c, v); + % temporal parameters tp = p2p_c.define_temporalparameters(); - -% set up plots -figure(1); hold on -p2p_c.plotretgrid(0, v, gray(64), 1, ['']); - -for ii=1:n_e_sample - disp(fprintf('Electrode %d of %d',ii,n_e_sample)); -v.e = all_v.e(ss(ii)); -c.e.radius = 0.25; -c = p2p_c.define_electrodes(c, v); -c = p2p_c.generate_ef(c); - amp = [50 100 150 200 250 500,750,1000,1500,2000,3000,4000]; -dur = 200*10^-3*ones(size(ampList)); -pw = .1 * 10^-3*ones(size(ampList)); -freq = 200 * 10^-3*ones(size(ampList)); - +dur = 200*10^-3*ones(size(amp)); +pw = .1 * 10^-3*ones(size(amp)); +freq = 200 * 10^-3*ones(size(amp)); Tsim = table(amp', pw', freq', dur'); Tsim.Properties.VariableNames = {'amp', 'pw','freq','dur'}; all_trl = p2p_c.loop_convolve_model(tp, Tsim); -figure(1) % show electrode location in visual space -plot([v.e.ecc].*cos([v.e.ang]),[v.e.ecc].*sin([v.e.ang]),'ko','MarkerFaceColor','w'); +% set up plots +figure(1); hold on +p2p_c.plotretgrid(0, v, gray(64), 1, ['';]); -v = p2p_c.generate_corticalelectricalresponse(c, v); +for ii=1:n_e_sample + disp(fprintf('Electrode %d of %d',ii,n_e_sample)); + figure(ii); clf; hold on + v.e = all_v.e(ss(ii)); + c.e.radius = 0.25; + c = p2p_c.define_electrodes(c, v); + c = p2p_c.generate_ef(c); + v = p2p_c.generate_corticalelectricalresponse(c, v); for tt=1:length(amp) trl = all_trl(tt); trl = p2p_c.generate_phosphene(v, tp, trl, all_trl(tt).resp); - img = mean(trl.maxphos, 3); % average over both eye-dominant cells - p2p_c.plotretgrid((img./max(img(:)))*256, v, gray(256), 2, ['';]); + if tt == length(amp) + img = mean(trl.maxphos, 3); % average over both eye-dominant cells + p2p_c.plotretgrid((img./max(img(:)))*256, v, gray(256), ii, ['';]); drawnow; hold on + p2p_c.draw_ellipse(trl, ii,[''], 1,'r'); + end trl.sim_radius= mean([trl.ellipse(1).sigma_x trl.ellipse(1).sigma_y]); trl.sim_diameter = 2 * trl.sim_radius; trl.sim_brightness = max(trl.maxphos(:)); trl.maxresp = max(trl.resp); sim_sizes(ii,tt) = trl.sim_diameter; - disp(['amp = ', num2str(ampList(tt)), ' size ', num2str(sim_sizes(ii,tt))]); - drawnow; + disp(['amp = ', num2str(amp(tt)), ' size ', num2str(sim_sizes(ii,tt))]); end end -% Plot normalized size as a function of Current (figure 3D) -max_size = max(sim_sizes')'; -y = sim_sizes./repmat(max_size,1,length(ampList)); - -figure(2); clf; hold on -plot(ampList(1:end)/1000,y(:,1:end),'k-','LineWidth',1); -plot(ampList(1:end)/1000,y(:,1:end),'s','Color','k','MarkerFaceColor',.75*[1,1,1],'MarkerSize',18,'LineWidth',2); -xlabel('Current (mA)'); -ylabel('Normalized phosphene size'); -set(gca,'XLim',[0,4.100]); -set(gca,'XTick',[0:4]); -set(gca,'YLIm',[0,1.1]); -set(gca,'FontSize',24); - -%% -% Plot phosphene size as a function of phosphene eccentricity (for +%% phosphene size as a function of phosphene eccentricity % stimulation at 1000 microamps) -figure(3); clf; hold on +figure(10); clf; hold on idx = find(amp == 1000); for ii=1:n_e_sample plot( all_v.e(ii).ecc, sim_sizes(ii,idx), 'ks', 'MarkerSize', 12,'MarkerFaceColor',[.75,.75,.75],'LineWidth',2); % Figure 4, Bosking 2017 -if ii>42 - plot(all_v.e(ii).ecc, sim_sizes(ii,idx), 'ks', 'MarkerSize', 12,'MarkerFaceColor',[1,.5,.5],'LineWidth',2); % Figure 4, Bosking 2017 -end + if ii>42 + plot(all_v.e(ii).ecc, sim_sizes(ii,idx), 'ks', 'MarkerSize', 12,'MarkerFaceColor',[1,.5,.5],'LineWidth',2); % Figure 4, Bosking 2017 + end end xlabel('Eccentricity (deg)'); ylabel('Phosphene size (deg)'); set(gca,'FontSize',20);set(gca, 'YLim', [0 6]) +%% Normalized size as a function of Current (figure 3D) +max_size = max(sim_sizes')'; +y = sim_sizes./repmat(max_size,1,length(amp)); +figure(2); clf; hold on +plot(amp(1:end)/1000,y(:,1:end),'k-','LineWidth',1); +plot(amp(1:end)/1000,y(:,1:end),'s','Color','k','MarkerFaceColor',.75*[1,1,1],'MarkerSize',18,'LineWidth',2); +xlabel('Current (mA)'); +ylabel('Normalized phosphene size'); +set(gca,'XLim',[0,4.100]); +set(gca,'XTick',[0:4]); +set(gca,'YLIm',[0,1.1]); +set(gca,'FontSize',24); function v = Bosking_getData(varargin) diff --git a/datasets/Winawer2016_data.xlsx b/datasets/Winawer2016_data.xlsx index 7016f5f..10f9036 100644 Binary files a/datasets/Winawer2016_data.xlsx and b/datasets/Winawer2016_data.xlsx differ diff --git a/fitFrequency.m b/fitFrequency.m index 21897d9..60bc5bd 100644 --- a/fitFrequency.m +++ b/fitFrequency.m @@ -1,15 +1,13 @@ %% fitFrequency.m % -% Fits a variety of frequency data +% Fits a variety of frequency data, examining how threshold varies with +% pulse width % % Written by gmb on feb 20, 2023 tp = p2p_c.define_temporalparameters(); -tp.tau2 = .15; -tp.refrac = 25; % rate of recovery from previous spike -tp.delta = .001; % minimum interspike interval (msec) % Calculate model response to a standard trial clear standard_trl; @@ -22,18 +20,16 @@ standard_trl= p2p_c.convolve_model(tp, standard_trl); tp.thresh_resp = standard_trl.maxresp;% use that as the threshold response -% paperList = {'Dobelle74'; 'Dobelle79'; 'Henderson79' ;'Girvin79'; 'Fernandez21';}; - -paperList = {'Dobelle74'; 'Henderson79' ;'Girvin79'}; % Henderson varied both frequency and duration :( -paperList = {'Dobelle74' ;'Girvin79'}; - -exList = {{'freq'},{'freq'},{'freq'}}; -electrodeList = {{'A1','A2','A3','A4','B2','B3','B4','C1','C2','C3','C4'},{'A'}}; - +% Dobelle 74 only had a limited number of freq/electrodes so including that data has a risk of overfitting +% paperList = {'Dobelle74'; 'Henderson79' ;'Girvin79'}; % Henderson varied both frequency and duration :( % legStr = {'Dobelle 1974, 0.5 msec pulse width, 1 sec duration',... % 'Henderson 1979, 0.5 msec pulse width, variable duration',... % 'Girvin 1979, 0.25 msec pulse width, 0.5 sec duration'}; + +paperList = {'Dobelle74' ;'Girvin79'}; +exList = {{'freq'},{'freq'},{'freq'}}; +electrodeList = {{'A1','A2','A3','A4','B2','B3','B4','C1','C2','C3','C4'},{'A'}}; legStr = {'Dobelle 1974, 0.5 msec pulse width, 1 sec duration',... 'Girvin 1979, 0.25 msec pulse width, 0.5 sec duration'}; @@ -43,9 +39,7 @@ for j=1:length(exList{i}) T = readtable(['datasets/', paperList{i}, '_data.xlsx']); count =count+1; - x{count} = []; - y1{count} = []; - y2{count} = []; + x{count} = []; y1{count} = []; y2{count} = []; for k=1:length(electrodeList{i}) disp(sprintf('Paper: %s, experiment: %s, electrode: %s'... ,paperList{i},exList{i}{j},electrodeList{i}{k})) @@ -79,9 +73,8 @@ end end -%% -% predict smooth freq curve using standard trial parameters. Should go -% through our standard trial's point: amp = 3 at pw = 1msec and freq = 50Hz +%% predict smooth freq curve using standard trial parameters. +% Should go through our standard trial's point: amp = 3 at pw = 1msec and freq = 50Hz % make fake data table freqList = exp(linspace(log(8), log(2000), 40)); @@ -93,26 +86,22 @@ % get thresholds [err, standard_thresh] = p2p_c.loop_find_threshold(tp,Tstandard); -%% -% Plotting (fast after running previous chunks) +%% plotting +% (fast after running previous chunks) -sz = 200; % symbol size +sz = 150; % symbol size xtick = [8,16,32,64,128,256,512,1024]; colList = {'b', 'r', 'c', 'm', 'g'}; alpha = .5; % transparency -sd = 0; % jitter - +sd = 0.15; % jitter fontSize = 12; - clear h -figure(1) -clf -hold on +figure(1); clf; hold on plot(log(freqList),standard_thresh,'k-','LineWidth',2); for i=1:length(x) - h(i)= scatter(x{i}+sd*randn(size(x{i})),y1{i},sz , 'MarkerFaceColor', colList{i},... - 'MarkerEdgeColor','k','MarkerFaceAlpha',alpha); + h(i)= scatter(x{i}+ sd*(rand(size(x{i}))-0.5),y1{i} + sd*(rand(size(x{i}))-0.5) ,sz , 'MarkerFaceColor', colList{i},... + 'MarkerEdgeColor','none','MarkerFaceAlpha',alpha); end set(gca,'XTick',log(xtick)); logx2raw(exp(1),0); ylabel('Threshold'); widen; grid; @@ -122,19 +111,15 @@ set(gca,'FontSize',fontSize); % Figure 2 holds the predicted thresholds for each experiment scaled -% to match the standard stimulus curve. The curves match perfectly - the -% shape - +% to match the standard stimulus curve. The curves match perfectly - the +% effect of frequency is independent of pulse width sd = 0; -figure(2) -clf -hold on +figure(2); clf; hold on plot(log(freqList),standard_thresh,'k-','LineWidth',2); for i=1:length(x) - h(i)= scatter(x{i}+sd*randn(size(x{i})),y2{i},sz , 'MarkerFaceColor', colList{i},... - 'MarkerEdgeColor','k','MarkerFaceAlpha',alpha); + h(i)= scatter(x{i}+sd*(rand(size(x{i}))-0.5),y2{i}+ sd*(rand(size(x{i}))-0.5),sz , 'MarkerFaceColor', colList{i},... + 'MarkerEdgeColor','none','MarkerFaceAlpha',alpha); end - set(gca,'XTick',log(xtick)); logx2raw(exp(1),0); ylabel('Threshold'); widen; grid; xlabel('Frequency (Hz)'); ylabel('Amplitude @ Threshold'); legend(h,legStr,'Location','NorthEast'); diff --git a/fitPulseWidth.m b/fitPulseWidth.m index 131c0dc..50b52df 100644 --- a/fitPulseWidth.m +++ b/fitPulseWidth.m @@ -3,11 +3,10 @@ % % Fits a variety of pulse width data, examining how threshold varies with % pulse width +% +% written gmb 2/2023 tp = p2p_c.define_temporalparameters(); -tp.refrac = 50; -tp.tau2 = .15; -tp.delta = .001; % Calculate model response to a standard trial clear standard_trl; @@ -22,47 +21,42 @@ % Now on to the real experiments... Note, Henderson79 has no pd % experiment, and Girvin79 has two. - +% Data from Brindley et al. 1968 are excluded because measured thresholds +% were unexpectedly non-monotonic as a function of frequency, suggesting +% an electronics issue paperList = {'Dobelle74'; 'Dobelle79'; 'Henderson79' ;'Girvin79'; 'Fernandez21';}; - exList = {{'pd'},{'pd'},{''},{'pd1','pd50'},{'pd'}}; - legStr = {'Dobelle 1974, 50Hz 1 sec duration',... 'Dobelle 1979, 50Hz 0.5 sec duration',... 'Girvin 1979, 1 pulse 0.5 sec duration','Girvin 1979 50Hz, 0.5 sec duration',... 'Fernandez 2021, 300Hz 0.167 sec duration'}; -% Data from Brindley et al. 1968 are excluded because measured thresholds -% were unexpectedly non-monotonic as a function of frequency, suggesting -% an electronics issue - - clear x y1 y2 count = 0; for i = 1:length(paperList) for j=1:length(exList{i}) - disp(sprintf('Paper: %s, experiment: %s',paperList{i},exList{i}{j})) + disp(sprintf('Paper: %s, experiment: %s \n',paperList{i},exList{i}{j})) T = readtable(['datasets/', paperList{i}, '_data.xlsx']); eid =strcmp(T.experiment,exList{i}{j}); if sum(eid) count = count+1; Texp = T(eid,:); [err, thresh] = p2p_c.loop_find_threshold(tp,Texp); - + Tsim = Texp; Tsim.freq = standard_trl.freq*ones(size(Tsim.freq)); Tsim.dur = standard_trl.dur*ones(size(Tsim.dur)); - + % get model response to standard stimulus at these pd's [err, standard_thresh] = p2p_c.loop_find_threshold(tp,Tsim); - + % regression to scale the actual data. scFac1 = standard_thresh'/Texp.amp'; - + % save results to plot later. x{count} = log(Texp.pw*1000); y1{count} = Texp.amp*scFac1; - + scFac2 = standard_thresh'/thresh'; y2{count} = thresh*scFac2; else @@ -71,9 +65,8 @@ end end -%% -% predict smooth pd curve using standard trial parameters. Should go -% through our standard trial's point: amp = 3 at pw = 1msec +%% predict smooth pd curve using standard trial parameters. +% Should go through our standard trial's point: amp = 3 at pw = 1msec % make fake data table pdList = exp(linspace(log(0.05), log(6.4), 40))/1000; @@ -83,28 +76,24 @@ Tstandard = table(pdList', freq', dur'); Tstandard.Properties.VariableNames = {'pw','freq','dur'}; % get thresholds -[err, standard_thresh] = p2p_c.loop_find_threshold(tp,Tstandard); +[err, standard_thresh] = p2p_c.loop_find_threshold(tp,Tstandard); -%% -% Plotting (fast after running previous chunks) +%% Plotting real data +% (fast after running previous chunks) -sz = 200; % symbol size +sz = 150; % symbol size xtick = [.05,.1,.2,.4,.8,1.6,3.2,6.4]; colList = {'b', 'r', 'c', 'm', 'g'}; alpha = .5; % transparency -sd = .1; % jitter +sd = .15; % jitter fontSize = 12; - -figure(1) -clf -hold on +figure(1); clf; hold on plot(log(pdList*1000),standard_thresh,'k-','LineWidth',2); for i=1:length(x) - h(i)= scatter(x{i}+sd*randn(size(x{i})),y1{i},sz , 'MarkerFaceColor', colList{i},... - 'MarkerEdgeColor','k','MarkerFaceAlpha',alpha); + h(i)= scatter(x{i}+sd*(rand(size(x{i}))-.5),y1{i},sz , 'MarkerFaceColor', colList{i},... + 'MarkerEdgeColor','none','MarkerFaceAlpha',alpha); end - set(gca,'XTick',log(xtick)); logx2raw; ylabel('Threshold'); widen; grid; xlabel('Pulse Width (msec)'); ylabel('Amplitude @ Threshold'); legend(h,legStr,'Location','NorthEast'); @@ -114,17 +103,13 @@ % Figure 2 holds the predicted thresholds for each experiment scaled % to match the standard stimulus curve. The curves match perfectly - the % shape - sd = 0; -figure(2) -clf -hold on +figure(2); clf; hold on plot(log(pdList*1000),standard_thresh,'k-','LineWidth',2); for i=1:length(x) - h(i)= scatter(x{i}+sd*randn(size(x{i})),y2{i},sz , 'MarkerFaceColor', colList{i},... - 'MarkerEdgeColor','k','MarkerFaceAlpha',alpha); + h(i)= scatter(x{i}+sd*(rand(size(x{i}))-.5),y2{i},sz , 'MarkerFaceColor', colList{i},... + 'MarkerEdgeColor','none','MarkerFaceAlpha',alpha); end - set(gca,'XTick',log(xtick)); logx2raw; ylabel('Threshold'); widen; grid; xlabel('Pulse Width (msec)'); ylabel('Amplitude @ Threshold'); legend(h,legStr,'Location','NorthEast'); @@ -134,5 +119,3 @@ - - diff --git a/p2p_c.m b/p2p_c.m index 384744a..78b2e77 100644 --- a/p2p_c.m +++ b/p2p_c.m @@ -45,6 +45,9 @@ on = mod(trl.t,1/trl.freq) <=(trl.pw*2); % turn it on on off = mod(trl.t-trl.pw,1/trl.freq) <=trl.pw; tmp = trl.amp.*(on-(2*off)); + if max(abs(tmp))>trl.amp + error('something wrong with generate_pt') + end else on = mod(trl.t,1/trl.freq) c.efthr * 255 - % x0 = c.v2c.X(pixNum); % x center - % y0 = c.v2c.Y(pixNum); % y center - % if strcmp(c.e(idx(ii)).hemi, 'lh') - % x0=-x0; y0 = -y0; - % end - % theta = pi-c.ORmap(pixNum); %orientation - % sigma_x = c.RFmap(pixNum) * c.ar; % major axis sd - % sigma_y = c.RFmap(pixNum); % minor axis sd - % - % % things you need to go from sigmas and theta to G - % aa = cos(theta)^2/(2*sigma_x^2) + sin(theta)^2/(2*sigma_y^2); - % bb = -sin(2*theta)/(4*sigma_x^2) + sin(2*theta)/(4*sigma_y^2); - % cc = sin(theta)^2/(2*sigma_x^2) + cos(theta)^2/(2*sigma_y^2); - % - % % oriented 2D Gaussian - % if strcmp(c.rftype, 'scoreboard') - % % scoreboard version - % tmp = double( c.e(idx(ii)).ef(pixNum))/255; - % rfmap(:, :, 1) = rfmap(:, :, 1) + (tmp * exp(-( (v.X-x0).^2/(0.01) + (v.Y-y0).^2/(.01)))); - % rfmap(:, :, 2) = rfmap(:, :, 1); - % else - % tmp = double(c.e(idx(ii)).ef(pixNum))/255; - % G = tmp * exp( - (aa*(v.X-x0).^2 + 2*bb*(v.X-x0).*(v.Y-y0) + cc*(v.Y-y0).^2)); - % rfmap(:, :, 1) = rfmap(:, :, 1) + c.ODmap(pixNum)*G; - % rfmap(:, :, 2) = rfmap(:, :, 2) + (1-c.ODmap(pixNum))*G; - % end - % end - % end - % if sum(rfmap(:)>0)<20 - % disp('WARNING! Too few pixels passed ef threshold.'); - % disp(' try lowering c.efthr, checking location of electrodes relative to cortical sheet & '); - % disp('checking the sampling resolution of cortex'); - % end - % - % v.e(idx(ii)).rfmap = rfmap./max(rfmap(:)); - % end - % end + function [trl,v] = generate_phosphene(v, tp, trl, varargin) if isnan(trl.freq) trl.maxphos = v.e(trl.e).rfmap; trl.resp = 1;% the scaling due to current integration else - if nargin==4 + if nargin==4 % short cut if already know the max response trl.resp = varargin{1}; else trl = p2p_c.convolve_model(tp, trl); end - end - trl.maxphos = v.e(trl.e).rfmap.*max(trl.resp); - % calculate the neural response over time across the spatial - % maps generated by generate_rfmap - trl.sim_area = (1/v.pixperdeg.^2) * sum(trl.maxphos(:) > v.drawthr)/2; % calculated area of phosphene based on mean of left and right eyes - - if ~isempty(trl.maxphos) - for i=1:2 % left and right eye - p = p2p_c.fit_ellipse_to_phosphene(trl.maxphos(:,:,i)>v.drawthr,v); - trl.ellipse(i).x = p.x0; trl.ellipse(i).y = p.y0; - trl.ellipse(i).sigma_x = p.sigma_x; trl.ellipse(i).sigma_y = p.sigma_y; - trl.ellipse(i).theta = p.theta; - end - % what rule to use to translate phosphene image to brightness? + trl.maxphos = v.e(trl.e).rfmap.*max(trl.resp); - beta = 6; % soft-max rule across pixels for both eyes - trl.sim_brightness = ((1/v.pixperdeg.^2) * sum(trl.maxphos(:).^beta)^(1/beta)); % IF CHECK - else - trl.sim_brightness = []; - end end - function tp = define_temporalparameters(varargin) - if nargin==0 - tp = []; - else - tp = varargin{1}; - end - if ~isfield(tp, 'dt'); tp.dt = .001 * 10^-3; end % 001 time sampling in ms, should be no larger than 1/10 of tau1 - if ~isfield(tp, 'tau1'); tp.tau1 =0.0003; end % fixed based on Nowak and Bullier, 1998 - if ~isfield(tp, 'refrac'); tp.refrac = 50 ; end % refractory period for spikes - if ~isfield(tp, 'delta'); tp.delta = 0.001 ; end % decay parameter for refractory period - - if ~isfield(tp, 'tSamp'); tp.tSamp = 1000;end % subsampling to speed things up - % fit based on Brindley, 1968, Tehovnk 2004 estimates 0.13-0.24 ms - if ~isfield(tp, 'tau2'); tp.tau2 = 0.15; end % 24-33 from retina - if ~isfield(tp, 'ncascades'); tp.ncascades = 3; end% number of cascades in the slow filter of the temporal convolution - if ~isfield(tp, 'gammaflag'); tp.gammaflag = 1; end % include second stage game - - % leak out of charge accumulation - % tp.flag_cl=0; % 1 if you want to charge to leak back out of the system - % tp.tau2_cl = tp.tau2_ca * 8; % used for the conv model, fe uses p.tau2_ca - - % nonlinearity parameters - if ~isfield(tp, 'model'); tp.model = 'compression'; end - if ~isfield(tp, 'power'); tp.power = 20.15; - elseif strcmp(tp.model, 'sigmoid') - disp('using sigmoid semisaturation constant') - tp.asymptote = 2000; - tp.e50 = 500; % electrical semisaturation constant - elseif strcmp(tp.model, 'normcdf') - disp('using normcdf semisaturation constant') - tp.asymptote = 1500; - tp.mean = 750; - tp.sigma = 175; - elseif strcmp(tp.model, 'weibull') - disp('using weibull semisaturation constant') - tp.asymptote = 1000; - tp.thresh = 600; - tp.beta = 3.5; + + % calculate the neural response over time across the spatial + % maps generated by generate_rfmap + trl.sim_area = (1/v.pixperdeg.^2) * sum(trl.maxphos(:) > v.drawthr)/2; % calculated area of phosphene based on mean of left and right eyes + + if ~isempty(trl.maxphos) + for i=1:2 % left and right eye + p = p2p_c.fit_ellipse_to_phosphene(trl.maxphos(:,:,i)>v.drawthr,v); + trl.ellipse(i).x = p.x0; trl.ellipse(i).y = p.y0; + trl.ellipse(i).sigma_x = p.sigma_x; trl.ellipse(i).sigma_y = p.sigma_y; + trl.ellipse(i).theta = p.theta; end + % what rule to use to translate phosphene image to brightness? + + beta = 6; % soft-max rule across pixels for both eyes + trl.sim_brightness = ((1/v.pixperdeg.^2) * sum(trl.maxphos(:).^beta)^(1/beta)); % IF CHECK + else + trl.sim_brightness = []; + end + end + function tp = define_temporalparameters(varargin) + if nargin==0 + tp = []; + else + tp = varargin{1}; + end + if ~isfield(tp, 'dt'); tp.dt = .001 * 10^-3; end % 001 time sampling in ms, should be no larger than 1/10 of tau1 + if ~isfield(tp, 'tau1'); tp.tau1 =0.0003; end % fixed based on Nowak and Bullier, 1998 + if ~isfield(tp, 'refrac'); tp.refrac = 50 ; end % refractory period for spikes + if ~isfield(tp, 'delta'); tp.delta = 0.001 ; end % decay parameter for refractory period + + if ~isfield(tp, 'tSamp'); tp.tSamp = 1000; end % subsampling to speed things up + % fit based on Brindley, 1968, Tehovnk 2004 estimates 0.13-0.24 ms + if ~isfield(tp, 'tau2'); tp.tau2 = 0.15; end % 24-33 from retina + if ~isfield(tp, 'ncascades'); tp.ncascades = 3; end% number of cascades in the slow filter of the temporal convolution + if ~isfield(tp, 'gammaflag'); tp.gammaflag = 1; end % include second stage gamma + if ~isfield(tp, 'simdur'); tp.simdur = 3; end + + % nonlinearity parameters + if ~isfield(tp, 'model'); tp.model = 'compression'; end + if ~isfield(tp, 'power'); tp.power = 20.15; + elseif strcmp(tp.model, 'sigmoid') + disp('using sigmoid semisaturation constant') + tp.asymptote = 2000; + tp.e50 = 500; % electrical semisaturation constant + elseif strcmp(tp.model, 'normcdf') + disp('using normcdf semisaturation constant') + tp.asymptote = 1500; + tp.mean = 750; + tp.sigma = 175; + elseif strcmp(tp.model, 'weibull') + disp('using weibull semisaturation constant') + tp.asymptote = 1000; + tp.thresh = 600; + tp.beta = 3.5; end + end - %% psychophysics - function [err, thresh] = loopall_find_threshold(tp,T) - % - % Runs the 'conv' model to get thresholds based on trials in the table 'T'. - % returns the SSE and thresholds. Table must contain fields holding the - % following parameters for each trial: - % pw pulse width (sec) - % dur trial duration (sec) - % freq pulse frequency (Hz) - % amp amplitude at detection threshold - if ~isfield(tp, 'nReps') - tp.nReps = 12; - end - thresh = NaN(1,size(T,1)); - for i=1:size(T,1) - % define trial parameters based on values in the table - clear trl; trl.pw = T.pw(i); trl.amp = 1; trl.dur = T.dur(i); trl.freq = T.freq(i); trl.simdur = 3; %sec - trl = p2p_c.define_trial(tp,trl); - - tp.scFacInd = 1; % separate 'scFac' for each experiment - if isfield(tp,'experimentList') % set tp_thresh accordingly for this trial - experimentNum = find(strcmp(T.experiment{i},tp.experimentList)); - if ~isempty(experimentNum) % set thresh_resp for this experiment. - tp.scFacInd = experimentNum; - end + %% psychophysics + function [err, thresh] = loopall_find_threshold(tp,T) + % + % Runs the 'conv' model to get thresholds based on trials in the table 'T'. + % returns the SSE and thresholds. Table must contain fields holding the + % following parameters for each trial: + % pw pulse width (sec) + % dur trial duration (sec) + % freq pulse frequency (Hz) + % amp amplitude at detection threshold + if ~isfield(tp, 'nReps') + tp.nReps = 12; + end + thresh = NaN(1,size(T,1)); + for i=1:size(T,1) + % define trial parameters based on values in the table + clear trl; trl.pw = T.pw(i); trl.amp = 1; trl.dur = T.dur(i); trl.freq = T.freq(i); trl.simdur = 3; %sec + trl = p2p_c.define_trial(tp,trl); + + tp.scFacInd = 1; % separate 'scFac' for each experiment + if isfield(tp,'experimentList') % set tp_thresh accordingly for this trial + experimentNum = find(strcmp(T.experiment{i},tp.experimentList)); + if ~isempty(experimentNum) % set thresh_resp for this experiment. + tp.scFacInd = experimentNum; end - thresh(i)= p2p_c.find_threshold(trl,tp); end + thresh(i)= p2p_c.find_threshold(trl,tp); + end - err = nansum((thresh-T.amp').^2); - % disp(sprintf('tau1 = %g, tau2 = %g, power = %5.2f err= %5.4f scFac= %5.4f', tp.tau1,tp.tau2,tp.power,err, tp.scFac)); - disp(fprintf('mean = %g, sigma = %g, err= %5.4f\n', tp.mean,tp.sigma,err)); - if isfield(tp,'experimentList') - for i = 1:length(tp.experimentList) - disp(fprintf('%10s: %g\n',tp.experimentList{i},tp.scFac(i))); - end + err = nansum((thresh-T.amp').^2); + % disp(sprintf('tau1 = %g, tau2 = %g, power = %5.2f err= %5.4f scFac= %5.4f', tp.tau1,tp.tau2,tp.power,err, tp.scFac)); + disp(fprintf('mean = %g, sigma = %g, err= %5.4f\n', tp.mean,tp.sigma,err)); + if isfield(tp,'experimentList') + for i = 1:length(tp.experimentList) + disp(fprintf('%10s: %g\n',tp.experimentList{i},tp.scFac(i))); end end - function [err, thresh] = loop_find_threshold(tp,T) - % - % Runs the 'conv' model to get thresholds based on trials in the table 'T'. - % returns the SSE and thresholds. Table must contain fields holding the - % following parameters for each trial: - % pw pulse width (sec) - % dur trial duration (sec) - % freq pulse frequency (Hz) - % amp amplitude at detection threshold - if ~isfield(tp, 'nReps') - tp.nReps = 12; - end - thresh = NaN(size(T,1), 1); - for i=1:size(T,1) - % define trial parameters based on values in the table - clear trl; trl.pw = T.pw(i); trl.amp = 1; trl.dur = T.dur(i); trl.freq = T.freq(i); trl.simdur = 3; %sec - trl = p2p_c.define_trial(tp,trl); - thresh(i)= p2p_c.find_threshold(trl,tp); - end + end + function [err, thresh] = loop_find_threshold(tp,T) + % + % Runs the 'conv' model to get thresholds based on trials in the table 'T'. + % returns the SSE and thresholds. Table must contain fields holding the + % following parameters for each trial: + % pw pulse width (sec) + % dur trial duration (sec) + % freq pulse frequency (Hz) + % amp amplitude at detection threshold + if ~isfield(tp, 'nReps') + tp.nReps = 12; + end + thresh = NaN(size(T,1), 1); + for i=1:size(T,1) + % define trial parameters based on values in the table + clear trl; trl.pw = T.pw(i); trl.amp = 1; trl.dur = T.dur(i); trl.freq = T.freq(i); trl.simdur = 3; %sec + trl = p2p_c.define_trial(tp,trl); + thresh(i)= p2p_c.find_threshold(trl,tp); + end + + if strmatch('amp',T.Properties.VariableNames) + err = nansum((thresh-T.amp').^2); + disp(sprintf('err= %5.4f', err)); + else + err = NaN; + end - if strmatch('amp',T.Properties.VariableNames) - err = nansum((thresh-T.amp').^2); - disp(sprintf('err= %5.4f', err)); - else - err = NaN; + if isfield(tp,'experimentList') + for i = 1:length(tp.experimentList) + disp(sprintf('%10s: %g',tp.experimentList{i},tp.scFac(i))); end + end + end + function err = fit_brightness(tp, T) + [loop_trl] = p2p_c.loop_convolve_model(tp,T); + y_est = [loop_trl.maxresp]; y = [T.brightness]; + y_est = reshape(y_est, length(y), 1); + y = reshape(y, length(y), 1); + ind = ~isnan(y_est) & ~isnan(y); + err = -corr(y(ind), y_est(ind)); + disp(sprintf('tau1 =%5.4f, tau2 =%5.4f, power =%5.4f, corr = %5.4f', tp.tau1, tp.tau2, tp.power,-err)); + end - if isfield(tp,'experimentList') - for i = 1:length(tp.experimentList) - disp(sprintf('%10s: %g',tp.experimentList{i},tp.scFac(i))); - end - end + function [loop_trl] = loop_convolve_model(tp,T) + % Runs the 'conv' model to get thresholds based on trials in the table 'T'. + % returns the SSE and thresholds. Table must contain fields holding the + % following parameters for each trial: + % pw pulse width (sec) + % dur trial duration (sec) + % freq pulse frequency (Hz) + % amp amplitude at detection threshold + + for i=1:size(T,1) + % define trial parameters based on values in the table + clear trl; trl.pw = T.pw(i); trl.amp = T.amp(i); trl.dur = T.dur(i); trl.freq = T.freq(i); trl.simdur = 3; %sec + trl = p2p_c.define_trial(tp,trl); + loop_trl(i) = p2p_c.convolve_model(tp, trl); end - function err = fit_brightness(tp, T) - [loop_trl] = p2p_c.loop_convolve_model(tp,T); - y_est = [loop_trl.maxresp]; y = [T.brightness]; - y_est = reshape(y_est, length(y), 1); - y = reshape(y, length(y), 1); - ind = ~isnan(y_est) & ~isnan(y); - err = -corr(y(ind), y_est(ind)); - disp(sprintf('tau1 =%5.4f, tau2 =%5.4f, power =%5.4f, corr = %5.4f', tp.tau1, tp.tau2, tp.power,-err)); - end - - function [loop_trl] = loop_convolve_model(tp,T) - % - % Runs the 'conv' model to get thresholds based on trials in the table 'T'. - % returns the SSE and thresholds. Table must contain fields holding the - % following parameters for each trial: - % pw pulse width (sec) - % dur trial duration (sec) - % freq pulse frequency (Hz) - % amp amplitude at detection threshold - - for i=1:size(T,1) - % define trial parameters based on values in the table - clear trl; trl.pw = T.pw(i); trl.amp = T.amp(i); trl.dur = T.dur(i); trl.freq = T.freq(i); trl.simdur = 3; %sec - trl = p2p_c.define_trial(tp,trl); - - % define impulse response - if isfield(tp,'tSamp') - if tp.tSamp~=1% down-sample the time-vectors - t = trl.t(1:tp.tSamp:end); - end - else - t = trl.t; - end - dt = t(2)-t(1); - h = p2p_c.gamma(tp.ncascades,tp.tau2,t); % Generate the n-cascade impulse response - tid = find(cumsum(h)*dt>.999,1,'first'); % Shorten the filter if needed to speed up the code. - if ~isempty(tid) - h = h(1:tid); - else - sprintf('Warning: gamma hdr might not have a long enough time vector'); - end - trl.imp_resp = h; % close enough to use h - loop_trl(i) = p2p_c.convolve_model(tp, trl); + end + function trl = define_impulse_response(tp, trl) + if isfield(tp,'tSamp') + if tp.tSamp~=1% down-sample the time-vectors + t = trl.t(1:tp.tSamp:end); end + else + t = trl.t; end - function amp = find_threshold(trl, tp) - % Find amplitudes at threshold with the convolve model. - % takes in trial, tp, and optional fitParams - % finds and returns the trl.amp for which the max output of the - % model for that trial, trial.resp, is equal to fitParams.thr - - if ~isfield(tp, 'nReps') - tp.nReps = 12; - end - % first find the lowest 'hi' response - hi = 1; - resp = 0; - while resp tp.thresh_resp - hi = trl.amp; - else - lo = trl.amp; - end + dt = t(2)-t(1); + h = p2p_c.gamma(tp.ncascades,tp.tau2,t); % Generate the n-cascade impulse response + tid = find(cumsum(h)*dt>.999,1,'first'); % Shorten the filter if needed to speed up the code. + if ~isempty(tid) + h = h(1:tid); + else + sprintf('Warning: gamma hdr might not have a long enough time vector'); + end + trl.imp_resp = h; % close enough to use h + trl.t = t; + end + function amp = find_threshold(trl, tp) + % Find amplitudes at threshold with the convolve model. + % takes in trial, tp, and optional fitParams + % finds and returns the trl.amp for which the max output of the + % model for that trial, trial.resp, is equal to fitParams.thr + + if ~isfield(tp, 'nReps') + tp.nReps = 12; + end + % first find the lowest 'hi' response + hi = 1; + resp = 0; + while resp tp.thresh_resp + hi = trl.amp; + else + lo = trl.amp; end - amp = (hi+lo)/2; end - function trl = convolve_model(tp, trl) + amp = (hi+lo)/2; + end + function trl = convolve_model(tp, trl) % Implements 'finite_element' using the closed-form solution to % the respose to a pluse Can be faster than 'finite_element'. % Assumes square pulse trains. @@ -686,6 +637,8 @@ ptid = find(diff(trl.pt))+1; % R will hold the values of R1 at the event times. + + Rtmp = zeros(1,length(ptid)); wasRising = 0; % Loop through the events, calculating R1 at the end of the event @@ -717,28 +670,27 @@ trl.spikes = R1(spikeId); ptid = ptid(spikeId); if tp.gammaflag - % if ~isfield(trl, 'imp_resp') % danger - if pre-computed, - % it wont change if tau2 changes. - dt = t(2)-t(1); - h = p2p_c.gamma(tp.ncascades,tp.tau2,t); % Generate the n-cascade impulse response - tid = find(cumsum(h)*dt>=.999,1,'first'); % Shorten the filter if needed to speed up the code. - if ~isempty(tid) - h = h(1:tid); - else - disp('Warning: gamma hdr might not have a long enough time vector'); - end - trl.imp_resp = h; % close enough to use h - % end + + dt = t(2)-t(1); + h = p2p_c.gamma(tp.ncascades,tp.tau2,t); % Generate the n-cascade impulse response + tid = find(cumsum(h)*dt>=.999,1,'first'); % Shorten the filter if needed to speed up the code. + if ~isempty(tid) + h = h(1:tid); + else + disp('Warning: gamma hdr might not have a long enough time vector'); + end + trl.imp_resp = h; % close enough to use h + impFrames = [0:(length(trl.imp_resp)-1)]; resp = zeros(1,length(t)+length(trl.imp_resp)); % zero stuff out - + %reduction in spikes by inter-spike intervals: interspike = [0,diff(trl.t(ptid))]; trl.spikes = trl.spikes.*(1-exp(-tp.refrac*(interspike+tp.delta))); for i=1:length(trl.spikes) - id = find(t>trl.t(ptid(i)),1,'first'); - resp(id+impFrames) = ... - resp(id+impFrames) + trl.spikes(i)*trl.imp_resp; + id = find(t>trl.t(ptid(i)),1,'first'); + resp(id+impFrames) = ... + resp(id+impFrames) + trl.spikes(i)*trl.imp_resp; end resp = p2p_c.nonlinearity(tp, resp); trl.maxresp = max(resp); % detection when maxresp goes above a threshold @@ -752,209 +704,209 @@ end - %% utilities - function out=chronaxie(p,pw) - out = p.amp./(p.tau*(1-exp(-pw/p.tau))); - end - - function y=gamma(n,k,t) - % y=gamma(n,k,t) - % returns a gamma function on vector t - % y=(t/k).^(n-1).*exp(-t/k)/(k*factorial(n-1)); - % which is the result of an n stage leaky integrator. - - % 6/27/95 Written by G.M. Boynton at Stanford University - % 4/19/09 Simplified it for Psychology 448/538 at U.W. - % - y = (t/k).^(n-1).*exp(-t/k)/(k*factorial(n-1)); - y(t<0) = 0; - end - function y = nonlinearity(tp,x) - if ~isfield(tp, 'scFac'); scFac = 1; - else scFac = tp.scFac; end - % some of our favorite static nonlinearities: - switch tp.model - case 'sigmoid' - y = scFac .* x.^tp.power./(x.^tp.power + tp.sigma.^2); - case 'normcdf' - y = normcdf(x, tp.mean, tp.sigma); - y(y<0) = 0; - case 'weibull' - y = scFac*p2p_c.weibull(tp,x); - case 'power' - y = scFac*x.^tp.power; - case 'exp' - y = scFac*x.^tp.k; - case 'compression' - y = tp.power.*tanh(x/tp.power); - case 'linear' - y = x; - end + %% utilities + function out=chronaxie(p,pw) + out = p.amp./(p.tau*(1-exp(-pw/p.tau))); + end + + function y=gamma(n,k,t) + % y=gamma(n,k,t) + % returns a gamma function on vector t + % y=(t/k).^(n-1).*exp(-t/k)/(k*factorial(n-1)); + % which is the result of an n stage leaky integrator. + + % 6/27/95 Written by G.M. Boynton at Stanford University + % 4/19/09 Simplified it for Psychology 448/538 at U.W. + % + y = (t/k).^(n-1).*exp(-t/k)/(k*factorial(n-1)); + y(t<0) = 0; + end + function y = nonlinearity(tp,x) + if ~isfield(tp, 'scFac'); scFac = 1; + else scFac = tp.scFac; end + % some of our favorite static nonlinearities: + switch tp.model + case 'sigmoid' + y = scFac .* x.^tp.power./(x.^tp.power + tp.sigma.^2); + case 'normcdf' + y = normcdf(x, tp.mean, tp.sigma); + y(y<0) = 0; + case 'weibull' + y = scFac*p2p_c.weibull(tp,x); + case 'power' + y = scFac*x.^tp.power; + case 'exp' + y = scFac*x.^tp.k; + case 'compression' + y = tp.power.*tanh(x/tp.power); + case 'linear' + y = x; end - function [p] = weibull(params, x) - % [p] = Weibull(params, x) - % - % The Weibull function based on this equation: - % - % k = (-log((1-e)/(1-g)))^(1/b) - % f(x) = 1 - ((1-g) * exp(-(k*x/t).^b)) - % - % Where g is performance expected at chance, e is performance level that - % defines the threshold, b is the slope of the Weibull function, and t is - % the threshold - % - % Inputs: - % params A structure containing the parameters of the Weibull - % function: - % b Slope - % t Stimulus intensity threshold as defined by 'params.e'. - % When x = 'params.t', then y = 'params.e' - % g Performance expected at chance, proportion - % e Threshold performance, proportion - % - % x Intensity values of the stimuli - % - % Output: - % p Output of the Weibull function as a function of the - % intensity values, x - - % Written by G.M. Boynton - 11/13/2007 - % Edited by Kelly Chang - February 13, 2017 - % Edited by Ione Fine - February 22, 2017 - - if ~isfield(params, 'g') - params.g = 0.5; - end - if ~isfield(params, 'e') - params.e = (0.5)^(1/3); - end - k = (-log((1-params.e)/(1-params.g)))^(1/params.b); - p = 1 - ((1-params.g) * exp(-(k*x/params.t).^params.b)); - end - - %% transforms - % transforms - function c2v_out = c2v(c, z) - % takes in imaginary numbers, and finds out where cortical values are in visual space (map) - c2v_out = c.k*log(z + c.a); - end - function v2c_out = v2c(c, z) - % takes in imaginary numbers, places visual values into the cortical grid (mapinv) - v2c_out = exp(z/c.k)-c.a; - end - - %% plotting functions - % plotting functions - function plotcortgrid(img, c, varargin) - % plotcortgrid(img, c) - % plotcortgrid(img, c, cmap,figNum, evalstr) - % takes as input: - % cortical image - % the structure c that defines the cortical surface - % optional arguments: - % colormap, figure number and a string to evaluate - % (e.g. ''title('''corticalsurface''')' or 'subplot(1, 2,1)'; - - if nargin<3 || isempty(varargin{1}); cmap = gray(256); else cmap = varargin{1}; end - if nargin<4 || isempty(varargin{2}); figNum = 1; else figNum = varargin{2}; end - if nargin<5 || isempty(varargin{3}); evalstr = ''; else evalstr = varargin{3}; end - - if isfield(c,'cropPix') - img(c.cropPix) = NaN; - img= img+2; - cmap = [0,0,0;cmap]; - end + end + function [p] = weibull(params, x) + % [p] = Weibull(params, x) + % + % The Weibull function based on this equation: + % + % k = (-log((1-e)/(1-g)))^(1/b) + % f(x) = 1 - ((1-g) * exp(-(k*x/t).^b)) + % + % Where g is performance expected at chance, e is performance level that + % defines the threshold, b is the slope of the Weibull function, and t is + % the threshold + % + % Inputs: + % params A structure containing the parameters of the Weibull + % function: + % b Slope + % t Stimulus intensity threshold as defined by 'params.e'. + % When x = 'params.t', then y = 'params.e' + % g Performance expected at chance, proportion + % e Threshold performance, proportion + % + % x Intensity values of the stimuli + % + % Output: + % p Output of the Weibull function as a function of the + % intensity values, x - fH=figure(figNum); - eval(evalstr); colormap(cmap); - if ~isempty(img) - image(c.x, c.y, img); hold on - end - xlabel('mm'); ylabel('mm') - set(gca,'YDir','normal'); - plot(c.v2c.gridAngZ, '-', 'Color', c.gridColor); - plot(c.v2c.gridEccZ, '-', 'Color', c.gridColor); - - axis equal; axis tight - set(gca,'XLim',[min(c.x(:)),max(c.x(:))]); - set(gca,'YLim',[min(c.y(:)),max(c.y(:))]); - drawnow; - end - function plotretgrid(img, v, varargin) - % plotretgrid(img, c) - % plotretgrid(img, c, cmap,figNum, evalstr) - % takes as input: - % retinal image - % the structure v that defines the retinal surface - % optional arguments: - % colormap, figure number and a string to evaluate - % (e.g. ''title('''corticalsurface''')' or 'subplot(1, 2,1)'; - - if nargin<3 || isempty(varargin{1}); cmap = gray(256); else; cmap = varargin{1}; end - if nargin<4 || isempty(varargin{2}); figNum = 1; else; figNum = varargin{2}; end - if nargin<5 || isempty(varargin{3}); evalstr = ''; else; evalstr = varargin{3}; end - - fH=figure(figNum); hold on - eval(evalstr); - image(v.x, v.y, img); hold on - - colormap(cmap); set(gca,'YDir','normal'); - - plot(v.zAng,'-','Color', v.gridColor); plot(v.zEcc,'-','Color', v.gridColor); - plot(-v.zAng,'-','Color', v.gridColor); plot(-v.zEcc,'-','Color', v.gridColor); - - axis equal; axis tight - xlabel('degrees'); ylabel('degrees') - set(gca,'XLim',[min(v.x(:)),max(v.x(:))]); - set(gca,'YLim',[min(v.y(:)),max(v.y(:))]); - drawnow; - end - function p = fit_ellipse_to_phosphene(img,v) - M00 = sum(sum(img)); - M10 = sum(sum(v.X.*img)); M01 = sum(sum(v.Y.*img)); M11 = sum(sum(v.X.*v.Y.*img)); - M20 = sum(sum(v.X.^2.*img)); M02 = sum(sum(v.Y.^2.*img)); - p.x0 = M10/M00; p.y0 = M01/M00; - mu20 = M20/M00 - p.x0^2; mu02 = M02/M00 - p.y0^2; mu11 = M11/M00 - p.x0*p.y0; - a = (mu20+mu02)/2; b = .5*sqrt(4*mu11^2+(mu20-mu02)^2); - lambda_1 = a+b; lambda_2 = a-b; - p.theta = -.5*atan2(2*mu11,mu20-mu02); - p.sigma_x = 2*sqrt(lambda_1); p.sigma_y = 2*sqrt(lambda_2); - end - function fillSymbols(h,colList) - if ~exist('h', 'var'); h = get(gca,'Children'); end - for i=1:length(h) - if ~exist('colList','var'); col = get(h(i),'Color'); - else - if iscell(colList); col = colList{i}; - else; col = colList(i,:); end - end - set(h(i),'MarkerFaceColor',col); - end + % Written by G.M. Boynton - 11/13/2007 + % Edited by Kelly Chang - February 13, 2017 + % Edited by Ione Fine - February 22, 2017 + + if ~isfield(params, 'g') + params.g = 0.5; end - function draw_ellipse(trl, figNum, spstr, varargin) - figure(figNum); hold on - eval(spstr); - if nargin >3 - eye = varargin{1}; - else - eye = 1; - end + if ~isfield(params, 'e') + params.e = (0.5)^(1/3); + end + k = (-log((1-params.e)/(1-params.g)))^(1/params.b); + p = 1 - ((1-params.g) * exp(-(k*x/params.t).^params.b)); + end + + %% transforms + % transforms + function c2v_out = c2v(c, z) + % takes in imaginary numbers, and finds out where cortical values are in visual space (map) + c2v_out = c.k*log(z + c.a); + end + function v2c_out = v2c(c, z) + % takes in imaginary numbers, places visual values into the cortical grid (mapinv) + v2c_out = exp(z/c.k)-c.a; + end - if nargin>4 - lineColor = varargin{2}; + %% plotting functions + % plotting functions + function plotcortgrid(img, c, varargin) + % plotcortgrid(img, c) + % plotcortgrid(img, c, cmap,figNum, evalstr) + % takes as input: + % cortical image + % the structure c that defines the cortical surface + % optional arguments: + % colormap, figure number and a string to evaluate + % (e.g. ''title('''corticalsurface''')' or 'subplot(1, 2,1)'; + + if nargin<3 || isempty(varargin{1}); cmap = gray(256); else cmap = varargin{1}; end + if nargin<4 || isempty(varargin{2}); figNum = 1; else figNum = varargin{2}; end + if nargin<5 || isempty(varargin{3}); evalstr = ''; else evalstr = varargin{3}; end + + if isfield(c,'cropPix') + img(c.cropPix) = NaN; + img= img+2; + cmap = [0,0,0;cmap]; + end + + fH=figure(figNum); + eval(evalstr); colormap(cmap); + if ~isempty(img) + image(c.x, c.y, img); hold on + end + xlabel('mm'); ylabel('mm') + set(gca,'YDir','normal'); + plot(c.v2c.gridAngZ, '-', 'Color', c.gridColor); + plot(c.v2c.gridEccZ, '-', 'Color', c.gridColor); + + axis equal; axis tight + set(gca,'XLim',[min(c.x(:)),max(c.x(:))]); + set(gca,'YLim',[min(c.y(:)),max(c.y(:))]); + drawnow; + end + function plotretgrid(img, v, varargin) + % plotretgrid(img, c) + % plotretgrid(img, c, cmap,figNum, evalstr) + % takes as input: + % retinal image + % the structure v that defines the retinal surface + % optional arguments: + % colormap, figure number and a string to evaluate + % (e.g. ''title('''corticalsurface''')' or 'subplot(1, 2,1)'; + + if nargin<3 || isempty(varargin{1}); cmap = gray(256); else; cmap = varargin{1}; end + if nargin<4 || isempty(varargin{2}); figNum = 1; else; figNum = varargin{2}; end + if nargin<5 || isempty(varargin{3}); evalstr = ''; else; evalstr = varargin{3}; end + + fH=figure(figNum); hold on + eval(evalstr); + image(v.x, v.y, img); hold on + + colormap(cmap); set(gca,'YDir','normal'); + + plot(v.zAng,'-','Color', v.gridColor); plot(v.zEcc,'-','Color', v.gridColor); + plot(-v.zAng,'-','Color', v.gridColor); plot(-v.zEcc,'-','Color', v.gridColor); + + axis equal; axis tight + xlabel('degrees'); ylabel('degrees') + set(gca,'XLim',[min(v.x(:)),max(v.x(:))]); + set(gca,'YLim',[min(v.y(:)),max(v.y(:))]); + drawnow; + end + function p = fit_ellipse_to_phosphene(img,v) + M00 = sum(sum(img)); + M10 = sum(sum(v.X.*img)); M01 = sum(sum(v.Y.*img)); M11 = sum(sum(v.X.*v.Y.*img)); + M20 = sum(sum(v.X.^2.*img)); M02 = sum(sum(v.Y.^2.*img)); + p.x0 = M10/M00; p.y0 = M01/M00; + mu20 = M20/M00 - p.x0^2; mu02 = M02/M00 - p.y0^2; mu11 = M11/M00 - p.x0*p.y0; + a = (mu20+mu02)/2; b = .5*sqrt(4*mu11^2+(mu20-mu02)^2); + lambda_1 = a+b; lambda_2 = a-b; + p.theta = -.5*atan2(2*mu11,mu20-mu02); + p.sigma_x = 2*sqrt(lambda_1); p.sigma_y = 2*sqrt(lambda_2); + end + function fillSymbols(h,colList) + if ~exist('h', 'var'); h = get(gca,'Children'); end + for i=1:length(h) + if ~exist('colList','var'); col = get(h(i),'Color'); else - lineColor = 'g'; - end - theta = linspace(-pi,pi,101); - - for e=1:length(eye) - r = sqrt( (trl.ellipse(eye(e)).sigma_x*trl.ellipse(eye(e)).sigma_y)^2./ ... - (trl.ellipse(eye(e)).sigma_y^2*cos(theta).^2 + trl.ellipse(eye(e)).sigma_x^2*sin(theta).^2)); - x = trl.ellipse(eye(e)).x+r.*cos(theta-trl.ellipse(eye(e)).theta); - y = trl.ellipse(eye(e)).y+r.*sin(theta-trl.ellipse(eye(e)).theta); - plot(x,y,'-','LineWidth',1,'Color',lineColor); + if iscell(colList); col = colList{i}; + else; col = colList(i,:); end end + set(h(i),'MarkerFaceColor',col); + end + end + function draw_ellipse(trl, figNum, spstr, varargin) + figure(figNum); hold on + eval(spstr); + if nargin >3 + eye = varargin{1}; + else + eye = 1; + end + + if nargin>4 + lineColor = varargin{2}; + else + lineColor = 'g'; + end + theta = linspace(-pi,pi,101); + for e=1:length(eye) + r = sqrt( (trl.ellipse(eye(e)).sigma_x*trl.ellipse(eye(e)).sigma_y)^2./ ... + (trl.ellipse(eye(e)).sigma_y^2*cos(theta).^2 + trl.ellipse(eye(e)).sigma_x^2*sin(theta).^2)); + x = trl.ellipse(eye(e)).x+r.*cos(theta-trl.ellipse(eye(e)).theta); + y = trl.ellipse(eye(e)).y+r.*sin(theta-trl.ellipse(eye(e)).theta); + plot(x,y,'-','LineWidth',1,'Color',lineColor); end + end end +end diff --git a/temporal_fit_Winawer_data.m b/temporal_fit_Winawer_data.m index 0b76453..ae96f80 100644 --- a/temporal_fit_Winawer_data.m +++ b/temporal_fit_Winawer_data.m @@ -1,6 +1,6 @@ %temporal_fit_Winawer_data.m % -% fits Winawer bightness data, corresponds to figure 2c +% fits Winawer bightness data, corresponds to figure 2c clear all; close all @@ -8,41 +8,62 @@ T.brightness( T.brightness==-1000) = NaN; % weird hack because otherwise brightness stuck as a cell array tp = p2p_c.define_temporalparameters(); - +tp.model = 'compression'; +tp.simdur = 3; rng(11171964) % fix the random number generator, used Geoff's birthday in paper -rval = randperm(size(T,1)); -rsamp(1).val = rval(1:ceil(length(rval)/2)); -rsamp(2).val = rval(ceil(length(rval)/2)+1:end); +gvals = find(T.amp>0); % only include trials where there was a simulus and a reported brightness +gvals = gvals(randperm(length(gvals))); +rsamp(1).val = gvals(1:ceil(length(gvals)/2)); +rsamp(2).val = gvals(ceil(length(gvals)/2)+1:end); titleStr = {'Training Data', 'Test Data'}; -FITFLAG = 0; +FITFLAG = 1; if FITFLAG Texp = T(rsamp(1).val,:); - freeParams = {'power', 'tau2'}; + freeParams = {'power'}; tp = fit('p2p_c.fit_brightness',tp,freeParams,Texp); end -colorList = [1 0 0; 0 1 0; 1 .7 0; .3 .3 1; 0 0 1 ]; % roughly match colors to Winawer paper -for r = 1:2 % test and train data - for site =1:5 +colorList = [1 0 0; 0 1 0; 1 .7 0; 1 .3 1; 0 0 1 ]; % roughly match colors to Winawer paper +alpha = 0.5; sz = 150; sd = .1; % jitter for plotting + +for site =1:5 + % plot model predictions vs. data + for r = 1:2 % test and train data eid = intersect(rsamp(r).val, find(T.electrode==site)); - Texp = T(eid, :); - trl= p2p_c.loop_convolve_model(tp, Texp); - subplot(1, 2, r) - x = [trl.maxresp]; y = [Texp.brightness] ; - x= reshape(x, length(x), 1); y = reshape(y, length(x), 1); - ind = ~isnan(x) & ~isnan(y); - if sum(ind)>0 - plot(x, y, 'o', 'MarkerFaceColor', colorList(site, :), 'MarkerEdgeColor', 'none'); hold on - corrval = corr(x(ind), y(ind)); + if ~isempty(eid) + Texp = T(eid, :); + trl= p2p_c.loop_convolve_model(tp, Texp); + subplot(1, 3, r) + x = [trl.maxresp]; y = [Texp.brightness] ; + x = reshape(x, length(x), 1); y = reshape(y, length(x), 1); + x = x + sd*(rand(size(x))-0.5); y = y + sd*(rand(size(x))-0.5); + h = scatter(x, y, sz, 'MarkerFaceColor', colorList(site, :), 'MarkerEdgeColor', 'none', 'MarkerFaceAlpha', alpha); + hold on + t = text(15, site/2, ['corr = ', num2str(round(corr(x, y), 3))]); + set(t, 'Color',colorList(site, :)); xlabel('Model Estimate'); - ylabel('Reported Brightness'); - set(gca, 'XLim', [0 25]); - set(gca, 'YLim', [0 11]); - t = text(15, site/2, ['corr = ', num2str(round(corrval, 3))]); - set(t, 'Color',colorList(site, :) ) - title(titleStr{r}); end + title(titleStr{r}); end + % look at correlations with total charge + subplot(1, 3, 3) + eid = find(T.electrode==site & ~isnan(T.brightness)); + if ~isempty(eid) + Texp = T(eid, :); + x = [Texp.totalcharge]; y = [Texp.brightness] ; + x = x + sd*(rand(size(x))-0.5); y = y + sd*(rand(size(x))-0.5); + h = scatter(x, y, sz, 'MarkerFaceColor', colorList(site, :), 'MarkerEdgeColor', 'none', 'MarkerFaceAlpha', alpha); + hold on + t = text(15, site/2, ['corr = ', num2str(round(corr(x, y), 3))]); + set(t, 'Color',colorList(site, :)); + xlabel('Total Charge'); + end +end +for i = 1:3 + subplot(1, 3, i) + ylabel('Reported Brightness'); +% set(gca, 'XLim', [0 25]); +% set(gca, 'YLim', [0 11]); end